1 Load & Inspect Data

This lab illustrates methods for analyzing data collected in experiments using a between-subjects, factorial design. The first example analyzed data from a 2(A) x 2(B) factorial design that is stored in the data frame df0 Load the data, inspect the data frame df0, and examine the design using xtabs. The data in each condition are plotted in Figure 1. Some R commands that are useful for summarizing data from factorial designs can be found here.

First, initialize R, load the data, and then use summary and xtabs to inspect the data and confirm that the design is balanced with 10 participants per condition.

options(contrasts=c('contr.sum','contr.poly'))
load(url("http://pnb.mcmaster.ca/bennett/psy710/datasets/rlab5-df0-2x2.rda"))
summary(df0) # A&B are factors; Y is numeric
#   A       B            Y         
#  a1:20   b1:20   Min.   : 60.92  
#  a2:20   b2:20   1st Qu.: 78.65  
#                  Median : 99.17  
#                  Mean   :101.07  
#                  3rd Qu.:123.15  
#                  Max.   :175.35
xtabs(~A+B,data=df0) # balanced; n=10 per condition
#     B
# A    b1 b2
#   a1 10 10
#   a2 10 10

Next, we are going to plot the data in each condition. Note how I use the formula Y~A*B in boxplot.

op <- par(no.readonly = T)
par(mar=c(5,5,5,2)+0.1, cex.axis=1.25, cex.lab=1.5,cex.main=1.5)
par(mfrow=c(1,1))
boxplot(Y~A*B,df0,main="A x B Cell Means")
Figure 1. Data from the 2 x 2 design.

Figure 1. Data from the 2 x 2 design.

2 Perform ANOVA

Now we perform a 2(A) x 2(B) ANOVA. The ANOVA table is shown below. The main effect of A is significant, but the main effect of B is not. The AxB interaction also is significant.

# 2 x 2 ANOVA:
options(contrasts=c("contr.sum","contr.poly")) # IMPORTANT!
aov.01 <- aov(Y~A*B,data=df0)
# aov.01 <- aov(Y ~ 1 + A + B + A:B,data=df0)
summary(aov.01)
#             Df Sum Sq Mean Sq F value  Pr(>F)   
# A            1   5096    5096  10.361 0.00273 **
# B            1    284     284   0.578 0.45207   
# A:B          1   2803    2803   5.699 0.02235 * 
# Residuals   36  17705     492                   
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

2.1 Questions

  1. What null hypotheses are being evaluated by each line in the ANOVA table?
# The null hypotheses for the main effects (A and B) are that the marginal means are equal. The null hypothesis for the interaction si that the effects of A are the same at all levels of B or, equivalently, that the effects of B are the same at all levels of A.
  1. Create a graph that illustrates the main effect of A.
# The main effects are related to the MARGINAL distributions
# of A ignoring B and of B ignoring A.
# The  main effect of A is plotted in Figure 2.

par(op)
par(mar=c(5,5,5,2)+0.1, cex.axis=1.125, cex.lab=1.5,cex.main=1.5)
par(mfrow=c(1,1))
boxplot(Y~A,df0,main="Main Effect of A") # ignore B
Figure 2. Marginal distribution of A.

Figure 2. Marginal distribution of A.

3 Check Assumptions

3.1 Normality

ANOVA assumes that the residuals are distributed normally. We can check this assumption with shapiro.test.

shapiro.test(residuals(aov.01)) # check normality
# 
#   Shapiro-Wilk normality test
# 
# data:  residuals(aov.01)
# W = 0.98818, p-value = 0.9454

We also can evaluate the normality assumption with a qq plot.

3.2 Constant Variance

ANOVA also assumes that the residuals have constant variance across conditions. We can check this assumption with bartlett.test. However, that command does not accept a factorial formula of the form Y~A+B+A:B.

# bartlett.test(Y~A*B,data=df0) # DOES NOT WORK!"

Because bartlett.test does not work with factorial model formulae, we have to use interaction create a new one-way grouping variable that has one level for each condition, and then re-submit the data to bartlett.test. The test is not significant, so we do not have any reason to reject the null hypothesis that the variance is constant across groups. (N.B. Inspection of Figure 1 leads to the same conclusion).

df0$group <- interaction(df0$A,df0$B) # create new factor
levels(df0$group)
# [1] "a1.b1" "a2.b1" "a1.b2" "a2.b2"
bartlett.test(Y~group,data=df0)
# 
#   Bartlett test of homogeneity of variances
# 
# data:  Y by group
# Bartlett's K-squared = 6.2391, df = 3, p-value = 0.1005

4 Analyzing a Main Effect

Our ANOVA found a main effect of A which leads us to reject the null hypothesis of no differences among the marginal means of A. To learn more about the pattern of differences among marginal means, we often perform linear contrasts or pairwise tests on the marginal means. In this example, such tests are not necessary because A only has two levels, and therefore main effect of A is unambiguous (i.e., the two levels differ significantly). Also, an analysis of the main effect might be premature and/or misleading because we found a significant AxB interaction, which implies that the main effect of A varies across the levels of B. In situations where the interaction is significant but we are interested in variation among levels of A, we normally would decompose the interaction into simple main effects and then look for patterns in the simple main effects. Despite these concerns, I perform the tests anyway for illustrative purposes.

4.1 Tukey HSD

Here I analyze main effect of A with pairwise Tukey tests. Remember: it is not necessary to obtain a significant omnibus \(F\) test before performing the Tukey pairwise test.

# perform Tukey HSD test on levels of factor A
TukeyHSD(aov.01,which="A") 
#   Tukey multiple comparisons of means
#     95% family-wise confidence level
# 
# Fit: aov(formula = Y ~ A * B, data = df0)
# 
# $A
#           diff      lwr      upr     p adj
# a2-a1 22.57333 8.350657 36.79601 0.0027251
# with emmeans
library(emmeans)
aov.emm.A <- emmeans(aov.01,specs="A")
contrast(aov.emm.A,method="pairwise",adjust="Tukey")
#  contrast estimate   SE df t.ratio p.value
#  a1 - a2     -22.6 7.01 36  -3.219  0.0027
# 
# Results are averaged over the levels of: B

4.2 Linear Contrasts

A significant main effect also can be analyzed with a linear contrast.

levels(df0$A)
# [1] "a1" "a2"
myC <- c(-1,1)
contrasts(df0$A) <- cbind(myC)
aov.01c <- aov(Y~A*B,data=df0) # note formula
summary(aov.01c,split=list(A=list(myC=1)))
#             Df Sum Sq Mean Sq F value  Pr(>F)   
# A            1   5096    5096  10.361 0.00273 **
#   A: myC     1   5096    5096  10.361 0.00273 **
# B            1    284     284   0.578 0.45207   
# A:B          1   2803    2803   5.699 0.02235 * 
#   A:B: myC   1   2803    2803   5.699 0.02235 * 
# Residuals   36  17705     492                   
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The 2nd line in the anova table indicates that my contrast was significant, so I reject the null hypothesis that the marginal means of a1 and a2 are equal. Furthermore, the 5th line in the anova table indicates that my contrast differed between the levels of B. Inspection of the anova table shows that the results for my contrast are identical to the results for the main effect of A, and the results for my contast x B interaction are the same as the results for the AxB interaction. The reason the results are identical is that factors A and B have two levels and one degree-of-freedom, and therefore my contrasts are equivalent to the omnibus F tests.

The next block of code illustrates how to perform the contrast with emmeans. Unlike the aov command, emmeans does not automatically test the hypothesis that the contrast varies across the levels of B.

library(emmeans)
myC <- c(-1,1)
aov.emm.A <- emmeans(aov.01,specs="A")
contrast(aov.emm.A,method=list(myC),adjust="none")
#  contrast estimate   SE df t.ratio p.value
#  c(-1, 1)     22.6 7.01 36   3.219  0.0027
# 
# Results are averaged over the levels of: B

5 Plotting Interactions

The following code illustrates a simple way to visualize the A x B interaction.

# see: https://rpubs.com/tf_peterson/interactionplotDemo
par(mar=c(5,5,5,2)+0.1, 
    cex.axis=1.125, 
    cex.lab=1.35,
    cex.main=1.5)
interaction.plot(x.factor=df0$A,
                 trace.factor=df0$B,
                 response=df0$Y,
                 type="b",
                 ylab="Y",
                 xlab="A",
                 trace.label="B",
                 main="AxB Interaction")
Figure 3. A x B interaction.

Figure 3. A x B interaction.

6 Simple Main Eeffects

There are a variety of ways to analyze an interaction. First, we will analyze the simple main effect of A at each level of B. The first block of code uses aov to perform a 1-way ANOVA on subsets of the data and then recalculates the \(F\) and \(p\) values using MS error and degrees-of-freedom from the main ANOVA. The simple main effect of \(A\) at \(b_1\) is significant, but the simple main effect of \(A\) at \(b_2\) is not.

# simple main effect using aov
# summary(aov.01)
MS.resid <- 492
df.resid <- 36
levels(df0$B)
# [1] "b1" "b2"
aov.A_at_b1 <- aov(Y~A,data=subset(df0,B=="b1"))
aov.A_at_b2 <- aov(Y~A,data=subset(df0,B=="b2"))
summary(aov.A_at_b1)
#             Df Sum Sq Mean Sq F value  Pr(>F)    
# A            1   7728    7728   19.43 0.00034 ***
# Residuals   18   7161     398                    
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(aov.A_at_b2)
#             Df Sum Sq Mean Sq F value Pr(>F)
# A            1    170   170.1    0.29  0.597
# Residuals   18  10544   585.8
# recalculate F & p values:
( F.recalc <- c(7728,170.1)/MS.resid ) # recalculate F's
# [1] 15.7073171  0.3457317
( p.recalc <- 1-pf(F.recalc,df1=1,df2=df.resid)) # recalculate p's
# [1] 0.0003355768 0.5602100635
( p.recalc < .05) # 1st simple main effect is significant
# [1]  TRUE FALSE

The same analysis can be performed with emmeans. Note that the \(F\) and \(p\) values correspond to the recalculated values above. You can read more about how emmeans analyzes interactions here.

# simple main effect using emmeans
# https://cran.r-project.org/web/packages/emmeans/vignettes/interactions.html
library(emmeans)
aov.emm <- emmeans(aov.01,specs="A",by="B")
#aov.emm <- emmeans(aov.01,specs=~A+B) # yields the same result
joint_tests(aov.emm,by="B")

6.1 Question

  1. Evaluate the simple main effect of B at each level of A.
# simple main effect using aov
# summary(aov.01)
MS.resid <- 492
df.resid <- 36
levels(df0$A)
# [1] "a1" "a2"
aov.B_at_a1 <- aov(Y~B,data=subset(df0,A=="a1"))
aov.B_at_a2 <- aov(Y~B,data=subset(df0,A=="a2"))
summary(aov.B_at_a1)
#             Df Sum Sq Mean Sq F value  Pr(>F)   
# B            1   2436  2435.9   8.415 0.00953 **
# Residuals   18   5211   289.5                   
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(aov.B_at_a2)
#             Df Sum Sq Mean Sq F value Pr(>F)
# B            1    651   650.9   0.938  0.346
# Residuals   18  12494   694.1
# recalculate F & p:
( F.recalc <- c(2435.9,650.9)/MS.resid )
# [1] 4.951016 1.322967
( p.recalc <- 1-pf(F.recalc,df1=1,df2=df.resid))
# [1] 0.03243052 0.25764364
( p.recalc < .05) # 1st simple main effect is significant
# [1]  TRUE FALSE


# library(emmeans)
aov.emm.B <- emmeans(aov.01,specs="B",by="A")
joint_tests(aov.emm,by="A")

7 Interaction Contrasts

This section uses data in the data frame biofeed. Load the data into R, inspect the variables, and perform a factorial ANOVA that evaluates the effects of biofeedback and drug on score. The interaction, illustrated in Figure 4, is significant

load(url("http://pnb.mcmaster.ca/bennett/psy710/datasets/biofeed.rda"))
summary(biofeed)
#  biofeedback    drug        score      
#  Bio :15     drugX:10   Min.   :160.0  
#  None:15     drugY:10   1st Qu.:181.5  
#              drugZ:10   Median :196.5  
#                         Mean   :195.0  
#                         3rd Qu.:206.0  
#                         Max.   :228.0
bio.aov.01 <- aov(score~drug*biofeedback,data=biofeed)
summary(bio.aov.01)
#                  Df Sum Sq Mean Sq F value  Pr(>F)   
# drug              2   1882   941.0   5.210 0.01321 * 
# biofeedback       1   1904  1904.0  10.542 0.00343 **
# drug:biofeedback  2   1248   624.0   3.455 0.04801 * 
# Residuals        24   4335   180.6                   
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
par(mar=c(5,5,5,2)+0.1, 
    cex.axis=1.125, 
    cex.lab=1.35,
    cex.main=1.5)
interaction.plot(x.factor=biofeed$drug,
                 trace.factor=biofeed$biofeedback,
                 response=biofeed$score,
                 type="b",
                 ylab="Score",
                 xlab="Drug",
                 trace.label="biofeedback",
                 main="Drug x Biofeedback")
Figure 4. Drug x Biofeedback Interaction.

Figure 4. Drug x Biofeedback Interaction.

In this section we analyze the interaction by evaluating an interaction contrast. Let us suppose that we we predicted that score would be lower in the drug X condition than the drug Y and Z conditions, and that this drug effect would be greater in the biofeedback-present condition. This prediction can be tested by evaluating a linear contrast that compares drug X to the average of the means in the drug Y and Z conditions, and then seeing if that contrast differs significantly between the two biofeedback conditions. We perform that contrast next.

# create contrast weights
levels(biofeed$drug)
# [1] "drugX" "drugY" "drugZ"
cDrug <- c(-1,1/2,1/2) # drug X vs (Y + Z)/2
levels(biofeed$biofeedback)
# [1] "Bio"  "None"
cBiofeed <- c(1,-1) # biofeedback vs none

# assign contrasts to factors:
contrasts(biofeed$drug) <- cbind(cDrug)
contrasts(biofeed$biofeedback) <- cbind(cBiofeed)

# perform ANOVA (again)
bio.aov.02 <- aov(score~drug*biofeedback,data=biofeed)
summary(bio.aov.02,split=list(drug=list(cDrug=1),biofeedback=list(cBiofeed=1)))
#                                    Df Sum Sq Mean Sq F value  Pr(>F)   
# drug                                2   1882   941.0   5.210 0.01321 * 
#   drug: cDrug                       1   1837  1837.1  10.171 0.00394 **
# biofeedback                         1   1904  1904.0  10.542 0.00343 **
#   biofeedback: cBiofeed             1   1904  1904.0  10.542 0.00343 **
# drug:biofeedback                    2   1248   624.0   3.455 0.04801 * 
#   drug:biofeedback: cDrug.cBiofeed  1    528   528.1   2.924 0.10019   
# Residuals                          24   4335   180.6                   
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The contrast under the main effect of drug is significant, which means that we reject the null hypothesis that the mean score in drugX equals the average of the mean scores in the drugY and drugZ conditions. The interaction contrast, listed under the drug:biofeedback interation, is not significant. So we do not reject the null hypothesis that the contrast between drugX and the average of the drugY and drugZ conditions is the same in the two biofeedback condition.

To perform the interaction contrast with emmeans we need to create a new factor that contains one level for each condition, and then apply appropriate contrast weights to the new factor. In other words, we need to perform the contrast as though we used a one-way design instead of a factorial design. First, we create the new factor and assign contrast weights.

# create new factor with 1 level per condition:
biofeed$condition <- with(biofeed,interaction(drug,biofeedback))
levels(biofeed$condition)
# [1] "drugX.Bio"  "drugY.Bio"  "drugZ.Bio"  "drugX.None" "drugY.None" "drugZ.None"
# create contrast weights
cDrug <- c(-1,1/2,1/2,-1,1/2,1/2) # X vs (Y + Z)/2
cBio <- c(1,1,1,-1,-1,-1) # Biofeedback vs None
cDrugXBio <- cDrug*cBio # interaction contrast
contrasts(biofeed$condition) <- cDrugXBio # assign contrast grouping factor
levels(biofeed$condition)
# [1] "drugX.Bio"  "drugY.Bio"  "drugZ.Bio"  "drugX.None" "drugY.None" "drugZ.None"
cDrugXBio # c(-1, 1/2, 1/2) **MINUS** c(-1, 1/2, 1/2)
# [1] -1.0  0.5  0.5  1.0 -0.5 -0.5

Note that the weights for the interaction contrast equal the weights for the comparison of drugX to drugY and drugZ in the biofeedback condition, c(-1, 1/2, 1/2), MINUS the drug contrast in the none condition, -1*c(-1, 1/2, 1/2) = c(1, -1/2, -1/2). Next we use those weights to perform the contrast in emmeans.

library(emmeans)
bio.aov.03 <- aov(score~condition,data=biofeed) # perform ANOVA
biofeed.emm <- emmeans(bio.aov.03,specs="condition") # create emmeans object
contrast(biofeed.emm,method=list(DxB=cDrugXBio),adjust="none") # evaluate contrast
#  contrast estimate   SE df t.ratio p.value
#  DxB          17.8 10.4 24   1.710  0.1002

Finally, we can also use aov to perform this contrast:

# perform ANOVA (again) note forumula!
#bio.aov.03 <- aov(score~condition,data=biofeed)
summary(bio.aov.03,split=list(condition=list(DrugXBio=1)))
#                       Df Sum Sq Mean Sq F value  Pr(>F)   
# condition              5   5034  1006.8   5.574 0.00152 **
#   condition: DrugXBio  1    528   528.1   2.924 0.10019   
# Residuals             24   4335   180.6                   
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

8 (2 X 3) factorial ANOVA

This section uses data in the data frame df1 contains data from a 2 (group) by 3 (task) factorial design. Load the data into R and inspect the variables.

load(url("http://pnb.mcmaster.ca/bennett/psy710/datasets/rlab5-df1.rda"))
sapply(df1,class)
#     group      task         Y 
#  "factor"  "factor" "numeric"
xtabs(~group+task,data=df1)
#          task
# group     t1 t2 t3
#   patient 10 10 10
#   control 10 10 10
summary(df1)
#      group    task          Y         
#  patient:30   t1:20   Min.   : 40.87  
#  control:30   t2:20   1st Qu.: 83.70  
#               t3:20   Median : 97.22  
#                       Mean   : 99.61  
#                       3rd Qu.:109.94  
#                       Max.   :157.90
  1. Plot the data with interaction.plot.
with(df1,interaction.plot(x.factor=task,
                          trace.factor=group,
                          response=Y) )

  1. Conduct an ANOVA to assess the main effects of A and B and the AxB interaction.
df1.aov.01 <- aov(Y~group*task,data=df1)
summary(df1.aov.01)
#             Df Sum Sq Mean Sq F value Pr(>F)   
# group        1   2210  2210.1   4.470 0.0391 * 
# task         2   5653  2826.3   5.716 0.0056 **
# group:task   2   3292  1646.2   3.329 0.0433 * 
# Residuals   54  26701   494.5                  
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  1. Evaluate the normality and constant-variance assumptions.
shapiro.test(residuals(df1.aov.01)) # normality
# 
#   Shapiro-Wilk normality test
# 
# data:  residuals(df1.aov.01)
# W = 0.99154, p-value = 0.953
df1$condition <- interaction(df1$group,df1$task)
bartlett.test(Y~condition,data=df1) # constant variance
# 
#   Bartlett test of homogeneity of variances
# 
# data:  Y by condition
# Bartlett's K-squared = 6.9568, df = 5, p-value = 0.2239
  1. Plot the marginal means for group and task.
op <- par(no.readonly = TRUE)
par(mfrow=c(1,2))
boxplot(Y~group,data=df1,main="Main Effect (Group)")
boxplot(Y~task,data=df1,main="Main Effect (Task)")

  1. Perform pairwise Tukey tests on the marginal means of task.
TukeyHSD(df1.aov.01,which="task")
#   Tukey multiple comparisons of means
#     95% family-wise confidence level
# 
# Fit: aov(formula = Y ~ group * task, data = df1)
# 
# $task
#             diff       lwr        upr     p adj
# t2-t1 -16.047329 -32.99383  0.8991681 0.0669637
# t3-t1 -23.215775 -40.16227 -6.2692786 0.0047896
# t3-t2  -7.168447 -24.11494  9.7780501 0.5679412
  1. Use aov to evaluate the hypotheses that i) performance in task 1 does not differ from average performance in the other two tasks; and ii) that this task contrast does not differ between groups.
# create contrast for  main effect:
levels(df1$group)
# [1] "patient" "control"
w.group <- c(-1,1) # patient vs control
contrasts(df1$group) <- cbind(w.group)
levels(df1$task)
# [1] "t1" "t2" "t3"
w.task <- c(-1,1/2,1/2) # task 1 vs avg(task2 & task3)
contrasts(df1$task) <- cbind(w.task) 
# perform anova:
df1.aov.02 <- aov(Y~group*task,data=df1)
summary(df1.aov.02,split=list(group=list(pVSc=1),task=list(t1vst23=1)))
#                            Df Sum Sq Mean Sq F value  Pr(>F)   
# group                       1   2210    2210   4.470 0.03913 * 
#   group: pVSc               1   2210    2210   4.470 0.03913 * 
# task                        2   5653    2826   5.716 0.00560 **
#   task: t1vst23             1   5139    5139  10.392 0.00215 **
# group:task                  2   3292    1646   3.329 0.04331 * 
#   group:task: pVSc.t1vst23  1   3234    3234   6.541 0.01338 * 
# Residuals                  54  26701     494                   
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# The `group` contrast is significant and is identical to the main effect. 
# The `task` contrast also is significant, so performance in task 1 differs 
# from average performance in tasks 2 and 3 (ignoring, or combining across, 
# groups). Finally, the group x contrast interaction is significant, indicating 
# that the **task contrast differed between groups**. Note that we could 
# now analyze simple main effects and see if the task contrast is significant
# in either group.