*! version 1.1  dgc/mh Mar 1997  STB-40 ssa10

program define stmc
   version 5.0
   st_is
   local varlist "req ex min(1)"
   local if "opt"
   local in "opt"
   local options "Compare(string) FCodes(string) BY(string)"
   local options "`options' ORigin(string) noMiss Level(integer $S_level)"
   parse "`*'"
   parse "`varlist'", parse(" ")
   local e "`1'"
   mac shift
   local str "`*'"
   if "`by'"!="" {
      confirm ex var `by'
   }
   if "`origin'"!="" {
      confirm ex var `origin'
   }
   di in gr "Mantel-Cox comparisons"
   st_show
   di
   if "`origin'"!="" {
      di in gr "     time origin:  " in ye "`origin'"
   }      
   if "`fcodes'"!="" {
      di
      di in gr "   failure codes:  " in ye "`fcodes'"
   }

   local d : char _dta[st_d]
   local tout : char _dta[st_t]
   local tin : char _dta[st_t0]
   local weight : char _dta[st_wt]
   local wt : char _dta[st_wv]
   local aw : char _dta[st_w]

   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'" != "" {
      di in gr "Mantel-Haenszel estimates 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 tests 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'"
   }
   di in gr _col(3) "controlling for: " in ye "time (by clicks)" _continue
   if "`str'"!="" {
      di in gr " and: " in ye "`str'"
   }
   else {
      di
   }
   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 id t rtype wx wx2 sw sx sx2 q r u v ef q2 r2 qr
   preserve
   quietly {
      mark `touse' `if' `in' 
      if "`compare'"!="" {
         recode `e' `grp1' = 1 `grp2' = 2  * = .
      }
      markout `touse' `d' `e' `tout' `tin' `origin' `wt'
      if "`miss'"!="" {
         markout `touse' `by' `str'
      }
      keep if `touse' 
      keep `d' `e' `tout' `tin' `by' `str' `wt' `origin' 
      if "`wt'"=="" {
         tempvar wt
         gen `wt' = 1
      }
      if "`fcodes'"=="" {
         replace `d' = (`d'!=0)
      }
      else {
         recode `d' `fcodes' = 1 * = 0
      }
      if "`tin'"=="" {
         local ttype : type `tout'
         tempvar tin
         gen `ttype' `tin' = 0
      }
      if "`origin'"!="" {
         replace `tin' = `tin' - `origin'
         replace `tout' = `tout' - `origin'
      } 
      if "`compare'" != "" {
         gen double `wx' = `wt'*(`e'==1)
         gen double `wx2' = `wx'
      }
      else {
         gen double `wx' = `wt'*`e'
         gen double `wx2' = `wt'*(`e'^2) 
      }     
      gen long `id'=_n
      expand 2
      sort `id'
      by `id': gen double `t' = cond(_n==1, `tin', `tout')
      by `id': gen int `rtype' = cond(_n==1, 2, `d'==0)
      by `id': replace `wt' = -`wt' if _n==2 
      by `id': replace `wx' = -`wx' if _n==2 
      by `id': replace `wx2' = -`wx2' if _n==2 
      sort `by' `strata' `t' `rtype'
      by `by' `str' `t' `rtype': replace `wt'=sum(`wt')
      by `by' `str' `t' `rtype': replace `wx'=sum(`wx')
      by `by' `str' `t' `rtype': replace `wx2'=sum(`wx2')
      by `by' `str' `t' `rtype': keep if _n == _N
      if "`by'`str'"!="" {
         by `by' `str': gen double `sw' = sum(`wt')
         by `by' `str': gen double `sx' = sum(`wx')
         by `by' `str': gen double `sx2' = sum(`wx2')
      }
      else {
         gen double `sw' = sum(`wt')
         gen double `sx' = sum(`wx')
         gen double `sx2' = sum(`wx2')
      }
      keep if `rtype'==0
      * Now there is one record for each failure time
      gen `u' = -`wx' + `wt'*(`sx'-`wx')/(`sw'-`wt')
      gen `v' = -`wt'*((`sx2'-`wx2')/(`sw'-`wt') - /*
              */ ((`sx'-`wx')/(`sw'-`wt'))^2)
   }
   /* Two group comparison */
   if "`compare'"!="" {
      quietly {
         gen `q' = -`wx'*(`sw'-`sx')/(`sw'-`wt')
         gen `r' = `sx'*(`wx'-`wt')/(`sw'-`wt')
         if "`by'"!="" {
            by `by': replace `u' = sum(`u')
            by `by': replace `v' = sum(`v')
            by `by': replace `q' = sum(`q')
            by `by': replace `r' = sum(`r')
            by `by': keep if _n==_N & `v'>0
            local hdf = _N - 1
            gen double `q2' = sum((`q'^2)/`v')
            gen double `r2' = sum((`r'^2)/`v')
            gen double `qr' = sum(`q'*`r'/`v')
            gen _RR = `q'/`r'
            gen `ef' = exp(invnorm(`level'*0.005 + 0.5) * sqrt(`v'/(`q'*`r')))
            gen _Lower = _RR/`ef' 
            gen _Upper = _RR*`ef'
            noi di in gr _n "RR estimate, and lower and upper `level'% "/*
               */ "confidence limits"
            format _RR _Lower _Upper %6.3f
            noi list `by' _RR _Lower _Upper, noobs
            drop `ef' _RR _Lower _Upper
         }
         else {
            local hdf = 0
         }
         replace `u' = sum(`u')
         replace `v' = sum(`v')
         replace `q' = sum(`q')
         replace `r' = sum(`r')
         keep if _n==_N
         gen _RR = `q'/`r'
         gen `ef' = exp(invnorm(`level'*0.005 + 0.5) * sqrt(`v'/(`q'*`r')))
         gen _Lower = _RR/`ef' 
         gen _Upper = _RR*`ef'
         gen _Chisq = `u'^2/`v'
         gen _p_value = chiprob(1, _Chisq)
      }
      di 
      di in gr "Overall Mantel-Haenszel estimate, controlling for: " _continue 
      if "`origin'"=="" {
         di in gr "time " _continue
      }
      else {
         di in gr "time-from-" in ye "`origin' " _continue 
      }
      di in ye "`str' `by'"
      format _RR _Lower _Upper _Chisq %6.3f
      format _p_value %8.5f
      list _RR _Lower _Upper _Chisq _p_value, noobs
      if `hdf'>0 {
         local het = `q2'[1]/_RR[1] + `r2'[1]*_RR[1] - 2*`qr'[1] 
         di in gr _n "Approx chisq for unequal RRs (effect modification)" /*
            */ in ye %8.2f =`het' " (" `hdf' " df, p = " /*
            */ %7.5f = chiprob(`hdf', `het') ")"
      }
   }
   else {
      quietly { 
         if "`by'"!="" {
            by `by': replace `u' = sum(`u')
            by `by': replace `v' = sum(`v')
            by `by': keep if _n==_N & `v'>0
            local hdf = _N - 1
            gen _RR = exp(`u'/`v')
            gen `ef' = exp(invnorm(`level'*0.005 + 0.5) * sqrt(1/`v'))
            gen _Lower = _RR/`ef' 
            gen _Upper = _RR*`ef'
            noi di in gr _n "Approximate RR estimate, and lower and upper "/*
               */ "`level'% confidence limits"
            format _RR _Lower _Upper %6.3f
            noi list `by' _RR _Lower _Upper, noobs
            drop `ef' _RR _Lower _Upper
         }
         else {
            local hdf = 0
         }
         replace `u' = sum(`u')
         replace `v' = sum(`v')
         gen _Chisq = sum(`u'^2/`v')
         local schi2 = _Chisq[_N]
         keep if _n==_N
         gen _RR = exp(`u'/`v')
         gen `ef' = exp(invnorm(`level'*0.005 + 0.5) * sqrt(1/`v'))
         gen _Lower = _RR/`ef' 
         gen _Upper = _RR*`ef'
         replace _Chisq = `u'^2/`v'
         gen _p_value = chiprob(1, _Chisq)
      } 
      format _RR _Lower _Upper _Chisq %6.3f
      format _p_value %8.5f
      di 
      di in gr "Overall Mantel-Haenszel estimate, controlling for: " _continue 
      if "`origin'"=="" {
         di in gr "time " _continue
      }
      else {
         di in gr "time-from-" in ye "`origin' " _continue 
      }
      di in ye "`str' `by'"
      list _RR _Lower _Upper _Chisq _p_value, noobs
      if `hdf'>0 {
         local het = `schi2' - _Chisq[1]
         di in gr _n "Approx chisq for unequal RRs (effect modification)" /*
            */ in ye %8.2f =`het' " (" `hdf' " df, p = " /*
            */ %7.5f = chiprob(`hdf', `het') ")"
      }
   }
   /* Save final rate ratio in S_1 */
   global S_1 = _RR[1]
end