rnegbinRw {bayesm} | R Documentation |
rnegbinRw
implements a Random Walk Metropolis Algorithm for the Negative
Binomial (NBD) regression model. beta | alpha and alpha | beta are drawn with two
different random walks.
rnegbinRw(Data, Prior, Mcmc)
Data |
list(y,X) |
Prior |
list(betabar,A,a,b) |
Mcmc |
list(R,keep,s_beta,s_alpha,beta0 |
Model: y ~ NBD(mean=lambda, over-dispersion=alpha).
lambda=exp(x'beta)
Prior: beta ~ N(betabar,A^{-1})
alpha ~ Gamma(a,b).
note: prior mean of alpha = a/b, variance = a/(b^2)
list arguments contain:
y
X
betabar
A
a
b
R
keep
s\_beta
s\_alpha
a list containing:
betadraw |
R/keep x nvar array of beta draws |
alphadraw |
R/keep vector of alpha draws |
llike |
R/keep vector of log-likelihood values evaluated at each draw |
acceptrbeta |
acceptance rate of the beta draws |
acceptralpha |
acceptance rate of the alpha draws |
The NBD regression encompasses Poisson regression in the sense that as alpha goes to
infinity the NBD distribution tends toward the Poisson.
For "small" values of alpha, the dependent variable can be extremely variable so that
a large number of observations may be required to obtain precise inferences.
Sridhar Narayanam & Peter Rossi, Graduate School of Business, University of Chicago, Peter.Rossi@ChicagoGsb.edu.
For further discussion, see Bayesian Statistics and Marketing
by Rossi, Allenby, McCulloch.
http://faculty.chicagogsb.edu/peter.rossi/research/bsm.html
## if(nchar(Sys.getenv("LONG_TEST")) != 0) {R=1000} else {R=10} set.seed(66) simnegbin = function(X, beta, alpha) { # Simulate from the Negative Binomial Regression lambda = exp(X %*% beta) y=NULL for (j in 1:length(lambda)) y = c(y,rnbinom(1,mu = lambda[j],size = alpha)) return(y) } nobs = 500 nvar=2 # Number of X variables alpha = 5 Vbeta = diag(nvar)*0.01 # Construct the regdata (containing X) simnegbindata = NULL beta = c(0.6,0.2) X = cbind(rep(1,nobs),rnorm(nobs,mean=2,sd=0.5)) simnegbindata = list(y=simnegbin(X,beta,alpha), X=X, beta=beta) Data1 = simnegbindata Mcmc1 = list(R=R) out = rnegbinRw(Data=Data1,Mcmc=Mcmc1) cat("Summary of alpha/beta draw",fill=TRUE) summary(out$alphadraw,tvalues=alpha) summary(out$betadraw,tvalues=beta) if(0){ ## plotting examples plot(out$betadraw) }