fastmcmc: Rapidly converging Markov chain Monte Carlo algorithm for...

In lmm: Linear Mixed Models

Rapidly converging Markov chain Monte Carlo algorithm for Bayesian inference in linear mixed models

Description

Simulates posterior draws of parameters in linear mixed models using the rapidly converging Markov chain Monte Carlo (MCMC) procedure described by Schafer (1998), which combines a Metropolis-Hastings algorithm with a modified Gibbs sampler.

Prior to the MCMC simulation, the posterior mode of the variance parameters is found using the algorithm of "fastmode". The results from a call to "fastmode" are returned along with the MCMC results.

For a description of the model and the prior distribution, see the "Details" section below.

Usage

fastmcmc(y, subj, pred, xcol, zcol, prior, seed, vmax,
   occ, start.mode, maxits=100, eps=0.0001, iter=1000, 
   start.mcmc, df=4)

Arguments

y	vector of responses. This is simply the individual yi vectors stacked upon one another. Each element of y represents the observed response for a particular subject-occasion, or for a particular unit within a cluster.
subj	vector of same length as y, giving the subject (or cluster) indicators i for the elements of y. For example, suppose that y is c(y1,y2,y3,y4) where length(y1)=2, length(y2)=3, length(y3)=2, and length(y4)=7. Then subj should be c(1,1,2,2,2,3,3,4,4,4,4,4,4,4).
pred	matrix of covariates used to predict y. The number of rows should be length(y). The first column will typically be constant (one), and the remaining columns correspond to other variables appearing in Xi and Zi.
xcol	vector of integers indicating which columns of pred will be used in Xi. That is, pred[,xcol] is the Xi matrices (stacked upon one another.
zcol	vector of integers indicating which columns of pred will be used in Zi. That is, pred[,zcol] is the Zi matrices (stacked upon one another).
prior	A list with four components specifying the hyperparameters of the prior distribution applied to sigma2 and psi. The components must be named "a", "b", "c", and "Dinv". All are scalars except for "Dinv", which is a matrix of dimension c(length(zcol),length(zcol)).
seed	Seed for random number generator. This should be a positive integer.
vmax	optional matrix of dimension c(max(occ),max(occ)) from which the Vi matrices will be extracted. In a longitudinal dataset, vmax would represent the Vi matrix for an individual with responses at all possible occasions 1,2,...,nmax=max(occ); for individuals with responses at only a subset of these occasions, the Vi will be obtained by extracting the rows and columns of vmax for those occasions. If no vmax is specified by the user, an identity matrix is used. In most applications of this model one will want to have Vi = identity, so most of the time this argument can be omitted.
occ	vector of same length as y indicating the "occasions" for the elements of y. This argument is relevant only if a non-identity vmax is specified. In a longitudinal dataset where each individual is measured on at most nmax distinct occasions, each element of y corresponds to one subject-occasion, and the elements of occ should be coded as 1,2,...,nmax to indicate these occasion labels. (You should label the occasions as 1,2,...,nmax even if they are not equally spaced in time; the actual times of measurement will be incorporated into the matrix "pred".)
start.mode	optional starting values of the parameters for the mode-finding procedure. If this argument is not given then the function chooses its own starting values. This argument should be a list of three elements named "beta", "psi", and "sigma2". Note that "beta" should be a vector of the same length as "xcol", "psi" should be a matrix of dimension c(length(zcol),length(zcol)), and "sigma2" should be a scalar.
maxits	maximum number of cycles of the mode-finding procedure. The algorithm runs to convergence or until "maxits" iterations, whichever comes first.
eps	convergence criterion for the mode-finding procedure. The algorithm is considered to have converged if the relative differences in all parameters from one iteration to the next are less than eps–that is, if all(abs(new-old)<eps*abs(old)).
iter	number of cycles of the MCMC procedure to be performed.
start.mcmc	optional starting values of the parameters for the MCMC procedure. If this argument is not given, then the procedure is started at the posterior mode.
df	degrees of freedom for the multivariate t approximation in the Metropolis-Hastings algorithm.

Details

The algorithm is described in Section 5 of Schafer (1998).

The model, which is typically applied to longitudinal or clustered responses, is

yi = Xi%*%beta + Zi%*%bi + ei , i=1,...,m,

where

yi = (ni x 1) response vector for subject or cluster i; Xi = (ni x p) matrix of covariates; Zi = (ni x q) matrix of covariates; beta = (p x 1) vector of coefficients common to the population (fixed effects); bi = (q x 1) vector of coefficients specific to subject or cluster i (random effects); and ei = (ni x 1) vector of residual errors.

The vector bi is assumed to be normally distributed with mean zero and unstructured covariance matrix psi,

bi \sim N(0,psi) independently for i=1,...,m.

The residual vector ei is assumed to be

ei \sim N(0,sigma2*Vi)

where Vi is a known (ni x ni) matrix. In most applications, Vi is the identity matrix.

The prior distribution applied to the within-unit residual variance is scaled inverted-chisquare,

sigma2 \sim a / chisq(b),

where chisq(b) denotes a chisquare random variable with b degrees of freedom, and a and b are user-defined hyperparameters. Values for the hyperparmeters may be chosen by regarding a/b as a rough prior guess for sigma2, and as the imaginary degrees of freedom on which this guess is based.

The prior distribution applied to the between-unit covariance matrix is inverted Wishart,

psiinv \sim W(c,D),

where psiinv is the inverse of the between-unit covariance matrix psi, and W(c,D) denotes a Wishart distribution with degrees of freedom c and scale matrix D. Values for the hyperparameters may be chosen by regarding Dinv/c (the inverse of D divided by c) as a rough prior guess for psi, and c as the imaginary degrees of freedom on which this guess is based.

An improper uniform prior density function is applied to the fixed effects beta.

Value

a list containing the following components.

beta	simulated value of coefficients beta after "iter" cycles of the MCMC algorithm. This is a vector of the same length as xcol.
sigma2	simulated value of the residual variance sigma2 after "iter" cycles of the MCMC algorithm.
psi	simulated value of the between-unit covariance matrix psi after "iter" cycles of the MCMC algorithm.
sigma2.series	vector of length "iter" containing the entire history of simulated values of sigma2. That is, sigma2.series[t] contains the value of sigma2 at cycle t.
psi.series	array of dimension c(length(zcol),length(zcol),iter) containing the entire history of simulated values of psi. That is, psi.series[,,t] contains the value of psi at cycle t.
ratios	vector of length "iter" containing the entire history of acceptance ratios from the Metropolis-Hastings algorithm. These ratios diagnose the quality of the multivariate t approximation. If the approximation were perfect, all of these ratios would be equal to one.
reject	logical vector of length "iter" indicating, for each cycle of the algorithm, whether the Metropolis-Hastings candidate was accepted (T) or rejected (F).
mode.list	a list containing the results of the mode-finding procedure. The contents of this list are identical to those produced by "fastmode". For more information, see the help file for "fastmode".

References

Schafer, J.L. (1998) Some improved procedures for linear mixed models. Submitted to Journal of the American Statistical Association.

See Also

ecmeml, ecmerml, fastml, fastrml, fastmode, mgibbs, example

Examples

## Not run: 
For a detailed example, see the file "example.R" distributed
with this library.


## End(Not run)

lmm documentation built on Aug. 20, 2023, 1:08 a.m.