Two-Way ANOVA Model, #2

Soybean Example

> soybean <− data.frame(yield=scan(),

+ density=rep( rep(c("regular","thick"),c(4,4)), 3 ),

+ fertilizer=factor( rep(1:3,c(8,8,8)), labels=c("low","medium",'high') )

+ )

1: 37.5 36.5 38.6 36.5 37.4 35.0 38.1 36.5

9: 48.1 48.3 48.6 46.4 36.7 36.4 39.3 37.5

17: 48.5 46.1 49.1 48.2 45.7 45.7 48.0 46.4

25:

Read 24 items

> soybean

   yield density fertilizer
1   37.5 regular        low
2   36.5 regular        low
3   38.6 regular        low
4   36.5 regular        low
5   37.4   thick        low
6   35.0   thick        low
7   38.1   thick        low
8   36.5   thick        low
9   48.1 regular     medium
10  48.3 regular     medium
11  48.6 regular     medium
12  46.4 regular     medium
13  36.7   thick     medium
14  36.4   thick     medium
15  39.3   thick     medium
16  37.5   thick     medium
17  48.5 regular       high
18  46.1 regular       high
19  49.1 regular       high
20  48.2 regular       high
21  45.7   thick       high
22  45.7   thick       high
23  48.0   thick       high
24  46.4   thick       high

> levels(soybean$density)

[1] "regular" "thick"  

> levels(soybean$fertilizer)

[1] "low"    "medium" "high"

The session below 'simulates' hand calculations (for demonstration purpose).

> print(Tij. <− xtabs(yield~.,data=soybean))

         fertilizer
density     low medium  high
  regular 149.1  191.4 191.9
  thick   147.0  149.9 185.8

> print(nij <− xtabs(rep(1,24)~density+fertilizer,data=soybean))

         fertilizer
density   low medium high
  regular   4      4    4
  thick     4      4    4

> Tij./nij # treatment means

         fertilizer
density      low medium   high
  regular 37.275 47.850 47.975
  thick   36.750 37.475 46.450

> print(Ti.. <− xtabs(yield~density,data=soybean))

density
regular   thick 
  532.4   482.7 

> Ti../xtabs(rep(1,24)~density,data=soybean) # means for density

density
 regular    thick 
44.36667 40.22500

> print(T.j. <− xtabs(yield~fertilizer,data=soybean))

fertilizer
   low medium   high 
 296.1  341.3  377.7

> T.j./xtabs(rep(1,24)~fertilizer,data=soybean) # means for fertilizer

fertilizer
    low  medium    high 
37.0125 42.6625 47.2125

The session below shows the analysis.

> soybean$comb <− factor(paste("(", unclass(soybean$density), ",",

+ unclass(soybean$fertilizer), ")", sep=""))

> soybean # with 'super' factor combo

   yield density fertilizer  comb
1   37.5 regular        low (1,1)
2   36.5 regular        low (1,1)
3   38.6 regular        low (1,1)
4   36.5 regular        low (1,1)
5   37.4   thick        low (2,1)
6   35.0   thick        low (2,1)
7   38.1   thick        low (2,1)
8   36.5   thick        low (2,1)
9   48.1 regular     medium (1,2)
10  48.3 regular     medium (1,2)
11  48.6 regular     medium (1,2)
12  46.4 regular     medium (1,2)
13  36.7   thick     medium (2,2)
14  36.4   thick     medium (2,2)
15  39.3   thick     medium (2,2)
16  37.5   thick     medium (2,2)
17  48.5 regular       high (1,3)
18  46.1 regular       high (1,3)
19  49.1 regular       high (1,3)
20  48.2 regular       high (1,3)
21  45.7   thick       high (2,3)
22  45.7   thick       high (2,3)
23  48.0   thick       high (2,3)
24  46.4   thick       high (2,3)

> # The following three methods of testing homogeneity

> # of variances give consistent conclusion (that

> # the variances are homogeneity).

> # Note that the second test is a (robust) nonparametric one.

> with(soybean, car::levene.test(yield, comb))

Levene's Test for Homogeneity of Variance
      Df F value Pr(>F)
group  5  0.1231 0.9854
      18

> fligner.test(yield~comb, data=soybean)

        Fligner-Killeen test of homogeneity of variances

data:  yield by comb 
Fligner-Killeen:med chi-squared = 2.4364, df = 5, p-value =
0.786

> bartlett.test(yield~comb, data=soybean)

        Bartlett test of homogeneity of variances

data:  yield by comb 
Bartlett's K-squared = 0.5172, df = 5, p-value = 0.9915

> oneway.test(yield~comb, data=soybean) # Welch F

        One-way analysis of means (not assuming equal variances)

data:  yield and comb 
F = 73.2282, num df = 5.000, denom df = 8.378, p-value =
1.131e-06

> oneway.test(yield~comb, data=soybean, var.equal=T) # F test

        One-way analysis of means

data:  yield and comb
F = 91.7897, num df = 5, denom df = 18, p-value = 3.638e-12

> soybean.aov <− aov(yield ~ density*fertilizer, data=soybean)

> summary(soybean.aov)

                   Df Sum Sq Mean Sq F value    Pr(>F)    
density             1 102.92  102.92  74.007 8.580e-08 ***
fertilizer          2 417.77  208.89 150.203 5.897e-12 ***
density:fertilizer  2 117.56   58.78  42.268 1.583e-07 ***
Residuals          18  25.03    1.39                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

> # The interaction is highly significant. This is evidenced by

> # the interaction plots below. See

> # Rfactory.html for explanation of

> # the commands below.

>

> cols <− paste("dark",c("magenta","green","blue"),sep="")

> windows(width=6, height=3, pointsize=10)

> par(mfrow=c(1,2), mar=c(4,4,1,2.5)+.1)

> par(oma=c(2,0,3,0)+.1, bty='l')

> with(soybean,

+ {interaction.plot(density,fertilizer,yield,type="b",

+ col=cols[1:3],pch=1:3)

+ box(which="figure",col="grey75")

+ interaction.plot(fertilizer,density,yield,type="b",

+ col=cols[1:2],pch=1:2)

+ box(which="figure",col="grey75")

+ }

+ )

> title(sub="Soybean Example",outer=T,line=0.5,adj=0)

> title(main="Interaction Plots",outer=T)

> box(which="outer",col="grey80")

> dev.off()

null device 
          1

(soybean01.png here)

The two main effects are also highly significant as evidenced by the main effect plots below. Note that the estimated experimental error (MSE) is quite small, and as consequences, it tends to be easy to detect small effects.

> ybar.A <− with(soybean,tapply(yield,density,mean))

> ybar.B <− with(soybean,tapply(yield,fertilizer,mean))

> ybar <− with(soybean, mean(yield))

> ybar.A # Factor A (density) treatment means

 regular    thick 
44.36667 40.22500 

> ybar.B # Factor B treatment means

    low  medium    high 
37.0125 42.6625 47.2125 

> ybar # grand mean

[1] 42.29583

> ylim <− range(pretty(range(ybar.A,ybar.B,ybar)))

> windows(width=6,height=3,pointsize=10)

> par(mfrow=c(1,2),mar=c(4,4,.5,.5)+.1, oma=c(2,0,3,0)+.1)

> plot(1:2,ybar.A,xlab="planting density",ylab="mean yield",

+ ylim=ylim,xaxt="n",type='l',lwd=2)

> abline(h=ybar, lty="dotted", col="darkorange")

> axis(1,at=1:2,lab=levels(soybean$density))

> box(which="figure", col="cyan",lwd=2)

> plot(1:3,ybar.B,xlab="fertilizer level",ylab="mean yield",

+ ylim=ylim,xaxt="n",type='l',lwd=2)

> axis(1,at=1:3,labels=levels(soybean$fertilizer))

> box(which="figure",col="cyan",lwd=2)

> abline(h=ybar, lty="dotted", col="darkorange")

> title(outer=T,

+ main="Main Effect Plots\nsoybean example")

> title(outer=T, line=0.5,

+ sub=bquote( bar(y)[...]==.(round(ybar,2)) )

+ )

> box(which="outer", col="blue", lwd=1.5)

> box(which="inner", col="darkorange2")

> dev.off()

null device 
          1

(soybean02.png here)

> ybar.AB <− tapply(soybean$yield, soybean[2:3], mean)

> ybar.AB # cell means

         fertilizer
density      low medium   high
  regular 37.275 47.850 47.975
  thick   36.750 37.475 46.450

> simple.A <− ybar.AB[1,]-ybar.AB[2,]

> first <− c(1,1,2); second <− c(2,3,3)

> simple.B <− ybar.AB[,c(1,1,2)]-ybar.AB[,c(2,3,3)]

> colnames(simple.B) <− paste(colnames(ybar.AB)[first],

+ colnames(ybar.AB)[second], sep=" - ")

> simple.B

         fertilizer
density   low - medium low - high medium - high
  regular      -10.575      -10.7        -0.125
  thick         -0.725       -9.7        -8.975

> MSE <− sum(residuals(soybean.aov)^2) / soybean.aov$df.resid

> MSE

[1] 1.390694

> stdev.simple <− sqrt(2*MSE/4) # n=4 replicates

> stdev.simple

[1] 0.8338748

> # Use protected LSD: reject each H0 simple effect

> # if |pt. est| > t.025,dfstdev.simple for

> # level-0.05 test.

> hw <− qt(.025, soybean.aov$df.residual, lower=F)*stdev.simple

> hw # half interval width

[1] 1.751906

> ifelse(abs(simple.A)>hw, "Reject H0", "Retain H0")

        low      medium        high 
"Retain H0" "Reject H0" "Retain H0"

> ifelse(abs(simple.B)>hw, "Reject H0", "Retain H0")

         fertilizer
density   low - medium low - high  medium - high
  regular "Reject H0"  "Reject H0" "Retain H0"  
  thick   "Retain H0"  "Reject H0" "Reject H0"

> # The P values are obtained by

> # 2 × P(tdf > |t0|), where

> # t0 = (pt. est) / stdev.simple

> t0.simple.A <− simple.A / stdev.simple

> t0.simple.A

       low     medium       high 
 0.6295909 12.4419154  1.8288117 

> t0.simple.B <− simple.B / stdev.simple

> t0.simple.B

         fertilizer
density   low - medium low - high medium - high
  regular   -12.681760  -12.83166    -0.1499026
  thick      -0.869435  -11.63244   -10.7630063

> P.simple.A <− 2 * pt(abs(t0.simple.A), soybean.aov$df.residual, lower=F)

> round(P.simple.A,5)

    low  medium    high 
0.53687 0.00000 0.08405

> P.simple.B <− 2 * pt(abs(t0.simple.B), soybean.aov$df.residual, lower=F)

> round(P.simple.B, 5)

         fertilizer
density   low - medium low - high medium - high
  regular      0.00000          0       0.88251
  thick        0.39606          0       0.00000

> ci.simple.A <− cbind(Estimate=simple.A,

+ lwr=simple.A-hw, upr=simple.A+hw)

> ci.simple.A

       Estimate       lwr       upr
low       0.525 -1.226906  2.276906
medium   10.375  8.623094 12.126906
high      1.525 -0.226906  3.276906

> print( simple.B <− t(simple.B) ) # transpose

               density
fertilizer      regular  thick
  low - medium  -10.575 -0.725
  low - high    -10.700 -9.700
  medium - high  -0.125 -8.975

> ci.simple.B <− cbind(Estimate=simple.B[1:6],

+ lwr=simple.B[1:6]-hw, upr=simple.B[1:6]+hw)

> FUN <− function(x,y) paste("(", x, ") within ", y, sep="")

> rownames(ci.simple.B) <− outer(rownames(simple.B),

+ colnames(simple.B), FUN)

> ci.simple.B

                               Estimate        lwr       upr
(low - medium) within regular   -10.575 -12.326906 -8.823094
(low - high) within regular     -10.700 -12.451906 -8.948094
(medium - high) within regular   -0.125  -1.876906  1.626906
(low - medium) within thick      -0.725  -2.476906  1.026906
(low - high) within thick        -9.700 -11.451906 -7.948094
(medium - high) within thick     -8.975 -10.726906 -7.223094

>

> # For Scheffe method, replace t.025,df above

> # by 5 × F0.05;5,df. That is, use

> hw <− qf(.05, 5,soybean.aov$df.residual, lower=F)*stdev.simple

> hw

[1] 2.312212

>

> # For Bonferroni method, adjust P values accordingly or adjust

> # t critical value.

> #

The session below shows the plot of simple A effects. It is built on the low-level graphics functions abline, axis, segments, points, title, and mtext. The high-level graphics function plot is used to set up the coordinate system. This is actually how such graph can be built from scratch. It is fairly close to creating a custom graphics function of your own (anybody would like to try?)

> windows(width=6, height=3, pointsize=10)

> plot(ci.simple.A[,"Estimate"], y <− 3:1,

+ xlab="regular-thick", xlim=range(pretty(ci.simple.A)),

+ yaxt="n", ylab="", ylim=c(y[3]-0.5,y[1]+0.5),

+ sub="Soybean Example",

+ type="n") # no points plotted, setup coordinates only

> abline(h=y, col="grey80") # horizontal grid lines

> axis(2, at=y, labels=rownames(ci.simple.A), las=1) # y axis

> segments(ci.simple.A[, "lwr"], y,

+ ci.simple.A[, "upr"], y,

+ lwd=2)

> points(ci.simple.A[, "Estimate"], y,

+ pch=21,

+ cex=2,

+ bg="green3")

> points(ci.simple.A[, "lwr"], y, pch="(")

> points(ci.simple.A[, "upr"], y, pch=")")

> abline(v=0, lty="dotted", col="darkorange4")

> title(main="95% Protected LSD Confidence Interval",

+ line=1.5)

> mtext("simple A (density) effects", side=3, at=0.5, adj=1)

> dev.off()

null device 
          1

(soybean03.png here)
To get plots for simple B effects, do the following changes to the codes above:
  1. use the command below following windows command
    par(mar=c(5,11,3,1)+.1)
    That is, leave enougn left margin (here 11 lines plus small offset 0.1). You need to experiment it. Or change the rownames of ci.simple.B to shorter ones (something like (L-M)wR and (M-H)wT).
  2. Change the assignment of y coordinates (within plot command) to
    y <− 6:1
  3. Change the assignment of ylim (within plot command) to
    ylim=c(y[6]-0.5,y[1]+0.5)
  4. Change the assignment of xlab (within plot command) to
    xlab="simple B effect"
  5. Change each occurrence of ci.simple.A to ci.simple.B.
  6. Omit mtext command.