bayesm (version 1.1-2)

rnegbinRw: 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

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

concept

• MCMC
• NBD regression
• Negative Binomial regression
• Poisson regression
• Metropolis algorithm
• bayes

Details

Model: $y$ $\sim$ $NBD(mean=lambda, over-dispersion=alpha)$. $lambda=exp(x'beta)$ Prior: $beta$ $\sim$ $N(betabar,A^{-1})$ $alpha$ $\sim$ $Gamma(a,b)$. note: prior mean of $alpha = a/b$, $variance = a/(b^2)$ list arguments contain:

• y

{ nobs vector of counts (0,1,2,...)} X{nobs x nvar matrix} betabar{ nvar x 1 prior mean (def: 0)} A{ nvar x nvar pds prior prec matrix (def: .01I)} a{ Gamma prior parm (def: .5)} b{ Gamma prior parm (def: .1)} R{ number of MCMC draws} keep{ MCMC thinning parm: keep every keepth draw (def: 1)} s_beta{ scaling for beta| alpha RW inc cov matrix (def: 2.93/sqrt(nvar)} s_alpha{ scaling for alpha | beta RW inc cov matrix (def: 2.93)}

References

For further discussion, see Bayesian Statistics and Marketing by Allenby, McCulloch, and Rossi. http://gsbwww.uchicago.edu/fac/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)
Data = simnegbindata
betabar = rep(0,nvar)
A = 0.01 * diag(nvar)
a = 0.5; b = 0.1
Prior = list(betabar=betabar, A=A, a=a, b=b)

keep =1
s_beta=2.93/sqrt(nvar); s_alpha=2.93
Mcmc = list(R=R, keep = keep, s_beta=s_beta, s_alpha=s_alpha)
out = rnegbinRw(Data, Prior, Mcmc)

cat("betadraws ",fill=TRUE)
mat=apply(out$betadraw,2,quantile,probs=c(.01,.05,.5,.95,.99))
mat=rbind(beta,mat); rownames(mat)[1]="beta"; print(mat)
cat("alphadraws ",fill=TRUE)
mat=apply(matrix(out$alphadraw),2,quantile,probs=c(.01,.05,.5,.95,.99))
mat=rbind(as.vector(alpha),mat); rownames(mat)[1]="alpha"; print(mat)