[R-sig-ME] Testing for significance of a single random-effect factor in lme4

[Claus Wilke]

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