[R] lmer applied to a wellknown (?) example

[Henrik Parn]

Dear all,

During my pre-R era I tried (yes, tried) to understand mixed models by 
working through the 'rat example' in Sokal and Rohlfs Biometry (2000) 
3ed p 288-292. The same example was later used by Crawley (2002) in his 
Statistical Computing p 363-373 and I have seen the same data being used 
elsewhere in the litterature.

Because this example is so thoroughly described, I thought it would be 
very interesting to analyse it also using lmer and to see how the 
different approaches and outputs differs - from the more or less manual 
old-school (?) approach in Sokal, aov in Crawley, and to mixed models by 
lmer.

In the example, three treatments (Treatment) with two rats (Rat) each 
(i.e six unique rats in total). Three liver preparations (Liver) are 
taken from each rat (i.e 18 unique liver preparations), and two glycogen 
readings (Glycogen) are taken from each liver preparation (36 readings).

We want to test if treatments has affected the glycogen levels. The 
readings are nested in preparation and the preparations nested in rats.

The data can be found here (or on p. 289 in Sokal):
http://www.bio.ic.ac.uk/research/mjcraw/statcomp/data/rats.txt
//
I was hoping to use the rat example as some kind of reference on my way 
to understand mixed models and using lmer. However, first I wish someone 
could check my suggested models!

My suggestions:

attach(rats)
rats$Glycogen <- as.numeric(Glycogen)
rats$Treatment <- as.factor(Treatment)
rats$Rat <- as.factor(Rat)
rats$Liver <- as.factor(Liver)
str(rats)

model1 <- lmer(Glycogen ~ Treatment + (1|Liver) + (1|Rat), data=rats)
summary(model1)

Was that it?

I also tried to make the 'liver-in-rat' nesting explicit  (as suggested 
in 'Examples from...')
 
model2 <- lmer(Glycogen ~ Treatment + (1|Rat:Liver) + (1|Rat), data=rats)
summary(model2)

but then the random effects differs from model1.

Does the non-unique coding of rats and preparations in the original data 
set mess things up? Do I need to recode the ids to unique levels...

rats$rat2 <- as.factor(rep(1:6, each=6))
rats$liver2 <- as.factor(rep(1:18, each=2))
str(rats)

...and then:

model3 <- lmer(Glycogen ~ Treatment + (1|liver2) + (1|rat2), data=rats)
# or maybe
model3 <- lmer(Glycogen ~ Treatment + (1|rat2:liver2) + (1|rat2), data=rats)
 

Can anyone help me to get this right! Thanks in advance!

P.S.
Thanks to all contributors to lme/lmer topic on the list (yes, I have 
searched for 'lmer nested'...) and also the examples provided by John 
Fox' 'Linear mixed models, Appendix to An R and S-PLUS companion...' and 
Douglas Bates' 'Examples from Multilevel Software...' and R-news 5/1. 
Very helpful, but as usually I bet I missed something...Sorry.

Regards,

Henrik 

-- 
************************
Henrik Pärn
Department of Biology
NTNU
7491 Trondheim
Norway

+47 735 96282 (office)
+47 909 89 255 (mobile)
+47 735 96100 (fax)