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

Gregory Warnes gregory.warnes at mac.com
Sat Mar 8 17:02:34 CET 2008


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.



More information about the R-help mailing list