[R] lme4 and incomplete block design

Emmanuel Charpentier emm.charpentier at free.fr
Sun Nov 8 10:34:04 CET 2009


Le dimanche 08 novembre 2009 à 00:05 +0100, Christian Lerch a écrit :
> Dear list members,
> 
> I try to simulate an incomplete block design in which every participants 
> receives 3 out of 4 possible treatment. The outcome in binary.
> 
> Assigning a binary outcome to the BIB or PBIB dataset of the package 
> SASmixed gives the appropriate output.
> With the code below, fixed treatment estimates are not given for each of 
> the 4 possible treatments, instead a kind of summary measure(?) for 
> 'treatment' is given.
> 
> block<-rep(1:24,each=3)
> treatment<-c(1, 2, 3, 2, 1, 4, 3, 4, 1, 4, 3, 2, 1, 3, 4, 2, 4, 3, 3, 1, 
> 2, 4, 2, 1, 1, 4, 2, 2, 3, 1, 3, 2, 4, 4, 1, 3, 1, 2, 3, 2, 1, 4, 3, 4, 
> 1, 4, 3, 2, 1, 3, 4, 2, 4, 3, 3, 1, 2, 4, 2, 1, 1, 4, 2, 2, 3, 1, 3, 2, 
> 4, 4, 1, 3)
> outcome<-c(0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 
> 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 
> 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 
> 0, 0, 1, 0)
> data<-data.frame(block,treatment,outcome)
> lmer(outcome~treatment +(1|block), family=binomial, data=data)
> 
> Is this a problem with the recovery of interblock estimates?

No...

>                                                               Are there 
> special rules how incomplete block designs should look like to enable 
> this recovery?

Neither...

Compare :

> library(lme4)
Le chargement a nécessité le package : Matrix
Le chargement a nécessité le package : lattice
> 
> summary(lmer(outcome~treatment +(1|block), family=binomial,
+              data=data.frame(block<-rep(1:24,each=3),
+                treatment<-c(1, 2, 3, 2, 1, 4, 3, 4, 1, 4, 3, 2, 1, 3, 4, 2, 4,
+                             3, 3, 1, 2, 4, 2, 1, 1, 4, 2, 2, 3, 1, 3, 2, 4, 4,
+                             1, 3, 1, 2, 3, 2, 1, 4, 3, 4, 1, 4, 3, 2, 1, 3, 4,
+                             2, 4, 3, 3, 1, 2, 4, 2, 1, 1, 4, 2, 2, 3, 1, 3, 2,
+                             4, 4, 1, 3),
+                outcome<-c(0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0,
+                           0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0,
+                           0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1,
+                           0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0))
+              ))
Generalized linear mixed model fit by the Laplace approximation 
Formula: outcome ~ treatment + (1 | block) 
   Data: data.frame(block <- rep(1:24, each = 3), treatment <- c(1, 2,      3, 2, 1, 4, 3, 4, 1, 4, 3, 2, 1, 3, 4, 2, 4, 3, 3, 1, 2,      4, 2, 1, 1, 4, 2, 2, 3, 1, 3, 2, 4, 4, 1, 3, 1, 2, 3, 2,      1, 4, 3, 4, 1, 4, 3, 2, 1, 3, 4, 2, 4, 3, 3, 1, 2, 4, 2,      1, 1, 4, 2, 2, 3, 1, 3, 2, 4, 4, 1, 3), outcome <- c(0, 0,      0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0,      0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0,      0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0,      0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0)) 
   AIC  BIC logLik deviance
 86.28 93.1 -40.14    80.28
Random effects:
 Groups Name        Variance Std.Dev.
 block  (Intercept) 0.60425  0.77734 
Number of obs: 72, groups: block, 24

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)  
(Intercept) -1.27783    0.71767  -1.780    0.075 .
treatment    0.01162    0.25571   0.045    0.964  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Correlation of Fixed Effects:
          (Intr)
treatment -0.892

with :

> summary(lmer(outcome~treatment +(1|block), family=binomial,
+              data=data.frame(block<-rep(1:24,each=3),
+                treatment<-factor(c(1, 2, 3, 2, 1, 4, 3, 4, 1, 4, 3, 2, 1, 3, 4,
+                                    2, 4, 3, 3, 1, 2, 4, 2, 1, 1, 4, 2, 2, 3, 1,
+                                    3, 2, 4, 4, 1, 3, 1, 2, 3, 2, 1, 4, 3, 4, 1,
+                                    4, 3, 2, 1, 3, 4, 2, 4, 3, 3, 1, 2, 4, 2, 1,
+                                    1, 4, 2, 2, 3, 1, 3, 2, 4, 4, 1, 3)),
+                outcome<-c(0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0,
+                           0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0,
+                           0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1,
+                           0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0))
+              ))
Generalized linear mixed model fit by the Laplace approximation 
Formula: outcome ~ treatment + (1 | block) 
   Data: data.frame(block <- rep(1:24, each = 3), treatment <- factor(c(1,      2, 3, 2, 1, 4, 3, 4, 1, 4, 3, 2, 1, 3, 4, 2, 4, 3, 3, 1,      2, 4, 2, 1, 1, 4, 2, 2, 3, 1, 3, 2, 4, 4, 1, 3, 1, 2, 3,      2, 1, 4, 3, 4, 1, 4, 3, 2, 1, 3, 4, 2, 4, 3, 3, 1, 2, 4,      2, 1, 1, 4, 2, 2, 3, 1, 3, 2, 4, 4, 1, 3)), outcome <- c(0,      0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0,      0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1,      0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0,      0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0)) 
   AIC   BIC logLik deviance
 87.33 98.72 -38.67    77.33
Random effects:
 Groups Name        Variance Std.Dev.
 block  (Intercept) 0.86138  0.9281  
Number of obs: 72, groups: block, 24

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept)  -1.9246     0.7117  -2.704  0.00684 **
treatment2    1.3910     0.8568   1.624  0.10446   
treatment3    0.4527     0.9163   0.494  0.62124   
treatment4    0.4526     0.9163   0.494  0.62131   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Correlation of Fixed Effects:
           (Intr) trtmn2 trtmn3
treatment2 -0.775              
treatment3 -0.721  0.598       
treatment4 -0.721  0.598  0.558


In the first case (your original "data"), "treatment" is interpreted as
a numeric (quantitative) variable , and whr lmre estimtes is a logistic
regression coefficient of the outcome n this numeric variable. Probbly
nonsensical, unless you hve reason to thin that your factor is ordered
and should be treated as numeric).

In the second case, "treatment" is a factor, so you get an estimate for
each treatment level except the first, to be interpreted as difference
of means with the first level.

I fell in that trap myself a few times, and took the habit to give evels
to my fctors tht cannot be interpreted as numbers (such as f<-paste("F",
as.character(v))).

> Any help is appreciated!

HTH,

					Emmanuel Charpentier



More information about the R-help mailing list