*! version 1.1.1 : 6-1-94  (c) 1994 J. Hilbe
* GLM: Log negative binomial with k estimation (Chi2 based convergence)
* expected IM initialization, observed thereafter
* Prior analytic and frequency weights allowed; if/in, noconstant, nolog allowed

program define negbin
  version 3.1
  local varlist "req ex"
  #delimit ;
  local options "`options' LTolera(real -1) ITerate(integer 0)
     Offset(string) EForm LEvel(integer $S_level) Resid Expec(real 1)
     NOLog noCONStant";
  # delimit cr
  local weight "fweight aweight"
  local in "opt"
  local if "opt"
  parse "`*'"
  parse "`varlist'",parse(" ")
  qui  {
 tempvar w wt mu dev n z eta oldev ofs chi cdisp ocdisp alpha ofs llog touse
  local y "`1'"
  mac shift

  if `ltolera'<0  { local ltolera 1e-8  }
  if `iterate'<=0 { local iterate 50 }

  if `level'<10 | `level'>99  {
     local level 95
   }
  if "`constan'"!="" { local cons "nocons" }

  if "`eform'"~="" { local eform "eform(IRR)" }

  if "`offset'"!="" { gen `ofs'=`offset'}
  else { gen `ofs' = 0 }

  gen byte `touse'=1 `if' `in'
  replace `touse'=0 if `y'==. | `touse'==.
  if "`exp'"~="" {
     gen float `wt' `exp' if `touse'
     replace `wt'=. if `wt'<=0
     replace `touse'=0 if `wt'==.
     if ("`weight'"=="aweight") {
        sum `wt'
        replace `wt'=`wt'/_result(3)
     }
  }
  else gen float `wt'=`touse' if `touse'
  qui reg `y' `*' if `touse'
  mac list
  local nobs=_result(1)

* POISSON INITIALIZATION
  gen `llog' = 1  if `touse'
  local i = 1
  local ddev = 1
  gen `z' = 0      if `touse'
  gen `dev' = 0     if `touse'
  gen `oldev' = 1    if `touse'
  egen `mu' = mean(`y')   if `touse'
  replace `mu' = (`y'+ `mu')/2
  gen `eta'= log(`mu')    if `touse'
  while ((abs(`ddev') > `ltolera' & `i'<`iterate') | (`i'==0))  {
         replace `z' = (`eta'+(`y'-`mu')/`mu') - `ofs'  if `touse'
         regress `z' `*' [iw=`mu'*`wt'], mse1  dep(`y') `cons'
         cap drop `eta'
         predict `eta'  if `touse'
         replace `eta' = `eta' + `ofs'
         replace `mu' = exp(`eta')
         replace `oldev' = `dev'
         replace `dev' = (`y'*log(`y'/`mu') - (`y'-`mu')) if `y'>0
         replace `dev' = `mu' if `y'==0
         replace `dev' = sum(`dev'*`wt')
         replace `dev' = 2*`dev'[_N]
         local ddev = `dev' - `oldev'
         local i = `i'+1
     }
* END POISSON LOOP
    local pdev=`dev'                    /* Poisson deviance   */
    gen `chi'=(`y'-`mu')^2/`mu'  if `touse'
    replace `chi'=sum(`chi'*`wt')
    replace `chi'= `chi'[_N]            /* Poisson chi2       */
    local pred = _result(3)
    local df = `nobs'-`pred'-("`cons'"=="") /* degrees of freedom */
    if ("`weight'"=="fweight") {
       tempvar fwt
       egen `fwt'= sum(`wt')
       local fobs =`fwt'
       local df = `fobs'-`pred'-1
       if "`constan'"!="" {
         local df = `fobs'-`pred'
       }
    }
    cap drop `fwt'
    gen `cdisp' = `chi'/`df'  if `touse'  /* chi2 dispersion    */
    local pdisp = `cdisp'
    gen `alpha' = 1/`cdisp'  if `touse'   /* initial alpha(k)   */


* MAIN NEGATIVE BINOMIAL LOOP
  local j = 1
  local ddisp = 1
  gen `ocdisp' = 1  if `touse'
  replace `dev'=1
  gen `w' = 1    if `touse'
  while ((abs(`ddisp') > .0001 & `j'<`iterate') | (`j'==0))  {
      replace `ocdisp'= `cdisp'
      local i = 1
      local ddev = 1
      replace `z' = 0
      replace `dev' = 0
      replace `oldev' = 1
      drop `mu'
      egen `mu' = mean(`y')  if `touse'
      replace `mu' = (`y'+ `mu')/2
      replace `eta'= log(`mu')  if `touse'        /* NB log link */
      while ((abs(`ddev') > `ltolera' & `i'<`iterate') | (`i'==0))  {
          if `i'<=`expec' {
             replace `w' = `mu'
             replace `z'=(`eta'+(`y'-`mu')/`w') -`ofs'
             replace `w'=`w'^2/(`mu'+`alpha'*`mu'*`mu')
          }
          else {
            #delimit ;
            replace `w' =(`mu'/(1+`alpha'*`mu')) + (`y'-`mu')*
               ( (`alpha'*`mu')/(1+(2*`alpha'*`mu')+(`alpha'^2*`mu'^2)));
            #delimit cr
            replace `z'=(`eta'+(`y'-`mu')/( (1+`alpha'*`mu')*`w')) - `ofs'
          }
          regress `z' `*' [iw=`w'*`wt'] , mse1 dep(`y') `cons'
          cap drop `eta'
          predict `eta'
          replace `eta'=`eta'+`ofs'
          replace `mu'= exp(`eta')
          replace `oldev'=`dev'
          #delimit ;
          replace `dev' = `y'*log(`y'/`mu')-(1+`alpha'*`y')/`alpha' *
                          log((1+`alpha'*`y')/(1+`alpha'*`mu')) if `y'>0;
          replace `dev' = log(1+`alpha'*`mu')/`alpha' if `y'==0;
          #delimit cr
          replace `dev' = sum(`dev'*`wt')
          replace `dev' = 2*`dev'[_N]
          local ddev = `dev'-`oldev'
          local i = `i'+1
       }
    replace `chi'=(`y'-`mu')^2/(`mu'+`alpha'*`mu'*`mu')
    replace `chi'=sum(`chi'*`wt')
    replace `chi'= `chi'[_N]
    local pred = _result(3)
    local df = `nobs'-`pred'-("`cons'"=="") /* degrees of freedom */
    if ("`weight'"=="fweight") {
       tempvar fwt
       egen `fwt'= sum(`wt')
       local fobs =`fwt'
       local df = `fobs'-`pred'-1
       cap drop `fwt'
       if "`constan'"!="" {
          local df = `fobs'-`pred'
       }
    }
    replace `cdisp' = `chi'/`df'
    if "`nolog'"=="" {
      noi di `j' in gr ":    Deviance: " in ye %9.4g  `dev' /*
         */ in gr "  Alpha (k): " in ye %9.4g `alpha' /*
         */ in gr "  Disp:  " in ye %9.4g `cdisp'
    }
    replace `alpha' = `cdisp'* `alpha'        /*  alpha dampening  */
    local ddisp = `cdisp' - `ocdisp'
    local j=`j'+1
}
replace `llog' = `y'*log(`alpha'*`mu')-(`y'+1/`alpha')* /*
    */ log(1+`alpha'*`mu')+lngamma(`y'+1/`alpha')-lngamma(`y'+1) - /*
    */ lngamma(1/`alpha')
replace `llog'=sum(`llog'*`wt')
replace `llog'=`llog'[_N]

noi di in gr _n  _col(55) "No obs.     = " in ye %9.0g `nobs'
noi di in gr "Poisson Dp = " in ye %9.4g `pdisp' in gr  _col(55) /*
    */ "Poisson Dev = " in ye %9.4g `pdev'
noi di in gr "Alpha (k)  = " in ye %9.4g `alpha' in gr _col(55) /*
    */ "Neg Bin Dev = " in ye %9.4g `dev'
noi di in gr _col(55) "Prob > Chi2 = " in ye %9.4g chiprob(1,`pdev'-`dev') _n
noi di in gr "Chi2       = " in ye %9.0g `cdisp'*`df' in gr _col(55) /*
    */ "Deviance    = " in ye %9.0g `dev'
noi di in gr "Prob>chi2  = " in ye %9.0g chiprob(`df',`cdisp'*`df') /*
    */ in gr _col(55) "Prob>chi2   = " in ye %9.0g chiprob(`df',`dev')
noi di in gr "Dispersion = " in ye %9.0g `cdisp' in gr _col(55) /*
    */ "Dispersion  = " in ye %9.0g `dev'/`df' _n
scalar S_al = `alpha'
noi di in gr  "Log negative binomial regression"
qui regress,  noheader `eform' level(`level')

tempname b V
mat `b' = get(_b)
mat `V' = get(VCE)
mat post `b' `V', depname(`y') obs(`nobs')
noi mat mlout, level(`level') `eform'
noi di in gr "Loglikelihood = " `llog'



if "`resid'"!="" {
   tempvar Presid Dresid stdp nbv Hat Spear Sdev delta
   cap drop _mu _eta _Presid _Dresid _hat _Spear _Sdev
   gen `Presid' = (`y'-`mu')/sqrt(`mu'+`alpha'*`mu'* `mu') * sqrt(`wt') /*
      */ if `touse'
   gen `Dresid' =  sqrt(2*(`y'*log(`y'/`mu')-(1+`alpha'*`y')/`alpha' * /*
      */ log((1+`alpha'*`y')/(1+`alpha'*`mu')))) if `y'>0
   replace `Dresid' =  sqrt(2*(log(1+`alpha'*`mu')/`alpha')) if `y'==0
   replace `Dresid' = -`Dresid' if (`y'-`mu')<0
   replace `Dresid' = `Dresid' if `touse'
   predict `stdp',stdp
   gen `nbv' = `mu'+`alpha'*`mu'*`mu'  if `touse'
   gen `Hat' = (`stdp'^2*`nbv')  if `touse'
   gen `Spear' = (`y'-`mu')/sqrt(`nbv' * (1-`Hat'))  if `touse'
   gen `Sdev' =  sqrt(2*(`y'*log(`y'/`mu')-(1+`alpha'*`y')/`alpha' * /*
      */ log((1+`alpha'*`y')/(1+`alpha'*`mu')))/(1-`Hat')) if `y'>0
   replace `Sdev' = sqrt(2*(log(1+`alpha'*`mu')/`alpha')/(1-`Hat')) if `y'==0
   replace `Sdev' = -`Sdev' if (`y'-`mu')<0
   replace `Sdev' = `Sdev' if `touse'
   gen _mu = `mu'  if `touse'
   gen _eta = `eta'  if `touse'
   gen _Presid = `Presid'  if `touse'
   gen _Dresid = `Dresid'  if `touse'
   gen _hat    = `Hat'   if `touse'
   gen _Spear  = `Spear'  if `touse'
   gen _Sdev   = `Sdev'  if `touse'
   lab var _mu "Fit"
   lab var _eta "Linear Predictor"
   lab var _Presid "Pearson residual"
   lab var _Dresid "Deviance residual"
   lab var _hat "Hat diagonal"
   lab var _Spear "St. Pearson"
   lab var _Sdev "St. deviance"
   noi di in gr _n "Created variables: " in ye /*
      */ "_mu, _eta, _Presid, _Dresid, _hat, _Spear, _Sdev"
   }

}
end