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 <[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.