* MAXIMUM LIKELIHOOD NEGATIVE BINOMIAL REGRESSION : Joseph Hilbe 8Feb1998
* OPTIONS: tolerance,log,iteration,offset, robust, cluster, score
program define nbinreg
 version 5.0
 local options "Level(integer $S_level) IRr"
 if "`*'"=="" | substr("`1'",1,1)=="," {
   if "$S_E_cmd"=="nbinreg" {
     error 301
   }
   parse "`*'"
   if `level'<10 | `level'> 99 {
     di in red "level() must be between 10 and 99"
     exit 198
   }
 }
 else {
   local varlist "req ex"
   local options "`options' LTolerance(real -1) noLOg ITerate OFfset(string) SCore(string) SCORE2 Robust CLuster(string) *"
   local in "opt"
   local if "opt"
   local weight "fweight aweight"
   parse "`*'"
   parse "`varlist'",parse(" ")
   if "`log'"!="" { local log "quietly" }
   global S_mloff "`offset'"

 tempvar touse mysamp
 tempname b f V
 mark `touse' `if' `in'
 markout `touse' `varlist' `offset'

  if "`weight'" != "" {
        tempvar wvar
        gen double `wvar' `exp'
  }
  else  local wvar   1

  if ("`weight'"=="aweight")  {
        qui sum `wvar'
        qui replace `wvar' = `wvar'/_result(3)
  }

  if "`offset'" != "" {
        local ovar "`offset'"
  }
  else  local ovar   0
  global S_cen "`censor'"
* ESTIMATION OF LL0
 quietly {
   qui poisson `1'
   mat coef=get(_b)
   eq `1': `1'
   eq lnalpha:
   ml begin
   ml function nbinlf		
   ml method lf
   ml model `b'= `1' lnalpha, from(coefs) depv(10)
   ml sample `mysamp' [`weight' `exp'] if `touse' 
   `log' ml maximize `f' `V', `options' 
   local lf0 = `f'
   drop `mysamp'
 } 

* ESTIMATION OF FULL MODEL
   qui poisson `varlist'
   local pll = $S_E_ll
   mat coef=get(_b)
   eq `1': `varlist'
   eq lnalpha:
   ml begin
   ml function nbinlf
   ml method lf
   ml model `b' = `1' lnalpha, from(coefs) depv(10)
   ml sample `mysamp' [`weight' `exp'] if `touse' 
   `log' ml maximize `f' `V', `options' 
   ml post nbin, title("Negative Binomial Estimates") lf0(`lf0') pr2
   local nbll = $S_E_ll
   global S_1 = 2*(`nbll' - `pll')

* ROBUST CALCULATIONS
   local y "`1'"
   mac shift
   global rhs `*' 
   tempname b V
   tempvar xb e1 e2 muvar alpha amuvar digam1 digam2

   mat `b' = get(_b)
   mat `V' = get(VCE)
   predict double `xb', index

* SCORE VARIABLES
   qui gen `muvar' = exp(`xb' + `ovar')
   qui gen `alpha' = exp([lnalpha]_b[_con])
   qui gen `amuvar' = `muvar'*`alpha'
   qui gen `digam1' = ( lngamma( (`y'+1/`alpha') +.0001)-lngamma((`y'+1/`alpha')-.0001) )/.0002
   qui gen `digam2' = ( lngamma(1/`alpha' +.0001)-lngamma(1/`alpha'-.0001) )/.0002

*SCORES
   qui gen double `e1' = (`y' - `muvar')/(1+ `amuvar') if `touse'
   qui gen double `e2' = (-1/`alpha')* ( (`alpha' * (`muvar'-`y'))/(1+`amuvar')) - ln(1+`amuvar') + `digam1' - `digam2' if `touse'

*  SCORE FUNCTIONS TO FILE
   if "`score'"!="" {
     gen `score' = `e1'
   }

   if "`score2'"!="" {
     gen score2 = `e2'
   }

* CLUSTER OPTIONS
   if "`robust'"!="" & "`cluster'"=="" {
      _robust `e1' `e2' if `touse'
      qui test $rhs, min
      local ch2 = _result(6)
   }
   if "`cluster'"!="" {
      _robust `e1' `e2' if `touse', cluster(`cluster')
      qui test $rhs, min
      local ch2 = _result(6)
   }
}
* OUTPUT
if "`irr'"!="" {
      local eform "eform(IRR)"
   }
ml mlout nbin, level(`level') `eform'
di in gr "   alpha    " in ye %9.0g exp([lnalpha]_b[_cons]) /*
   */ in gr "   [lnalpha]_cons = ln(alpha)"
di _col(21) in gr "(LR test against Poisson, chi2(1) = " in ye %9.0g /*
   */ $S_1 /*
   */ in gr " P = ", in ye %6.4f chiprob(1,$S_1) in gr ")"

end