[R] Example of cell means model

Douglas Bates bates at stat.wisc.edu
Wed Oct 15 15:40:33 CEST 2003


This is an example from chapter 11 of the 6th edition of Devore's
engineering statistics text.  It happens to be a balanced data set in
two factors but the calculations will also work for unbalanced data.

I create a factor called 'cell' from the text representation of the
Variety level and the Density level using '/' as the separator
character. The coefficients for the linear model 'Yield ~ cell - 1'
are exactly the cell means.

> library(Devore6)
> data(xmp11.07)
> str(xmp11.07)
`data.frame':	36 obs. of  3 variables:
 $ Yield  : num  10.5 9.2 7.9 12.8 11.2 13.3 12.1 12.6 14 10.8 ...
 $ Variety: Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
 $ Density: Factor w/ 4 levels "1","2","3","4": 1 1 1 2 2 2 3 3 3 4 ...
> xmp11.07$cell = with(xmp11.07, factor(paste(Variety, Density, sep = '/')))
> xtabs(~ Variety + Density, xmp11.07)
       Density
Variety 1 2 3 4
      1 3 3 3 3
      2 3 3 3 3
      3 3 3 3 3
> means = xtabs(Yield ~ Variety+Density, xmp11.07)/xtabs(~ Variety + Density, xmp11.07)
> means
       Density
Variety 1         2         3         4        
      1  9.200000 12.433333 12.900000 10.800000
      2  8.933333 12.633333 14.500000 12.766667
      3 16.300000 18.100000 19.933333 18.166667
> coef(fm1 <- lm(Yield ~ cell - 1, xmp11.07))
  cell1/1   cell1/2   cell1/3   cell1/4   cell2/1   cell2/2   cell2/3   cell2/4 
 9.200000 12.433333 12.900000 10.800000  8.933333 12.633333 14.500000 12.766667 
  cell3/1   cell3/2   cell3/3   cell3/4 
16.300000 18.100000 19.933333 18.166667 

The residual sum of squares for this model is the same as that for the
model with crossed terms

> deviance(fm1)
[1] 38.04
> deviance(fm2 <- lm(Yield ~ Variety * Density, xmp11.07))
[1] 38.04

but the coefficients are quite different because they represent a
different parameterization.

> coef(fm2)
      (Intercept)          Variety2          Variety3          Density2 
       9.20000000       -0.26666667        7.10000000        3.23333333 
         Density3          Density4 Variety2:Density2 Variety3:Density2 
       3.70000000        1.60000000        0.46666667       -1.43333333 
Variety2:Density3 Variety3:Density3 Variety2:Density4 Variety3:Density4 
       1.86666667       -0.06666667        2.23333333        0.26666667 

I hope this answers your question.  Sorry for the delay in
responding to you.

"Francisco Vergara" <gerifalte28 at hotmail.com> writes:

> Many thanks for your reply.  An example would be very helpful to make
> sure that I understand correctly what you described.




More information about the R-help mailing list