[R] Problem With Model.Tables Function

[John Maindonald]

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