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

From: Claus Wilke <cwilke_at_mail.utexas.edu>
Date: Fri, 04 Apr 2008 15:32:32 -0500


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

______________________________________________
R-help_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Received on Fri 04 Apr 2008 - 20:38:48 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Fri 04 Apr 2008 - 21:30:26 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive