[R] how to constrast with factorial experiment
(Ted Harding)
Ted.Harding at nessie.mcc.ac.uk
Thu Aug 24 22:44:01 CEST 2006
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 ------------------------------
More information about the R-help
mailing list