Dear BUGS list members:
I wish to acquaint myself with fitting generalized linear models in
BUGS. I am therefore attempting to replicate a simple negative
binomial regression model that I have estimated via max-likelihood and
know well. In the model, the count variable v is expressed as
v[n] ~ negative-binomial (mean=mu.v[n], overdispersion=omega)
and I wish to estimate effects of covariates X on mu.v
mu.v[n] <- exp(beta*X[n]).
I also wish the overdispersion estimate omega---the point of all this
will be to embed omega in a hierarchical setting later on.
Setting X=constant+var1+var2, I get virtually identical
maximum-likelihood estimates in stata with
nbreg v var1 var2, exposure(constant)
and in R with
m1 <- glm.nb(v ~ var1+var2).
I checked Gelman+Hill (p. 336), some example from Ntzoufras, and the
WinBUGS FAQs (http://www.mrc-bsu.cam.ac.uk/bugs/faqs/contents.shtml#q5)
for reference and have three questions. My code appears below.
(1) I am pretty much confused about how BUGS (in my case, OpenBUGS)
parameterizes the dnegbin() distribution. Lines with comments ##2 and
##3 in my code are alternative specifications of mu.v from the
references I checked. Which, if any, is right?
(2) How is omega recovered? Line with comment ##5 is inspired from
Gelman+Hill. Does BUGS invert this element as it does with the normal
distribution (precision=1/variance)?
(3) The Bayesian model, as such, is un-identified, requiring
constraints of some sort over beta parameters. What are good
alternatives? Lines with comments ##6 and ##7 suggest that beta[1] is
determined by other betas, but I clearly have no clue here.
Suggestions and/or code on some comparable problem will be much
appreciated. I will append interesting/useful responses them and send
this back to the list.
Thank you very much in advance.
model {
for (n in 1:N){ ##1 loop over observations
# stochastic component
v[n] ~ dnegbin(A[n] , B);
A[n] <- B / (B + mu.v[n]); ##2 Ntzoufras version
A[n] <- B * mu.v[n]; ##3 As in Gelman+Hill
# link and linear predictor
log(v[n]) <- beta[1]
+beta[2]*var1[n] ##4 Exposure
+beta[3]*var2[n];
}
omega <- 1 + 1/B ##5 As in Gelman+Hill
###################################################
## CONSTRAINT ON BETA 1 FOR IDENTIFICATION HERE? ##
###################################################
##6 What should I use here???
############################
## NON-INFORMATIVE PRIORS ##
############################
for (k in 2:3){ ##7 k=1 omitted bec constrained?
beta[k] ~ dnorm(0, .001);
}
B ~ dgamma( 0.001, 0.001 )
}
------------
Eric Magar
ITAM, Mexico City
-------------------------------------------------------------------
This list is for discussion of modelling issues and the BUGS software.
For help with crashes and error messages, first mail [log in to unmask]
To mail the BUGS list, mail to [log in to unmask]
Before mailing, please check the archive at www.jiscmail.ac.uk/lists/bugs.html
Please do not mail attachments to the list.
To leave the BUGS list, send LEAVE BUGS to [log in to unmask]
If this fails, mail [log in to unmask], NOT the whole list
|