rhierNegbinRw {bayesm} R Documentation

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

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

### 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) ### 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 ### 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 Narayanan (Stanford GSB) and Peter Rossi (Anderson School, UCLA), perossichi@gmail.com. ### 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))
}


[Package bayesm version 3.1-6 Index]