Stata The Stata listserver
[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]

Re: st: Truncated negative binomial with endogenous stratification


From   Cmeisner@worldbank.org
To   statalist@hsphsun2.harvard.edu
Subject   Re: st: Truncated negative binomial with endogenous stratification
Date   Thu, 28 Apr 2005 13:11:53 -0400

Thank you for the quick response!  Yes, I would like to persue this further if
possible. If you want to switch to another channel, just let me know.  I include
the code below for the nbinlf.ado likelihood  function (I think it's right), and
the trnbin0.ado code.  I renamed the file with an 'end' suffix for endogeneity.

Likelihood function:

* 0_truncated negative binomial log likelihood function
* 17Feb1998 Joseph Hilbe   STB-47 sg102

program define nbin0lfend
 tempvar mu a
 qui {
     local lnf    "`1'"
     local I      "`2'"
     local alpha  "`3'"
     if "$S_mloff" != "" {
       tempvar Io
       qui gen double `Io' = `I' + $S_mloff
     }
     else  local Io "`I'"
     gen double `a' = exp(`alpha')
     gen double `mu' = exp(`Io') * `a'
     replace `lnf' = ln($S_mldepn) + ($S_mldepn - 1) * ln(`mu'/(1+`mu'))   - /*
     */  ln(1+`mu')/`a' + lngamma($S_mldepn + 1/`a') - /*
     */  lngamma($S_mldepn + 1) - lngamma(1/`a')  -    /*
     */  ln(1-(1+`mu')^(-1/`a'))
  }
end

Estimation routine:

* 0-TRUNCATED NEGATIVE BINOMIAL REGRESSION : Joseph Hilbe 17Feb1998 STB-47 sg102
* OPTIONS: tolerance,log,iteration,offset, robust, cluster, score
program define trnbin0end
 version 5.0
 local options "Level(integer $S_level) IRr"
 if "`*'"=="" | substr("`1'",1,1)=="," {
   if "$S_E_cmd"=="trnbin0end" {
     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"
   #delimit ;
   local options "`options' LTolerance(real -1) noLOg ITerate
   OFfset(string) SCore(string) SCORE2 Robust CLuster(string) *";
   #delimit cr
   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 nbin0lfend
   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 nbin0lfend
   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'
   #delimit ;
   ml post nbin0, title("0-Truncated Negative Binomial w/ End. Strat.
Estimates")
   lf0(`lf0') pr2;
   #delimit cr
   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'
   #delimit ;
   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;
   #delimit cr

*SCORES
   tempvar tr0 tra ialpha
   qui gen double `ialpha' = -1/`alpha' if `touse'
   #delimit ;
   qui gen double `tr0' = ( `muvar' * -(1+`amuvar')^(-(1+`alpha')/
   `alpha') ) / (-1+(1+`amuvar')^(`ialpha')) if `touse';
   qui gen double `tra' = - (1 +`amuvar')^ (`ialpha')  *
   ( 1/(`alpha'^2) * ln(1+`amuvar') - (`muvar'/ (`alpha'*
   (1+`amuvar')))) / ( 1-(1+`amuvar')^(`ialpha')) if `touse';
   qui gen double `e1' = (`y' - `muvar')/(1+ `amuvar') - `tr0'
   if `touse';
   qui gen double `e2' = `ialpha'* ( (`alpha' * (`muvar'-`y'))/
   (1+`amuvar')) - ln(1+`amuvar') + `digam1' - `digam2' -
   `tra' if `touse';
#delimit cr

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

* 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 nbin0, 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



                                                                                                                                                               
                      Jhilbe@aol.com                                                                                                                           
                      Sent by:                         To:       statalist@hsphsun2.harvard.edu                                                                
                      owner-statalist@hsphsun2.        cc:                                                                                                     
                      harvard.edu                      Subject:  Re: st: Truncated negative binomial with endogenous stratification                            
                                                                                                                                                               
                                                                                                                                                               
                      04/28/2005 12:57 PM                                                                                                                      
                      Please respond to                                                                                                                        
                      statalist                                                                                                                                
                                                                                                                                                               
                                                                                                                                                               
                                                                                                                                                               




I am interested in your problem. I  suspect that the log-likelihood formula
used in nbinlf.ado for estimating the  zero truncated negative binomial is OK.
It's not a difficult formula and it was  checked many times prior to
distribution. I am wondering if the problem you are  experiencing relates to the
older
code in which the program was written. I  believe it was version 5 or 6.  I
had to entirely re-write my censored  Poisson regression algorithm using version

8 code in order to add certain survey  capabilities to it, including
pweights. I renamed it cepois, to distinguish it  from the 1997 version,
cenpois, and
later posted it to the Boston SSC  site.  I'll be happy to help with your
project of modifying the code to  allow endogenous stratification. I suspect,
though, that it will entail a  re-write to version 8 code. Or we can see what
capabilities version 9 has to  make the project easier.

If you are interesting, contact me privately  and we'll go from there.

Joe  Hilbe
====================================
I attempted to modify Joseph  Hilbe's 0-truncated negative binomial routine
(nbinlf.ado likelihood  function) to accomodate for endogenous stratification
(e.g. on-site  sampling), however I am running into problems during the
maximization  process; probably inherent in the trnbin0.ado file which
performs
the  estimation.  Has anyone else attempted to do this?  Would anyone be
able to
help?

Cheers,
Craig  Meisner

-------------------------------------------------------------
Craig  Meisner
The World Bank
Development Economics Research  Group
Infrastructure and Environment
MSN MC2-205, 1818 H Street,  NW
Washington DC 20433
Tel: 202-473-6852, Fax: 202-522-3230
E-mail:  cmeisner@worldbank.org

*
*   For searches and help try:
*   http://www.stata.com/support/faqs/res/findit.html
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/



*
*   For searches and help try:
*   http://www.stata.com/support/faqs/res/findit.html
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/



© Copyright 1996–2014 StataCorp LP   |   Terms of use   |   Privacy   |   Contact us   |   What's new   |   Site index