cplm (version 0.4-1)

bcpglmm: Bayesian Compound Poisson Generalized Linear Mixed Models

Description

This function fits Bayesian compound Poisson generalized linear mixed models using MCMC methods.

Usage

bcpglmm(formula, link = "log", data, inits = NULL, weights, 
    offset, subset, na.action, contrasts = NULL, n.chains = 3, 
    n.iter = 2000, n.burnin = floor(n.iter/2), 
    n.thin = max(1, floor(n.chains * (n.iter - n.burnin)/n.sims)), 
    n.sims = 1000, n.report = 1000, prior.beta.mean = NULL, 
    prior.beta.var = NULL, bound.phi = 100, bound.p = c(1.01, 1.99), 
    prior.Sigma = NULL,  tune.iter=4000, n.tune=10, tune.weight=0.25, ...)

Arguments

formula

a two-sided linear formula object describing the fixed-effects part of the model, with the response on the left of a ~ operator and the terms, separated by + operators, on the right. The vertical bar character "|" separates an expression for a model matri

link

a specification for the model link function. This can be either a literal character string or a numeric number. If it is a character string, it must be one of "log", "identity", "sqrt" or "inverse". If it is numeric, it is the same as the link.power

data

an optional data frame, list or environment (or object coercible by as.data.frame to a data frame) containing the variables in the model.

subset, weights, na.action, offset, contrasts

further model specification arguments as in glm; see there for details.

inits

a list of initial values to be used for each chain. It must be of length n.chains. For each element, it is a named list with five components 'beta' (fixed effects), 'phi' (dispersion), 'p' (index parameter), 'u' (random effects) and 'Sigma' (

n.chains, n.iter, n.burnin, n.thin, n.sims, n.report

parameters that control the number of chains, iterations, burnins, thinning, simulations to keep and report intervals. See bcpglm for details.

prior.beta.mean

a vector of prior means for the model coefficients. Default is a vector of zeros.

prior.beta.var

a vector of prior variance for the model coefficients. Default is a vector of 10000's.

bound.phi

upper bound of the uniform prior for the dispersion parameter phi. The default is 100.

bound.p

a vector of lower and upper bound for the index parameter $p$. The default is c(1.01,1.99).

tune.iter

number of iterations used for tuning the variance for the Normal proposal distribution. These iterations will not be used in the final output. Default is 4000, and set it to be zero if the tuning process is not desired.

n.tune

if positive, the tune.iter iterations will be divided into n.tune loops. Default is 10.

tune.weight

the weight to be given to the old covariance matrix as opposed to the covariance matrix estimated from simulated samples. Default is 0.25.

prior.Sigma

a list that specifies the prior for each of the variance components.

...

not used

Value

• bcpglmm returns an object of class "bcpglmm". See bcpglmm-class for details of the return values as well as various methods available for this class.

Details

This function implements the MCMC algorithm for the Bayesian compound Poisson generalized linear mixed model. It is similar to the bcpglm function, but only the direct density evaluation method is available currently. It employs the block-wise Metropolis-Hasting for the random effects ($u$) and direct simulation for the variance components ($\Sigma$) due to conjugacy. The defualt prior of the variance component for one group is either inverse-Gamma (igamma(scale=0.001, shape=0.001)) if there is one variance component parameter, or inverse-Wishart (iwish(df=3,scale=diag(1,df))) if there are more than one parameters. These can be changed in the prior.Sigma argument.

See Also

The users are recommended to see the documentation for bcpglmm-class, bcpglm, cpglmm, mcmc, and tweedie for related information.

Examples

fit1 <- bcpglmm(RLD~Zone*Stock+(1|Plant), data = fineroot,  
                   n.chains=3, n.iter=25000, n.burnin=5000,
                   n.sims=2000, n.report=5000)
gelman.diag(fit1$sims.list)                   
summary(fit1)

# now change the prior distribution for the random component
fit2 <- bcpglmm(RLD~Zone*Stock+(1|Plant), data = fineroot,  
                   n.chains=3, n.iter=25000, n.burnin=5000,
                   n.sims=2000, n.report=5000, tune.iter=6000, n.tune=20,                  
                   prior.Sigma = list(igamma(scale=0.01,shape=0.01)))
summary(fit2)