bayesm (version 3.1-6)

rhierNegbinRw: MCMC Algorithm for Hierarchical Negative Binomial Regression

Description

rhierNegbinRw implements an MCMC algorithm 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)

Value

A list containing:

llike

\(R/keep x 1\) vector of values of log-likelihood

betadraw

\(nreg x nvar x R/keep\) array of beta draws

alphadraw

\(R/keep x 1\) vector of alpha draws

acceptrbeta

acceptance rate of the beta draws

acceptralpha

acceptance rate of the alpha draws

Arguments

Data

list(regdata, Z)

Prior

list(Deltabar, Adelta, nu, V, a, b)

Mcmc

list(R, keep, nprint, s_beta, s_alpha, alpha, c, Vbeta0, Delta0)

Author

Sridhar Narayanan (Stanford GSB) and Peter Rossi (Anderson School, UCLA), .

Details

Model and Priors

\(y_i\) \(\sim\) NBD(mean=\(\lambda\), over-dispersion=alpha)
\(\lambda = exp(X_i\beta_i)\)

\(\beta_i\) \(\sim\) \(N(\Delta'z_i,Vbeta)\)

\(vec(\Delta|Vbeta)\) \(\sim\) \(N(vec(Deltabar), Vbeta (x) Adelta)\)
\(Vbeta\) \(\sim\) \(IW(nu, V)\)
\(alpha\) \(\sim\) \(Gamma(a, b)\) (unless Mcmc$alpha specified)
Note: prior mean of \(alpha = a/b\), variance \(= a/(b^2)\)

Argument Details

Data = list(regdata, Z) [Z optional]

regdata:	list of lists with data on each of nreg units
regdata[[i]]$X:	\(nobs_i x nvar\) matrix of \(X\) variables
regdata[[i]]$y:	\(nobs_i x 1\) vector of count responses
Z:	\(nreg x nz\) matrix of unit characteristics (def: vector of ones)

Prior = list(Deltabar, Adelta, nu, V, a, b) [optional]

Deltabar:	\(nz x nvar\) prior mean matrix (def: 0)
Adelta:	\(nz x nz\) PDS prior precision matrix (def: 0.01*I)
nu:	d.f. parameter for Inverted Wishart prior (def: nvar+3)
V:	location matrix of Inverted Wishart prior (def: nu*I)
a:	Gamma prior parameter (def: 0.5)
b:	Gamma prior parameter (def: 0.1)

Mcmc = list(R, keep, nprint, s_beta, s_alpha, alpha, c, Vbeta0, Delta0) [only R required]

R:	number of MCMC draws
keep:	MCMC thinning parameter -- keep every keepth draw (def: 1)
nprint:	print the estimated time remaining for every nprint'th draw (def: 100, set to 0 for no print)
s_beta:	scaling for beta | alpha RW inc cov (def: 2.93/sqrt(nvar))
s_alpha:	scaling for alpha | beta RW inc cov (def: 2.93)
alpha:	over-dispersion parameter (def: alpha ~ Gamma(a,b))
c:	fractional likelihood weighting parm (def: 2)
Vbeta0:	starting value for Vbeta (def: I)
Delta0:	starting value for Delta (def: 0)

References

For further discussion, see Chapter 5, Bayesian Statistics and Marketing by Rossi, Allenby, and McCulloch.

See Also

rnegbinRw

Examples

if(nchar(Sys.getenv("LONG_TEST")) != 0) {R=2000} else {R=10}
set.seed(66)

# Simulate from the Negative Binomial Regression
simnegbin = function(X, beta, alpha) {
  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(4,2), c(0.1,-1))
alpha = 5
Vbeta = rbind(c(2,1), c(1,2))

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

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