Dale Steele
2008-Mar-07 02:07 UTC
[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;
Thompson, David (MNR)
2008-Mar-07 21:20 UTC
[R] Finding Interaction and main effects contrasts for two-way ANOVA
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 secondavailable 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;