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.