*! version 1.0  dgc/mh Nov 1996  STB-40 ssa10

program define dymh
   version 5.0
   local varlist "req ex min(2)"
   local if "opt"
   local in "opt"
   local weight "iweight pweight fweight"
   local options "Compare(string) BY(string) Exposure(string)"
   local options "`options' noMiss Level(integer $S_level)"
   parse "`*'"
   parse "`varlist'", parse(" ")
   local d "`1'"
   mac shift
   local e "`1'"
   mac shift
   local str "`*'"
   if "`by'"!="" {
      confirm ex var `by'
   }
   local y "`exposure'"
   if "`y'"=="" {
      di in re "must specify exposure()"
      exit
   }
   else {
      confirm var `y'
   }
   if "`weight'" != "" {
      tempvar wt
      gen `wt' `exp'
   }
   preserve
   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'
   }
   di
   if "`compare'" != "" {
      if "`str'" != "" {
         di in gr "Mantel-Haenszel estimate of the rate ratio"
      }
      else {
         di in gr "Maximum likelihood estimate of the rate ratio"
      }
      di in gr _col(3) "comparing: " in ye "`e'==`grp1'" in gr " vs " /*
         */ in ye "`e'==`grp2'" 
   }
   else {
      di in gr "Score test for trend of rates with: " in ye "`e'" 
      di in gr _col(3) "with an " in bl "approximate " in gr /*
         */ "estimate of the " _n /*
         */ _col(3) "rate ratio for one unit increase in " in ye "`e'"
   }
   if "`str'" != "" {
      di in gr _col(3) "controlling for: " in ye "`str'"
   } 
   if "`by'" != "" {
      di in gr _col(3) "by: " in ye "`by'"
   }
   if "`weight'"=="iweight" | "`weight'"=="pweight" {
      di in bl "Warning: `weight's used - "/*
         */ "confidence intervals & p-values may be wrong"
   }

   tempvar touse strata d1 d2 y1 y2 dt yt q r u v rr ef cl cu ch pv sm
   quietly {
      mark `touse' `if' `in' `aw'
      if "`compare'"!="" {
         recode `e' `grp1' = 1 `grp2' = 2  * = .
      }
      markout `touse' `d' `e' `y' `wt'
      if "`wt'"=="" {
         local wt=1
      }
      if "`miss'"!="" {
         markout `touse' `by' `str'
      }
      keep if `touse' 
      * Two group comparison
      if "`compare'"!="" {
         if "`str'`by'"!="" {
            sort `str' `by'
            by `str' `by': gen `d1' = sum(`wt'*cond(`e'==1,`d',0)) 
            by `str' `by': gen `d2' = sum(`wt'*cond(`e'==2,`d',0)) 
            by `str' `by': gen `y1' = sum(`wt'*cond(`e'==1,`y',0))
            by `str' `by': gen `y2' = sum(`wt'*cond(`e'==2,`y',0)) 
            by `str' `by': keep if _n == _N
         }
         else {
            gen `d1' = sum(`wt'*cond(`e'==1,`d',0)) 
            gen `d2' = sum(`wt'*cond(`e'==2,`d',0)) 
            gen `y1' = sum(`wt'*cond(`e'==1,`y',0))
            gen `y2' = sum(`wt'*cond(`e'==2,`y',0)) 
            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'
         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 & `v'>0
            local hdf = _N - 1
            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'
            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
            local hdf = 0
         }
         gen `rr' = `q'/`r'
         gen `ef' = exp(invnorm(`level'*0.005 + 0.5) * sqrt(`v'/(`q'*`r')))
      }
      else {
         if "`str'`by'"!="" {
            sort `str' `by'
            by `str' `by': gen `dt' = sum(`wt'*`d') 
            by `str' `by': gen `yt' = sum(`wt'*`y') 
            by `str' `by': gen `d1' = sum(`wt'*`d'*`e') 
            by `str' `by': gen `y1' = sum(`wt'*`y'*`e')
            by `str' `by': gen `y2' = sum(`wt'*`y'*`e'*`e')
            by `str' `by': keep if _n == _N
         }
         else {
            gen `dt' = sum(`wt'*`d') 
            gen `yt' = sum(`wt'*`y') 
            gen `d1' = sum(`wt'*`d'*`e') 
            gen `y1' = sum(`wt'*`y'*`e')
            gen `y2' = sum(`wt'*`y'*`e'*`e')
            keep if _n == _N
         }
         gen `u' = `d1' - `dt'*`y1'/`yt'
         drop `d1'
         gen `v' = `dt'*(`y2'-`y1'*`y1'/`yt')/`yt'
         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
   di in gr _n "RR estimate, and lower and upper `level'% confidence limits"
   format RR Lower Upper Chisq %6.3f
   format p_value %8.5f
   if "`by'" != "" {
      list `by' RR Lower Upper if _n<_N, noobs nodisplay
      di in gr _n "Overall 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 = " /*
            */ %7.5f = 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
   }
   * Save overall RR in S_1
   global S_1 = RR[_N]
   end