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

Christoph Buser buser at stat.math.ethz.ch
Thu Aug 31 12:28:05 CEST 2006


Dear Henrik

There is an article in the R-News "Fitting linear mixed models
in R" in which you can find some examples for the syntax of
nested and non-nested design.

http://cran.r-project.org/doc/Rnews/Rnews_2005-1.pdf

Hope this helps

Christoph 

--------------------------------------------------------------
Christoph Buser <buser at stat.math.ethz.ch>
Seminar fuer Statistik, LEO C13
ETH Zurich	8092 Zurich	 SWITZERLAND
phone: x-41-44-632-4673		fax: 632-1228
http://stat.ethz.ch/~buser/
--------------------------------------------------------------

Henrik Parn writes:
 > 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 at 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.



More information about the R-help mailing list