bayesm (version 3.0-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

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,…)

• Xnobs 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)

• nprint print the estimated time remaining for every nprint'th draw (def: 100)

• 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 Rossi, Allenby, McCulloch. http://www.perossi.org/home/bsm-1

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=Mcmc1)

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

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