*! version 2.0  dgc/mh 8 Sept 1995
program define mhodds
version 3.1
local varlist "req ex min(2)"
local if "opt"
local in "opt"
local weight "fweight"
local options "Binomial(string) Compare(string) BY(string) Level(integer $S_level)"
parse "`*'"
parse "`varlist'", parse(" ")

local d `1'
local e `2'

preserve

* checks that weight and binomial are not both present

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

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

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

* puts variables to be controlled for in local macro str

local str=""
while "`3'" != "" {
   local str="`str' `3'"
   mac shift
}

* if compare blank and two levels for xvar set compare = "compare"
* and use the two levels as default: high/low

if "`compare'" == "" {
   quietly inspect `e'
   if _result(7) == 2 {
      local compare "compare"
      quietly summarize `e'
      local grp1 = _result(6)
      local grp2 = _result(5)
   }
}

* if compare not blank, use the levels specified.

else {
   parse "`compare'", parse(",")
   local grp1 `1'
   local grp2 `3'
}

* sets up titles

if "`compare'" != "" {

   if "`str'" != "" {
      di in gr _n "Mantel-Haenszel estimate of the odds ratio"
      di in gr "Comparing `e'==`grp1' vs `e'==`grp2', controlling for `str'"
   }
   else {
      di in gr _n "Maximum likelihood estimate of the odds ratio"
      di in gr "Comparing `e'==`grp1' vs `e'==`grp2'"
   }
}
else {

   di in gr _n "Score test for trend of odds with `e'"
   if "`str'" != "" {
      di in gr "controlling for `str'"
   } 
}
if "`by'" != "" {
   di in gr "by `by'"
}

tempvar touse strata h d1 d2 h1 h2 p1 p2 dt ht pt q r u v rr ef cl cu ch pv wt sm

* sets up the variable h for controls

if "`binomial'" == "" {
   gen `h' = 1 - `d'
}
else {
   gen `h' = `binomial' - `d'
}
if "`weight'" == "" {
   gen `wt' = 1
}
else {
   gen `wt' `exp'
}
quietly {
   egen int `touse' = rmiss(`varlist' `by') `if' `in'
   replace `touse' = cond(`touse'>0,.,1) 
   

* sets up variable to index strata which is uniformly 1 for no strata   

   if "`str'" != "" {
      egen int `strata' = group(`str')
   }
   else {
      gen int `strata' = 1
   }
   
* If compare not blank then use odds ratios

   if "`compare'" != "" {
      keep if `touse'==1 & (`e' == `grp1' | `e' == `grp2')
      sort `strata' `by'
      by `strata' `by': gen `d1' = sum(cond(`e'==`grp1',`d',0)*`wt')
      by `strata' `by': gen `d2' = sum(cond(`e'==`grp2',`d',0)*`wt')
      by `strata' `by': gen `h1' = sum(cond(`e'==`grp1',`h',0)*`wt')
      by `strata' `by': gen `h2' = sum(cond(`e'==`grp2',`h',0)*`wt')
      by `strata' `by': keep if _n == _N
      gen `dt' = `d1' + `d2'
      gen `ht' = `h1' + `h2'
      gen `p1' = `h1' + `d1'
      gen `p2' = `h2' + `d2' 
      gen `pt' = `p1' + `p2'
      gen `q' = `d1'*`h2'/`pt'
      gen `r' = `d2'*`h1'/`pt'
      gen `u' = `q' - `r'
      gen `v' = `p1'*`p2'*`dt'*`ht'/((`pt'-1)*(`pt'^2))
      qui count if `v'>0&`v'!=.
      if _result(1)<_N {
	 noi di in blu "WARNING: only " _result(1) " of the " _N " strata " /*
	 */ "formed in this analysis" _n "contribute information about the "/*
	 */"effect of the explanatory variable"
      }
      
* if by not blank then do heterogeneity test
      
      if "`by'" != "" {
	 sort `by'
	 by `by': replace `q' = sum(`q')
	 by `by': replace `r' = sum(`r')
	 by `by': replace `u' = sum(`u')
	 by `by': replace `v' = sum(`v')
	 by `by': keep if _n == _N
	 local N1 = _N + 1
	 set obs `N1'
	 gen `sm' = sum(`q')
	 replace `q' = `sm' if _n == _N
	 replace `sm' = sum(`r') 
	 replace `r' = `sm' if _n == _N
	 replace `sm' = sum(`u') 
	 replace `u' = `sm' if _n == _N
	 replace `sm' = sum(`v') 
	 replace `v' = `sm' if _n == _N
	 local Q = `q'[_N]
	 local R = `r'[_N]
	 replace `sm' = ((`q'*`R' - `r'*`Q')^2)/(`Q'*`R'*`v') if _n < _N
	 replace `sm' = 0 if _n == _N | `v' == 0
	 inspect `sm'
	 local hdf = _result(4)
	 replace `sm' = sum(`sm')
	 local het = `sm'[_N]
      }
      else {
	 replace `q' = sum(`q')
	 replace `r' = sum(`r')
	 replace `u' = sum(`u')
	 replace `v' = sum(`v')
	 keep if _n == _N
      }
      drop `h1' `h2' `d1' `d2' `ht' `dt' `p1' `p2' `pt'
      gen `rr' =`q'/`r'
      gen `ef' = exp(invnorm(`level'*0.005 + 0.5) * sqrt(`v'/(`q'*`r')))
   }
   
* if compare is blank then use trend with xvar

   else {
      keep if `touse'==1
      sort `strata' `by'
      by `strata' `by': gen `dt' = sum(`d'*`wt')
      by `strata' `by': gen `pt' = sum((`d'+`h')*`wt')
      by `strata' `by': gen `d1' = sum(`d'*`e'*`wt')
      by `strata' `by': gen `p1' = sum((`d'+`h')*`e'*`wt')
      by `strata' `by': gen `p2' = sum((`d'+`h')*`e'*`e'*`wt')
      by `strata' `by': keep if _n == _N
      gen `u' = `d1' - (`dt'*`p1'/`pt')
      gen `v' = `dt'*(`pt'-`dt')*(`p2'-`p1'*`p1'/`pt')/(`pt'*(`pt'-1))
      qui count if `v'>0&`v'!=.
      if _result(1)<_N {
	 noi di in blu "WARNING: only " _result(1) " of the " _N " strata " /*
	 */ "formed in this analysis" _n "contribute information about the "/*
	 */"effect of the explanatory variable"
      }
      
* if by is not blank do heterogeneity test      
      
      if "`by'" != "" {
	 sort `by'
	 by `by': replace `u' = sum(`u')
	 by `by': replace `v' = sum(`v')
	 by `by': keep if _n == _N
	 local N1 = _N + 1
	 set obs `N1'
	 gen `sm' = sum(`u')
	 replace `u' = `sm' if _n == _N
	 replace `sm' = sum(`v') 
	 replace `v' = `sm' if _n == _N
	 replace `sm' = `u'^2/`v'
	 replace `sm' = - `sm' if _n == _N
	 inspect `sm'
	 local hdf = _result(4)
	 replace `sm' = sum(`sm')
	 local het = `sm'[_N]
      }
      else {
	 replace `u' = sum(`u')
	 replace `v' = sum(`v')
	 keep if _n == _N
      }
      drop `d1' `dt' `p1' `p2' `pt'
      gen `rr' = exp(`u'/`v')
      gen `ef' = exp(invnorm(`level'*0.005 + 0.5)/sqrt(`v'))
   }
   
* Calculate the OR or trend with confidence limits

   gen `cl' = `rr'/`ef'
   gen `cu' = `rr'*`ef'
   gen `ch' = (`u'^2)/`v'
   gen `pv' = chiprob(1,`ch')

* end of quietly

}
rename `rr' OR
rename `cl' Lower
rename `cu' Upper
rename `ch' Chisq
rename `pv' p_value
display ""
di in gr "OR estimate, lower and upper `level'% confidence limits, and"
di in gr "chi-squared test for OR=1 (1 degree of freedom)"
if "`compare'" == "" {
   di in gr "(The OR estimate is an approximation to the odds ratio "
   di in gr "for one unit increase in `e')"
}
format OR Lower Upper Chisq p_value %6.3f

if "`by'" != "" {
   list `by' OR Lower Upper Chisq p_value if _n<_N , noobs nodisplay

   di in gr _n "Mantel-Haenszel estimate controlling for:" in ye " `str' `by'" 
   list OR Lower Upper Chisq p_value if _n==_N , noobs nodisplay
   if `hdf' > 1 {
      di in gr _n "Approx chisq for unequal ORs (effect modification)" /*
      */ in ye %8.2f =`het' " (" `hdf'-1 " df, p = " %5.3f = chiprob(`hdf'-1, `het') ")"
   }
   else {
      di in gr _n "Too few informative strata for test for unequal OR's"
   }
}
else {
list OR Lower Upper Chisq p_value , noobs nodisplay
}

end