[R] how to constrast with factorial experiment
szhan at uoguelph.ca
szhan at uoguelph.ca
Thu Aug 24 21:41:03 CEST 2006
Hello, R users,
I have two factors (treat, section) anova design experiment where
there are 3 replicates. The objective of the experiment is to test if
there is significant difference of yield between top (section 9 to 11)
and bottom (section 9 to 11) of the fruit tree under treatment. I
found that there are interaction between two factors. I wonder if I
can contrast means from levels of one factor (section) under another
factor (treat)? if so, how to do it in R and how to interpret the
output?
Here is the data and commands I used to test the differece between
section 1 to 8 and 9 to 11 under treatment. But I don't know if I was
right, how to interpret the out and whether there are significant
difference between section 1 to 8 and section 9 to 11 under treatment.
yield replicate treat section
35.55 1 Ctl 1
53.70 1 Ctl 2
42.79 1 Ctl 3
434.81 1 Ctl 4
705.96 1 Ctl 5
25.91 1 Ctl 6
57.53 1 Ctl 7
41.45 1 Ctl 8
85.54 1 Ctl 9
51.23 1 Ctl 10
188.24 1 Ctl 11
35.71 2 Ctl 1
45.15 2 Ctl 2
40.10 2 Ctl 3
312.76 2 Ctl 4
804.05 2 Ctl 5
28.22 2 Ctl 6
68.51 2 Ctl 7
46.15 2 Ctl 8
123.14 2 Ctl 9
33.78 2 Ctl 10
121.28 2 Ctl 11
30.96 3 Ctl 1
36.10 3 Ctl 2
47.19 3 Ctl 3
345.80 3 Ctl 4
644.61 3 Ctl 5
27.73 3 Ctl 6
56.63 3 Ctl 7
42.63 3 Ctl 8
61.25 3 Ctl 9
59.43 3 Ctl 10
109.87 3 Ctl 11
143.50 1 Trt 1
82.76 1 Trt 2
125.03 1 Trt 3
493.76 1 Trt 4
868.48 1 Trt 5
45.09 1 Trt 6
249.43 1 Trt 7
167.28 1 Trt 8
274.72 1 Trt 9
176.40 1 Trt 10
393.10 1 Trt 11
93.75 2 Trt 1
63.83 2 Trt 2
117.50 2 Trt 3
362.68 2 Trt 4
659.40 2 Trt 5
62.10 2 Trt 6
218.24 2 Trt 7
210.98 2 Trt 8
291.48 2 Trt 9
209.36 2 Trt 10
454.68 2 Trt 11
119.62 3 Trt 1
66.50 3 Trt 2
87.37 3 Trt 3
414.01 3 Trt 4
707.70 3 Trt 5
44.40 3 Trt 6
142.59 3 Trt 7
137.37 3 Trt 8
181.03 3 Trt 9
131.65 3 Trt 10
310.18 3 Trt 11
> dat1<-read.delim("c:/testcontr.txt", header=T)
> dat1$treat<-as.factor(dat1$treat)
> dat1$replicate<-as.factor(dat1$replicate)
> dat1$section<-as.factor(dat1$section)
> attach(dat1)
> obj<-lm(log2(yield)~treat*section)
> anova(obj)
Analysis of Variance Table
Response: log2(yield)
Df Sum Sq Mean Sq F value Pr(>F)
treat 1 24.608 24.608 297.8649 < 2.2e-16 ***
section 10 99.761 9.976 120.7565 < 2.2e-16 ***
treat:section 10 6.708 0.671 8.1197 2.972e-07 ***
Residuals 44 3.635 0.083
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
> contrasts(section)<-c(3,3,3,3,3,3,3,3,-8,-8,-8)
> objnew<-lm(log2(yield)~treat*section)
> summary(objnew)
Call:
lm(formula = log2(yield) ~ treat * section)
Residuals:
Min 1Q Median 3Q Max
-0.49647 -0.14913 -0.01521 0.17471 0.51105
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.288403 0.050034 125.682 < 2e-16 ***
treatTrt 1.221219 0.070759 17.259 < 2e-16 ***
section1 -0.008502 0.010213 -0.832 0.409675
section2 -0.491175 0.165945 -2.960 0.004942 **
section3 2.569427 0.165945 15.484 < 2e-16 ***
section4 3.556067 0.165945 21.429 < 2e-16 ***
section5 -1.157069 0.165945 -6.973 1.25e-08 ***
section6 -0.003562 0.165945 -0.021 0.982971
section7 -0.487770 0.165945 -2.939 0.005223 **
section8 0.106181 0.165945 0.640 0.525585
section9 -0.776882 0.165945 -4.682 2.74e-05 ***
section10 0.759168 0.165945 4.575 3.87e-05 ***
treatTrt:section1 -0.049000 0.014444 -3.392 0.001474 **
treatTrt:section2 0.160825 0.234682 0.685 0.496757
treatTrt:section3 -0.949101 0.234682 -4.044 0.000208 ***
treatTrt:section4 -1.118870 0.234682 -4.768 2.07e-05 ***
treatTrt:section5 -0.295937 0.234682 -1.261 0.213950
treatTrt:section6 0.538638 0.234682 2.295 0.026549 *
treatTrt:section7 0.796518 0.234682 3.394 0.001468 **
treatTrt:section8 -0.548744 0.234682 -2.338 0.023984 *
treatTrt:section9 -0.191029 0.234682 -0.814 0.420033
treatTrt:section10 -0.556642 0.234682 -2.372 0.022137 *
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
Residual standard error: 0.2874 on 44 degrees of freedom
Multiple R-Squared: 0.973, Adjusted R-squared: 0.9601
F-statistic: 75.55 on 21 and 44 DF, p-value: < 2.2e-16
Thanks,
Joshua
More information about the R-help
mailing list