Fitting Generalized Linear Mixed Models using MCML

Description

This function fits generalized linear mixed models (GLMMs) by approximating the likelihood with ordinary Monte Carlo, then maximizing the approximated likelihood.

Usage

glmm(fixed, random, varcomps.names, data, family.glmm, m, 
varcomps.equal, doPQL = TRUE,debug=FALSE, p1=1/3,p2=1/3, p3=1/3,
rmax=1000,iterlim=1000, par.init, zeta=5)

Arguments

Details

Let beta be a vector of fixed effects and let u be a vector of random effects. Let X and Z be design matrices for the fixed and random effects, respectively. The random effects are assumed to be normally distributed with mean 0 and variance matrix D, where D is diagonal with entries from the unknown vector nu. Letting g be the link function, g(mu)=X beta+ ZU. For example, if the response is Bernoulli, then the logit function is the link.

Models for glmm are specified symbolically. A typical fixed effects model has the form response ~ terms where response is the (numeric) response vector and terms is a series of terms which specifies a linear predictor for response. A terms specification of the form first + second indicates all the terms in first together with all the terms in second with duplicates removed.

A specification of the form first:second indicates the set of terms obtained by taking the interactions of all terms in first with all terms in second. The specification first*second indicates the cross of first and second. This is the same as first + second + first:second.

The terms in the formula will be re-ordered so that main effects come first, followed by the interactions, all second-order, all third-order and so on: to avoid this, pass a terms object as the formula.

The random effects for glmm are also specified symbolically. The random effects model specification is typically a list. Each element of the random list has the form response ~ 0 + term. The 0 centers the random effects at 0. If you want your random effects to have a nonzero mean, then include that term in the fixed effects. Each variance component must have its own formula in the list.

To set some variance components equal to one another, use the varcomps.equal argument. The argument varcomps.equal should be a vector whose length is equal to the length of the random effects list. The vector should contain positive integers, and the first element of the varcomps.equal should be 1. To set variance components equal to one another, assign the same integer to the corresponding elements of varcomps.equal. For example, to set the first and second variance components equal to each other, the first two elements of varcomps.equal should be 1. If varcomps.equal is omitted, then the variance components are assumed to be distinct.

Each distinct variance component should have a name. The length of varcomps.names should be equal to the number of distinct variance components. If varcomps.equal is omitted, then the length of varcomps.names should be equal to the length of random.

Monte Carlo likelihood approximation relies on an importance sampling distribution. Though infinitely many importance sampling distributions should yield the correct MCMLEs eventually, the importance sampling distribution used in this package was chosen to reduce the computation cost. When doPQL is TRUE, the importance sampling distribution relies on PQL estimates (as calculated in this package). When doPQL is FALSE, the random effect estimates in the distribution are taken to be 0, the fixed effect estimates are taken to be 0, and the variance component estimates are taken to be 1.

This package's importance sampling distribution is a mixture of three distributions: a t centered at 0 with scale matrix determined by the PQL estimates of the variance components and with zeta degrees of freedom, a normal distribution centered at the PQL estimates of the random effects and with a variance matrix containing the PQL estimates of the variance components, and a normal distribution centered at the PQL estimates of the random effects and with a variance matrix based on the Hessian of the penalized log likelihood. The first component is included to guarantee the gradient of the MCLA has a central limit theorem. The second component is included to mirror our best guess of the distribution of the random effects. The third component is included so that the numerator and the denominator are similar when calculating the MCLA value.

The Monte Carlo sample size m should be chosen as large as possible. You may want to run the model a couple times to begin to understand the variability inherent to Monte Carlo. There are no hard and fast rules for choosing m, and more research is needed on this area. For a general idea, I believe the BoothHobert model produces stable enough estimates at m=10^3 and the salamander model produces stable enough estimates at m=10^5, as long as doPQL is TRUE.

To see the summary of the model, use summary().

Value

glmm returns an object of class glmm is a list containing at least the following components:

The function summary (i.e., summary.glmm) can be used to obtain or print a summary of the results. The generic accessor function coef (i.e., coef.glmm) can be used to extract the coefficients.

Author(s)

Christina Knudson

References

Geyer, C. J. (1994). On the convergence of Monte Carlo maximum likelihood calculations. Journal of the Royal Statistical Society, Series B, 61:261–274.

Geyer, C. J. and Thompson, E. (1992). Constrained Monte Carlo maximum likelihood for dependent data. Journal of the Royal Statistical Society, Series B, 54:657–699.

Sung, Y. J. and Geyer, C. J. (2007). Monte Carlo likelihood inference for missing data models. Annals of Statistics, 35:990–1011.

Examples

#First, using the basic Booth and Hobert dataset
#to fit a glmm with a logistic link, one variance component,
#one fixed effect, and an intercept of 0. The Monte Carlo
#sample size is 100 to save time.
library(glmm)
data(BoothHobert)
set.seed(1234)
mod.mcml1<-glmm(y~0+x1,list(y~0+z1),varcomps.names=c("z1"),data=BoothHobert,
family.glmm=bernoulli.glmm,m=100,doPQL=TRUE)
mod.mcml1$beta
mod.mcml1$nu
summary(mod.mcml1)
coef(mod.mcml1)

#Next, a model setting two variance components equal.
data(Booth2)
set.seed(1234)
mod.mcml3<-glmm(y~0+x1,list(y~0+z1,~0+z2),varcomps.names=c("z"), 
varcomps.equal=c(1,1), data=Booth2,family.glmm=bernoulli.glmm,
m=100,doPQL=FALSE)
mod.mcml3$beta
mod.mcml3$nu
summary(mod.mcml3)

#Now, a model with crossed random effects. There are two distinct
#variance components. To get more accurate answers for this model,
#use a larger Monte Carlo sample size, such as m=10^4 or 10^5
#and doPQL=TRUE.
set.seed(1234)
data(salamander)
m<-10
sal<-glmm(Mate~0+Cross,random=list(~0+Female,~0+Male),varcomps.names=c("F","M"),
data=salamander,family.glmm=bernoulli.glmm,m=m,debug=TRUE,doPQL=FALSE)
summary(sal)

Results