[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