Statalist The Stata Listserver


[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]

Subject: st: How can I use this V9 module (hnblogit) with Version 8.0?


From   [email protected]
To   [email protected]
Subject   Subject: st: How can I use this V9 module (hnblogit) with Version 8.0?
Date   Tue, 16 Jan 2007 13:26:27 EST

Peter et al:

I wrote these programs for the second edition of  "Generalized Linear Models 
and Extensions" (Stata Press) which is currently at  the printer, and for my 
Cambridge Univ Press book entitled "Negative Binomial  Regression", for which I 
am currently going over the XML typescript pages for  indexing. I have made 
changes to the Poisson and negative binomial hurdle models  that are posted on 
the SSC site, --- the enhanced versions should be available  at each text's 
web site. I'll also repost them to SSC, but I'm swamped now  trying to finish up 
the indexing. Foremost changes relate to the differential  exponentiation of 
estimation tables and in calculating the correct post  estimation predict 
command. I'll also likewise enhance the geometric hurdle  models.  I am pasting 
the new version into this posting that supercedes the  older code found at the 
SSC site.

I wish I could help Peter out on his  request, but I don't have the time now. 
It should not be difficult to convert  the code he has pasted in below for my 
earlier version into Stata version 8  code. As it is, I used Stata's version 
9 ML capabilities that allow passthrough  features, thereby necessitating 
substantially less code than otherwise needed.  Perhaps someone can help him with 
his request. I suggest using the enhanced  version I have posted here rather 
than the version Peter pasted into his request  of yesterday. 

The NEW version is below my name. The NB-logit hurdle  predict command, 
hnblogit_p, is not shown here. I believe that the LL function,  jhnb_logit_ll, is 
the same as before. 

Joe Hilbe

=======================================================================
*!  Version 1.1
*  NEGATIVE BINOMIAL-LOGIT HURDLE REGRESSION :   
*  Joseph  Hilbe  :  7Oct2005  eform:9Jan2006 alpha  & hnblogit_p:26Sep2006

program hnblogit, eclass  byable(recall)  sortpreserve  ///
properties(swml svyr  svyj  svyb or rrr irr eform)

version  9.1
syntax [varlist]  [if] [in] [fw aw pw iw] [,   ///
CENsor(string)  Level(cilevel)    ///
OFFset(passthru)  EXposure(passthru)  ///
CLuster(passthru)  Robust    ///
EForm IRr OR ORIrr     ///
noLOG FROM(string asis) *  ]    

gettoken  lhs rhs : varlist

// Syntax  check
local expopts "`eform' `irr' `or' `orirr'"
if `:word  count `expopts'' > 1 {
di as error "only one of {opt eform},  {opt irr}, " _c
di as error "{opt or}, or {opt orirr} may be  specified"
exit 198
}

mlopts  mlopts, `options'

if ("`weight'" != "")   local  weight  "[`weight'   `exp']"
if (`"`from'"' !=  `""')  local initopt `"init(`from')"'

ml model lf  jhnb_logit_ll      ///
(logit: `lhs' =   `rhs',  `offset' `exposure')    ///
(negbinomial:  `lhs' =  `rhs', `offset' `exposure')     ///
/lnalpha       ///
`if' `in'   `weight',      ///
`mlopts' `robust'   `cluster'    ///
title("Negative Binomial-Logit   Hurdle Regression") ///
maximize `log'   `initopt'

ereturn scalar k_aux = 1
ereturn scalar   k_eform = 2
if "`orirr'" == "" {  
ml display,  level(`level') `eform' `irr' `or' plus
}
else  {    
ml display, level(`level') or plus  neq(1)

// Eq 2 
// Do the  header
di as text _col(14) "{c |}" _col(23) "IRR    Std. Err." _c
di as text _col(44) "z    P>|z|"  _c
local levstr "[`=strsubdp("`level'")'% Conf.  Interval]"
local len = length("`levstr'")
di as  text _col(`=79-`len'') "`levstr'"
di as text "{hline 13}{c  +}{hline 64}"
di as result "negbinomial" _col(14) as text "{c  |}"

// Coefficients
tempname  beta
matrix `beta' = e(b)
// Get 2nd eq  matrix
matrix `beta' = `beta'[1,"negbinomial:"]
local  bnames `:colnames `beta''
local cv =  invnorm((100+`level')/200)
local i = 1
foreach b of local  bnames {
if "`b'" == "_cons" {
continue //  No IRR for constant
}

local  bstr = abbrev("`b'", 12)

local bval = `beta'[1,  `i']
local seb = [negbinomial]_se[`b']
local z =  `bval'/`seb' 
local pval = 2*normal(-abs(`z'))
local lb = `bval' - `cv'*`seb'
local ub = `bval' +  `cv'*`seb'

di as text %12s "`bstr'" _col(14) "{c |}"  _c
di as result _col(17) %9.0g exp(`bval') _c
di  as result _col(28) %9.0g `seb'*exp(`bval') _c
di as result  _col(38) %8.2f `z' _c
di as result _col(49) %5.3f `pval'  _c
di as result _col(58) %9.0g exp(`lb') _c
di  as result _col(70) %9.0g exp(`ub')

local i =  `i' + 1

}
di as text "{hline 13}{c  +}{hline 64}"
_diparm lnalpha, level(`level')
}
_diparm lnalpha, level(`level') exp label("alpha")
di as  text "{hline 13}{c BT}{hline 64}"

qui {
*  AIC
tempname aic
local nobs e(N)
local npred  e(df_m)
local df = e(N) - e(df_m) -1
local llike    e(ll)
scalar `aic' = ((-2*`llike') +    2*(`npred'+1))/`nobs'
}

di in gr _col(1) "AIC Statistic =  " in ye %9.3f `aic'

ereturn local cmd "hnblogit"
ereturn  local predict  "hnblogit_p"

end
======================================================================


000000000000000000000000000000000000000000000000000000000000000000000000000000
000000
Dear  experienced Stata users,

After a couple of years of putting it off, I  have recently started using
Stata. Joseph Hilbe has writted a module  "hnblogit" which would be very
useful for my PhD research. It is listed as a  Version 9 module.
Unfortunately, the Australian Bureau of Statistics holds  data which
currently I can only access remotely using Stata V8.0. I would be  very
grateful if somebody could tell me if this module would work in V8.0 if  I
simply changed the "Version 9.1" line to "Version 8.0". If not, is the  code
easily modifiable to make it work?

Thanks in  advance,
Peter.

Peter Siminski
PhD Student
Social Policy  Research Centre
University of New South Wales,  Australia
[email protected]

P.S. Here are the two  ado-files which make up the hnblogit module (sourced
from  http://econpapers.repec.org/software/bocbocode/s456401.htm ):

*! Version  1.0.0
* NEGATIVE BINOMIAL-LOGIT HURDLE REGRESSION :  Joseph Hilbe   :  7Oct2005
program hnblogit, eclass properties(svyb svyj   svyr)
version 9.1
syntax [varlist]  [if] [in] [fweight pweight aweight  iweight] [,  ///
CENsor(string)   Level(cilevel)                 ///
OFFset(passthru)   EXposure(passthru)             ///
CLuster(passthru) IRr Robust  noLOG FROM(string asis) *  ]
gettoken  lhs rhs :  varlist

mlopts mlopts,  `options'

if ("`weight'" != "")   local  weight "[`weight'   `exp']"
if (`"`from'"' !=  `""') local initopt `"init(`from')"'

ml model lf  jhnb_logit_ll (logit: `lhs' =  `rhs', `offset'  `exposure')
///
(negbinomial: `lhs' = `rhs', `offset' `exposure')    /lnalpha
///
`if' `in'  `weight',                      ///
`mlopts' `robust'  `cluster'              ///
title("Negative Binomial-Logit  Hurdle Regression")            ///
maximize `log'  `initopt'                 ///


ereturn scalar k_aux = 1

ml display, level(`level') `irr'

qui {
* AIC
tempvar aic
local nobs e(N)
local  npred e(df_m)
local df = e(N) - e(df_m)  -1
local llike  e(ll)
gen `aic'  = ((-2*`llike') +  2*(`npred'+1))/`nobs'


}
di in gr _col(1)  "AIC Statistic = " in ye %9.3f `aic'
end


*! version 1.0.0   30Sep2005
* Negative Binomial-Logit Hurdle: log likelihood function :Joseph  Hilbe
program define jhnb_logit_ll
version 9
args lnf beta1 beta2  alpha
tempvar pi mu a
qui gen double `a'  =  exp(`alpha')
qui gen double `pi' = exp(`beta1')
qui gen  double `mu' = exp(`beta2') * `a'
qui replace `lnf' = cond($ML_y1==0,  ln(1/(1+`pi')), ln(`pi'/(1+`pi')) +
///
$ML_y1 *  ln(`mu'/(1+`mu'))   -   ///
ln(1+`mu')/`a' +  lngamma($ML_y1 +  1/`a') -  ///
lngamma($ML_y1 + 1) -   lngamma(1/`a')   -     ///
ln(1-(1+`mu')^(-1/`a') ) )
end
 
*
*   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–2024 StataCorp LLC   |   Terms of use   |   Privacy   |   Contact us   |   What's new   |   Site index