[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