rnegbinRw {bayesm}R Documentation

MCMC Algorithm for Negative Binomial Regression

Description

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.

Usage

rnegbinRw(Data, Prior, Mcmc)

Arguments

Data list(y,X)
Prior list(betabar,A,a,b)
Mcmc list(R,keep,s_beta,s_alpha,beta0

Details

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:

Value

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

Note

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.

Author(s)

Sridhar Narayanam & Peter Rossi, Graduate School of Business, University of Chicago, Peter.Rossi@ChicagoGsb.edu.

References

For further discussion, see Bayesian Statistics and Marketing by Rossi, Allenby, McCulloch.
http://faculty.chicagogsb.edu/peter.rossi/research/bsm.html

See Also

rhierNegbinRw

Examples

##
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)
}


[Package bayesm version 2.2-2 Index]