[R] lme4 and incomplete block design
Christian Lerch
t.c.l at gmx.net
Sun Nov 8 20:21:49 CET 2009
Many thanks, Bill and Emmanuel!
Christian
Emmanuel Charpentier schrieb:
> 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
>
> ______________________________________________
> R-help at r-project.org 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