Multi-Way ANOVA Model

Bottler Example

> y <- scan()

1:   7  9   10 11   15 14

7:   9 10   12 11   17 16

13:  9 10   12 13   17 19

19: 11 11   16 15   20 21

25:

Read 24 items

> bottler <- data.frame(volume=y,

+   expand.grid(replicate=1:2, carbon=seq(from=10,to=14,by=2),

+     speed=c(200,250), pressure=c(25,30)))

> rm(y) # remove y

> # Command below shows that the data entry is correctly done.

> xtabs(volume~carbon+speed+replicate+pressure, data=bottler)

, , replicate = 1, pressure = 25

      speed
carbon 200 250
    10   7   9
    12  10  12
    14  15  17

, , replicate = 2, pressure = 25

      speed
carbon 200 250
    10   9  10
    12  11  11
    14  14  16

, , replicate = 1, pressure = 30

      speed
carbon 200 250
    10   9  11
    12  12  16
    14  17  20

, , replicate = 2, pressure = 30

      speed
carbon 200 250
    10  10  11
    12  13  15
    14  19  21

> #

> # For the example below, we will avoid long names by creating another

> # data frame object with shorter names and force the

> # quantitative predictors to become factors.

> colnames(bottler)

[1] "volume"      "replicate"   "carbon" "speed"      
[5] "pressure"   

> bottler2 <- bottler[,c(1,3,5,4)]

> colnames(bottler2)

[1] "volume"      "carbon" "pressure"    "speed"      

> colnames(bottler2)[-1] <- LETTERS[1:3]

> colnames(bottler2)

[1] "volume" "A"      "B"      "C"     

> for(i in 2:4)bottler2[[i]] <- as.factor(bottler2[[i]])

> #

> # Now fit full model.

> bottler.lm <- lm(volume ~ .^3, data=bottler2)

> anova(bottler.lm)

Analysis of Variance Table

Response: volume
          Df  Sum Sq Mean Sq  F value    Pr(>F)    
A          2 252.750 126.375 178.4118 1.186e-09 ***
B          1  45.375  45.375  64.0588 3.742e-06 ***
C          1  22.042  22.042  31.1176 0.0001202 ***
A:B        2   5.250   2.625   3.7059 0.0558081 .  
A:C        2   0.583   0.292   0.4118 0.6714939    
B:C        1   1.042   1.042   1.4706 0.2485867    
A:B:C      2   1.083   0.542   0.7647 0.4868711    
Residuals 12   8.500   0.708                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

> # It appears that a reduced model on all main effects and

> # A:B (carbon × pressure) can be fit. Since all

> # factors are quantitative, orthogonal polynomial contrasts

> # are considered.

> bottler2.lm <- lm(volume ~ A*B + C,

+   contrast=list(A="contr.poly", B="contr.poly", C="contr.poly"),

+   data=bottler2)

> anova(bottler2.lm)

Analysis of Variance Table

Response: volume
          Df  Sum Sq Mean Sq  F value    Pr(>F)    
A          2 252.750 126.375 191.6766 2.178e-12 ***
B          1  45.375  45.375  68.8216 2.218e-07 ***
C          1  22.042  22.042  33.4312 2.210e-05 ***
A:B        2   5.250   2.625   3.9814   0.03818 *  
Residuals 17  11.208   0.659                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

> summary(bottler2.lm)

Call:
lm(formula = volume ~ A * B + C, data = bottler2, contrasts =
   list(A = "contr.poly", B = "contr.poly", C = "contr.poly"))

Residuals:
     Min       1Q   Median       3Q      Max 
-1.29167 -0.47917 -0.04167  0.58333  1.20833 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  13.1250     0.1657  79.188  < 2e-16 ***
A.L           5.5685     0.2871  19.397 4.93e-13 ***
A.Q           0.7655     0.2871   2.666   0.0163 *  
B.L           1.9445     0.2344   8.296 2.22e-07 ***
C.L           1.3553     0.2344   5.782 2.21e-05 ***
A.L:B.L       1.1250     0.4060   2.771   0.0131 *  
A.Q:B.L      -0.2165     0.4060  -0.533   0.6007    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Residual standard error: 0.812 on 17 degrees of freedom
Multiple R-Squared: 0.9667,     Adjusted R-squared: 0.955 
F-statistic: 82.26 on 6 and 17 DF,  p-value: 1.294e-11

> #

> # It seems to be reasonable to fit the following polynomial

> # regression model by treating the quantitative predictors

> # as is. That is, the model

> # volume = const + carbon + carbon^2 + pressure + speed

> # + carbon*pressure

> # is contemplated.

> # The model formula below is equivalent to

> # volume ~ carbon + I(carbon^2) + pressure + speed

> # + I(carbon*pressure)

> bottler3.lm <- lm(data=bottler,

+   formula=volume~poly(carbon,2,raw=T)+pressure+speed

+                  +carbon:pressure)

> summary(bottler3.lm)

Call:
lm(formula = volume ~ poly(carbon, 2, raw = T) + pressure + 
    speed + carbon:pressure, data = bottler)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.3542 -0.4167 -0.1250  0.6146  1.2708 

Coefficients:
                           Estimate Std. Error t value Pr(>|t|)
(Intercept)               36.000000  18.092541   1.990   0.0620
poly(carbon, 2, raw = T)1 -6.750000   2.341001  -2.883   0.0099
poly(carbon, 2, raw = T)2  0.234375   0.086135   2.721   0.0140
pressure                  -0.800000   0.481806  -1.660   0.1141
speed                      0.038333   0.006497   5.900 1.38e-05
pressure:carbon            0.112500   0.039784   2.828   0.0112
                                  
(Intercept)               .  
poly(carbon, 2, raw = T)1 ** 
poly(carbon, 2, raw = T)2 *  
pressure
speed                     ***
pressure:carbon           *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Residual standard error: 0.7957 on 18 degrees of freedom
Multiple R-Squared: 0.9661,     Adjusted R-squared: 0.9567 
F-statistic: 102.7 on 5 and 18 DF,  p-value: 1.375e-12

> # Residual analysis on this fit should be carried out so that

> # the model can be used for predication.

> #

> # The session below demonstrate a view of a 3-way interaction

> # plot (in fact, a two-way interaction plot was produced).

> cols <- rainbow(7)

> windows(width=5, height=4, pointsize=10)

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

> with(bottler,

+      {trace <- paste("(", carbon, ",", pressure, ")", sep="")

+       interaction.plot(speed, trace, volume, type="b",

+             trace.label="(carbon,pressure)",

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

+      }

+     )

> title(main="Three-Way Interaction Plot\nan example view")

> title(sub=" Bottler Example", adj=0)

> dev.off()

null device 
          1

(bottler01.png here)

The profiles appear to be strikingly parallel indicating no 3-way interaction.