How to test random effects in a multilevel model in R

Question

I have been reading a good book called Applied Longitudinal Data Analysis: Modeling Change and Event Occurrence by Judith Singer and John Willet. The book shows that by modeling in 2 levels, we can model the individual change in level 1 and in level 2 model for systematic interindividual differences in change.

The R codes for the examples only show how to use lme() to estimate the fixed and random effects. However, the text suggested that we should test the variance components to determine whether the random effects are significant or not.

For example, one of the codes does only the following:

library(nlme)

model.a <- lme(alcuse~ 1, alcohol1, random= ~1 |id)
summary(model.a)

Linear mixed-effects model fit by REML
 Data: alcohol1 
       AIC      BIC    logLik
  679.0049 689.5087 -336.5025

Random effects:
 Formula: ~1 | id
        (Intercept)  Residual
StdDev:   0.7570578 0.7494974

Fixed effects: alcuse ~ 1 
                Value  Std.Error  DF  t-value p-value
(Intercept) 0.9219549 0.09629638 164 9.574139       0

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-1.8892070 -0.3079143 -0.3029178  0.6110925  2.8562135 

Number of Observations: 246
Number of Groups: 82

But the text lists the following:

• fixed effect: 0.922*** (s = 0.096) -> available in the output
• within person variance: 0.562*** (s = 0.062) -> can be obtained from the output (random effect residual std. dev squared)
• between person variance: 0.564*** (s = 0.119)

My work involves a lot of analysis for longitudinal data so I really need to understand this idea. Your help is very much appreciated.

Answer 1

First point:

You need to be careful if you want to test whether the variance of a random effect is 0: The standard $\chi^2$-asymptotics for the LR that are used in anova() do not apply for LRTs or restricted likelihood ratio tests on variance components since the null hypothesis is on the edge of the parameter space. Using this wrong reference distribution will yield ridiculously conservative tests.

Sources: Self & Liang (1987), Crainiceanu & Ruppert (2003), Greven, Crainiceanu, Küchenhoff (2008).

For exact LR tests on variances in linear mixed models with uncorrelated random effects you can use my package RLRsim. For generalized linear mixed models or models with correlated random effects, I would strongly recommend a parametric bootstrap to approximate the correct p-value if the LRT or RLRT (the latter has more power, see sources above) are somewhere in the vicinity of the critical value for the (wrong) standard reference distribution. Code for lme4 at the end of the post.

Second point:

It may be dangerous to use standard confidence intervals for variance parameter estimates, since their distribution is usually very skewed and so using a symmetric interval is a little simplistic- Douglas Bates has some material on this here.

---

Parametric bootstrap code for lme4-models:

#m0 is the lmer model under the null hypothesis (i.e. the smaller model)
#mA is the lmer model under the alternative
bootstrapAnova <- function(mA, m0, B=1000){
     oneBootstrap <- function(m0, mA){
         d <- drop(simulate(m0))
         m2 <-refit(mA, newresp=d)
         m1 <-refit(m0, newresp=d)
         return(anova(m2,m1)$Chisq[2])
     }  
     nulldist <-  if(!require(multicore)){
         replicate(B, oneBootstrap(m0, mA))
     } else {
         unlist(mclapply(1:B, function(x) oneBootstrap(m0, mA)))
     }   
     ret <- anova(mA, m0)
     ret$"Pr(>Chisq)"[2] <- mean(ret$Chisq[2] < nulldist)
     names(ret)[7] <- "Pr_boot(>Chisq)"
     attr(ret, "heading") <- c(attr(ret, "heading")[1], 
          paste("Parametric bootstrap with", B,"samples."),
          attr(ret, "heading")[-1])
     attr(ret, "nulldist") <- nulldist
     return(ret)
}
#use like this (and increase B if you want reviewers to believe you):
(bRLRT <- bootstrapAnova(mA=<BIG MODEL>, m0=<SMALLER MODEL>))

Answer 2

(I would post this as a comment to the previous answer by 'chi', but don't see that that is possible here.)

Be very careful with intervals() and anova() from the nlme package, these have exactly the flaws that @fabians points out above -- they rely on the standard chi-squared distribution applied to a quadratic approximation of the shape of the likelihood profile for the variances. The newer lme4a package (on r-forge) allows the creation of likelihood profiles for the variance parameters, which takes care of the quadratic approximation part (although not the distributional part).

Also, it is fairly hotly debated whether dropping non-significant variance components is a good idea or not (don't have an immediate reference, but this has been discussed on r-sig-mixed-models).

Answer 3

The intervals() function should provide you with $100(1-\alpha)$ confidence intervals for the random effects in your model, see help(intervals.lme) for more information. You can also test if any of the variance components can be droped from the model by using anova() (which amounts to do an LRT between two nested models).

Answer 4

One of the best resources on multilevel analysis in R is John Fox's web appendix to the text "An R and S-PLUS Companion to Applied Regression". It provides a great overview of the method and means to calculate some of the more familiar measures from the R NLME output. The appendix is available on CRAN.

Here is the link:

http://cran.r-project.org/doc/contrib/Fox-Companion/appendix-mixed-models.pdf