*! version 3.0 Nov 1996     STB-46 sg98

program define rpoisson
  version 4.0
  di
  di in gr "Poisson regression with frailty"
  parse "`*'", parse(" ,")
  local options "IRr Level(integer $S_level) Robust"
  if !("`1'"=="," | "`*'"=="") {
    local varlist "req ex min(1)"
    local if "opt"
    local in "opt"
    local options "`options' CLuster(string) Exposure(string) FVar(real 0)"
    local options "`options' noLog Method(string)"
    parse "`*'"
    local rden `exposure'
    local ifgrp
    if "`cluster'"!="" {
      confirm var `cluster'
      local ifgrp "cl(`cluster')"
    }
    local ifexp
    if "`rden'"!="" {
      confirm var `rden'
      local ifexp "rd(`rden')"
    }
    parse "`varlist'", parse(" ")
    local yv `1'
    mac shift
    local xvars "`*'"
    if "`method'"=="" {
      local method "ps"
    }
    if "`method'"=="ps" {
      local mdesc "maximum pseudolikelihood"
    }
    else if "`method'"=="pe" {
      local mdesc "posterior expectation of variance of frailties"
    }
    else if "`method'"=="ml" {
      local mdesc "maximum likelihood"
    }
    else if "`method'"=="fi" {
      local mdesc "fixed value, `fvar'"
    }
    else {
      di in re "Unknown method option"
      exit
    }
    di in gr _col(3) "frailty variance estimation: " in ye "`mdesc'"
    di
    preserve
    quietly {
      tempvar touse
      mark `touse' `if' `in'
      markout `touse' `varlist' `rden' `cluster'
      count if `touse'
      local nused = _result(1)
      keep if `touse'
    }
    _rpois_1 `varlist', `ifexp' `ifgrp' `log' fv(`fvar') m(`method')
    if _rc!=0 {
        exit _rc
    }
    local ngrp "$S_1"
    local dev "$S_2"
    local fvar "$S_3"
    local mdf "$S_4"
    local logl = "$S_5"
  }
  else {
    if "$S_E_cmd" != "rpoisson" { 
      di in red  "rpoisson was not the last command"
      exit 301  
    }
    local nused "$S_E_nobs"
    local if "$S_E_if"
    local in "$S_E_in"
    local yv "$S_E_depv"
    local rden "$S_E_off"
    local cluster "$S_E_cl"
    local ngrp "$S_E_ncl"
    local xvars "$S_E_xv"
    local logl "$S_E_ll"
    local fvar "$S_E_fvar"
    local mdf "$S_E_mdf"
    local dev "$S_E_dev"
    parse "`*'"
    if "`robust'"!="" {
      tempvar touse
      mark `touse' `if' `in'
      markout `touse' `yv' `xvars' `rden' `cluster'
      if "`cluster'"!="" {
        sort `cluster'
      }
    }
  }
  /* Display the results */
  di
  di in gr "Number of observations = " in ye `nused'
  if "`cluster'"!="" {
    di in gr "Frailty varies by " in ye "`cluster'" " (" `ngrp' " clusters)"
  }
  else {
    di in gr "Frailty operates at the unit level"
  }
  di in gr "Estimated frailty variance = " in ye %8.4f `fvar'
  local rdf = `nused' - `mdf' - 1 
  di in gr "Deviance = " in ye %8.3f `dev' " (" %4.0f `rdf' " df)"
  di
  if "`irr'"!="" {
    local ef "ef(Rate ratio)"
  }
  matrix mlout ,`ef' level(`level')
  if "$S_E_se"=="robust" {
    di 
    di in gr "(Robust standard errors are estimated using Huber's formula)"
  }
  /* robust SEs */
  if "`robust'"!="" & "$S_E_se"!="robust" {
    di
    di in gr "Estimates with robust standard errors:"
    quietly {
      tempvar lp mu wy wt
      matrix beta = get(_b)
      matrix score double `lp' = beta
      gen double `mu' = exp(`lp')
      if "`rden'"!="" {
        replace `mu' = `mu'*`rden'
      }
      gen double `wy' = `lp' + (`yv' - `mu')/`mu'
      if "`cluster'"=="" {
        gen double `wt' = `mu'/(1+ `fvar' * `mu')
        _rpois_2 `wy' `xvars' [iweight=`wt'] if `touse', /*
          */ dep(`yv') robust
      }
      else {
        _rpois_2 `wy' `xvars' [iw = `mu'] if `touse', /*
           */ dep(`yv') robust cl(`cluster') ra(`fvar')
      }
      global S_E_se "robust"
    }
    di
    matrix mlout ,`ef' level(`level')
  }
  global S_E_cmd "rpoisson"
  global S_E_depv "`yv'"
  global S_E_off "`rden'"
  global S_E_cl "`cluster'"
  global S_E_ncl "`ngrp'"
  global S_E_xv "`xvars'"
  global S_E_nobs "`nused'"
  global S_E_fvar "`fvar'"
  global S_E_mdf "`mdf'"
  global S_E_ll "`logl'"
  global S_E_dev "`dev'"
  if "`if'"!="" {
    global S_E_if "`if'"
  }
  if "`if'"!="" {
    global S_E_in "`in'"
  }
end

program define _rpois_1
  version 4.0
  local varlist "req ex min(1)"
  local options "CLuster(string) RDen(string) Fvar(real 0) noLog"
  local options "`options' Method(string)"
  parse "`*'"
  parse "`varlist'", parse(" ")
  local yv `1'
  mac shift 
  local xvars "`*'"
  tempvar mu lp wy wt cn dv so se lt
  quietly {
    gen double `wt' = .
    if "`cluster'" == "" {
      gen double `so' = `yv'
    }
    else {
      sort `cluster'
      by `cluster': gen double `so' = sum(`yv')
      by `cluster': gen double `lt' = (_n==_N)
      replace `so' = 0 if !`lt'
      count if `lt'
      local ngrp = _result(1)
    }
    gen double `se' = 0
    local it = 0
    local jt = 0
    local kt = 0
    local couter 1
    while `couter'  {
      local it=`it'+1
      if `it'==1 {
        gen double `mu' = `yv' + 0.5
        gen double `wy' = log(`mu')
        if  "`rden'"!="" {
          replace `wy' = `wy' - log(`rden')
        }
        regr `wy' `xvars'
        local mdf = _result(3)
        local rdf = _result(5)
        predict double `lp', xb
      }
      /* Iteratively reweighted least squares */
      if "`log'"!="nolog" {
        noi di "Cycle " `it' "; Estimation of fixed effects"
      }
      local jt = 0
      local cont = 1
      while  `cont'!=0 {
        local jt = `jt'+1
        replace `mu' = exp(`lp')
        if "`rden'"!="" {
          replace `mu' = `mu' * `rden'
        }
        replace `wy' = `lp' + (`yv' - `mu')/`mu'
        if "`cluster'"=="" {
          /* Calculate deviance */
          gen double `dv' = -2*`yv'*log(`mu')
          replace `dv' = `dv' + 2*`yv'*log(`yv') if `yv'!=0
          if `fvar'>0 {
            replace `dv' = `dv' + 2*(`yv'+1/`fvar')* /*
            */  log((1+`fvar'*`mu')/(1+`fvar'*`yv'))
          }
          else {
            replace `dv' = `dv' + 2*(`mu' - `yv')
          }  
          replace `se' = `mu'
        }
        else {
          /* Calculate deviance */
          gen double `dv' = -2*`yv'*log(`mu')
          replace `dv' = `dv' + 2*`yv'*log(`yv') if `yv'!=0
          by `cluster': replace `se' = sum(`mu')
          replace `se' = 0 if !`lt'
          if `fvar'>0 {
            replace `dv' = `dv' + 2*(`so'+1/`fvar')* /*
            */  log((1+`fvar'*`se')/(1+`fvar'*`so')) if `lt'
          }
          else {
            replace `dv' = `dv' + 2*(`se' - `so') if `lt'
          }  
        }
        replace `dv' = sum(`dv')
        local dev = `dv'[_N]
        drop `dv'
        if "`log'"!="nolog" {
          noi di "   Iteration " `jt' "; deviance = " %8.3f `dev'
        }
        /* Convergence Test */
        if `jt'>1 {
          if (`dev' - `last')>0.005 {
            noi di in blue "Warning - deviance not reduced"
            local cont = 0
          }
          else if (`last'-`dev')<0.001 { 
            local cont = 0 
          }
        }
        local last=`dev'
        if `cont'!=0 {
          /* New regression estimate */
          if "`cluster'"=="" {
            replace `wt' = `mu'/(1+`fvar'*`mu')
            _rpois_2 `wy' `xvars' [iw = `wt'], dep(`yv')
            matrix beta = get(_b)
            drop `lp'
            matrix score double `lp' = beta
          }
          else {
            _rpois_2 `wy' `xvars' [iw = `mu'], /*
            */ cl(`cluster') ra(`fvar') dep(`yv')
            matrix beta = get(_b)
            drop `lp'
            matrix score double `lp' = beta
          }
        }
      }
      /* Frailty variance estimation */
      if "`method'"!="fi" {
        if "`log'"!="nolog" {
          noi di "Cycle " `it' "; Estimation of frailty variance"
        }
        local kt = 0
        local cont = 1
        while `cont'!=0 {
          local kt = `kt'+1
          local last = `fvar'
          if "`method'" == "ml" {
            _rpois_3 `so' `se', fv(`fvar')
            local fstep = ${S_1}/${S_2}          
            if `fvar'>`fstep' {
              local fvar = `fvar' - `fstep'
            } 
            else {
              local fvar = 0
            }
          }
          else {
            replace `wt'=1/(1+`fvar'*`se')^2
            if "`method'" == "ps" {
              replace `wy' = sum(`wt'*((`so'-`se')^2-`se'))
            }
            else { /* method == "pe" */
              replace `wy' = sum(`wt'*((`so'-`se')^2+`so'-2*`se'))
            }
            replace `wt' = `wt'*(`se'^2)
            replace `wt' = sum(`wt')
            local fvar = `wy'[_N]/`wt'[_N]
            if `fvar' < 0 { 
              local fvar = 0 
            }     
          }
          if "`log'"!="nolog" {
            noi di "   Iteration " `kt' "; frailty variance = " %6.3f `fvar' 
          }
          /* Convergence test */
          if (abs(`fvar'-`last')<0.001) {
             local cont = 0
          }
        }
        local couter = (`kt'!=1)
      }
      else {
        local couter 0
      }
    }
    /* Calculate log likelihood 
    if `fvar'>0 {
      replace `so' = lngamma(`so' + 1/`fvar') - lngamma(1/`fvar') + /*
                         */ `so'*log(`fvar')
      replace `so' = sum(`so')
      local logl = `so'[_N] - `dev'/2
    }
    */
  }
  /* Return results to calling program */
  global S_1 "`ngrp'" 
  global S_2 "`dev'"
  global S_3 "`fvar'"
  global S_4 "`mdf'"
  global S_5 "`logl'"
  capture
end

program define _rpois_2
  version 4.0 
  local varlist "req ex min(1)"
  local if "opt"
  local in "opt"
  local options "CLuster(string)  RAtio(real 0)  Dep(string) RObust" 
  local weight "iweight"
  parse "`*'"
  preserve
  tempvar k w sw sx cn wy lp touse
  quietly {
    if "`weight'"=="" {
      gen double `w' = 1
    }
    else {
      gen double `w' `exp'
    }
  mark `touse' `if' `in'
  markout `touse' `varlist' `cluster' `w' 
  quietly {
    keep if `touse'
    if "`cluster'"=="" {
      parse "`varlist'", parse(" ")
      local wy "`1'"
      mac shift
      local xvars "`*'"
      matrix accum A = `wy' `xvars' [iw=`w']
      matrix XWX = A[2...,2...]
      matrix XWY = A[1,2...]
      matrix XWXi = syminv(XWX)
      matrix beta = XWY*XWXi
      matrix score double `lp' = beta
      if "`robust'"=="" {
        matrix post beta XWXi, dep(`dep')
      }
      else {
        replace `w' = (`w'*(`wy' - `lp'))^2
        matrix accum A = `xvars' [iw = `w']
        matrix A = XWXi*A
        matrix A = A*XWXi
        matrix beta = get(_b)         /* Restore old coefficients */
        matrix post beta A, dep(`dep')
      }
    }
    else {
      by `cluster': gen double `sw' = sum(`w')
      gen double `k' = 1 - sqrt(1/(1+`ratio'*`sw'))
      by `cluster': gen double `cn' = 1 - `k'[_N]
      parse "`varlist'", parse(" ")
      by `cluster': gen double `sx' = sum(`w'*`1')
      by `cluster': gen double `wy' = `1' - `k'[_N]*`sx'[_N]/`sw'[_N]
      mac shift
      local xvars ""
      local xlist ""
      while "`1'"!="" {
        by `cluster': replace `sx' = sum(`w'*`1')
        tempvar xv
        by `cluster': gen double `xv' = `1' - `k'[_N]*`sx'[_N]/`sw'[_N]
        local xvars "`xvars' `1'"
        local xlist "`xlist' `xv'"
        mac shift
      }
      matrix accum A = `wy' `xlist' `cn' [iw=`w'], noconstant 
      matrix XWX = A[2...,2...]
      matrix XWY = A[1,2...]
      matrix XWXi = syminv(XWX)
      matrix beta = XWY*XWXi
      matrix colnames beta = `xvars' _cons
      matrix colnames XWXi = `xvars' _cons
      matrix rownames XWXi = `xvars' _cons
      matrix score double `lp' = beta
      if "`robust'"=="" {
        matrix post beta XWXi, dep(`dep')
      }
      else {
        by `cluster': replace `sx' = sum(`w'*`lp')
        by `cluster': replace `wy' = `wy' - `lp' + `k'[_N]*`sx'[_N]/`sw'[_N]
        replace `cn' = `cn'*`w'*`wy'
        by `cluster': replace `cn' = sum(`cn') 
        parse "`xlist'", parse(" ")
        while "`1'" != "" {
          replace `1' = `1'*`w'*`wy'
          by `cluster': replace `1' = sum(`1')
          mac shift
        }
        by `cluster': replace `k' = (_n==_N)
        matrix accum A = `xlist' `cn' if `k', noconstant
        matrix A = XWXi*A
        matrix A = A*XWXi
        matrix beta = get(_b)  /* Restore old coefficients */
        matrix post beta A, dep(`dep')
      }
    }
  }
end

/*
  Calculate derivatives of log likelihood 
  with respect to frailty variance gamma. Arguments are observed 
  and expected frequencies
*/

program define _rpois_3
  version 4.0
  local varlist "req ex min(2) max(2)"
  local options "FVar(real 0) "
  parse "`*'"
  parse "`varlist'", parse(" ")
  tempvar ll d1 d2 td te used
  quietly {
    local D "`1'"
    local E "`2'"
    gen byte `used' = (`E'>0)
    if `fvar'==0 {
      gen double `d1' = - sum((`D'-`E')^2 - `D') 
      gen double `d2' = sum(`D')
      global S_1 = `d1'[_N]
      global S_2 = `d2'[_N]
    }
    else {
      local tau = 1/`fvar'
      local ltau = log(`tau')
      _lng double `ll' double `d1' double `d2', d(`D') t(`tau')
      gen double `td' = `tau' + `D'
      gen double `te' = `tau' + `E'
      replace `ll' = `ll' + `tau'*`ltau' - `td'*log(`te') if `used'
      replace `d1' = `d1' + 1 + `ltau' - log(`te') - `td'/`te' if `used'
      replace `d2' = `d2' + 1/`tau' + (`td' - 2*`te')/(`te'^2) if `used'
      replace `ll' = sum(`ll')
      replace `d1' = sum(`d1')
      replace `d2' = sum(`d2')
      local D1 = `d1'[_N]
      local D2 = `d2'[_N]
      global S_1 = -`D1'*(`tau'^2)
      global S_2 = (`D2'*`tau' + 2*`D1')*(`tau'^3)
    }
  }
  end

/*
   Log gamma function (logG) and its derivatives

   [logG(d + t) - logG(t)] and first two derivatives with respect to t

   approx() gives the threshold on d above which approximation cuts in
*/

program define _lng 
   version 5.0
   local varlist "req new min(3) max(3)"
   local options "d(string) t(real 0) approx(int 10)"
   parse "`*'"
   confirm var `d'
   parse "`varlist'", parse(" ")
   local lga `1'
   local dig `2'
   local tri `3'
   qui {
      replace `lga' = 0
      replace `dig' = 0
      replace `tri' = 0
   }
   tempvar Z update
   gen byte `update' = (`d' > `approx')
   gen double `Z' = `d' + `t'
   qui count if `update'
   local ldc = _result(1)
   * If large values of Z, reduce to managable values using approximation
   if `ldc'>0 {
      local z =   `approx'+`t'
      local lgaa = (`z'-0.5)*log(`z') - `z' + (1 - 1/(30*`z'))/(12*`z')
      local diga = log(`z') - (1+(1-1/(10*`z'^2))/(6*`z'))/(2*`z')
      local tria = (1 + (1 + (1 - 1/(5*`z'^2))/(3*`z'))/(2*`z'))/`z'
      qui {
         #delim ;
         replace `lga'= (`Z'-0.5)*log(`Z') - `Z' + (1 - 1/(30*`Z'))/(12*`Z')
            - `lgaa' if `update';
         replace `dig'= log(`Z') - (1+(1-1/(10*`Z'^2))/(6*`Z'))/(2*`Z')
            - `diga' if `update';
         replace `tri'= (1 + (1 + (1 - 1/(5*`Z'^2))/(3*`Z'))/(2*`Z'))/`Z'
            - `tria' if `update';
         replace `Z' = `z' if `update';
         #delim cr 
      }
   }
   local anyleft = 1
   local t1 = `t' + 0.001
   while `anyleft' {
      qui {
         replace `update' = `Z'>`t1'
         replace `Z' = `Z' - 1 if `update'
         replace `lga' = `lga' + log(`Z') if `update'
         replace `dig' = `dig' + 1/`Z' if `update'
         replace `tri' = `tri' - 1/`Z'^2 if `update'
         count if `update'
         local anyleft = _result(1)
      }
   }
   end