Finding Interaction and main effects contrasts for two-wayANOVA

Thanks to those who have replied to my original query.  However, I'm
still confused on how obtain estimates, standard error and F-tests for
main effect and interaction contrasts which agree with the SAS code
with output appended below.

for example,
 ## Given the dataset (from Montgomery)
twoway <- read.table("http://dsteele.veryspeedy.net/sta501/twoway.txt",
col.names=c('material', 'temp','voltage'),colClasses=c('factor',
'factor', 'numeric'))

 ## the model
fit <- aov(voltage ~ material*temp, data=twoway)

material.means <- tapply(twoway$voltage, twoway$material, mean)
temp.means <- tapply(twoway$voltage, twoway$temp, mean)
cell.means <- tapply(twoway$voltage, twoway[,1:2], mean)

Contrasts of Interest ....

cell.means[2,1] - cell.means[2,2] - cell.means[3,1] + cell.means[3,2]

[1] 37.75

material.means[1] - material.means[2]

        1
-25.16667

temp.means[1] - temp.means[3]

      50
80.66667


I expected the following code to provide the estimates above  for
(material 1 - material 2) and (temp1 - temp3), but get unexpected
results...

library(gmodels)
fit.contrast(fit, "material", rbind("(1 - 2)" =c(1, -1, 0) ))

                Estimate Std. Error   t value  Pr(>|t|)
material(1 - 2)      -21   18.37407 -1.142915 0.2631074

fit.contrast(fit, "temp", rbind("50 - 80" =c(1, 0, -1) ))

            Estimate Std. Error  t value     Pr(>|t|)
temp50 - 80    77.25   18.37407 4.204294 0.0002572756

Thanks.  --Dale

/* SAS code */
proc glm data=twoway;
class material temp;
model voltage = material temp material*temp;
contrast '21-22-31+32' material*temp 0 0 0 1 -1 0 -1 1 0;
estimate '21-22-31+32' material*temp 0 0 0 1 -1 0 -1 1 0;
contrast 'material1-material2' material 1 -1 0;
estimate 'material1-material2' material 1 -1 0;
contrast 'temp50 - temp80' temp 1 0 -1;
estimate 'temp50 - temp80' temp 1 0 -1;
run;

SAS output

 Contrast                   DF    Contrast SS    Mean Square   F Value   Pr > F

 21-22-31+32                 1     1425.06250     1425.06250      2.11   0.1578
 material1-material2         1     3800.16667     3800.16667      5.63   0.0251
 temp50 - temp80             1    39042.66667    39042.66667     57.82   <.0001


                                              Standard
  Parameter                   Estimate           Error    t Value    Pr > |t|

  21-22-31+32               37.7500000      25.9848603       1.45      0.1578
  material1-material2      -25.1666667      10.6082748      -2.37      0.0251
  temp50 - temp80           80.6666667      10.6082748       7.60      <.0001

On Sat, Mar 8, 2008 at 11:02 AM, Gregory Warnes <gregory.warnes at mac.com> wrote:

 Dale,

 You might find it fruitful to look at the help pages for fit.contrast
 () and estimble() functions in the gmodels package, and the contrast
 () functions in the Hmisc package.

 -Greg

 On Mar 7, 2008, at 4:20PM , Thompson, David ((MNR)) wrote:

 >  Dale,
 >
 > Other than the first SAS contrast, does the following demonstrate what
 > your asking for?

 >> summary(twoway)

 >  material temp       voltage
 >  1:12     50:12   Min.   : 20
 >  2:12     65:12   1st Qu.: 70
 >  3:12     80:12   Median :108
 >                   Mean   :106
 >                   3rd Qu.:142
 >                   Max.   :188

 >> contrasts(twoway$material)

 >   2 3
 > 1 0 0
 > 2 1 0
 > 3 0 1

 >> contrasts(twoway$temp)

 >    65 80
 > 50  0  0
 > 65  1  0
 > 80  0  1

 >> fit <- aov(voltage ~ material*temp, data=twoway)
 >> summary.aov(fit)

 >               Df Sum Sq Mean Sq F value  Pr(>F)
 > material       2  10684    5342    7.91  0.0020 **
 > temp           2  39119   19559   28.97 1.9e-07 ***
 > material:temp  4   9614    2403    3.56  0.0186 *
 > Residuals     27  18231     675
 > ---
 > Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
 >
 >
 > # setting (partial) contrasts

 >> contrasts(twoway$material) <- c(1,-1,0) # ignoring the second

 > available df

 >> contrasts(twoway$temp) <- c(0,1,-1) # ignoring the second
 >> available df
 >> contrasts(twoway$material)

 >   [,1]  [,2]
 > 1    1 -0.41
 > 2   -1 -0.41
 > 3    0  0.82

 >> contrasts(twoway$temp)

 >    [,1]  [,2]
 > 50    0 -0.82
 > 65    1  0.41
 > 80   -1  0.41
 >

 >> summary.aov(fit, split=list(material=list('m1-m2'=1), temp=list
 >> ('t50 -

 > t80'=1)))
 >                                  Df Sum Sq Mean Sq F value  Pr(>F)
 > material                          2  10684    5342    7.91 0.00198 **
 >   material: m1-m2                 1   3800    3800    5.63 0.02506 *
 > temp                              2  39119   19559   28.97 1.9e-07 ***
 >   temp: t50 - t80                 1  11310   11310   16.75 0.00035 ***
 > material:temp                     4   9614    2403    3.56 0.01861 *
 >   material:temp: m1-m2.t50 - t80  1   4970    4970    7.36 0.01146 *
 > Residuals                        27  18231     675
 > ---
 > Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
 >
 > # other examples of setting contrasts
 > # compare m1 vs m2 and m2 vs m3

 >> contrasts(twoway$material) <- matrix(c(1,-1,0,1,1,-2), nrow=3)
 >> contrasts(twoway$material)

 >   [,1] [,2]
 > 1    1    0
 > 2   -1    1
 > 3    0   -1
 > # compare m1 vs m2 and m1+m2 vs m3

 >> contrasts(twoway$material) <- matrix(c(1,-1,0,1,1,-2), nrow=3)
 >> contrasts(twoway$material)

 >   [,1] [,2]
 > 1    1    1
 > 2   -1    1
 > 3    0   -2
 >
 > I'm not sure if 'summary.aov' is the only lm-family summary method
 > with
 > the split argument.
 >
 > DaveT.
 > *************************************
 > Silviculture Data Analyst
 > Ontario Forest Research Institute
 > Ontario Ministry of Natural Resources
 > david.john.thompson at ontario.ca
 > http://ofri.mnr.gov.on.ca
 > *************************************

 >> -----Original Message-----
 >> From: Steele [mailto:dale.w.steele at gmail.com]
 >> Sent: March 6, 2008 09:08 PM
 >> To: r-help at stat.math.ethz.ch
 >> Subject: [R] Finding Interaction and main effects contrasts
 >> for two-way ANOVA
 >>
 >> I've tried  without success to calculate interaction and main effects
 >> contrasts using R.  I've found the functions C(), contrasts(),
 >> se.contrasts() and fit.contrasts() in package gmodels.  Given the url
 >> for a small dataset and the two-way anova model below, I'd like to
 >> reproduce the results from appended SAS code.  Thanks.  --Dale.
 >>
 >>  ## the dataset (from Montgomery)
 >> twoway <- read.table("http://dsteele.veryspeedy.net/sta501/
 >> twoway.txt",
 >> col.names=c('material', 'temp','voltage'),colClasses=c('factor',
 >> 'factor', 'numeric'))
 >>
 >>  ## the model
 >> fit <- aov(voltage ~ material*temp, data=twoway)
 >>
 >> /* SAS code */
 >> proc glm data=twoway;
 >> class material temp;
 >> model voltage = material temp material*temp;
 >> contrast '21-22-31+32' material*temp 0 0 0 1 -1 0 -1 1 0;
 >> estimate '21-22-31+32' material*temp 0 0 0 1 -1 0 -1 1 0;
 >> contrast 'material1-material2' material 1 -1 0;
 >> estimate 'material1-material2' material 1 -1 0;
 >> contrast 'temp50 - temp80' temp 1 0 -1;
 >> estimate 'temp50 - temp80' temp 1 0 -1;
 >> run;