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@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.