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? Are there special rules how incomplete block designs should look like to enable this recovery? Any help is appreciated! Regards, Christian
lme4 and incomplete block design
3 messages · Christian Lerch, Emmanuel Charpentier
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
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.