[R] lme4: How to specify nested factors, meaning of : and %in%

[Claus Wilke]

Hello list,

I'm trying to figure out how exactly the specification of nested random 
effects works in the lmer function of lme4. To give a concrete example, 
consider the rat-liver dataset from the R book (rats.txt from: 
http://www.bio.ic.ac.uk/research/mjcraw/therbook/data/ ).

Crawley suggests to analyze this data in the following way:
	library(lme4)
	attach(rats)
	Treatment <- factor(Treatment)
	Rat <-factor(Rat)
	Liver<-factor(Liver)
	m1<-lmer(Glycogen~Treatment+(1|Treatment/Rat/Liver))
The problem that I have with this analysis is that Treatment gets both a fixed 
and a random effect [1]. So I would like to remove the fixed effect from 
Treatment, but I'm not sure how to do this. Is the following correct?
	m2<-lmer(Glycogen~Treatment+(1|Treatment:Rat)+(1|Treatment:Rat:Liver))

From playing around with the lmer function, it seems that
	(1|Treatment/Rat/Liver)
is the same as
	(1|Treatment)+(1|Treatment:Rat)+(1|Treatment:Rat:Liver)
The two formulas give exactly the same fit. If everything I have said so far 
is correct, then I'm wondering what the %in% operator does. Naively, I would 
think that the following is again the same as (1|Treatment/Rat/Liver):
	(1|Treatment)+(1|Rat %in% Treatment)+(1|Liver %in% (Rat %in% Treatment))
Yet fitting this formula to the data gives different estimates for the random 
effects than fitting (1|Treatment/Rat/Liver). Can anybody tell me what's 
going on?

Below follows the R output for the various fits. Thanks a lot for your help,
  Claus Wilke

[1] The more I think about this example, the more I belive that the formula 
should actually be Glycogen~Treatment+(1|Rat/Treatment/Liver).  However, for 
the sake of the argument, please assume that rats are nested within 
treatments, because that corresponds to the case I actually want to analyze.

> (m1<-lmer(Glycogen~Treatment+(1|Treatment/Rat/Liver)))
Linear mixed-effects model fit by REML
Formula: Glycogen ~ Treatment + (1 | Treatment/Rat/Liver)
   AIC   BIC logLik MLdeviance REMLdeviance
 231.6 241.1 -109.8      234.9        219.6
Random effects:
 Groups                Name        Variance Std.Dev.
 Liver:(Rat:Treatment) (Intercept) 14.1617  3.7632
 Rat:Treatment         (Intercept) 36.0843  6.0070
 Treatment             (Intercept)  4.7039  2.1689
 Residual                          21.1678  4.6008
number of obs: 36, groups: Liver:(Rat:Treatment), 18; Rat:Treatment, 6; 
Treatment, 3

Fixed effects:
            Estimate Std. Error t value
(Intercept)  140.500      5.184  27.104
Treatment2    10.500      7.331   1.432
Treatment3    -5.333      7.331  -0.728

Correlation of Fixed Effects:
           (Intr) Trtmn2
Treatment2 -0.707
Treatment3 -0.707  0.500
	
> (m1a<-lmer(Glycogen~Treatment+(1|Treatment)+(1|Treatment:Rat)+(1|
Treatment:Rat:Liver)))
Linear mixed-effects model fit by REML
Formula: Glycogen ~ Treatment + (1 | Treatment) + (1 | Treatment:Rat) +      
(1 | Treatment:Rat:Liver)
   AIC   BIC logLik MLdeviance REMLdeviance
 231.6 241.1 -109.8      234.9        219.6
Random effects:
 Groups              Name        Variance Std.Dev.
 Treatment:Rat:Liver (Intercept) 14.1617  3.7632
 Treatment:Rat       (Intercept) 36.0843  6.0070
 Treatment           (Intercept)  4.7039  2.1689
 Residual                        21.1678  4.6008
number of obs: 36, groups: Treatment:Rat:Liver, 18; Treatment:Rat, 6; 
Treatment, 3

Fixed effects:
            Estimate Std. Error t value
(Intercept)  140.500      5.184  27.104
Treatment2    10.500      7.331   1.432
Treatment3    -5.333      7.331  -0.728

Correlation of Fixed Effects:
           (Intr) Trtmn2
Treatment2 -0.707
Treatment3 -0.707  0.500

> (m2<-lmer(Glycogen~Treatment+(1|Treatment:Rat)+(1|Treatment:Rat:Liver)))
Linear mixed-effects model fit by REML
Formula: Glycogen ~ Treatment + (1 | Treatment:Rat) + (1 | 
Treatment:Rat:Liver)
   AIC   BIC logLik MLdeviance REMLdeviance
 229.6 237.5 -109.8      234.3        219.6
Random effects:
 Groups              Name        Variance Std.Dev.
 Treatment:Rat:Liver (Intercept) 14.162   3.7632
 Treatment:Rat       (Intercept) 36.084   6.0070
 Residual                        21.168   4.6008
number of obs: 36, groups: Treatment:Rat:Liver, 18; Treatment:Rat, 6

Fixed effects:
            Estimate Std. Error t value
(Intercept)  140.500      4.708  29.842
Treatment2    10.500      6.658   1.577
Treatment3    -5.333      6.658  -0.801

Correlation of Fixed Effects:
           (Intr) Trtmn2
Treatment2 -0.707
Treatment3 -0.707  0.500

> (m3<-lmer(Glycogen~Treatment+(1|Treatment)+(1|Rat %in% Treatment)+(1|
Liver %in% (Rat %in% Treatment))))
Linear mixed-effects model fit by REML
Formula: Glycogen ~ Treatment + (1 | Treatment) + (1 | Rat %in% Treatment) +      
(1 | Liver %in% (Rat %in% Treatment))
   AIC   BIC logLik MLdeviance REMLdeviance
 244.6 254.1 -116.3      247.2        232.6
Random effects:
 Groups                          Name        Variance Std.Dev.
 Treatment                       (Intercept) 11.9371  3.4550
 Rat %in% Treatment              (Intercept)  3.9790  1.9948
 Liver %in% (Rat %in% Treatment) (Intercept)  3.9790  1.9948
 Residual                                    53.7172  7.3292
number of obs: 36, groups: Treatment, 3; Rat %in% Treatment, 1; Liver %in% 
(Rat %in% Treatment), 1

Fixed effects:
            Estimate Std. Error t value
(Intercept)  140.500      4.937  28.460
Treatment2    10.500      5.729   1.833
Treatment3    -5.333      5.729  -0.931

Correlation of Fixed Effects:
           (Intr) Trtmn2
Treatment2 -0.580
Treatment3 -0.580  0.500

-- 
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