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() function¹.

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 :

adjust	method
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
                            adjust = "bonferroni" # bonferroni adjustment
                            ) 
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:

within subject error bars (based on the Cosineau-Morey-O’Brien method)
between subject error bars where data is first aggregated per participant
mean error bars based on the cell means (ignoring any repeated measures)
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
)

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).↩︎

Analysis of Variance with afex