Here is my assessment: The following seems pertinent to making sense of the result from model.tables(,type="means"). As I see it, users should be warned off using model.tables() with data from designs in which treatments are not balanced over blocking factors. df <- data.frame(bl=factor(c(1:4,1,2,4,1,2,3)), tr <- factor(rep(1:3,c(4,3,3))), y=c(10,12,9,11,13,15,16,18,22,17)) options(digits=4) # First fit the block effect, unadjusted for treatments blocks.aov <- aov(y~bl, data=df) # Take residuals from these block effects res <- residuals(blocks.aov) # Fit treatment effects to these residuals res.aov <- aov(res~df$tr) b <- summary.lm(res.aov)$coef # EITHER (1) # Add overall mean to the fitted treatment effects # NB: The residuals above summed to 1. So the average of # the fitted treatment effects is zero.
mean(df$y) + b[1]+c(0,b[2:3]) # with options(digits=4)
df$tr2 df$tr3 10.68 14.47 18.97
model.tables(df.aov, type="means")
Tables of means
Grand mean
14.3
bl
1 2 3 4
13.67 16.33 13 13.5
rep 3.00 3.00 2 2.0
tr
1 2 3
10.68 14.47 18.97
rep 4.00 3.00 3.00
# OR (2)
b <- summary.lm(res.aov)$coef b[1]+c(0, b[2:3])
df$tr2 df$tr3 -3.6250 0.1667 4.6667
# Equivalent to model.tables(df.aov)
---------------------------------------------------------- The objection to the treatment effects from model.tables() is that they are based on deviations from block means that have not been adjusted to allow for the different relative numbers of the different treatments in those blocks. The effect can be severe when different blocks have greatly different relative numbers of different treatments. Thus model.tables() should not be used even for balanced incomplete block designs. # Choices that are defendable include:
#1 Add treatment effects from above least squares fit,
# constrained to sum to zero, to overall mean.
options(contrasts=c("contr.sum","contr.poly"))
b <- summary.lm(aov(y~bl+tr, data=df))$coef
mean(df$y)+ c(b[5:6], -sum(b[5:6]))
tr1 tr2 10.16 13.67 19.07
#2 Add treatment effects from least squares fit, # constrained to sum to zero, to mean of block means. b[1] + c(b[5:6],-sum(b[5:6]))
tr1 tr2 10.50 14.01 19.41
#3 Use estimates from lme, with blocks as random effects
Brian Ripley, responding to Gary Whysong, wrote:
For the record, these R results agree with S-PLUS on your examples. But they are not much documented in the unbalanced case, so perhaps `for experts only'.
. .
On Sat, 10 Mar 2001, Gary Whysong wrote: [...]
blocks2<-factor(c(1,2,3,4,1,2,4,1,2,3)) trtmts2<-factor(c(1,1,1,1,2,2,2,3,3,3,)) data2<-c(10,12,9,11,13,15,16,18,22,17) unbalanced<-aov(data2~blocks2+trtmts2) summary(unbalanced)
Df Sum Sq Mean Sq F value Pr(>F) blocks2 3 18.267 6.089 7.4341 0.0410993 * trtmts2 2 126.557 63.279 77.2587 0.0006367 *** Residuals 4 3.276 0.819 --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
model.tables(unbalanced,"means")
Tables of means
Grand mean
14.3
blocks2
1 2 3 4
13.67 16.33 13 13.5
rep 3.00 3.00 2 2.0
trtmts2
1 2 3
10.68 14.47 18.97
rep 4.00 3.00 3.00
We find that the treatment means (trtmts2) are incorrect although the
number of replications indicated are correct. Block means (blocks2) are
correct.
The treatment means should be: 10.5, 14.67, and 19.0, respectively.
Why? These are *model*.tables not *data*.tables. You have to adjust for block effects, and they are unbalanced. Blocks 3 and 4 have lower responses than 1 and 2, and they are missing for treatments 2 and 3. Seems to adjust correctly to me. Note the order of terms matters in unbalanced models.
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html Send "info", "help", or "[un]subscribe" (in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._