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

From: Henrik Parn <henrik.parn_at_bio.ntnu.no>
Date: Thu 31 Aug 2006 - 03:02:07 EST


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)
______________________________________________ R-help@stat.math.ethz.ch 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 Thu Aug 31 06:17:16 2006

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.1.8, at Thu 31 Aug 2006 - 22:28:34 EST.

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