bayesm (version 3.1-6)

rnegbinRw: MCMC Algorithm for Negative Binomial Regression

Description

rnegbinRw implements a Random Walk Metropolis Algorithm for the Negative Binomial (NBD) regression model where \(\beta|\alpha\) and \(\alpha|\beta\) are drawn with two different random walks.

Usage

rnegbinRw(Data, Prior, Mcmc)

Value

A list containing:

betadraw

\(R/keep x k\) matrix of beta draws

alphadraw

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

llike

\(R/keep x 1\) vector of log-likelihood values evaluated at each draw

acceptrbeta

acceptance rate of the beta draws

acceptralpha

acceptance rate of the alpha draws

Arguments

Data

list(y, X)

Prior

list(betabar, A, a, b)

Mcmc

list(R, keep, nprint, s_beta, s_alpha, beta0, alpha)

Author

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

Details

Model and Priors

\(y\) \(\sim\) \(NBD(mean=\lambda, over-dispersion=alpha)\)
\(\lambda = exp(x'\beta)\)

\(\beta\) \(\sim\) \(N(betabar, A^{-1})\)
\(alpha\) \(\sim\) \(Gamma(a, b)\) (unless Mcmc$alpha specified)
Note: prior mean of \(alpha = a/b\), \(variance = a/(b^2)\)

Argument Details

Data = list(y, X)

y:	\(n x 1\) vector of counts (\(0,1,2,\ldots\))
X:	\(n x k\) design matrix

Prior = list(betabar, A, a, b) [optional]

betabar:	\(k x 1\) prior mean (def: 0)
A:	\(k x k\) PDS prior precision matrix (def: 0.01*I)
a:	Gamma prior parameter (not used if Mcmc$alpha specified) (def: 0.5)
b:	Gamma prior parameter (not used if Mcmc$alpha specified) (def: 0.1)

Mcmc = list(R, keep, nprint, s_beta, s_alpha, beta0, alpha) [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 matrix (def: 2.93/sqrt(k))
s_alpha:	scaling for alpha | beta RW inc cov matrix (def: 2.93)
alpha:	over-dispersion parameter (def: alpha ~ Gamma(a,b))

References

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

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=list(R=R))

cat("Summary of alpha/beta draw", fill=TRUE)
summary(out$alphadraw, tvalues=alpha)
summary(out$betadraw, tvalues=beta)

## plotting examples
if(0){plot(out$betadraw)}