R: Markov Chain Monte Carlo for the Hierarchical Poisson Linear...
MCMChpoisson
R Documentation
Markov Chain Monte Carlo for the Hierarchical Poisson Linear
Regression Model using the log link function
Description
MCMChpoisson generates a sample from the posterior
distribution of a Hierarchical Poisson Linear Regression Model using
the log link function and Algorithm 2 of Chib and Carlin
(1999). This model uses a multivariate Normal prior for the fixed
effects parameters, an Inverse-Wishart prior on the random effects
variance matrix, and an Inverse-Gamma prior on the variance modelling
over-dispersion. The user supplies data and priors, and a sample from
the posterior distribution is returned as an mcmc object, which can be
subsequently analyzed with functions provided in the coda package.
A two-sided linear formula of the form 'y~x1+...+xp'
describing the fixed-effects part of the model, with the response on
the left of a '~' operator and the p fixed terms, separated by '+'
operators, on the right. Response variable y must be 0 or 1
(Binomial process).
random
A one-sided formula of the form '~x1+...+xq' specifying
the model for the random effects part of the model, with the q
random terms, separated by '+' operators.
group
String indicating the name of the grouping variable
in data, defining the hierarchical structure of the model.
data
A data frame containing the variables in the model.
burnin
The number of burnin iterations for the sampler.
mcmc
The number of Gibbs iterations for the sampler. Total
number of Gibbs iterations is equal to
burnin+mcmc. burnin+mcmc must be divisible by 10 and
superior or equal to 100 so that the progress bar can be displayed.
thin
The thinning interval used in the simulation. The number of
mcmc iterations must be divisible by this value.
seed
The seed for the random number generator. If NA, the
Mersenne Twister generator is used with default seed 12345; if an
integer is passed it is used to seed the Mersenne twister.
verbose
A switch (0,1) which determines whether or not the
progress of the sampler is printed to the screen. Default is 1: a
progress bar is printed, indicating the step (in %) reached by the
Gibbs sampler.
beta.start
The starting values for the beta
vector. This can either be a scalar or a p-length vector. The
default value of NA will use the OLS beta estimate of
the corresponding Gaussian Linear Regression without random
effects. If this is a scalar, that value will serve as the starting
value mean for all of the betas.
sigma2.start
Scalar for the starting value of the residual
error variance. The default value of NA will use the OLS estimates
of the corresponding Gaussian Linear Regression without random
effects.
Vb.start
The starting value for variance matrix of the random
effects. This must be a square q-dimension matrix. Default value of NA
uses an identity matrix.
mubeta
The prior mean of beta. This can either be
a scalar or a p-length vector. If this takes a scalar value, then
that value will serve as the prior mean for all of the betas. The
default value of 0 will use a vector of zeros for an uninformative
prior.
Vbeta
The prior variance of beta. This can either
be a scalar or a square p-dimension matrix. If this takes a scalar value,
then that value times an identity matrix serves as the prior
variance of beta. Default value of 1.0E6 will use a diagonal matrix
with very large variance for an uninformative flat prior.
r
The shape parameter for the Inverse-Wishart prior on variance
matrix for the random effects. r must be superior or equal to
q. Set r=q for an uninformative prior. See the NOTE for more details
R
The scale matrix for the Inverse-Wishart prior on variance matrix
for the random effects. This must be a square q-dimension
matrix. Use plausible variance regarding random effects for the
diagonal of R. See the NOTE for more details
nu
The shape parameter for the Inverse-Gamma prior on the
residual error variance. Default value is nu=delta=0.001 for
uninformative prior.
delta
The rate (1/scale) parameter for the Inverse-Gamma prior on the
residual error variance. Default value is nu=delta=0.001 for
uninformative prior.
FixOD
A switch (0,1) which determines whether or not the
variance for over-dispersion (sigma2) should be fixed (1) or not
(0). Default is 0, parameter sigma2 is estimated. If FixOD=1, sigma2
is fixed to the value provided for sigma2.start.
...
further arguments to be passed
Details
MCMChpoisson simulates from the posterior distribution sample
using the blocked Gibbs sampler of Chib and Carlin (1999), Algorithm
2. The simulation is done in compiled C++ code to maximize
efficiency. Please consult the coda documentation for a comprehensive
list of functions that can be used to analyze the posterior sample.
The model takes the following form:
y_i ~ Poisson(lambda_i)
With latent variables phi(lambda),
phi being the log link function:
NOTE: We do not provide default parameters for the priors on
the precision matrix for the random effects. When fitting one of
these models, it is of utmost importance to choose a prior that
reflects your prior beliefs about the random effects. Using the
dwish and rwish functions might be useful in choosing
these values.
Value
mcmc
An mcmc object that contains the posterior sample. This
object can be summarized by functions provided by the coda
package. The posterior sample of the deviance D, with
D=-2log(prod_i P(y_i|lambda_i)), is also provided.
lambda.pred
Predictive posterior mean for the exponential of
the latent variables. The approximation of Diggle et
al. (2004) is used to marginalized with respect to over-dispersion terms:
Siddhartha Chib and Bradley P. Carlin. 1999. “On MCMC Sampling in
Hierarchical Longitudinal Models.” Statistics and Computing. 9:
17-26.
Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007.
Scythe Statistical Library 1.0.http://scythe.wustl.edu.
Andrew D. Martin and Kyle L. Saunders. 2002. “Bayesian Inference for
Political Science Panel Data.” Paper presented at the 2002 Annual Meeting
of the American Political Science Association.