options(contrasts=c("contr.sum","contr.poly"))

3-Way ANOVA

load(url("http://pnb.mcmaster.ca/bennett/psy710/datasets/aov-3way-data.rda") )
summary(myDf09)
##      score        A       B       C          groupID  
##  Min.   :17.23   a0:60   b0:60   c0:60   a0.b0.c0:15  
##  1st Qu.:19.06   a1:60   b1:60   c1:60   a1.b0.c0:15  
##  Median :19.79                           a0.b1.c0:15  
##  Mean   :19.82                           a1.b1.c0:15  
##  3rd Qu.:20.48                           a0.b0.c1:15  
##  Max.   :22.21                           a1.b0.c1:15  
##                                          (Other) :30
aov.01 <- aov(score~A*B*C,data=myDf09)
anova(aov.01)
## Analysis of Variance Table
## 
## Response: score
##            Df  Sum Sq Mean Sq F value    Pr(>F)    
## A           1   9.395  9.3952  9.3291 0.0028187 ** 
## B           1   4.935  4.9352  4.9005 0.0288756 *  
## C           1   0.261  0.2612  0.2594 0.6115684    
## A:B         1  14.578 14.5778 14.4752 0.0002319 ***
## A:C         1   0.014  0.0143  0.0142 0.9053737    
## B:C         1   0.733  0.7333  0.7281 0.3953070    
## A:B:C       1   0.000  0.0002  0.0002 0.9877120    
## Residuals 112 112.794  1.0071                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

analyze 2-way interactions

with(myDf09,interaction.plot(x.factor=A,trace.factor=B,response=score,type="b",ylab="score",pch=c(19,21)))
mtext(side=3,text="AxB interaction (ignoring C)",line=0.5,cex=1.25)

# simple main effect of A at B0 (ignoring C)
curDf <- subset(myDf09,B=="b0")
a.b0 <- aov(score~A,data=curDf)
summary(a.b0)
##             Df Sum Sq Mean Sq F value Pr(>F)
## A            1   0.28  0.2834   0.258  0.613
## Residuals   58  63.61  1.0967
# simple main effect of A at B1 (ignoring C)
curDf <- subset(myDf09,B=="b1")
a.b1 <- aov(score~A,data=curDf)
summary(a.b1)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## A            1  23.69  23.690   27.38 2.41e-06 ***
## Residuals   58  50.19   0.865                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

analyze 3-way interaction

summary(myDf10)
##      score         A       B       C          groupID  
##  Min.   : 5.562   a1:60   b1:60   c1:60   a1.b1.c1:15  
##  1st Qu.:15.407   a2:60   b2:60   c2:60   a2.b1.c1:15  
##  Median :19.652                           a1.b2.c1:15  
##  Mean   :19.841                           a2.b2.c1:15  
##  3rd Qu.:24.495                           a1.b1.c2:15  
##  Max.   :30.584                           a2.b1.c2:15  
##                                           (Other) :30
aov.100 <- aov(score~A*B*C,data=myDf10)
anova(aov.100)
## Analysis of Variance Table
## 
## Response: score
##            Df  Sum Sq Mean Sq F value    Pr(>F)    
## A           1  326.32  326.32 12.9611 0.0004753 ***
## B           1    8.36    8.36  0.3322 0.5655501    
## C           1  488.06  488.06 19.3852 2.453e-05 ***
## A:B         1  118.26  118.26  4.6970 0.0323332 *  
## A:C         1  171.49  171.49  6.8111 0.0102971 *  
## B:C         1   88.57   88.57  3.5181 0.0633076 .  
## A:B:C       1  185.38  185.38  7.3632 0.0077097 ** 
## Residuals 112 2819.84   25.18                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
MS.resid.full <- 25.18
df.resid.full <- 112
par(mfrow=c(1,2),cex=1.25)
with(subset(myDf10,C=="c1"),
     interaction.plot(x.factor=A,
                      trace.factor=B,
                      ylab="score",
                      response=score,
                      type="b",
                      pch=c(19,21),
                      main=bquote(C[1]),ylim=c(10,30)))

with(subset(myDf10,C=="c2"),
     interaction.plot(x.factor=A,
                      trace.factor=B,
                      ylab="score",
                      response=score,
                      type="b",
                      pch=c(19,21),
                      main=bquote(C[2]),
                      ylim=c(10,30)))

A x B interaction at c1

# AxB interaction at c1:
aov.101.c1 <- aov(score~A*B,data=subset(myDf10,C=="c1"))
anova(aov.101.c1)
## Analysis of Variance Table
## 
## Response: score
##           Df  Sum Sq Mean Sq F value    Pr(>F)    
## A          1  485.46  485.46 17.3717 0.0001076 ***
## B          1   75.68   75.68  2.7083 0.1054294    
## A:B        1    3.76    3.76  0.1344 0.7152685    
## Residuals 56 1564.95   27.95                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# recalculate F and p values:
# simple AxB interaction at c1
F.AxB.c1 <- 3.76/MS.resid.full
1-pf(F.AxB.c1,df1=1,df2=df.resid.full)
## [1] 0.6999145
# simple main effect of A at c1
F.A.c1 <- 485.46/MS.resid.full
1-pf(F.A.c1,df1=1,df2=df.resid.full)
## [1] 2.571595e-05
# simple main effect of B at c1
F.B.c1 <- 75.68/MS.resid.full
1-pf(F.B.c1,df1=1,df2=df.resid.full)
## [1] 0.08573115

A x B interaction at c2

# AxB interaction at c2:
aov.101.c2 <- aov(score~A*B,data=subset(myDf10,C=="c2"))
anova(aov.101.c2)
## Analysis of Variance Table
## 
## Response: score
##           Df  Sum Sq Mean Sq F value    Pr(>F)    
## A          1   12.35  12.347  0.5510 0.4610248    
## B          1   21.25  21.252  0.9484 0.3343153    
## A:B        1  299.89 299.885 13.3825 0.0005628 ***
## Residuals 56 1254.89  22.409                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# recalculate AxB F and p values:
F.AxB.c2 <- 299.885/MS.resid.full
1-pf(F.AxB.c2,df1=1,df2=df.resid.full)
## [1] 0.0007885303

simple simple main effects of A

# evaluate simple simple main effect of A at C2 & B1:
aov.101.c2.b1 <- aov(score~A,data=subset(myDf10,B=="b1"&C=="c2"))
anova(aov.101.c2.b1)
## Analysis of Variance Table
## 
## Response: score
##           Df Sum Sq Mean Sq F value Pr(>F)  
## A          1  95.27  95.267   3.695 0.0648 .
## Residuals 28 721.91  25.783                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(F.A.c2.b1 <- 95.267/MS.resid.full)
## [1] 3.783439
1-pf(F.A.c2.b1,df1=1,df2=df.resid.full)
## [1] 0.05426882
# evaluate simple simple main effect of A at C2 & B2:
aov.101.c2.b2 <- aov(score~A,data=subset(myDf10,B=="b2"&C=="c2"))
anova(aov.101.c2.b2)
## Analysis of Variance Table
## 
## Response: score
##           Df Sum Sq Mean Sq F value   Pr(>F)   
## A          1 216.96 216.964  11.398 0.002172 **
## Residuals 28 532.98  19.035                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(F.A.c2.b2 <- 216.964/MS.resid.full)
## [1] 8.616521
1-pf(F.A.c2.b2,df1=1,df2=df.resid.full)
## [1] 0.004042999

simple simple main effects of B

# evaluate simple simple main effect of B at C2 & A1:
aov.101.c2.a1 <- aov(score~B,data=subset(myDf10,A=="a1"&C=="c2"))
anova(aov.101.c2.a1)
## Analysis of Variance Table
## 
## Response: score
##           Df Sum Sq Mean Sq F value   Pr(>F)   
## B          1 240.40 240.402   12.18 0.001618 **
## Residuals 28 552.66  19.738                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(F.B.c2.a1 <- 240.402/MS.resid.full)
## [1] 9.547339
1-pf(F.B.c2.a1,df1=1,df2=df.resid.full)
## [1] 0.00252606
# evaluate simple simple main effect of B at C2 & A2:
aov.101.c2.a2 <- aov(score~B,data=subset(myDf10,A=="a2"&C=="c2"))
anova(aov.101.c2.a2)
## Analysis of Variance Table
## 
## Response: score
##           Df Sum Sq Mean Sq F value  Pr(>F)  
## B          1  80.74  80.736  3.2192 0.08359 .
## Residuals 28 702.23  25.080                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(F.B.c2.a2 <- 80.736/MS.resid.full)
## [1] 3.206354
1-pf(F.B.c2.a2,df1=1,df2=df.resid.full)
## [1] 0.07605422