Re: [R] Finding Interaction and main effects contrasts for two-wayANOVA

From: Dale Steele <Dale_Steele_at_brown.EDU>
Date: Sat, 08 Mar 2008 15:05:55 -0500

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;
> >
>
>
> > ______________________________________________
> > R-help_at_r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide http://www.R-project.org/posting-
> > guide.html
> > and provide commented, minimal, self-contained, reproducible code.
>
>



R-help_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Sat 08 Mar 2008 - 20:10:33 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Sun 09 Mar 2008 - 13:30:20 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive