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

From: Douglas Bates <dmbates_at_gmail.com>
Date: Wed 13 Jul 2005 - 23:05:57 EST

On 7/12/05, Eduardo.Garcia@uv.es <Eduardo.Garcia@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)

> 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@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help