[R] fixed and random effects in lme

[Douglas Bates]

Federico Calboli <f.calboli at ucl.ac.uk> writes:

> Hi All,
> 
> I would like to ask a question on fixed and random effecti in lme. I am
> fiddlying around Mick Crawley dataset "rats" :
> 
> http://www.bio.ic.ac.uk/research/mjcraw/statcomp/data/
>  
> The advantage is that most work is already done in Crawley's book (page 361
> onwards) so I can check what I am doing.
> 
> I am tryg to reproduce the nested analysis on page 368:
> 
> model<-aov(Glycogen~Treatment/Rat/Liver + Error(Treatment/Rat/Liver), rats)
> 
> using lme.
> The code: 
> 
> model1<- lme(Glycogen~Treatment, random = ~1|Rat/Liver, data=rats)
> VarCorr(model1)
> 
>             Variance     StdDev  
> Rat =       pdLogChol(1)         
> (Intercept) 20.6019981   4.538942
> Liver =     pdLogChol(1)         
> (Intercept)  0.0540623   0.232513
> Residual    42.4362241   6.514309
> 
> Does NOT give me the same variance componets I find in Crawley's book (page
> 371 onwards).
> The code:
> 
> model2<- lme(Glycogen~Treatment, random = ~1|Treatment/Rat/Liver, data=rats)
> VarCorr(model2)
> 
>  	Variance     StdDev  
> Treatment = pdLogChol(1)         
> (Intercept) 12.54061     3.541272
> Rat =       pdLogChol(1)         
> (Intercept) 36.07900     6.006580
> Liver =     pdLogChol(1)         
> (Intercept) 14.17434     3.764883
> Residual    21.16227     4.600247
> 
> 
> gets me very similar results (I would guess the differences are due to
> rounding and the fact I am using R 1.6.2 and Crawley used S+).
> 
> My problem is: as *Treatment* is a fixed factor, why should I put it in
> both the fixed term side and random terms side of my code to get the right
> numbers? I fail to get this at all. Any elucidation would be appreciated.

Rat is nested within Treatment and Liver is nested within Rat within
Treatment.  Although level 1 of Rat in Treatment 1 is unrelated to
level 1 of Rat in Treatment 2, lme won't recognize that when you
specify the random effects as ~ 1 | Rat/Liver.  There are two ways out
of this: you can specify the nesting as you did
   ~ 1|Treatment/Rat/Liver 
but, as you point out, that doesn't always make sense, or you can
define a new set of groups as shown below

> library(nlme)
Loading required package: nls 
Loading required package: lattice 
> Rats = read.table('/tmp/rats.txt', header = TRUE)
> Rats$Trat = getGroups(Rats, form = ~ 1|Treatment/Rat, level = 2)
> str(Rats)
`data.frame':	36 obs. of  5 variables:
 $ Glycogen : int  131 130 131 125 136 142 150 148 140 143 ...
 $ Treatment: int  1 1 1 1 1 1 1 1 1 1 ...
 $ Rat      : int  1 1 1 1 1 1 2 2 2 2 ...
 $ Liver    : int  1 1 2 2 3 3 1 1 2 2 ...
 $ Trat     : Factor w/ 6 levels "1/1","1/2","2/1",..: 1 1 1 1 1 1 2 2 2 2 ...
> fm = lme(Glycogen ~ Treatment, data = Rats, random=~1|Trat/Liver)
> VarCorr(fm)
            Variance     StdDev  
Trat =      pdLogChol(1)         
(Intercept) 82.71680     9.094878
Liver =     pdLogChol(1)         
(Intercept) 14.17466     3.764925
Residual    21.16529     4.600575

Is this result comparable to Crawley's? 

Note that it can be meaningful to have a factor appear as both a fixed
effect and a random effect if, as a random effect, it is nested within
a random effect.