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 <[EMAIL PROTECTED]> 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
>> > [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.
--
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@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.
