Two-Way ANOVA Model, #2 (addendum)

Soybean Example (continued)

This addendum demonstrates the use of multcomp package to perform simple effects analyses, and residual analysis plots.

> # To make the names of the combined treatment levels

> # easier to use, the following is done.

> soybean$comb <− with(soybean,

+   factor(paste(substr(density,start=1,stop=1),

+                substr(fertilizer,1,1),

+                sep=""

+               )

+         )

+   )

> soybean

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

>

> # The ANOVA result below is for testing cell means.

> # The error term is the same as that of fitting

> # the full model with the two factors.

> soybean.aov <− aov(yield ~ comb, data=soybean)

> summary(soybean.aov)

            Df Sum Sq Mean Sq F value    Pr(>F)    
comb         5 638.26  127.65   91.79 3.638e-12 ***
Residuals   18  25.03    1.39                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

> summary(aov(yield~density*fertilizer,data=soybean))

                   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

>

> # Since the interaction is highly significant, we do simple

> # effects analyses by making use of the multcomp package.

> require(multcomp)

Loading required package: multcomp
Loading required package: mvtnorm
[1] TRUE

>

> ### Simple A effects analyses.

> # Use symbolic specification for linear function hypothesis.

> # For instance, the test H0rl−μtl=0 is

> # specified by "rl-tl=0".

> simple.A <− glht(soybean.aov,

+   linfct=mcp(comb=c("rl-tl=0","rm-tm=0","rh-th=0")))

>

> # After obtaining simple A effects object, do summary

> # for hypothesis tests, confint for confidence

> # intervals, and plot for confidence interval plot.

>

> ## Inferences without multiplicity adjustment.

> summary(simple.A, test=adjusted("none"))

         Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: User-defined Contrasts

Fit: aov(formula = yield ~ comb, data = soybean)

Linear Hypotheses:
             Estimate Std. Error t value  p value    
rl - tl == 0   0.5250     0.8339   0.630    0.537    
rm - tm == 0  10.3750     0.8339  12.442 2.81e-10 ***
rh - th == 0   1.5250     0.8339   1.829    0.084 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
(Adjusted p values reported -- none method)

> # Use t.025,18 for protected LSD (i.e., no adjustment)

> confint(simple.A, calpha=qt(.025,18,lower=F))

         Simultaneous Confidence Intervals for General Linear Hypotheses

Multiple Comparisons of Means: User-defined Contrasts


Fit: aov(formula = yield ~ comb, data = soybean)

Estimated Quantile = 2.1009

Linear Hypotheses:
             Estimate lwr     upr    
rl - tl == 0  0.5250  -1.2269  2.2769
rm - tm == 0 10.3750   8.6231 12.1269
rh - th == 0  1.5250  -0.2269  3.2769

Error in if (attr(x, "type") == "adjusted") { : 
        argument is of length zero

> # The last two error message, I think, is a bug. It can

> # be 'tricked' to avoid producing such error message (an

> # example is given later).

>

> ## Inferences with multiplicity adjustment (i.e., simultaneous

> # inferences).

> summary(simple.A) # hypothesis tests for simple A effects

 Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: User-defined Contrasts

Fit: aov(formula = yield ~ comb, data = soybean)

Linear Hypotheses:
             Estimate Std. Error t value p value    
rl - tl == 0   0.5250     0.8339   0.630   0.894    
rm - tm == 0  10.3750     0.8339  12.442  <1e-04 ***
rh - th == 0   1.5250     0.8339   1.829   0.224    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
(Adjusted p values reported)

> confint(simple.A) # confidence intervals

 Simultaneous Confidence Intervals for General Linear Hypotheses

Multiple Comparisons of Means: User-defined Contrasts

Fit: aov(formula = yield ~ comb, data = soybean)

Estimated Quantile = 2.6187

Linear Hypotheses:
             Estimate lwr     upr    
rl - tl == 0  0.5250  -1.6587  2.7087
rm - tm == 0 10.3750   8.1913 12.5587
rh - th == 0  1.5250  -0.6587  3.7087

95% family-wise confidence level

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

> plot(simple.A, sub="Soybean Example")

> title(main="(adjusted quantile=2.6187)", line=0.5,cex.main=1)

> dev.off()

null device 
          1

(soybean04.png here)

>

> ### Inferences on simple B effects

> simple.B <− glht(soybean.aov,

+   linfct=mcp(comb=c("rl-rm=0","rl-rh=0","rm-rh=0",

+                     "tl-tm=0","tl-th=0","tm-th=0")

+             )

+   )

> summary(simple.B, test=adjusted("none"))

 Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: User-defined Contrasts

Fit: aov(formula = yield ~ comb, data = soybean)

Linear Hypotheses:
             Estimate Std. Error t value  p value    
rl - rm == 0 -10.5750     0.8339 -12.682 2.06e-10 ***
rl - rh == 0 -10.7000     0.8339 -12.832 1.70e-10 ***
rm - rh == 0  -0.1250     0.8339  -0.150    0.883    
tl - tm == 0  -0.7250     0.8339  -0.869    0.396    
tl - th == 0  -9.7000     0.8339 -11.632 8.31e-10 ***
tm - th == 0  -8.9750     0.8339 -10.763 2.85e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
(Adjusted p values reported -- none method)

>

> # Now, here is the trick to 'fool' confint

> # The print command is not really needed, it

> # is used here only to show the quantile for individual

> # inferences.

> print( crit <− qt(0.025, 18, lower=FALSE) )

[1] 2.100922

> # The trick! follows.

> attr(crit, "type") <− "adjusted"

> confint(simple.B, calpha=crit)

 Simultaneous Confidence Intervals for General Linear Hypotheses

Multiple Comparisons of Means: User-defined Contrasts

Fit: aov(formula = yield ~ comb, data = soybean)

Estimated Quantile = 2.1009

Linear Hypotheses:
             Estimate lwr      upr     
rl - rm == 0 -10.5750 -12.3269  -8.8231
rl - rh == 0 -10.7000 -12.4519  -8.9481
rm - rh == 0  -0.1250  -1.8769   1.6269
tl - tm == 0  -0.7250  -2.4769   1.0269
tl - th == 0  -9.7000 -11.4519  -7.9481
tm - th == 0  -8.9750 -10.7269  -7.2231

95% family-wise confidence level

>

> # Now, do confidence interval plot (non-adjusted). It

> # makes a call to confint explicitly with user-specified

> # critical value.

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

> plot(confint(simple.B, calpha=crit), sub="Soybean Example")

> title(line=0.5, cex.main=1,

+   main="Simple B Effects (protected LSD)")

> # The next is added (for fun!) with a plotmath (LaTeX) style

> # text string on the plot. Use help(plotmath)

> # for 'how to'.

> title(line=3, adj=0.1,

+   sub=bquote("critical="*t[.025,18]==.(round(qt(.025,18,low=F),4))))

> dev.off()

null device 
          1

(soybean05.png here)

>

> ## Simultaneous inferences for simple B effects

> summary(simple.B)

 Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: User-defined Contrasts

Fit: aov(formula = yield ~ comb, data = soybean)

Linear Hypotheses:
             Estimate Std. Error t value p value    
rl - rm == 0 -10.5750     0.8339 -12.682  <1e-04 ***
rl - rh == 0 -10.7000     0.8339 -12.832  <1e-04 ***
rm - rh == 0  -0.1250     0.8339  -0.150    1.00    
tl - tm == 0  -0.7250     0.8339  -0.869    0.88    
tl - th == 0  -9.7000     0.8339 -11.632  <1e-04 ***
tm - th == 0  -8.9750     0.8339 -10.763  <1e-04 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
(Adjusted p values reported)

> confint(simple.B)

 Simultaneous Confidence Intervals for General Linear Hypotheses

Multiple Comparisons of Means: User-defined Contrasts

Fit: aov(formula = yield ~ comb, data = soybean)

Estimated Quantile = 2.875

Linear Hypotheses:
             Estimate lwr      upr     
rl - rm == 0 -10.5750 -12.9724  -8.1776
rl - rh == 0 -10.7000 -13.0974  -8.3026
rm - rh == 0  -0.1250  -2.5224   2.2724
tl - tm == 0  -0.7250  -3.1224   1.6724
tl - th == 0  -9.7000 -12.0974  -7.3026
tm - th == 0  -8.9750 -11.3724  -6.5776

95% family-wise confidence level

>

> # The session below (not really needed to do this way)

> # gives a plotmath (LaTeX) style y-axis tick mark labels

> # (for the simple B effects). The trick is to turn off

> # the plotting of y axis with the graphical parameter

> # yaxt set to "n" (i.e., no axis). And

> # then use the low-level graphics function axis

> # to plot the custom labels. The confidence intervals

> # are produced top down at y-axis values (6 down to 1

> # in this case) for the 6 intervals. What makes

> # the math-style labels is the call to expression.

> # In the call below, it produces a vector of length 6 of

> # mode expression suitable to produce math-style

> # symbols.

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

> plot(simple.B, sub="Soybean Example", yaxt="n")

> title(main="(adjusted quantile=2.875)",line=0.5,cex.main=1)

> labs <− expression(mu[rl]-mu[rm], mu[rl]-mu[rh], mu[rm]-mu[rh],

+                    mu[tl]-mu[tm], mu[tl]-mu[th], mu[tm]-mu[th])

> axis(2, 6:1, labels=labs, las=1)

> # Again, R uses the four margins consistently, 1=bottom, 2=left,

> # 3=top, and 4=right. So the call to axis above

> # produced y-axis on the left (and yes, you can place y-axis

> # on the right as long as you leave enough room for R on the

> # right margin!) On some occasions, you want to produce two

> # axes for a dimension (say y). Try this:

> # par(mar=c(5,4,4,5)+.1)

> # qqnorm(rnorm(50), datax=TRUE, yaxt="n")

> # axis(2, at=(x <− seq(-2,2,.5)), las=1)

> # axis(4, at=x, las=1, labels=round(pnorm(x),3))

> # mtext("Normal Probability", side=4, line=4)

> dev.off()

null device 
          1

(soybean06.png here)

>

> ### Residual Analysis

> # The following produces two residual plots: residuals versus

> # fitted values, and the normal quantile plots of the residuals.

> # Please consult with the help documentation using

> # help(plot.lm) for details.

> # The graphics parameter oma used below ensures

> # enough room for overall title (placed in outer margin).

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

> par(mfrow=c(1,2), oma=c(0,0,2,0)+.1)

> plot(aov(yield~density*fertilizer,data=soybean),

+      which=1:2, add.smooth=FALSE)

> dev.off()

null device 
          1

(soybean07.png here)

> # To get individual plots, for instances,

> # qqnorm(resid(soybean.aov))

> # for residuals Q-Q plot

> # plot(fitted(soybean.aov), residuals(soybean.aov))

> # for plot of residuals versus fitted values.

> # Note that you can use resid for short

> # (for residuals).