×

The following packages are used in this section: dplyr, tibble, ggplot2, afex, emmeans, ggbeeswarm and teachingtools. dplyr, tibble and ggplot2 are part of the tidyverse family.

To make sure that the packages used in this section are installed and loaded, run the following code in your console.

xfun::pkg_attach2("afex","ggbeeswarm","emmeans","teachingtools","tidyverse")

In this session we’ll learn how to perform ANOVA in R. There are many different ways to perform the humble ANOVA in R, and there are lots of alternatives to ANOVA available in R, but we won’t be discussing those today. Instead, the approach that I’m going to take is to teach you how to something (almost) identical to the kind of ANOVA you might already be performing in SPSS. For this reason, we’ll focus on the afex package. It provides a one-stop shop for most of the things you’ll want to do with ANOVA.

We’ll cover all the basic kinds of ANOVA including within-subjects, between-subjects, and mixed ANOVA.

# Types of ANVOA

## Between-subjects one-way ANVOA

We’ll start super simple with a between subjects one-way ANOVA. We’ll use the data below. We just have our factor group, our dependent variable anxiety, and a column with the participant ids called id.

aov_data_oneway
# A tibble: 60 x 3
id    group   anxiety
<chr> <chr>     <dbl>
1 s001  Group 1      17
2 s002  Group 2      41
3 s003  Group 3      16
4 s004  Group 1      18
5 s005  Group 2      44
6 s006  Group 3       9
7 s007  Group 1      16
8 s008  Group 2      39
9 s009  Group 3      23
10 s010  Group 1      13
# … with 50 more rows

Now let’s try to do a one-way ANOVA. We’re going to compare anxiety across the different groups. We’ll use the aov_ez() function from the afex package. The structure of the basic function call is as follows.

afex::aov_ez(id = "id", # the name of our ID variable
dv = "anxiety",  # the name of our DV
between = "group", # the name of our between groups factor
data = aov_data_oneway # our data table
)
Converting to factor: group
Contrasts set to contr.sum for the following variables: group
Anova Table (Type 3 tests)

Response: anxiety
Effect    df   MSE          F  ges p.value
1  group 2, 57 26.75 141.52 *** .832   <.001
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1

When we run this we get a print out of an ANOVA table. Just having a print out to console might not be super useful. So we’ll save the output as an object which will allow us to do some work with it.

aov_oneway <- afex::aov_ez(id = "id", # the name of our ID variable
dv = "anxiety",  # the name of our DV
between = "group", # the name of our between groups factor
data = aov_data_oneway # our data table
)
Converting to factor: group
Contrasts set to contr.sum for the following variables: group

By saving the output as an object we can make some changes to the information that afex give us. For example, by default afex gives us an effect size called generalised eta-squared ($$\eta_{\textrm{G}}^2$$) rather than the partial eta-squared ($$\eta_{\textrm{p}}^2$$) that we get from SPSS. But we can ask afex to give us this instead. To do this, we’ll use the anova() function. And we’ll set an argument es (for effect size) to pes for $$\eta_{\textrm{p}}^2$$.

anova(object = aov_oneway, # our anova object
es = "pes" # pes for partial eta squared, or ges for generalised eta squared
)
Anova Table (Type 3 tests)

Response: anxiety
num Df den Df    MSE      F     pes    Pr(>F)
group      2     57 26.752 141.52 0.83237 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## Factorial ANOVA

To extend our example to factorial ANOVA is just a trivial extension of what we’ve already done. All we need is a dataset with another factor. If we take a look at the dataset named aov_data_factorial we’ll see that it has an additional factor called gender.

aov_data_factorial
# A tibble: 60 x 4
id    group   anxiety gender
<chr> <chr>     <dbl> <chr>
1 s001  Group 1      17 F
2 s002  Group 2      41 M
3 s003  Group 3      16 M
4 s004  Group 1      18 M
5 s005  Group 2      44 M
6 s006  Group 3       9 M
7 s007  Group 1      16 M
8 s008  Group 2      39 F
9 s009  Group 3      23 M
10 s010  Group 1      13 M
# … with 50 more rows

We can just add in gender in to the between argument in our call to aov_ez().

aov_factorial <- afex::aov_ez(id = "id", # the name of our ID variable
dv = "anxiety",  # the name of our DV
between = c("group","gender"), # the names of our factors as a vector
data = aov_data_factorial # our factorial data table
)
aov_factorial
Anova Table (Type 3 tests)

Response: anxiety
Effect    df   MSE          F  ges p.value
1        group 2, 54 27.44 134.66 *** .833   <.001
2       gender 1, 54 27.44       0.21 .004    .648
3 group:gender 2, 54 27.44       0.69 .025    .504
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1

## Repeated measures ANOVA

The extension to repeated measures ANOVA is again simply a matter of specifying our within subject factor(s). If you’re coming from the SPSS world you’ll be used to your data being in wide format. But for R we’ll have to have it in long format.

We’ll use the data aov_data_repeated

aov_data_repeated
# A tibble: 480 x 5
subjectID    RT emotion congruency group
<chr>     <dbl> <fct>   <fct>      <chr>
1 s001        506 happy   congruent  Group 1
2 s002        479 happy   congruent  Group 1
3 s003        459 happy   congruent  Group 1
4 s004        493 happy   congruent  Group 1
5 s005        489 happy   congruent  Group 1
6 s006        492 happy   congruent  Group 1
7 s007        507 happy   congruent  Group 1
8 s008        502 happy   congruent  Group 1
9 s009        515 happy   congruent  Group 1
10 s010        518 happy   congruent  Group 1
# … with 470 more rows

Now we can just specify the ANOVA using aov_ez() as before, but we’ll specify emotion and congruency in the within argument.

aov_repeated <- afex::aov_ez(
id = "subjectID", # the name of our ID variable
dv = "RT", # the name of our DV
within = c("emotion","congruency"), # the names of our factors
between = "group", # our between subjects factor
data = aov_data_repeated # our repeated measures data table
)
aov_repeated
Anova Table (Type 3 tests)

Response: RT
Effect           df    MSE    F  ges p.value
1                    group        1, 58 272.93 0.56 .001    .456
2                  emotion 2.75, 159.44 292.17 0.40 .003    .738
3            group:emotion 2.75, 159.44 292.17 0.50 .003    .666
4               congruency        1, 58 259.74 2.69 .006    .106
5         group:congruency        1, 58 259.74 1.75 .004    .191
6       emotion:congruency 2.77, 160.79 276.99 1.79 .011    .156
7 group:emotion:congruency 2.77, 160.79 276.99 0.19 .001    .886
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1

Sphericity correction method: GG 

One of the things you’ll notice with this output is that aov_ez() will automatically apply a correction for violations of sphericity if it’s needed. By default, it applies the Greenhouse-Geisser correct, and the table that it prints out includes the corrected DFs and the corresponding p value. If you instead want the Huynh—Feldt correction then you can specify this with the anova() function using the correction argument and setting it to HF.

anova(object = aov_repeated, # our anova object
correction = "HF" # HF for Huynh—Feldt, GG for Greenhouse-Geisser
)
Anova Table (Type 3 tests)

Response: RT
num Df den Df    MSE      F       ges Pr(>F)
group                    1.0000  58.00 272.93 0.5643 0.0012608 0.4556
emotion                  2.8994 168.16 277.02 0.3966 0.0026042 0.7488
group:emotion            2.8994 168.16 277.02 0.5016 0.0032907 0.6754
congruency               1.0000  58.00 259.74 2.6944 0.0057030 0.1061
group:congruency         1.0000  58.00 259.74 1.7484 0.0037082 0.1913
emotion:congruency       2.9253 169.67 262.49 1.7854 0.0111110 0.1532
group:emotion:congruency 2.9253 169.67 262.49 0.1949 0.0012249 0.8956

Note, however, that we don’t get the $$\epsilon$$ values or the uncorrected DFs. And, depending on your personal reporting preferences, you might want these. To get all the information, including the $$\epsilon$$ values, and the output from the Mauchly Tests for Sphericity, we can just use the summary() function1.

summary(aov_repeated)

Univariate Type III Repeated-Measures ANOVA Assuming Sphericity

Sum Sq num Df Error SS den Df    F value Pr(>F)
(Intercept)              66296515      1    15830     58 2.4291e+05 <2e-16 ***
group                         154      1    15830     58 5.6430e-01 0.4556
emotion                       319      3    46584    174 3.9660e-01 0.7556
group:emotion                 403      3    46584    174 5.0160e-01 0.6817
congruency                    700      1    15065     58 2.6944e+00 0.1061
group:congruency              454      1    15065     58 1.7484e+00 0.1913
emotion:congruency           1371      3    44536    174 1.7854e+00 0.1517
group:emotion:congruency      150      3    44536    174 1.9490e-01 0.8998
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Mauchly Tests for Sphericity

Test statistic p-value
emotion                         0.85789 0.12195
group:emotion                   0.85789 0.12195
emotion:congruency              0.87256 0.17166
group:emotion:congruency        0.87256 0.17166

Greenhouse-Geisser and Huynh-Feldt Corrections
for Departure from Sphericity

GG eps Pr(>F[GG])
emotion                  0.91635     0.7382
group:emotion            0.91635     0.6655
emotion:congruency       0.92405     0.1564
group:emotion:congruency 0.92405     0.8864

HF eps Pr(>F[HF])
emotion                  0.9664578  0.7488042
group:emotion            0.9664578  0.6753869
emotion:congruency       0.9750973  0.1532343
group:emotion:congruency 0.9750973  0.8956113

# Follow-up tests

Once you’ve run an ANOVA the next thing you’ll probably want to do is run some follow-up tests. Usually the first thing we’ll do when doing follow up test is calculate the estimated marginal means. To do this, we’ll use the emmeans::emmeans function from the emmeans.

For our simple one-way ANOVA (aov_oneway) all we need to do is specify the the name of our one predictor for which we want our estimated marginal means (i.e., group)

oneway_emm <- emmeans::emmeans(object = aov_oneway, # our ANOVA object
specs =  "group", # the name of our factor
)
oneway_emm
 group   emmean   SE df lower.CL upper.CL
Group 1   14.2 1.16 57     11.8     16.5
Group 2   40.8 1.16 57     38.4     43.1
Group 3   21.4 1.16 57     19.0     23.7

Confidence level used: 0.95 

When we have a more complex ANOVA like our factorial ANOVA, then we can specify the predictors for which we want estimated marginal means, but we can also specify any predictors we want to condition on using the by argument.

factorial_emm <- emmeans::emmeans(
object = aov_factorial, # out anova object
specs =  "group", # the predictor we want EMMs for
by = "gender" # the predictor to condition on
)

factorial_emm 
gender = F:
group   emmean   SE df lower.CL upper.CL
Group 1   13.9 1.98 54     9.89     17.8
Group 2   41.4 1.75 54    37.94     44.9
Group 3   19.6 1.98 54    15.60     23.5

gender = M:
group   emmean   SE df lower.CL upper.CL
Group 1   14.3 1.45 54    11.39     17.2
Group 2   40.2 1.58 54    37.02     43.3
Group 3   22.3 1.45 54    19.39     25.2

Confidence level used: 0.95 

We can also do the same for our repeated measures ANOVA (aov_repeated). And we can even collapse across predictors by not including them.

For example, to drop group:

repeated_emm <- emmeans::emmeans(
object = aov_repeated, # our anova object
spec = "emotion", # the preditor we want EMMs for
by = "congruency" # the predictor to condition on
)
repeated_emm 
congruency = congruent:
emotion emmean   SE  df lower.CL upper.CL
happy      498 2.73 454      493      504
sad        499 2.73 454      493      504
digust     501 2.73 454      495      506
anger      492 2.73 454      487      498

congruency = incongruent:
emotion emmean   SE  df lower.CL upper.CL
happy      498 2.73 454      493      504
sad        500 2.73 454      495      506
digust     501 2.73 454      495      506
anger      503 2.73 454      498      509

Results are averaged over the levels of: group
Warning: EMMs are biased unless design is perfectly balanced
Confidence level used: 0.95 

Or if we want to include group:

mixed_emm <- emmeans::emmeans(
object = aov_repeated, # our anova object
spec = c("congruency","emotion"), # the preditor we want EMMs for
by = "group" # the predictors to condition on
)
mixed_emm 
group = Group 1:
congruency  emotion emmean   SE  df lower.CL upper.CL
congruent   happy      497 4.95 449      487      506
incongruent happy      500 4.95 449      490      510
congruent   sad        494 4.95 449      485      504
incongruent sad        499 4.95 449      490      509
congruent   digust     500 4.95 449      490      510
incongruent digust     500 4.95 449      490      510
congruent   anger      491 4.95 449      481      501
incongruent anger      506 4.95 449      496      516

group = Group 2:
congruency  emotion emmean   SE  df lower.CL upper.CL
congruent   happy      500 2.48 399      495      505
incongruent happy      497 2.48 399      492      502
congruent   sad        503 2.48 399      498      508
incongruent sad        502 2.48 399      497      506
congruent   digust     501 2.48 399      497      506
incongruent digust     501 2.48 399      496      506
congruent   anger      494 2.48 399      489      499
incongruent anger      501 2.48 399      496      506

Warning: EMMs are biased unless design is perfectly balanced
Confidence level used: 0.95 

A note about replicating the SPSS output for Estimated Marginal Means

The afex and emmeans calculate the SD for estimated marginal means in a manner that is different to SPSS does it. Specifically, emmeans will give you the same SD/SE for the means and this will be calculated off the average SD. This is not what SPSS does. But it is possible to replicate the SPSS behaviour if you want. To do this you can use the ez package. ez is less flexible than afex, but it’s main goal is just replicating the numbers you get from SPSS. To calculate estimated marginal means in the style of SPSS you can just use the ez::ezStats() function. The inputs are almost identical to the inputs for afex::aov_ez(), expect for: 1) we leave out the " around the column names, 2) we replace the c with a . for concatenating our vectors, and 3) the input for the id column is called wid and not id.

# our old call just for reference
# aov_repeated <- afex::aov_ez(id = "subjectID", # the name of our ID variable
#                             dv = "RT", # the name of our DV
#                             within = c("emotion","congruency"), # the names of our factors
#                             between = "group", # our between subjects factor
#                             data = aov_data_repeated # our repeated measures data table
#                             )

xfun::pkg_attach2("ez") # load the ezANOVA package

ez::ezStats(wid = subjectID, # the name of our ID variable
dv = RT, # the name of our DV
within = .(emotion,congruency), # the names of our factors
between = group,
data = aov_data_repeated # our repeated measures data table
)
Warning: Converting "subjectID" to factor for ANOVA.
Warning: Converting "group" to factor for ANOVA.
Warning: Data is unbalanced (unequal N per group). Make sure you specified a
well-considered value for the type argument to ezANOVA().
Warning: Mixed within-and-between-Ss effect requested; FLSD is only appropriate
for within-Ss comparisons (see warning in ?ezStats or ?ezPlot).
Warning in ez::ezStats(wid = subjectID, dv = RT, within = .(emotion,
congruency), : Unbalanced groups. Mean N will be used in computation of FLSD
     group emotion  congruency  N   Mean       SD     FLSD
1  Group 1   happy   congruent 10 496.00 17.74511 8.152932
2  Group 1   happy incongruent 10 499.30 10.61498 8.152932
3  Group 1     sad   congruent 10 493.80 14.77836 8.152932
4  Group 1     sad incongruent 10 498.90 15.80752 8.152932
5  Group 1  digust   congruent 10 499.40 15.88290 8.152932
6  Group 1  digust incongruent 10 499.50 20.32650 8.152932
7  Group 1   anger   congruent 10 490.50 20.84999 8.152932
8  Group 1   anger incongruent 10 505.40  8.66923 8.152932
9  Group 2   happy   congruent 50 499.38 16.56760 8.152932
10 Group 2   happy incongruent 50 496.44 15.83739 8.152932
11 Group 2     sad   congruent 50 502.66 16.09488 8.152932
12 Group 2     sad incongruent 50 501.02 15.79782 8.152932
13 Group 2  digust   congruent 50 500.88 14.64915 8.152932
14 Group 2  digust incongruent 50 500.86 17.71084 8.152932
15 Group 2   anger   congruent 50 493.30 17.48264 8.152932
16 Group 2   anger incongruent 50 500.42 15.55122 8.152932

Once we’ve created our estimated marginal means using emmeans() we can compute various kinds of contrasts and comparisons. The most straightforward are simple pairwise comparisons. To do this we can use the pairs() function (pairs(), like summary() is another generic)

pairs(oneway_emm)
 contrast          estimate   SE df t.ratio p.value
Group 1 - Group 2    -26.6 1.64 57 -16.263 <.0001
Group 1 - Group 3     -7.2 1.64 57  -4.402 0.0001
Group 2 - Group 3     19.4 1.64 57  11.861 <.0001

P value adjustment: tukey method for comparing a family of 3 estimates 
pairs(repeated_emm)
congruency = congruent:
contrast       estimate   SE  df t.ratio p.value
happy - sad       -0.54 3.96 348 -0.136  0.9991
happy - digust    -2.45 3.96 348 -0.618  0.9263
happy - anger      5.79 3.96 348  1.461  0.4624
sad - digust      -1.91 3.96 348 -0.482  0.9631
sad - anger        6.33 3.96 348  1.597  0.3816
digust - anger     8.24 3.96 348  2.079  0.1620

congruency = incongruent:
contrast       estimate   SE  df t.ratio p.value
happy - sad       -2.09 3.96 348 -0.527  0.9524
happy - digust    -2.31 3.96 348 -0.583  0.9372
happy - anger     -5.04 3.96 348 -1.272  0.5817
sad - digust      -0.22 3.96 348 -0.056  0.9999
sad - anger       -2.95 3.96 348 -0.744  0.8791
digust - anger    -2.73 3.96 348 -0.689  0.9014

Results are averaged over the levels of: group
P value adjustment: tukey method for comparing a family of 4 estimates 
pairs(factorial_emm)
gender = F:
contrast          estimate   SE df t.ratio p.value
Group 1 - Group 2   -27.59 2.64 54 -10.450 <.0001
Group 1 - Group 3    -5.71 2.80 54  -2.041 0.1123
Group 2 - Group 3    21.87 2.64 54   8.285 <.0001

gender = M:
contrast          estimate   SE df t.ratio p.value
Group 1 - Group 2   -25.87 2.15 54 -12.056 <.0001
Group 1 - Group 3    -8.00 2.05 54  -3.893 0.0008
Group 2 - Group 3    17.87 2.15 54   8.328 <.0001

P value adjustment: tukey method for comparing a family of 3 estimates 

By default pairs() uses Tukey for adjusting p values. However, you can can select a few other options by using the adjust argument.

pairs(factorial_emm, adjust = "tukey") # the default
gender = F:
contrast          estimate   SE df t.ratio p.value
Group 1 - Group 2   -27.59 2.64 54 -10.450 <.0001
Group 1 - Group 3    -5.71 2.80 54  -2.041 0.1123
Group 2 - Group 3    21.87 2.64 54   8.285 <.0001

gender = M:
contrast          estimate   SE df t.ratio p.value
Group 1 - Group 2   -25.87 2.15 54 -12.056 <.0001
Group 1 - Group 3    -8.00 2.05 54  -3.893 0.0008
Group 2 - Group 3    17.87 2.15 54   8.328 <.0001

P value adjustment: tukey method for comparing a family of 3 estimates 

But we can set other values using the adjust argument:

pairs(factorial_emm, adjust = "sidak") 
gender = F:
contrast          estimate   SE df t.ratio p.value
Group 1 - Group 2   -27.59 2.64 54 -10.450 <.0001
Group 1 - Group 3    -5.71 2.80 54  -2.041 0.1323
Group 2 - Group 3    21.87 2.64 54   8.285 <.0001

gender = M:
contrast          estimate   SE df t.ratio p.value
Group 1 - Group 2   -25.87 2.15 54 -12.056 <.0001
Group 1 - Group 3    -8.00 2.05 54  -3.893 0.0008
Group 2 - Group 3    17.87 2.15 54   8.328 <.0001

P value adjustment: sidak method for 3 tests 

The full list is :

sidak Sidak’s
tukey Tukey’s
scheffe Scheffe’s
dunnettx Dunnett’s
mvt Multivariate t distribution
holm Holm’s
hochberg Hochberg’s
hommel Hommel’s
bonferroni Bonferroni
fdr Benjamini & Hochberg
none none

The output of the pairs() function will also default to giving p values. But if you don’t want p values and would prefer confidence intervals then you can just wrap the call to pairs() in a call to confint().

factorial_pairwise <- pairs(factorial_emm, # the EMMs from our factorial ANOVA
)
confint(factorial_pairwise, # our factorial pairwise comparisons
level = .90 # the confdience level we want / default is 95%
)
gender = F:
contrast          estimate   SE df lower.CL upper.CL
Group 1 - Group 2   -27.59 2.64 54    -33.4  -21.822
Group 1 - Group 3    -5.71 2.80 54    -11.8    0.401
Group 2 - Group 3    21.87 2.64 54     16.1   27.639

gender = M:
contrast          estimate   SE df lower.CL upper.CL
Group 1 - Group 2   -25.87 2.15 54    -30.6  -21.187
Group 1 - Group 3    -8.00 2.05 54    -12.5   -3.513
Group 2 - Group 3    17.87 2.15 54     13.2   22.561

Confidence level used: 0.9
Conf-level adjustment: bonferroni method for 3 estimates 

Or we could %>% the output from pairs() into confint()

factorial_emm %>% pairs() %>% confint()
gender = F:
contrast          estimate   SE df lower.CL upper.CL
Group 1 - Group 2   -27.59 2.64 54    -33.9   -21.22
Group 1 - Group 3    -5.71 2.80 54    -12.5     1.03
Group 2 - Group 3    21.87 2.64 54     15.5    28.24

gender = M:
contrast          estimate   SE df lower.CL upper.CL
Group 1 - Group 2   -25.87 2.15 54    -31.0   -20.70
Group 1 - Group 3    -8.00 2.05 54    -13.0    -3.05
Group 2 - Group 3    17.87 2.15 54     12.7    23.05

Confidence level used: 0.95
Conf-level adjustment: tukey method for comparing a family of 3 estimates 

You’ll notice that in our comparisons for factorial_emm, we have the pairwise comparisons for each group separately for each gender. However, we might for example be interested in comparing the gender difference across the different groups. For example, is the M - F difference different between Group 1 and Group 2, Group 2 and Group 3, and Group 1 and Group 3. To do this, we can just add another call to the pairs() function and use the simple argument.

Let’s do that, but we’ll do it in one step with %>% pipes:

aov_factorial %>%  # our ANOVA object
emmeans(spec = "gender", by = "group") %>% # condition on group
pairs() %>% # pair wise comparisons
pairs(simple= "group") # simple contrasts for group
contrast = F - M:
contrast1         estimate   SE df t.ratio p.value
Group 1 - Group 2    -1.71 3.40 54 -0.504  0.8699
Group 1 - Group 3     2.29 3.47 54  0.658  0.7886
Group 2 - Group 3     4.00 3.40 54  1.175  0.4729

P value adjustment: tukey method for comparing a family of 3 estimates 

# Plotting ANOVAs

Usually the next thing you’d be interested in doing is looking at some plots. The afex package provides an easy interface for plotting ANOVAs. There’s lot of options you can change for plotting, so if you want to know more there’s a lot of info on the afex website, but I’ll just cover a few basics.

The main workhorse for plotting ANOVAs with afex is the afex_plot() function. To generate a plot all you need to do is feed it the output from a call to aov_ez(). We’ll start with something simple, our one-way ANOVA which we have saved in the object aov_oneway.

The only thing we need to tell afex_plot() apart from the name of our ANOVA object is the name of the factor we want on the x-axis. For our one-way ANOVA this will just be group.

afex:::afex_plot(object = aov_oneway, # our anova object
x = "group" # what we want on the x axis
)

With factorial ANOVAs there’s lot of ways in which we can plot them, so afex_plot() gives us this flexibility. Here are just a few of the possible options.

Gender along the x-axis with group in different panels (panel argument).

afex_plot(object = aov_factorial, # our anova object
x = "gender", # what we want on the x-axis
panel = "group" # what we want in different panels
)

Group along the x-axis with lines joining the gender categories (trace argument).

afex_plot(object = aov_factorial, # the name of our anova object
x = "group", # what we want on the x-axis
trace = "gender" # the factor we want joined by lines
)

Gender along the x-axis with lines joining the group categories.

afex_plot(object = aov_factorial,
x = "gender",
trace = "group"
)

We can also make adjustments to our confidence intervals. For example, our repeated measures ANOVA (aov_repeated) has within subject factors and a between subjects factor. So we have to make a decision between:

1. within subject error bars (based on the Cosineau-Morey-O’Brien method)

2. between subject error bars where data is first aggregated per participant

3. mean error bars based on the cell means (ignoring any repeated measures)

4. none for removing the error bars

afex_plot(object = aov_repeated,
x = "emotion",
panel = "congruency",
trace = "group",
error = "within")
Warning: Panel(s) show a mixed within-between-design.
Error bars do not allow comparisons across all means.
Suppress error bars with: error = "none"

## Styling our plots

There are also many options of styling options for our plots.

By default, the raw data is shown on the plot, but we can turn that off using the data_plot argument.

afex_plot(object = aov_factorial,
x = "gender",
trace = "group",
data_plot = FALSE)

The x-axis points are dodged by default, but that can also be turned off using the dodge argument

afex_plot(object = aov_factorial,
x = "gender",
trace = "group",
data_plot = FALSE,
dodge = FALSE)

And the data are shown as points, but we could replace it with a box plot using the data_geom argument.

afex_plot(object = aov_factorial,
x = "gender",
trace = "group",
data_plot = FALSE,
dodge = FALSE,
data_geom = geom_boxplot
)

1. The summary() function is what is known as a generic. This means that it works with lots of different object not only the object produced by aov_ez(). We’ll see that we can also use it for objects produced from the psych::mediate() function (for mediation analysis) or the stats::lm() function (for regression analysis).↩︎