×

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.

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 `group`

s. 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
```

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
```

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^{1}.

`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
```

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
```

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
```

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"
```

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