[R-sig-ME] Error message in lmer
Douglas Bates
bates at stat.wisc.edu
Mon Apr 20 21:25:20 CEST 2009
I would go about things in a different way than is described in that book.
This is an example of what happens when the data are organized with
"implicitly nested" factors. My recollection of the data is that
there are 6 rats and 3 liver sections from each rat, creating a total
of 18 liver sections in which the glycogen was measured. It may have
made sense in the past to designate the rats as 1 and 2 with the
implicit assumption that rat 1 in treatment 1 is different from rat 1
in treatment 2 or rat 1 in treatment 3 but I don't think this is a
need any more for that potentially confusing designation.
If instead you think of the three treatments as a factor that we will
model with fixed effects paramaters and the 6 rats and 18 liver
sections as defining random effects then the analysis is
straightforward.
> rats <- read.delim("/home/bates/Desktop/rats.txt")
> str(rats)
'data.frame': 36 obs. of 4 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 ...
> rats$Treatment <- factor(rats$Treatment, labels = LETTERS[1:3])
> rats$rr <- with(rats, Treatment:factor(Rat))
> rats$ll <- with(rats, Treatment:factor(Rat):factor(Liver))
> str(rats)
'data.frame': 36 obs. of 6 variables:
$ Glycogen : int 131 130 131 125 136 142 150 148 140 143 ...
$ Treatment: Factor w/ 3 levels "A","B","C": 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 ...
$ rr : Factor w/ 6 levels "A:1","A:2","B:1",..: 1 1 1 1 1 1 2 2 2 2 ...
$ ll : Factor w/ 18 levels "A:1:1","A:1:2",..: 1 1 2 2 3 3 4 4 5 5 ...
> rats
Glycogen Treatment Rat Liver rr ll
1 131 A 1 1 A:1 A:1:1
2 130 A 1 1 A:1 A:1:1
3 131 A 1 2 A:1 A:1:2
4 125 A 1 2 A:1 A:1:2
5 136 A 1 3 A:1 A:1:3
6 142 A 1 3 A:1 A:1:3
7 150 A 2 1 A:2 A:2:1
8 148 A 2 1 A:2 A:2:1
9 140 A 2 2 A:2 A:2:2
10 143 A 2 2 A:2 A:2:2
11 160 A 2 3 A:2 A:2:3
12 150 A 2 3 A:2 A:2:3
13 157 B 1 1 B:1 B:1:1
14 145 B 1 1 B:1 B:1:1
15 154 B 1 2 B:1 B:1:2
16 142 B 1 2 B:1 B:1:2
17 147 B 1 3 B:1 B:1:3
18 153 B 1 3 B:1 B:1:3
19 151 B 2 1 B:2 B:2:1
20 155 B 2 1 B:2 B:2:1
21 147 B 2 2 B:2 B:2:2
22 147 B 2 2 B:2 B:2:2
23 162 B 2 3 B:2 B:2:3
24 152 B 2 3 B:2 B:2:3
25 134 C 1 1 C:1 C:1:1
26 125 C 1 1 C:1 C:1:1
27 138 C 1 2 C:1 C:1:2
28 138 C 1 2 C:1 C:1:2
29 135 C 1 3 C:1 C:1:3
30 136 C 1 3 C:1 C:1:3
31 138 C 2 1 C:2 C:2:1
32 140 C 2 1 C:2 C:2:1
33 139 C 2 2 C:2 C:2:2
34 138 C 2 2 C:2 C:2:2
35 134 C 2 3 C:2 C:2:3
36 127 C 2 3 C:2 C:2:3
> library(lme4)
Loading required package: Matrix
Attaching package: 'Matrix'
The following object(s) are masked from package:stats :
xtabs
The following object(s) are masked from package:base :
rcond
> (fm1 <- lmer(Glycogen ~ Treatment + (1|rr) + (1|ll), rats))
Linear mixed model fit by REML
Formula: Glycogen ~ Treatment + (1 | rr) + (1 | ll)
Data: rats
AIC BIC logLik deviance REMLdev
231.6 241.1 -109.8 234.3 219.6
Random effects:
Groups Name Variance Std.Dev.
ll (Intercept) 14.167 3.7639
rr (Intercept) 36.065 6.0054
Residual 21.167 4.6007
Number of obs: 36, groups: ll, 18; rr, 6
Fixed effects:
Estimate Std. Error t value
(Intercept) 140.500 4.707 29.851
TreatmentB 10.500 6.656 1.577
TreatmentC -5.333 6.656 -0.801
Correlation of Fixed Effects:
(Intr) TrtmnB
TreatmentB -0.707
TreatmentC -0.707 0.500
Also, if you look at the enclosed plot you will see that the main
reason for not having a significant difference due to treatment is
because the two rats who got treatment A had very different levels of
glycogen. Furthermore there is considerable section to section
variability within rat and even assay to assay variability within the
same section (look at rat B:1's data).
On Thu, Apr 16, 2009 at 7:40 PM, Luciano La Sala
<lucianolasala at yahoo.com.ar> wrote:
>
> Dear people,
>
> I am doing some exercises from The R Book by Crawley. When trying to fit a mixed model for the “Rat” data (page 649 on this book), after loading the data and specifying the model I get the following error messages:
>
> model <- lmer(Glycogen~Treatment+(1|Treatment/Rat/Liver)
>
> Error: Matrices must have same number of columns in rbind2(..1, r)
> In addition: Warning messages:
> 1: In Rat:Treatment :
> numerical expression has 36 elements: only the first used
> 2: In Rat:Treatment :
> numerical expression has 36 elements: only the first used
> 3: In Liver:(Rat:Treatment) :
> numerical expression has 36 elements: only the first used
> 4: In Rat:Treatment :
> numerical expression has 36 elements: only the first used
> 5: In Rat:Treatment :
> numerical expression has 36 elements: only the first used
>
> What's all this mean? Any tips as to where I may be going wrong?
>
> Thank you in advance!
>
> Luciano
>
>
>
>
>
> Yahoo! Cocina
> Recetas prácticas y comida saludable
> http://ar.mujer.yahoo.com/cocina/
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
More information about the R-sig-mixed-models
mailing list