*! version 2.0  dgc/mh 8 Sept 1995              (STB-27: ssa7)
program define mhrate
version 3.1
local varlist "req ex min(2)"
local if "opt"
local in "opt"
local options "Compare(string) BY(string) Exposure(string) Level(integer $S_level)"
parse "`*'"
parse "`varlist'", parse(" ")

local d `1'
local e `2'

preserve

if "`exposure'"==""{
	di in re "Needs person-years or expected numbers"
	exit
}


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


if "`compare'" == "" {
   quietly inspect `e'
   if _result(7) == 2 {
      local compare "compare"
      quietly summarize `e'
      local grp1 = _result(6)
      local grp2 = _result(5)
   }
}
else {
   parse "`compare'", parse(",")
   local grp1 `1'
   local grp2 `3'
}
if "`compare'" != "" {
   di ""
   if "`compare'"!="" {
	if "`str'" != "" {
		di in gr "Mantel-Haenszel estimate of the rate ratio"
		di in gr "Comparing `e'==`grp1' vs `e'==`grp2', controlling for `str'"
	}
	else {
		di in gr "Maximum likelihood estimate of the rate ratio"
		di in gr "Comparing `e'==`grp1' vs `e'==`grp2'"
	}
   }
}
else {
   di ""
   di in gr "Score test for trend of rates with `e'"
   if "`str'" != "" {
      di in gr "controlling for `str'"
   } 
}
if "`by'" != "" {
   di in gr "by `by'"
}

tempvar touse strata d1 d2 y1 y2 dt yt q r u v rr ef cl cu ch pv sm
quietly {

   egen int `touse' = rmiss(`varlist' `by') `if' `in'
   replace `touse' = cond(`touse'>0,.,1) 

   if "`str'" != "" {
      egen int `strata' = group(`str')
   }
   else {
      gen int `strata' = 1
   }
   if "`compare'" != "" {
      keep if `touse'==1 & (`e' == `grp1' | `e' == `grp2')
      sort `strata' `by'
      by `strata' `by': gen `d1' = sum(cond(`e'==`grp1',`d',0)) 
      by `strata' `by': gen `d2' = sum(cond(`e'==`grp2',`d',0)) 
      by `strata' `by': gen `y1' = sum(cond(`e'==`grp1',`exposure',0))
      by `strata' `by': gen `y2' = sum(cond(`e'==`grp2',`exposure',0)) 
      by `strata' `by': keep if _n == _N
      gen `dt' = `d1' + `d2'
      gen `yt' = `y1' + `y2'
      gen `q' = `d1'*`y2'/`yt'
      gen `r' = `d2'*`y1'/`yt'
      drop `d1'  `d2'
      gen `u' = `q' - `r'
      gen `v' = `dt'*`y1'*`y2'/(`yt'^2)
      drop `y1' `y2' `dt' `yt'
      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'" != "" {
	 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
      }
      
      gen `rr' = `q'/`r'
      gen `ef' = exp(invnorm(`level'*0.005 + 0.5) * sqrt(`v'/(`q'*`r')))
   }
   else {
      keep if `touse'==1 
      sort `strata' `by'
      by `strata' `by': gen `dt' = sum(`d') 
      by `strata' `by': gen `yt' = sum(`exposure') 
      by `strata' `by': gen `d1' = sum(`d'*`e') 
      by `strata' `by': gen `y1' = sum(`exposure'*`e')
      by `strata' `by': gen `y2' = sum(`exposure'*`e'*`e')
      by `strata' `by': keep if _n == _N
      gen `u' = `d1' - `dt'*`y1'/`yt'
      drop `d1'
      gen `v' = `dt'*(`y2'-`y1'*`y1'/`yt')/`yt'
      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"
      }
      
      drop `y1' `y2' `yt' `dt' 
      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
      }
      
      gen `rr' = exp(`u'/`v')
      gen `ef' = exp(invnorm(`level'*0.005 + 0.5)/sqrt(`v'))

   }
   gen `cl' = `rr'/`ef'
   gen `cu' = `rr'*`ef'
   gen `ch' = (`u'^2)/`v'
   gen `pv' = chiprob(1,`ch')
}
rename `rr' RR
rename `cl' Lower
rename `cu' Upper
rename `ch' Chisq
rename `pv' p_value
if "`compare'" != "" {
   di in gr _n "RR estimate, lower and upper `level'% confidence limits, and"
   di in gr "chi-squared test for RR=1 (1 degree of freedom)"
}
else {
   di in gr _n "RR estimate, lower and upper `level'% confidence limits, and"
   di in gr "chi-squared test for trend (1 degree of freedom)"_n
   di in blu "The RR estimate is an approximate estimate of the " _n /*
	*/"rate ratio for one unit increase in `e'"
}
format RR Lower Upper Chisq p_value %6.3f

if "`by'" != "" {
   list `by' RR Lower Upper Chisq p_value if _n<_N, noobs nodisplay
   di in gr _n "Mantel-Haenszel estimate controlling for:"in ye " `str' `by'"
   list RR Lower Upper Chisq p_value if _n==_N, noobs nodisplay   
   if `hdf' > 1 {
      di in gr _n "Approx chisq for unequal RRs (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 RR Lower Upper Chisq p_value , noobs nodisplay
}
end