# required packages:
::p_load(tidyverse,
pacman
cowplot, emmeans)
44 Dealing with Higher Order Factorial ANOVA
This walk-through is a continuation of our previous work with factorial ANOVA. If you haven’t already, please check the factorial ANOVA and Interaction vignettes. In truth, nothing terribly new is being introduced here; we are just ramping up the complexity of our ANOVA models. Whereas before we considered a scenario with just 2 factors, here we will consider a 3 factor ANOVA. For all practical purposes you typically do not want to go higher than a 3 or 4 factor ANOVA, simply because of the exponential increase in complexity. For example, consider below where each letter is an IV and the resulting tests in the full factorial omnibus ANOVA:
- 2 factor: a, b, a×b
- 3 factor: a, b, c, a×b, a×c, b×c, a×b×c
- 4 factor: a, b, c, d, a×b, a×c, a×d, b×c, b×d, c×d, a×b×c, a×b×d, a×c×d, b×c×d, a×b×c×d
- 5 factor: ummm, no… just don’t (unless you really want to)
Increasing the number of factors not only increases the practical difficulty of analysis, but more importantly, makes it more difficult to interpret your results. The primary culprit here is the potential presence of multiple interactions.
When dealing with interactions in a higher order factorial design you always start off with the highest order interaction. For example in the 3 factor case, if you have a 3-way interaction, that supersedes any other effects and interactions, and must be dealt with first. If there is no three-way interaction, then you can co about the business of addressing any two way interactions that are present.
44.1 Packages and data
Let’s load in our necessary packages, scripts and some data and try some examples. This write-up requires the following packages:
And now let’s re-acquaint ourselves with our familar example dataset, but with one new twist…
Background: Given the ease of access for new technologies and increasing enrollment demands, many grade schools are advocating that teachers switch over to E-courses, where students view an online, pre-recorded lecture on the course topic in lieu of sitting in a classroom in a live lecture. Given this push, a critical question remains regarding the impact of E-courses on primary school student outcomes. More it may be the case that certain subject content more readily lends itself to E-course presentations than other subjects. To address this question we tested students performance on a test one week after listening to the related lecture. Lectures were either experienced via online (Computer) or in a live classroom (Standard). In addition, the lecture content varied in topic (Physical science, Social science, History). Finally, to assess whether age had an impact on expected outcomes we assessed students in 5th and 8th grades.
Looking at the above we have a 2 (Grade: 5th, 8th) × 2 (Presentation: Computer, Standard) × 3 (Lecture: Physical, Social, History) design. Our omnibus ANOVA therefore will test for the following effects:
- Main Effects: Grade, Presentation, Lecture
- 2-way interactions: Grade × Presentation, Grade × Lecture, Lecture × Presentation
- 3-way interaction: Grade × Presentation × Lecture
44.2 EXAMPLE 1: no three-way interaction, single two way interaction
Let’s load in this data:
<- read_delim("http://tehrandav.is/courses/statistics/practice_datasets/ANOVA5_higherEx1.txt", delim = "\t") higherEx1
Rows: 72 Columns: 5
── Column specification ────────────────────────────────────────────────────────
Delimiter: "\t"
chr (3): Grade, Lecture, Presentation
dbl (2): Subject, Score
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
show(higherEx1)
# A tibble: 72 × 5
Subject Grade Lecture Presentation Score
<dbl> <chr> <chr> <chr> <dbl>
1 1 5th Physical Computer 53
2 2 5th Physical Computer 49
3 3 5th Physical Computer 47
4 4 5th Physical Computer 42
5 5 5th Physical Computer 51
6 6 5th Physical Computer 34
7 7 5th Physical Standard 44
8 8 5th Physical Standard 48
9 9 5th Physical Standard 35
10 10 5th Physical Standard 18
# ℹ 62 more rows
It’s already in long format, but I should tell R
that Grade
, Lecture
, and Presentation
are factors:
$Grade <- factor(higherEx1$Grade)
higherEx1$Lecture <- factor(higherEx1$Lecture)
higherEx1$Presentation <- factor(higherEx1$Presentation) higherEx1
44.2.1 First plot: 3 way interaction plot
One thing you may notice is since we have a more complex ANOVA, we have a more complex design with three factors we also have to construct a more complex plot. Whereas in the 2 factor case we could place one factor on the x-axis and group the other factor by line or shape or color; in the three factor case things become slightly more tricky.
One thing to consider is the logic of what we are doing with our plots. When we have two factors, say A & B, we are distinguishing between each level of B (by lines, shapes, colors) at each level of A (on the x-axis). In the 3 factor case where we have A,B, & C, we are distinguishing the B×C interaction (through the combination of shapes, lines, colors) on each level of A.
In the following plot, let’s place Lecture type along the x-axis. From here, we can create 2 groups of lines (each with 2 levels) representing the Grade × Presentation interaction. We note this by specifying the interaction
in group
of our baseline ggplot
call:
ggplot(higherEx1,mapping = aes(x = Lecture,y = Score, group=interaction(Presentation,Grade)))
From here, we can proceed as before. Note, however, depending on the complexity of your plot you may also need to specify the interaction in each stat_summary::aes()
as necessary.
Below I am changing the shape
by Presentation
and the linetype
by Grade
. The entire call, including the line above is:
ggplot(higherEx1,mapping = aes(x = Lecture,y = Score, group=interaction(Presentation,Grade))) +
stat_summary(geom = "pointrange",fun.data = "mean_se",position = position_dodge(.15), aes(shape=Presentation)) +
stat_summary(geom = "line", fun="mean", position = position_dodge(.15), aes(linetype=Grade)) +
theme_cowplot() +
# renaming axies
xlab("Lecture") +
ylab ("Score") +
# adding whitespace around plot (optional)
theme(plot.margin=unit(c(.25,.25,.25,.25),"in")) +
# adjusting the legend
theme(
legend.title = element_text(size = 12, face = "bold"),
legend.position = c(.1,.25))
OK, not exactly the best plot. Let’s see if we can create some more space to work with by forcing the y-axis to go to 0, instead of cutting off at 10. There’s also a trick that we can use to stack our legend keys horizontally, rather than vertically
ggplot(higherEx1,mapping = aes(x = Lecture,y = Score, group=interaction(Presentation,Grade))) +
stat_summary(geom = "pointrange",fun.data = "mean_se",position = position_dodge(.15), aes(shape=Presentation)) +
stat_summary(geom = "line", fun="mean", position = position_dodge(.15), aes(linetype=Grade)) +
theme_cowplot() +
# renaming axies
xlab("Lecture") +
ylab ("Score") +
# adding whitespace around plot (optional)
theme(plot.margin=unit(c(.25,.25,.25,.25),"in")) +
# new additions / changes
coord_cartesian(ylim=c(0,50)) +
# direction to stack legend keys
theme(legend.direction = "horizontal",
legend.position = c(.1,.25))
alternatively, we could also just stack the legend boxes side by side:
ggplot(higherEx1,mapping = aes(x = Lecture,y = Score, group=interaction(Presentation,Grade))) +
stat_summary(geom = "pointrange",fun.data = "mean_se",position = position_dodge(.15), aes(shape=Presentation)) +
stat_summary(geom = "line", fun="mean", position = position_dodge(.15), aes(linetype=Grade)) +
theme_cowplot() +
# renaming axies
xlab("Lecture") +
ylab ("Score") +
# adding whitespace around plot (optional)
theme(plot.margin=unit(c(.25,.25,.25,.25),"in")) +
# new additions / changes
coord_cartesian(ylim=c(0,50)) +
# direction to stack legend boxes instead of keys
theme(legend.box = "horizontal",
legend.position = c(.1,.25))
NOTE: For what its work, we could also obviate the need for this complex of a plot, and just facet_wrap()
by our third factor. For practice try taking the code above and faceting by ~Grade.
44.2.2 running the ANOVA:
From here we run the omnibus as usual:
<- aov(Score~Lecture*Presentation*Grade, data = higherEx1)
omnibus_aov
%>% rstatix::anova_test(effect.size="pes") omnibus_aov
ANOVA Table (type II tests)
Effect DFn DFd F p p<.05 pes
1 Lecture 2 60 21.333000 1.00e-07 * 4.16e-01
2 Presentation 1 60 80.029000 1.22e-12 * 5.72e-01
3 Grade 1 60 0.006000 9.37e-01 1.06e-04
4 Lecture:Presentation 2 60 13.366000 1.58e-05 * 3.08e-01
5 Lecture:Grade 2 60 0.003000 9.97e-01 1.10e-04
6 Presentation:Grade 1 60 0.002000 9.62e-01 3.81e-05
7 Lecture:Presentation:Grade 2 60 0.000763 9.99e-01 2.54e-05
In this case we only have Main effects for Lecture and Presentation, and a Lecture × Presentation interaction. Given this we can disregard Grade as its not contributing to any effects.
44.2.3 Replotting the 2-way interaction
In some cases it is advisable to re-plot the data factoring out Grade by removing Grade from baseline group=
as well as any summary_stat::aes()
, (essentially treating this if it was a 2 factorial design like last week). For the sake of your sanity, and your readers’ sanity) I would also suggest being consistent with how you present your conditions. For example, in the original plot, Presentation was differentiated by shape
. To be consistent, we need to do the same here).
(You’ll note in the example below I’m also customizing font size and weight purely for personal aesthetics):
ggplot(higherEx1,mapping = aes(x = Lecture,y = Score, group=Presentation)) +
stat_summary(geom = "pointrange",fun.data = "mean_se",position = position_dodge(.15), aes(shape=Presentation)) +
stat_summary(geom = "line", fun="mean", position = position_dodge(.15), aes(shape=Presentation)) + theme_cowplot() + theme(
axis.title = element_text(size = 16, face = "bold", lineheight = .55),
axis.text = element_text(size = 12),
legend.title = element_text(size = 12, face = "bold"),
legend.position = c(.25,.25)) +
xlab("Lecture") +
ylab ("Score") +
theme(plot.margin=unit(c(.25,.25,.25,.25),"in")) +
coord_cartesian(ylim=c(0,50)) +
# stack legend boxes horizontally:
theme(legend.box = "horizontal")
Warning in stat_summary(geom = "line", fun = "mean", position =
position_dodge(0.15), : Ignoring unknown aesthetics: shape
From here I would deal with the follow-ups as we did in our 2×3 example from the last write-up: Run the simple effects and any necessary posthocs, being sure to correct using the error terms from the omnibus ANOVA.
44.3 EXAMPLE 2: No three-way interaction, multiple main effects, multiple 2 way interactions:
Let’s load in a data set that is a little more complicated:
<- read_delim("http://tehrandav.is/courses/statistics/practice_datasets/ANOVA5_higherEx2.txt", delim = "\t") higherEx2
Rows: 72 Columns: 5
── Column specification ────────────────────────────────────────────────────────
Delimiter: "\t"
chr (3): Grade, Lecture, Presentation
dbl (2): Subject, Score
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
44.3.1 plotting the data
ggplot(higherEx2,mapping = aes(x = Lecture,y = Score, group=interaction(Presentation,Grade))) +
stat_summary(geom = "pointrange",fun.data = "mean_se",position = position_dodge(.15), aes(shape=Presentation)) +
stat_summary(geom = "line", fun="mean", position = position_dodge(.15), aes(linetype=Grade)) + theme_classic() + theme(
axis.title = element_text(size = 16, face = "bold", lineheight = .55),
axis.text = element_text(size = 12),
legend.title = element_text(size = 12, face = "bold"),
legend.position = c(.25,.25)) +
xlab("Lecture") +
ylab ("Score") +
theme(plot.margin=unit(c(.25,.25,.25,.25),"in")) +
coord_cartesian(ylim=c(0,60)) +
# stack legend boxes horizontally:
theme(legend.box = "horizontal")
Looking at this plot its quite apparent that the pattern of results is different from Example 1. Perhaps most simply, there is a general shared pattern of effects for each of the shapes (Presentation) with the exception of the 8th Grade-Computer condition. For the remainder there is a dip for Social lectures, but in this one condition Social actually increases. Let’s run our ANOVA.
44.3.2 running the omnibus ANOVA:
<- aov(Score~Lecture*Presentation*Grade, data = higherEx2)
omnibus_aov
%>% rstatix::anova_test(effect.size="pes") omnibus_aov
ANOVA Table (type II tests)
Effect DFn DFd F p p<.05 pes
1 Lecture 2 60 7.089000 2.00e-03 * 1.91e-01
2 Presentation 1 60 80.029000 1.22e-12 * 5.72e-01
3 Grade 1 60 5.645000 2.10e-02 * 8.60e-02
4 Lecture:Presentation 2 60 13.366000 1.58e-05 * 3.08e-01
5 Lecture:Grade 2 60 5.532000 6.00e-03 * 1.56e-01
6 Presentation:Grade 1 60 0.002000 9.62e-01 3.81e-05
7 Lecture:Presentation:Grade 2 60 0.000763 9.99e-01 2.54e-05
Our ANOVA reveals an abundance of results. There are main effects for each of our IVs. In addition there are two interaction effects: Lecture:Presentation
and Lecture:Grade
. Unless there is a strong theoretical reason not to (which I see none) we will need to examine each of these interactions in further detail.
44.3.3 testing the Lecture:Presentation
interaction:
44.3.3.1 2-way interaction plot
We can begin by recreating our interaction plot, this time only focusing on the IVs that are of interest (those that are interacting). In this first case, they are Lecture
and Presentation
.
ggplot(higherEx2,mapping = aes(x = Lecture,y = Score, group=Presentation)) +
stat_summary(geom = "pointrange",fun.data = "mean_cl_normal",position = position_dodge(.15), aes(shape=Presentation)) +
stat_summary(geom = "line", fun="mean", position = position_dodge(.15), aes(shape=Presentation)) + theme_cowplot() + theme(
axis.title = element_text(size = 16, face = "bold", lineheight = .55),
axis.text = element_text(size = 12),
legend.title = element_text(size = 12, face = "bold"),
legend.position = c(.25,.25)) +
xlab("Lecture") +
ylab ("Score") +
theme(plot.margin=unit(c(.25,.25,.25,.25),"in")) +
coord_cartesian(ylim=c(0,60)) +
# stack legend boxes horizontally:
theme(legend.box = "horizontal")
Warning in stat_summary(geom = "line", fun = "mean", position =
position_dodge(0.15), : Ignoring unknown aesthetics: shape
44.3.3.2 simple effect ANOVAs:
Provided what we see on this plot, it may make sense to first run a simple effects ANOVA for each Presentation and then run the appropriate follow-ups. I say this because the more obvious, and potentially easily interpret-able effects occur moving across the line series (i.e., lines trend up or down).
When testing for simple effects of Lecture on each Presentation type, we still need to include Grade in our analysis—this is to account for Grade’s contribution to our original model. Essentially we conduct a test for Lecture, Grade, and the Lecture by Grade interaction on each level of Presentation type. Using the template from the previous walkthroughs, this can be accomplished by simply telling R
to run joint_tests
separating the original model by
Presentation.
Using emmeans::joint_tests()
:
::joint_tests(omnibus_aov, by = "Presentation") emmeans
Presentation = Computer:
model term df1 df2 F.ratio p.value
Lecture 2 60 4.739 0.0123
Grade 1 60 2.710 0.1049
Lecture:Grade 2 60 2.823 0.0673
Presentation = Standard:
model term df1 df2 F.ratio p.value
Lecture 2 60 15.715 <.0001
Grade 1 60 2.937 0.0917
Lecture:Grade 2 60 2.710 0.0747
Both simple effects ANOVA of Lecture within Presentation were significant, suggesting that we should perform a post-hoc analyses of each:
emmeans(omnibus_aov, by = "Presentation", specs="Lecture") %>%
pairs(adjust="tukey")
NOTE: Results may be misleading due to involvement in interactions
Presentation = Computer:
contrast estimate SE df t.ratio p.value
History - Physical -7.92 3.02 60 -2.624 0.0292
History - Social -8.17 3.02 60 -2.707 0.0237
Physical - Social -0.25 3.02 60 -0.083 0.9962
Presentation = Standard:
contrast estimate SE df t.ratio p.value
History - Physical -3.00 3.02 60 -0.994 0.5833
History - Social 12.92 3.02 60 4.281 0.0002
Physical - Social 15.92 3.02 60 5.275 <.0001
Results are averaged over the levels of: Grade
P value adjustment: tukey method for comparing a family of 3 estimates
The results of the Tukey post-hoc suggests that History scores are significantly less than the other two groups when considering only the Computer presentation. For the standard presentation, scores in the History and Physical lecture were greater than the Social lecture.
44.3.4 testing the Lecture:Grade
interaction:
Now we do the same for the Lecture:Grade
interaction. This time removing Presentation
from our consideration.
44.3.4.1 2-way interaction plot
You’ll notice that in this interaction plot, and the plot above, I am using the shape/linetype conventions that I established in my original plot. This will make things much easier if I wind up comparing these plots to my 3-way plot. I know I’ve said this before, but ONE MORE TIME WITH FEELING!
ggplot(higherEx2,mapping = aes(x = Lecture,y = Score, group=Grade)) +
stat_summary(geom = "pointrange",fun.data = "mean_se",position = position_dodge(.15), aes(shape=Presentation)) +
stat_summary(geom = "line", fun="mean", position = position_dodge(.15), aes(linetype=Grade)) + theme_classic() + theme(
axis.title = element_text(size = 16, face = "bold", lineheight = .55),
axis.text = element_text(size = 12),
legend.title = element_text(size = 12, face = "bold"),
legend.position = c(.25,.25)) +
xlab("Lecture") +
ylab ("Score") +
theme(plot.margin=unit(c(.25,.25,.25,.25),"in")) +
coord_cartesian(ylim=c(0,60)) +
# stack legend boxes horizontally:
theme(legend.box = "horizontal")
44.3.5 simple effect ANOVAs:
Provided what we see on this plot, it may make the best sense to examine the pairwise differences in each Grade group by Lecture. Whereas in the previous interaction, the compelling changes occurred across lines, here what is more compelling is the presence/absence of gaps between the lines. In particular, while there are practically no differences due to Grade in the History or Physical Lectures, there is a larger gap in the Social condition.
Running these simple effects ANOVAs using emmeans:joint_tests()
. Keep in mind, since we are parsing the Lecture:Grade
interaction, be should only take a look at the effects of Grade
on each of our Presentation
levels. We can ignore the other output.
::joint_tests(omnibus_aov,by = "Lecture") emmeans
Lecture = History:
model term df1 df2 F.ratio p.value
Presentation 1 60 5.255 0.0254
Grade 1 60 0.001 0.9781
Presentation:Grade 1 60 0.001 0.9781
Lecture = Physical:
model term df1 df2 F.ratio p.value
Presentation 1 60 15.382 0.0002
Grade 1 60 0.000 1.0000
Presentation:Grade 1 60 0.003 0.9561
Lecture = Social:
model term df1 df2 F.ratio p.value
Presentation 1 60 86.123 <.0001
Grade 1 60 16.710 0.0001
Presentation:Grade 1 60 0.000 1.0000
As anticipated from looking at the plot, there was a significant difference between our Grades
in the Social group, but not the other two. Since there are only two levels of Grade
we don’t deen to run any other post hocs (the \(F\) test is our comparison between means).
44.3.6 Main takeaways / write-up
It’s typically best to go back and look at all of your plots.
the main effect for
Presentation
is meaningful: Regardless of what plays out on the other factors, students tended to perform better in the computer presentation. Best to report that.the
Lecture:Presentation
interaction resulted from the following:- in the Computer presentations scores were lowest in the History lecture (indifferent in other two)
- in the Standard presentations scores were lowest in the Social lecture (indifferent in the other two)
the
Lecture:Grade
interaction resulted from the 5th graders performing worse than the eight graders in Social lectures; while performance was not different in the other two lectures.
With this in mind, an example write-up:
Our results revealed main effects for each of our factors, in addition to Lecture × Presentation \([F(2,60)=13.37, p<.001, \eta_p^2=.31]\) and Lecture × Grade interactions \([F(2,60)=5.53, p=.006, \eta_p^2=.16]\).
Considering the former interaction, students presented material via computer format tended to perform poorest from the History Lecture, compared to Physical and Social (Tukey HSD, p<.05). However, when material was presented in the standard format, students performed poorest in the Social condition (p<.05), while the other two were not different.
Considering the latter interaction, a detailed analysis of performance as a function of lecture revaled that 8th Graders performed better than their 5th grade counterparts only in the Social lecture, \(F(1,60) = 16.71, p < .001, \eta_p^2=0.21\).
In all conditions, students performed better when presented the material via computer compared to the standard presentation, \([F(2,60)=7.09, p=.002, \eta_p^2=.19]\) (see Figure 1).
44.4 EXAMPLE 3: OMG, multiple two way interactions and a nasty three-way!!!
Finally let’s take a look at some data that is all over the place:
<- read_delim("http://tehrandav.is/courses/statistics/practice_datasets/ANOVA5_higherEx3.txt", delim = "\t") higherEx3
Rows: 72 Columns: 5
── Column specification ────────────────────────────────────────────────────────
Delimiter: "\t"
chr (3): Grade, Lecture, Presentation
dbl (2): Subject, Score
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
44.4.1 Plotting:
ggplot(higherEx3,mapping = aes(x = Lecture,y = Score, group=interaction(Presentation,Grade))) +
stat_summary(geom = "pointrange",fun.data = "mean_se",position = position_dodge(.5), aes(shape=Presentation)) +
stat_summary(geom = "line", fun="mean", position = position_dodge(.5), aes(linetype=Grade)) +
theme_classic() + theme(
axis.title = element_text(size = 16, face = "bold", lineheight = .55),
axis.text = element_text(size = 12),
legend.title = element_text(size = 12, face = "bold"),
legend.position = c(.25,.15)) +
xlab("Lecture") +
ylab ("Score") +
theme(plot.margin=unit(c(.25,.25,.25,.25),"in")) +
coord_cartesian(ylim=c(0,50)) +
# stack legend boxes horizontally:
theme(legend.box = "horizontal")
44.4.2 running the ANOVA:
<- aov(Score~Lecture*Presentation*Grade, data = higherEx3)
omnibus_aov
%>% rstatix::anova_test(effect.size="pes") omnibus_aov
ANOVA Table (type II tests)
Effect DFn DFd F p p<.05 pes
1 Lecture 2 60 11.644 5.34e-05 * 0.280
2 Presentation 1 60 36.149 1.17e-07 * 0.376
3 Grade 1 60 65.619 3.30e-11 * 0.522
4 Lecture:Presentation 2 60 4.761 1.20e-02 * 0.137
5 Lecture:Grade 2 60 3.784 2.80e-02 * 0.112
6 Presentation:Grade 1 60 16.259 1.58e-04 * 0.213
7 Lecture:Presentation:Grade 2 60 3.874 2.60e-02 * 0.114
O.M.F.G., everything is significant!!! What to do!!!
Remember, the first thing that we do is tackle the highest-order interaction. In this case we should tease-out the three-way. Again, we should look to our plot to help guide us in what to do. Looking at the plot it appears that different things are happening by Grade. The Presentation:Lecture
lines tend to stay parallel for the 5th graders, but not as much for the 8th graders. At the same time the gaps between the two lines are different by grade. So, my advice would be to attack the 3-way interaction by looking at the individual Lecture:Presentation
interactions on each Grade
.
44.4.3 Fifth Graders, Lecture:Presentation interaction.
::joint_tests(omnibus_aov, by="Grade") emmeans
Grade = 5th:
model term df1 df2 F.ratio p.value
Lecture 2 60 2.167 0.1235
Presentation 1 60 1.960 0.1666
Lecture:Presentation 2 60 0.024 0.9763
Grade = 8th:
model term df1 df2 F.ratio p.value
Lecture 2 60 13.261 <.0001
Presentation 1 60 50.448 <.0001
Lecture:Presentation 2 60 8.611 0.0005
You’ll note that our results include simple effects for Lecture
and Presentation
at each Grade
for as their interaction. That said, there are no effects when looking at the 5th
graders. However, the detailed Lecture:Presentation
analysis for the 8th graders revealed both main effects and their interaction. We need to follow-up on this interaction in the 8th-grade data. In this case we would look at the effect of Lecture on each level of Presentation for the eighth graders.
To do this we first subset the data further, by both Grade
and Presentation
:
::joint_tests(omnibus_aov, by=c("Grade", "Presentation")) emmeans
Grade = 5th, Presentation = Computer:
model term df1 df2 F.ratio p.value
Lecture 2 60 1.100 0.3393
Grade = 8th, Presentation = Computer:
model term df1 df2 F.ratio p.value
Lecture 2 60 2.264 0.1127
Grade = 5th, Presentation = Standard:
model term df1 df2 F.ratio p.value
Lecture 2 60 1.090 0.3428
Grade = 8th, Presentation = Standard:
model term df1 df2 F.ratio p.value
Lecture 2 60 19.607 <.0001
Because there were no effects for 5th
graders in the previous analysis, we can ignore that output. Focusing on the 8th
grade output, we see that there is no effect for lecture type for 8th
graders using the Computer
presentation. However, the Standard
presentation was significant. Following up on that:
emmeans(omnibus_aov, by=c("Grade", "Presentation"), specs = "Lecture") %>% pairs(adjust="tukey")
Grade = 5th, Presentation = Computer:
contrast estimate SE df t.ratio p.value
History - Physical -5.33 3.78 60 -1.411 0.3417
History - Social -1.17 3.78 60 -0.309 0.9489
Physical - Social 4.17 3.78 60 1.102 0.5164
Grade = 8th, Presentation = Computer:
contrast estimate SE df t.ratio p.value
History - Physical -7.83 3.78 60 -2.072 0.1042
History - Social -2.33 3.78 60 -0.617 0.8112
Physical - Social 5.50 3.78 60 1.455 0.3197
Grade = 5th, Presentation = Standard:
contrast estimate SE df t.ratio p.value
History - Physical -4.83 3.78 60 -1.279 0.4126
History - Social 0.00 3.78 60 0.000 1.0000
Physical - Social 4.83 3.78 60 1.279 0.4126
Grade = 8th, Presentation = Standard:
contrast estimate SE df t.ratio p.value
History - Physical -3.00 3.78 60 -0.794 0.7083
History - Social 18.83 3.78 60 4.983 <.0001
Physical - Social 21.83 3.78 60 5.776 <.0001
P value adjustment: tukey method for comparing a family of 3 estimates
Again, this output give us a lot of info, BUT all we should focus on is the section of interest, in this case 8th Grade, Standard Presentation
or the last section. Here, scores for the Social
test were less than the Physical
and History.
44.4.4 Constructing a narrative:
Given these results, what narrative do you think you can tell? Remember use the plots to unpack what the results are telling you. Our results focused on differences in the Lecture:Presentation
analysis for each grade. I would recommend that your story start there. Something to the effect of Lecture and Presentation interact for 8th graders but not 5th graders.
To remind ourselves:
- Main effects for Lecture, Presentation, and Grade
- Two way interactions
- Three-way interaction that we focus on
- The results of the 3-way interaction, by Grade:
- 5th graders = no effects
- 8th graders = Lecture × Presentation interaction, where there was a null effect for Lecture on the Computer presentation, but an effect for Lecture on the Standard Presentation.
Now the narrative:
… to answer this question we conducted a 2 (Grade) × 2 (Presentation) × 3 (Lecture) between effects ANOVA. Our ANOVA revealed a three-way, Grade × Presentation × Lecture type interaction, \(F\)(2, 60) = 3.87, \(p\) = .03, \(\eta_p^2\) = .11. To address this interaction we conducted separate Presentation × Lecture ANOVA for each Grade. Our results for 5th graders showed no effects nor an interaction (\(p\) > .05). For 8th graders, we found an Presentation × Lecture interaction, \(F\)(2, 60) = 8.61, \(p\) < .001. While there were no differences due to Lecture type for 8th graders presented information via computer (\(p\) > .05), 8th graders in the Standard presentation, \(F\)(2, 60) = 19.61, \(p\) < .001 tended to perform worse for the Social lecture (\(M\) ± \(SD\) = 12.33 ± 3.72) than for the History (31.17 ± 9.87) and Physical (34.17 ± 10.98) lectures (Tukey HSD, \(p\) < .05).
Now imagine trying to deal with a FOUR-WAY INTERACTION!!!! SHEESH!!!
On to the next one!