[R] testing for significance in random-effect factors using lmer

[Douglas Bates]

On 7/12/05, Eduardo.Garcia at uv.es <Eduardo.Garcia at uv.es> wrote:
> Hi, I would like to know whether it is possible to obtain a value of
> significance for random effects when aplying the lme or related
> functions. The default output in R is just a variance and standard
> deviation measurement.
> 
> I feel it would be possible to obtain the significance of these random
> effects by comparing models with and without these effects. However,
> I'm not used to perform this in R and I would thank any easy guide or
> example.

It is possible to do a likelihood ratio test on two fitted lmer models
with different specifications of the random effects.  The p-value for
such a test is calculated using the chi-squared distribution from the
asymptotic theory which does not apply in most such comparisons
because the parameter for the null hypothesis is on the boundary of
the parameter region.  The p-value shown will be conservative (that
is, it is an upper bound on the true p-value).

For example

> library(mlmRev)
Loading required package: lme4
Loading required package: Matrix
Loading required package: lattice
> options(show.signif.stars = FALSE)
> (fm1 <- lmer(normexam ~ standLRT + sex + type + (1|school), Exam))
Linear mixed-effects model fit by REML
Formula: normexam ~ standLRT + sex + type + (1 | school) 
   Data: Exam 
      AIC      BIC    logLik MLdeviance REMLdeviance
 9357.384 9395.237 -4672.692   9325.485     9345.384
Random effects:
 Groups   Name        Variance Std.Dev.
 school   (Intercept) 0.084367 0.29046 
 Residual             0.562529 0.75002 
# of obs: 4059, groups: school, 65

Fixed effects:
               Estimate  Std. Error   DF t value  Pr(>|t|)
(Intercept) -1.7233e-03  5.4982e-02 4055 -0.0313   0.97500
standLRT     5.5983e-01  1.2448e-02 4055 44.9725 < 2.2e-16
sexM        -1.6596e-01  3.2812e-02 4055 -5.0579 4.426e-07
typeSngl     1.6546e-01  7.7428e-02 4055  2.1369   0.03266
> (fm2 <- lmer(normexam ~ standLRT + sex + type + (standLRT|school), Exam))
Linear mixed-effects model fit by REML
Formula: normexam ~ standLRT + sex + type + (standLRT | school) 
   Data: Exam 
      AIC      BIC    logLik MLdeviance REMLdeviance
 9316.573 9367.043 -4650.287    9281.17     9300.573
Random effects:
 Groups   Name        Variance Std.Dev. Corr  
 school   (Intercept) 0.082477 0.28719        
          standLRT    0.015081 0.12280  0.579 
 Residual             0.550289 0.74181        
# of obs: 4059, groups: school, 65

Fixed effects:
               Estimate  Std. Error   DF t value  Pr(>|t|)
(Intercept)   -0.020727    0.052548 4055 -0.3944   0.69327
standLRT       0.554101    0.020117 4055 27.5433 < 2.2e-16
sexM          -0.167971    0.032281 4055 -5.2034 2.054e-07
typeSngl       0.176390    0.069587 4055  2.5348   0.01129
> anova(fm2, fm1)
Data: Exam
Models:
fm1: normexam ~ standLRT + sex + type + (1 | school)
fm2: normexam ~ standLRT + sex + type + (standLRT | school)
    Df     AIC     BIC  logLik  Chisq Chi Df Pr(>Chisq)
fm1  6  9357.4  9395.2 -4672.7                         
fm2  8  9316.6  9367.0 -4650.3 44.811      2  1.859e-10

At present the anova method for lmer objects does not allow comparison
with models that have no fixed effects.  Writing that code is on my
ToDo list but not currently at the top.  It is possible to use anova
to compare models fit by lme with models fit by lm (with the same
caveat about the calculated p-value being conservative).

An interesting alternative approach is to use Metropolis-Hastings
sampling for a MCMC chain based on the fitted model and create HPD
intervals from such a sample.  I have a prototype function to do this
for generalized linear mixed models in versions 0.97-3 and later of
the Matrix  package (currently hidden in the namespace and not
documented but the interested user can look at Matrix:::glmmMCMC).  It
happens that I developed the generalized linear version of this before
developing a version for linear mixed models but the lmm version will
be forthcoming.


> 
> Thanks.
> --
> ********************************
> Eduardo Moisés García Roger
> 
> Institut Cavanilles de Biodiversitat i Biologia
> Evolutiva - ICBIBE.
> Tel. +34963543664
> Fax  +34963543670
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
>