View source: R/rnegbinrw_rcpp.r
rnegbinRw | R Documentation |
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.
rnegbinRw(Data, Prior, Mcmc)
Data |
list(y, X) |
Prior |
list(betabar, A, a, b) |
Mcmc |
list(R, keep, nprint, s_beta, s_alpha, beta0, alpha) |
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)
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 keep th 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)) |
A list containing:
betadraw |
|
alphadraw |
|
llike |
|
acceptrbeta |
acceptance rate of the beta draws |
acceptralpha |
acceptance rate of the alpha draws |
The NBD regression encompasses Poisson regression in the sense that as alpha goes to infinity the NBD distribution tends toward 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.
Sridhar Narayanan (Stanford GSB) and Peter Rossi (Anderson School, UCLA), perossichi@gmail.com.
For further discussion, see Bayesian Statistics and Marketing by Rossi, Allenby, McCulloch.
rhierNegbinRw
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)}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.