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

From: Chuck Cleland <ccleland_at_optonline.net>
Date: Sun, 09 Mar 2008 06:45:04 -0400

On 3/8/2008 3:05 PM, Dale Steele wrote:
> 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

   Here is one way to reproduce the SAS contrasts:

   ## 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'))

library(gmodels)

fm <- lm(voltage ~ material:temp + 0, data = twoway)

cm <- rbind(

  '21-22-31+32        ' = c(  0,   1, -1,  0,  -1,1,   0,   0,   0),
  'material1-material2' = c(1/3,-1/3,  0,1/3,-1/3,0, 1/3,-1/3,   0),
  'temp50-temp80      ' = c(1/3, 1/3,1/3,  0,   0,0,-1/3,-1/3,-1/3))

estimable(fm, cm)
                      Estimate Std. Error   t value DF     Pr(>|t|)
21-22-31+32          37.75000   25.98486  1.452769 27 1.578118e-01
material1-material2 -25.16667   10.60827 -2.372362 27 2.505884e-02
temp50-temp80        80.66667   10.60827  7.604127 27 3.525248e-08

   This formulates the model so that each coefficient corresponds to one of the 9 cell means. For me, that makes specifying the contrasts much easier.

> /* 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.
-- 
Chuck Cleland, Ph.D.
NDRI, Inc.
71 West 23rd Street, 8th floor
New York, NY 10010
tel: (212) 845-4495 (Tu, Th)
tel: (732) 512-0171 (M, W, F)
fax: (917) 438-0894

______________________________________________
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 Sun 09 Mar 2008 - 10:48:59 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 - 11: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