[R] how to contrast with factorial experiment
szhan at uoguelph.ca
szhan at uoguelph.ca
Tue Aug 29 18:09:33 CEST 2006
Hello, R experts,
If I understand Ted's anwser correctly, then I can not contrast the
mean yields between sections 1-8 and 9-11 under "Trt" but I can
contrast mean yields for sections 1-3 and 6-11 because there exists
significant interaction between two factors (Trt:section4,
Trt:section5). Could I use the commands below to test
the difference between sections 1-3 and 6-11 ?
> contrasts(section)<-c(-2,-2,-2,0,0,1,1,1,1,1,1)
> newobj<-lm(log2(yield)~treat*section)
How can I infer that there is significant difference between sections
1-3 and sections 6-11 for the "Trt" from the output below?
> summary(newobj)
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.28840 0.05003 125.682 < 2e-16 ***
treatTrt 1.22122 0.07076 17.259 < 2e-16 ***
section1 0.17831 0.03911 4.559 4.08e-05 ***
section2 -0.23102 0.16595 -1.392 0.17087
section3 2.38170 0.16595 14.352 < 2e-16 ***
section4 3.36834 0.16595 20.298 < 2e-16 ***
section5 -1.56873 0.16595 -9.453 3.67e-12 ***
section6 -0.41522 0.16595 -2.502 0.01613 *
section7 -0.89943 0.16595 -5.420 2.38e-06 ***
section8 0.09522 0.16595 0.574 0.56901
section9 -0.78784 0.16595 -4.748 2.21e-05 ***
section10 0.74821 0.16595 4.509 4.79e-05 ***
treatTrt:section1 0.10101 0.05532 1.826 0.07461 .
treatTrt:section2 0.27270 0.23468 1.162 0.25151
treatTrt:section3 -1.22210 0.23468 -5.207 4.85e-06 ***
treatTrt:section4 -1.39187 0.23468 -5.931 4.26e-07 ***
treatTrt:section5 -0.76137 0.23468 -3.244 0.00225 **
treatTrt:section6 0.07320 0.23468 0.312 0.75658
treatTrt:section7 0.33108 0.23468 1.411 0.16535
treatTrt:section8 -0.13686 0.23468 -0.583 0.56276
treatTrt:section9 0.22086 0.23468 0.941 0.35180
treatTrt:section10 -0.14476 0.23468 -0.617 0.54054
---
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
Joshua
Quoting "(Ted Harding)" <Ted.Harding at nessie.mcc.ac.uk>:
> On 24-Aug-06 szhan at uoguelph.ca wrote:
>> 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)
> [I think you mean sections 1 to 8]
>
>> 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?
>
> I think you would be well advised to look at a plot of the data.
> For example, let Y stand for yield, R for replicate, T for treat
> and S for section.
>
> ix<-(T=="Trt");plot(S[ix],Y[ix],col="red",ylim=c(0,1000))
> ix<-(T=="Ctl");points(S[ix],Y[ix],col="blue")
>
>> From this it is clear that sections 4 and 5 are in a class of
> their own. Also, in sections 1-3 and 6-11 the "Ctl" yields
> are not only lower, but have smaller (in some cases hardly any)
> variance, compared with the "Trt" yields. The variances for
> sections 7,8,9,10,11 are greater than for 1,2,3,6 without
> great change in mean value.
>
> While there is an evident difference between "Trt" yields and
> "Ctrl" yields for sections 1-3 and 6-11, this is not so for
> sections 4 and 5.
>
> This sort of behaviour no doubt provides some reasons for the
> interaction you observed. You seem to have a quite complex
> phenomenon here!
>
> To some extent the problems with variance can be diminished by
> working with logarithms. Compare the previous plot with
>
> ix<-(T=="Trt");plot(S[ix],log10(Y[ix]),col="red",ylim=c(0,3))
> ix<-(T=="Ctl");points(S[ix],log10(Y[ix]),col="blue")
>
> (you have used log2() in your commands). The above observations
> can be seen reflected in R if you look at the output of
>
> summary(obj)
>
> where in particular:
>
> treatTrt:section2 -1.11691 0.33189 -3.365 0.001595 **
> treatTrt:section3 -0.45634 0.33189 -1.375 0.176099
> treatTrt:section4 -1.56627 0.33189 -4.719 2.42e-05 ***
> treatTrt:section5 -1.73604 0.33189 -5.231 4.48e-06 ***
> treatTrt:section6 -0.91311 0.33189 -2.751 0.008588 **
> treatTrt:section7 -0.07853 0.33189 -0.237 0.814055
> treatTrt:section8 0.17935 0.33189 0.540 0.591654
> treatTrt:section9 -0.28859 0.33189 -0.870 0.389277
> treatTrt:section10 0.06913 0.33189 0.208 0.835972
> treatTrt:section11 -0.29649 0.33189 -0.893 0.376543
>
> which, precisely, "contrasts means from levels of one factor
> (section) under another factor (treat)", and shows that most
> of the "interaction" arises in sections 4 and 5.
>
> Since sections 4 and 5 (in the middle of sections 1 to 8) are
> so exceptional, they will have strong influence on your comparison
> between sections 1-8 and sections 9-11. You need to think about
> what to do with sections 4 and 5!
>
>> 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
>>
>> ______________________________________________
>> R-help at stat.math.ethz.ch 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.
>
> --------------------------------------------------------------------
> E-Mail: (Ted Harding) <Ted.Harding at nessie.mcc.ac.uk>
> Fax-to-email: +44 (0)870 094 0861
> Date: 24-Aug-06 Time: 21:43:57
> ------------------------------ XFMail ------------------------------
>
> ______________________________________________
> R-help at stat.math.ethz.ch 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 at stat.math.ethz.ch 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.
More information about the R-help
mailing list