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
> [EMAIL PROTECTED]
> http://ofri.mnr.gov.on.ca
> *************************************
>> -----Original Message-----
>> From: Steele [mailto:[EMAIL PROTECTED]
>> Sent: March 6, 2008 09:08 PM
>> To: [EMAIL PROTECTED]
>> 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@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@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.
