*! version 2.0 dgc/mh Sep 1997   STB-41 ssa10.1

program define dyrate
    version 5.0
    local varlist "optional ex min(1)"
    local if "opt"
    local in "opt"
    local using "opt noprefix"
    local weight "fweight iweight pweight"
    local options "Jack CLuster(string) Exposure(string) Graph noLIst noMiss"
    local options "`options' SCale(real 1) SMR(string) L1(string)"
    local options "`options' level(real $S_level) B2(string) noWhisker"
    local options "`options' REPlace *"
    parse "`*'"
    parse "`varlist'", parse(" ")
    local d "`1'"
    mac shift
    local ev "`*'"
    local y "`exposure'"
    if "`exposure'"!="" {
       confirm var `exposure'
    }
    else {
       di in re "must specify exposure()
    }
    tempvar y
    gen `y'=`exposure'/`scale'
    if "`cluster'" != "" {
       local jack "jack"
       confirm var `cluster'
    }            
    if "`smr'"!="" {
       confirm var `smr'
    }
    if "`list'"=="nolist" & "`graph'"=="" & "`using'"=="" {
       di in re "this command has no output"
       exit
    }
    if "`weight'"!="" {
       tempvar W
       gen `W' `exp'
    }
    preserve
    if "`weight'"=="iweight"|"`weight'"=="pweight" {
       local jack "jack"
    }
    tempvar touse D Y wd wy rate jm jv jn cll clh
    mark  `touse' `if' `in' 
    if "`miss'"=="nomiss" {
       markout `touse' `W' `smr' `d' `y' `ev' 
    }
    else {
       markout `touse' `W' `smr' `d' `y'
    }
    if "`W'"=="" {
       local W = 1
    }
    qui {
       drop if !`touse'
       if "`smr'"!="" {
          replace `y' = `smr'*`y'
       }
       local inc = _N
       gen `wd' = `W'*`d'
       gen `wy' = `W'*`y'
       if "`ev'`cluster'"!="" {
          sort `ev' `cluster'
       }
       if "`jack'"!="" {
          if "`cluster'"!="" {
             by `ev' `cluster': replace `wd' = sum(`wd')
             by `ev' `cluster': replace `wy' = sum(`wy')
             by `ev' `cluster': drop if _n!=_N
             local cin = _N
          }
          if "`ev'"!="" {
             by `ev': gen `D' = sum(`wd')
             by `ev': gen `Y' = sum(`wy')
             by `ev': gen `jn' = _n
             by `ev': gen double `jm' = (`D'[_N]-(`wd'))/(`Y'[_N]-(`wy')) 
             by `ev': replace `jm' = _N*log(`D'[_N]/`Y'[_N]) /*
                */ - (_N-1)*log(`jm')
             gen double `jv' = `jm'^2
             by `ev': replace `jm' = sum(`jm')
             by `ev': replace `jv' = sum(`jv')
             by `ev': keep if _n == _N
          }
          else {
             gen `D' = sum(`wd')
             gen `Y' = sum(`wy')
             gen `jn' = _n
             gen double `jm' = (`D'[_N]-(`wd'))/(`Y'[_N]-(`wy')) 
             replace `jm' = _N*log(`D'[_N]/`Y'[_N]) /*
                */ - (_N-1)*log(`jm')
             gen double `jv' = `jm'^2
             replace `jm' = sum(`jm')
             replace `jv' = sum(`jv')
             keep if _n == _N
          }
       }
       else {
          if "`ev'"!="" {
             by `ev': gen `D' = sum(`wd')
             by `ev': gen `Y' = sum(`wy')
             by `ev': keep if _n == _N
          }
          else {
             gen `D' = sum(`wd')
             gen `Y' = sum(`wy')
             keep if _n == _N
          }
       } 
       gen `rate' = `D'/`Y'
       if "`jack'"!="" {
          replace `jm' = `jm'/`jn'
          replace `jv' = sqrt((`jv' - `jn'*(`jm'^2))/(`jn'*(`jn' - 1)))
       }
       else {
          gen `jm' = log(`D'/`Y')
          gen `jv' = sqrt(1/`D')
       }
       local z = invnorm(0.5 + `level'/200)
       gen `cll' = cond(`D'==0, ., exp(`jm' - `z'*`jv'))
       gen `clh' = cond(`D'==0, ., exp(`jm' + `z'*`jv'))
    }
    di
    if "`smr'"!="" {
       di in gr "Estimated SMRs and lower/upper bounds" /*
          */ " of `level'% confidence intervals"
    }
    else {
       di in gr "Estimated rates " _continue 
       if `scale'!=1 {
          di in gr "(per `scale') " _continue
       }
       di in gr "and lower/upper bounds of `level'% confidence intervals"
    }
    di in ye `inc' in gr _col(8) "records included in the analysis "
    if "`jack'"!="" {
       di in gr "The jackknife was used to calculate confidence intervals"
       if "`cluster'"!="" {
          di in ye `cin' in gr _col(8) "clusters, " /*
             */  "identified by the variable: " in ye "`cluster'"
       }
    }
    if "`ev'"=="" & "`graph'"!="" {
       di in bl "Warning: graphical output suppressed - no groups to plot"
       local graph
    }
    if "`graph'"!="" {
       parse "`ev'", parse(" ,")
       if "`2'"!="" {
          tempvar g
          qui egen `g'=group(`ev')
          local b2 "Groups by: `ev'"
       }
       else {
          local g "`ev'"
          local b2 "`ev'"
          local listg "noobs"
       }
    }
    else {
       local listg "noobs"
    }
    rename `D' _D
    rename `cll' _Lower
    rename `clh' _Upper
    format _D %4.0f 
    if "`smr'"!="" {
       rename `Y' _E
       format _E %8.2f
       rename `rate' _SMR
       _setfmt _SMR _Lower _Upper, w(8) s(4)
       if "`list'"!="nolist" {
          list `ev' _D _E _SMR _Lower _Upper, `listg' nodisp
       }
       if "`graph'"!="" {
          if "`l1'"=="" {
             local l1 "SMR"
          }
          if "`whisker'"=="nowhisker" {
             graph _SMR `g', s(o) c(.) /*
                */ l1("`l1'") b2("`b2'") `options' 
          }
          else {
          graph _SMR _Lower _Upper `g', s(oii) c(.||) /*
             */ l1("`l1'") b2("`b2'") `options' 
          }
       }
    }
    else {
       rename `Y' _Y
       _setfmt _Y, w(8) s(5)
       rename `rate' _Rate
       _setfmt _Rate _Lower _Upper, w(8) s(4)
       if "`list'"!="nolist" {
          list `ev' _D _Y _Rate _Lower _Upper, `listg' nodisp
       }
       if "`graph'"!="" {
          if "`l1'"=="" {
             if `scale'!=1 {
                local l1 "Rate (per `scale')"
             }
             else {
                local l1 "Rate"
          }
          if "`whisker'"=="nowhisker" {
             graph _Rate `g', s(o) c(.) /*
                */ l1("`l1'") b2("`b2'") `options' 
          }
          else {
             graph _Rate _Lower _Upper `g', s(oii) c(.||) /*
                */ l1("`l1'") b2("`b2'") `options' 
          }
       }
    }
    di
    if "`using'"!="" {
       label var _D "# Failures"
       label var _Lower "Lower `level'% confidence limit"
       label var _Upper "Upper `level'% confidence limit"
       if "`smr'"=="" {
          keep `ev' _D _Y _Rate _Lower _Upper
          label var _Y "Person-time observation"
          label var _Rate "Rate estimate" 
       }
       else {
          keep `ev' _D _E _SMR _Lower _Upper
          label var _E "Expected # failures"
          label var _Rate "SMR estimate" 
       } 
       desc
       save `using', `replace'
    }
    end

* Program to compute output format

program define _setfmt
    version 5.0
    local varlist "req ex"
    local options "Width(int 8) Sigfig(int 3)"
    parse "`*'"
    confirm ex var `varlist'
    parse "`varlist'", parse(" ")
    tempvar lv 
    qui gen double `lv' = .
    local after = 0
    local minus = 0
    while "`1'"!="" {
       local v "`1'"
       mac shift
       quietly {
          count if `v'<0
          local minus = (_result(1)>0) | `minus'
          replace `lv' = log10(abs(`v')) if `v'!=0
          summ `lv' 
          local count = _result(1)
          if `count'>0 {
             local aneed = int(`sigfig' - _result(5))            
             if `aneed' > `after' {
                local after "`aneed'"
             }
          }
          local lname = length("`v'")
          if `lname'>`width' {
             local width  "`lname'"
          }
       }
    }
    local amax = `width' - `minus' - 1
    if `after' > `amax' {
       local after = `amax'
    }
    format `varlist' %`width'.`after'f
    end