*! version 2.0  dgc/mh 8 Sept 1995
program define tabodds
version 3.1
local varlist "req ex min(1) max(2)"
local if "opt"
local in "opt"
local weight "fweight"

#delimit ;
local options "Binomial(string) Trend Graph
Level(int $S_level) *";
#delimit cr

parse "`*'"
parse "`varlist'", parse(" ")
local d `1'
recast float `d'

tempname md mh nd nh nt v het chitr ptrend one
tempvar touse odds h ci_low ci_high efac W n e

if "`2'"!=""{
local xp `2'
}
else {
gen `one'=1
local xp `one'
}

preserve

local level = `level'*0.005 +0.50

* checks that weight and binomial are not both present

if "`binomial'"!="" & "`weight'"!="" {
        di in re "Weight not allowed with binomial frequency records"
        exit
        }

/* checks whether the response variable is coded 0/1 for individual 
   or frequency records*/

if "`binomial'" == "" {
    qui count if `d' != float(0) & `d' != float(1) & `d' != .
    if _result(1)>0 {
        di in re "Response `d' not coded 0/1"
        exit
    }
}

if "`binomial'"=="" {
   gen `h'= 1-`d'
}
else {
   gen `h' = `binomial' - `d' 
}

* sets frequency weights

if "`weight'"!="" {
        gen `W' `exp'
}
else {
        gen `W' = 1
}


quietly {

drop if `d'==. | `xp'==.
if "`in'"!="" | "`if'"!="" {
keep `in' `if'
}

replace `d' = `d' * `W'
replace `h' = `h' * `W'

sort  `xp'
by  `xp': replace `d' = sum(`d')
by  `xp': replace `h' = sum(`h')
by  `xp': keep if _n==_N

gen `odds'=`d'/`h'
gen `efac'  = exp(invnorm(`level')*sqrt(1/`d'+1/`h'))
gen `ci_low'=`odds'/`efac'
gen `ci_high'=`odds'*`efac'

        summ `xp' [fweight = `d']
        scalar `md' = _result(3)
        scalar `nd' = _result(1)
        summ `xp' [fweight = `h']
        scalar `mh' = _result(3)
        scalar `nh' = _result(1)
        scalar `nt' = `nd' + `nh'
        gen `n' = `d' + `h'        

if "`trend'"!=""{
        summ `xp' [fweight = `n']
        scalar `v' = _result(4)*`nt'/(`nd'*`nh')
        scalar `chitr' = (`md' - `mh')^2/`v'
        scalar `ptrend' = chiprob(1,`chitr')
}
else {
        gen `e'=`nd'*`n'/`nt'
        replace `e'=(`d'-`e')*(`d'-`e')/`n'
        replace `e'=sum(`e')
        gen `het' = `e'[_N] * `nt'*(`nt'-1)/(`nd'*`nh')
}
        

        rename `d' _D
        rename `h' _H
        rename `odds' _odds
        rename `ci_low' ci_low
        rename `ci_high' ci_high

        format _odds ci_low ci_high %8.3f
}

        if "`graph'"!=""&"`2'"!=""{
        di "`options'"
                graph _odds `xp', c(l) `options'
        }

di in blue _n "table of cases (D), controls (H), and odds (D/H)"
if "`2'"!=""{
  list  `xp' _D _H _odds ci_low ci_high, noobs
}
else {
list  _D _H _odds ci_low ci_high, noobs
}

if "`2'"!=""{
  if "`trend'"!=""{
        di in gr _n(2) "Chisq test for trend = " /*
        */ in ye %8.4f  `chitr' " ( 1 df, p = " %6.3f `ptrend' " )"
  }
 
  else {
        di in gr _n "Chisq test for unequal odds = " /*
        */ in ye %8.2f  `het' " (" _N-1 " df, p = " %6.3f = chiprob(_N-1, `het') " )"
  }
}
end