Claus Wilke
2008-Apr-04 20:32 UTC
[R] lme4: How to specify nested factors, meaning of : and %in%
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
Claus Wilke
2008-Apr-04 20:45 UTC
[R] lme4: How to specify nested factors, meaning of : and %in%
A quick correction to my earlier email:> 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?I meant: "I would like to remove the random effect from Treatment". 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