rhierNegbinRw {bayesm}R Documentation

MCMC Algorithm for Negative Binomial Regression

Description

rhierNegbinRw implements an MCMC strategy for the hierarchical Negative Binomial (NBD) regression model. Metropolis steps for each unit level set of regression parameters are automatically tuned by optimization. Over-dispersion parameter (alpha) is common across units.

Usage

rhierNegbinRw(Data, Prior, Mcmc)

Arguments

Data list(regdata,Z)
Prior list(Deltabar,Adelta,nu,V,a,b)
Mcmc list(R,keep,s_beta,s_alpha,c,Vbeta0,Delta0)

Details

Model: y_i ~ NBD(mean=lambda, over-dispersion=alpha).
lambda=exp(X_ibeta_i)

Prior: beta_i ~ N(Delta'z_i,Vbeta).

vec(Delta|Vbeta) ~ N(vec(Deltabar),Vbeta (x) Adelta).
Vbeta ~ IW(nu,V).
alpha ~ Gamma(a,b).
note: prior mean of alpha = a/b, variance = a/(b^2)

list arguments contain:

Value

a list containing:

llike R/keep vector of values of log-likelihood
betadraw nreg x nvar x R/keep array of beta draws
alphadraw R/keep vector of alpha draws
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 to 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.

For ease of interpretation, we recommend demeaning Z variables.

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 and McCulloch, Chapter 5.
http://faculty.chicagogsb.edu/peter.rossi/research/bsm.html

See Also

rnegbinRw

Examples

##
if(nchar(Sys.getenv("LONG_TEST")) != 0) {R=2000} 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)
}

nreg = 100        # Number of cross sectional units
T = 50            # Number of observations per unit
nobs = nreg*T
nvar=2            # Number of X variables
nz=2              # Number of Z variables
              
# Construct the Z matrix
Z = cbind(rep(1,nreg),rnorm(nreg,mean=1,sd=0.125))

Delta = cbind(c(0.4,0.2), c(0.1,0.05))
alpha = 5
Vbeta = rbind(c(0.1,0),c(0,0.1))

# Construct the regdata (containing X)
simnegbindata = NULL
for (i in 1:nreg) {
    betai = as.vector(Z[i,]%*%Delta) + chol(Vbeta)%*%rnorm(nvar)
    X = cbind(rep(1,T),rnorm(T,mean=2,sd=0.25))
    simnegbindata[[i]] = list(y=simnegbin(X,betai,alpha), X=X,beta=betai)
}

Beta = NULL
for (i in 1:nreg) {Beta=rbind(Beta,matrix(simnegbindata[[i]]$beta,nrow=1))}
    
Data1 = list(regdata=simnegbindata, Z=Z)
Mcmc1 = list(R=R)

out = rhierNegbinRw(Data=Data1, Mcmc=Mcmc1)

cat("Summary of Delta draws",fill=TRUE)
summary(out$Deltadraw,tvalues=as.vector(Delta))
cat("Summary of Vbeta draws",fill=TRUE)
summary(out$Vbetadraw,tvalues=as.vector(Vbeta[upper.tri(Vbeta,diag=TRUE)]))
cat("Summary of alpha draws",fill=TRUE)
summary(out$alpha,tvalues=alpha)

if(0){
## plotting examples
plot(out$betadraw)
plot(out$alpha,tvalues=alpha)
plot(out$Deltadraw,tvalues=as.vector(Delta))
}

[Package bayesm version 2.1-2 Index]