Thank you to everyone who responded to my post. This is my first use of the
group and it is encouraging that people are so supportive. I will take all
of your comments on board: I will look into the V8 syntax for ML, remove
properties(svyb svyj svyr), encourage the ABS administrators to update to
8.2 (If not V9!) and do all this in renamed ado files. Any aditional advice
(especially on the ML) would be appreciated. And if I can make it work I
will post the code to the list.
Thanks,
Peter.
-----Original Message-----
From: [email protected]
[mailto:[email protected]]On Behalf Of [email protected]
Sent: Wednesday, 17 January 2007 5:26 AM
To: [email protected]
Subject: Subject: st: How can I use this V9 module (hnblogit) with
Version 8.0?
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
======================================================================
0000000000000000000000000000000000000000000000000000000000000000000000000000
00
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/
*
* 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/
