[R-sig-ME] Testing for significance of a single random-effect factor in lme4
Claus Wilke
cwilke at mail.utexas.edu
Fri May 1 22:22:35 CEST 2009
I'm wondering how to test for the significance of a single random-effect
factor in lme4. For example, in the famous Rail example, I might want to know
whether the identity of the rail makes any difference. The naive approach
would be:
> (a <- lmer( travel ~ (1|Rail) + 1, Rail ) )
> (b <- lmer( travel ~ 1, Rail ) )
> anova(a,b)
but this doesn't work because the second line causes an error:
> (b <- lmer( travel ~ 1, Rail ) )
Error in lmerFactorList(formula, fr, 0L, 0L) :
No random effects terms specified in formula
An alternative would be to use
> (b <- lm( travel ~ 1, Rail ) )
but then the anova function fails:
> anova(a,b)
Error in FUN(X[[1L]], ...) :
no slot of name "call" for this object of class "lm"
(This case was already discussed in an earlier email to this list:
https://stat.ethz.ch/pipermail/r-sig-mixed-models/2007q2/000197.html
)
So I'm wondering whether I can just add a dummy factor which is the same for
all data points and thus cannot explain any variance:
> Rail$dummy <- factor( rep( "X", nrow(Rail) ) ) # add a dummy factor
> (a <- lmer( travel ~ (1|Rail) + (1|dummy), Rail ) )
Linear mixed model fit by REML
Formula: travel ~ (1 | Rail) + (1 | dummy)
Data: Rail
AIC BIC logLik deviance REMLdev
130.2 133.7 -61.09 128.6 122.2
Random effects:
Groups Name Variance Std.Dev.
Rail (Intercept) 615.3111 24.8055
dummy (Intercept) 2.3951 1.5476
Residual 16.1667 4.0208
Number of obs: 18, groups: Rail, 6; dummy, 1
Fixed effects:
Estimate Std. Error t value
(Intercept) 66.50 10.29 6.464
> (b <- lmer( travel ~ (1|dummy), Rail ) )
Linear mixed model fit by REML
Formula: travel ~ (1 | dummy)
Data: Rail
AIC BIC logLik deviance REMLdev
164.7 167.4 -79.34 165.2 158.7
Random effects:
Groups Name Variance Std.Dev.
dummy (Intercept) 82.828 9.101
Residual 559.088 23.645
Number of obs: 18, groups: dummy, 1
Fixed effects:
Estimate Std. Error t value
(Intercept) 66.50 10.67 6.231
> anova(a,b)
Data: Rail
Models:
b: travel ~ (1 | dummy)
a: travel ~ (1 | Rail) + (1 | dummy)
Df AIC BIC logLik Chisq Chi Df Pr(>Chisq)
b 3 171.226 173.897 -82.613
a 4 136.648 140.209 -64.324 36.578 1 1.467e-09 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
It seems to me that most of the estimates (likelihood, residuals, estimates
for the Rail factor) are the same whether I use
lmer( travel ~ (1|Rail) + (1|dummy), Rail )
or
lmer( travel ~ (1|Rail) + 1, Rail )
but on the other hand AIC and BIC are somewhat different.
Thanks for your help,
Claus Wilke
--
Claus Wilke
Section of Integrative Biology
and Center for Computational Biology and Bioinformatics
University of Texas at Austin
1 University Station C0930
Austin, TX 78712
cwilke at mail.utexas.edu
512 471 6028
More information about the R-sig-mixed-models
mailing list