[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