aov_car {afex} | R Documentation |
Convenient ANOVA estimation for factorial designs
Description
These functions allow convenient specification of any type of ANOVAs (i.e.,
purely within-subjects ANOVAs, purely between-subjects ANOVAs, and mixed
between-within or split-plot ANOVAs) for data in the long format
(i.e., one observation per row). If the data has more than one observation
per individual and cell of the design (e.g., multiple responses per
condition), the data will be automatically aggregated. The default settings
reproduce results from commercial statistical packages such as SPSS or SAS.
aov_ez
is called specifying the factors as character vectors,
aov_car
is called using a formula similar to aov
specifying an error strata for the within-subject factor(s), and aov_4
is called with a lme4-like formula (all ANOVA functions return
identical results). The returned object can be passed to e.g., emmeans
for further analysis (e.g., follow-up tests, contrasts, plotting, etc.).
These functions employ Anova
(from the car package)
to provide test of effects avoiding the somewhat unhandy format of
car::Anova
.
Usage
aov_car(
formula,
data,
fun_aggregate = NULL,
type = afex_options("type"),
factorize = afex_options("factorize"),
observed = NULL,
anova_table = list(),
include_aov = afex_options("include_aov"),
return = afex_options("return_aov"),
...
)
aov_4(
formula,
data,
observed = NULL,
fun_aggregate = NULL,
type = afex_options("type"),
factorize = afex_options("factorize"),
return = afex_options("return_aov"),
anova_table = list(),
include_aov = afex_options("include_aov"),
...,
print.formula = FALSE
)
aov_ez(
id,
dv,
data,
between = NULL,
within = NULL,
covariate = NULL,
observed = NULL,
fun_aggregate = NULL,
transformation,
type = afex_options("type"),
factorize = afex_options("factorize"),
return = afex_options("return_aov"),
anova_table = list(),
include_aov = afex_options("include_aov"),
...,
print.formula = FALSE
)
Arguments
formula |
A formula specifying the ANOVA model similar to
|
data |
A |
fun_aggregate |
The function for aggregating the data before running the
ANOVA if there is more than one observation per individual and cell of the
design. The default |
type |
The type of sums of squares for the ANOVA. The default is given
by |
factorize |
logical. Should between subject factors be factorized (with
note) before running the analysis. The default is given by
|
observed |
|
anova_table |
|
include_aov |
Boolean. Allows suppressing the calculation of the aov
object. If TRUE the aov model is part of the returned |
return |
What should be returned? The default is given by
|
... |
Further arguments passed to |
print.formula |
|
id |
|
dv |
|
between |
|
within |
|
covariate |
|
transformation |
In |
Details
Details of ANOVA Specification
aov_ez
will concatenate
all between-subject factors using *
(i.e., producing all main effects
and interactions) and all covariates by +
(i.e., adding only the main
effects to the existing between-subject factors). The within-subject factors
do fully interact with all between-subject factors and covariates. This is
essentially identical to the behavior of SPSS's glm
function.
The formula
s for aov_car
or aov_4
must contain a single
Error
term specifying the ID
column and potential
within-subject factors (you can use mixed
for running
mixed-effects models with multiple error terms). Factors outside the
Error
term are treated as between-subject factors (the within-subject
factors specified in the Error
term are ignored outside the
Error
term; in other words, it is not necessary to specify them
outside the Error
term, see Examples).
Suppressing the intercept
(i.e, via 0 +
or - 1
) is ignored. Specific specifications of
effects (e.g., excluding terms with -
or using ^
) could be okay
but is not tested. Using the I
or poly
function
within the formula is not tested and not supported!
To run an ANCOVA you need to set factorize = FALSE
and make sure that
all variables have the correct type (i.e., factors are factors and numeric
variables are numeric and centered).
Note that the default behavior is to include calculation of the effect size
generalized eta-squared for which all non-manipluated (i.e.,
observed) variables need to be specified via the observed
argument to
obtain correct results. When changing the effect size to "pes"
(partial eta-squared) or "none"
via anova_table
this becomes
unnecessary.
Factor contrasts will be set to "contr.sum"
for all between-subject
factors if default contrasts are not equal to "contr.sum"
or
attrib(factor, "contrasts") != "contr.sum"
. (within-subject factors
are hard-coded "contr.sum"
.)
Statistical Issues
Type 3 sums of squares are default in afex. While some authors argue that so-called type 3 sums of squares are dangerous and/or problematic (most notably Venables, 2000), they are the default in many commercial statistical application such as SPSS or SAS. Furthermore, statisticians with an applied perspective recommend type 3 tests (e.g., Maxwell and Delaney, 2004). Consequently, they are the default for the ANOVA functions described here. For some more discussion on this issue see here.
Note that lower order effects (e.g., main effects) in type 3 ANOVAs are only
meaningful with
effects
coding. Therefore, contrasts are set to contr.sum
which
ensures meaningful results. For a discussion of the other (non-recommended)
coding schemes see
here.
Follow-Up Contrasts and Post-Hoc Tests
The S3 object returned
per default can be directly passed to emmeans::emmeans
for further
analysis. This allows to test any type of contrasts that might be of interest
independent of whether or not this contrast involves between-subject
variables, within-subject variables, or a combination thereof. The general
procedure to run those contrasts is the following (see Examples for a full
example):
Estimate an
afex_aov
object with the function returned here. For example:x <- aov_car(dv ~ a*b + (id/c), d)
Obtain a
emmGrid-class
object by runningemmeans
on theafex_aov
object from step 1 using the factors involved in the contrast. For example:r <- emmeans(x, ~a:c)
Create a list containing the desired contrasts on the reference grid object from step 2. For example:
con1 <- list(a_x = c(-1, 1, 0, 0, 0, 0), b_x = c(0, 0, -0.5, -0.5, 0, 1))
Test the contrast on the reference grid using
contrast
. For example:contrast(r, con1)
To control for multiple testing p-value adjustments can be specified. For example the Bonferroni-Holm correction:
contrast(r, con1, adjust = "holm")
Note that emmeans allows for a variety of advanced settings and
simplifications, for example: all pairwise comparison of a single factor
using one command (e.g., emmeans(x, "a", contr = "pairwise")
) or
advanced control for multiple testing by passing objects to multcomp.
A comprehensive overview of the functionality is provided in the
accompanying vignettes (see
here).
Since version 1.0, afex per default uses the multivariate
model
(i.e., the lm
slot of the afex_aov
object) for follow-up tests
with emmeans. Compared to the univariate
model (i.e., the
aov
slot), this can handle unbalanced data and addresses sphericity
better. To use the older (and not recommended) model = "univariate"
make sure to set include_aov = TRUE
when estimating the ANOVA.
Starting with afex version 0.22, emmeans is not
loaded/attached automatically when loading afex. Therefore,
emmeans now needs to be loaded by the user via
library("emmeans")
or require("emmeans")
.
Methods for afex_aov
Objects
A full overview over the
methods provided for afex_aov
objects is provided in the corresponding
help page: afex_aov-methods
. The probably most important ones
for end-users are summary
, anova
, and nice
.
The summary
method returns, for ANOVAs containing within-subject
(repeated-measures) factors with more than two levels, the complete
univariate analysis: Results without df-correction, the Greenhouse-Geisser
corrected results, the Hyunh-Feldt corrected results, and the results of the
Mauchly test for sphericity.
The anova
method returns a data.frame
of class "anova"
containing the ANOVA table in numeric form (i.e., the one in slot
anova_table
of a afex_aov
). This method has arguments such as
correction
and es
and can be used to obtain an ANOVA table with
different correction than the one initially specified.
The nice
method also returns a data.frame
, but rounds
most values and transforms them into characters for nice printing. Also has
arguments like correction
and es
which can be used to obtain an
ANOVA table with different correction than the one initially specified.
Value
aov_car
, aov_4
, and aov_ez
are wrappers for
Anova
and aov
, the return value is
dependent on the return
argument. Per default, an S3 object of class
"afex_aov"
is returned containing the following slots:
"anova_table"
An ANOVA table of class
c("anova", "data.frame")
."aov"
aov
object returned fromaov
(should not be used to evaluate significance of effects, but can be passed toemmeans
for post-hoc tests)."Anova"
object returned from
Anova
, an object of class"Anova.mlm"
(if within-subjects factors are present) or of classc("anova", "data.frame")
."lm"
the object fitted with
lm
and passed toAnova
(i.e., an object of class"lm"
or"mlm"
). Also returned ifreturn = "lm"
."data"
a list containing: (1)
long
(the possibly aggregated data in long format used foraov
),wide
(the data used to fit thelm
object), andidata
(if within-subject factors are present, theidata
argument passed tocar::Anova
). Also returned ifreturn = "data"
.
In addition, the object has the following attributes: "dv"
,
"id"
, "within"
, "between"
, and "type"
.
The print method for afex_aov
objects
(invisibly) returns (and prints) the same as if return
is
"nice"
: a nice ANOVA table (produced by nice
) with the
following columns: Effect
, df
, MSE
(mean-squared
errors), F
(potentially with significant symbols), ges
(generalized eta-squared), p
.
Functions
-
aov_4()
: Allows definition of ANOVA-model usinglme4::lmer
-like Syntax (but still fits a standard ANOVA). -
aov_ez()
: Allows definition of ANOVA-model using character strings.
Note
Calculation of ANOVA models via aov
(which is done per default)
can be comparatively slow and produce comparatively large objects for
ANOVAs with many within-subjects factors or levels. To avoid this
calculation set include_aov = FALSE
. You can also disable this
globally with: afex_options(include_aov = FALSE)
The id variable and variables entered as within-subjects (i.e.,
repeated-measures) factors are silently converted to factors. Levels of
within-subject factors are converted to valid variable names using
make.names(...,unique=TRUE)
. Unused factor levels are
silently dropped on all variables.
Contrasts attached to a factor as an attribute are probably not preserved and not supported.
The workhorse is aov_car
. aov_4
and aov_ez
only
construe and pass an appropriate formula to aov_car
. Use
print.formula = TRUE
to view this formula.
In contrast to aov
aov_car
assumes that all factors to
the right of /
in the Error
term are belonging together.
Consequently, Error(id/(a*b))
and Error(id/a*b)
are identical
(which is not true for aov
).
Author(s)
Henrik Singmann
The design of these functions was influenced by ezANOVA
from package ez.
References
Cramer, A. O. J., van Ravenzwaaij, D., Matzke, D., Steingroever, H., Wetzels, R., Grasman, R. P. P. P., ... Wagenmakers, E.-J. (2015). Hidden multiplicity in exploratory multiway ANOVA: Prevalence and remedies. Psychonomic Bulletin & Review, 1-8. doi:10.3758/s13423-015-0913-5
Maxwell, S. E., & Delaney, H. D. (2004). Designing Experiments and Analyzing Data: A Model-Comparisons Perspective. Mahwah, N.J.: Lawrence Erlbaum Associates.
Venables, W.N. (2000). Exegeses on linear models. Paper presented to the S-Plus User's Conference, Washington DC, 8-9 October 1998, Washington, DC. Available from: http://www.stats.ox.ac.uk/pub/MASS3/Exegeses.pdf
See Also
Various methods for objects of class afex_aov
are available:
afex_aov-methods
nice
creates the nice ANOVA tables which is by default printed.
See also there for a slightly longer discussion of the available effect
sizes.
mixed
provides a (formula) interface for obtaining p-values for
mixed-models via lme4. The functions presented here do not estimate
mixed models.
Examples
##########################
## 1: Specifying ANOVAs ##
##########################
# Example using a purely within-subjects design
# (Maxwell & Delaney, 2004, Chapter 12, Table 12.5, p. 578):
data(md_12.1)
aov_ez("id", "rt", md_12.1, within = c("angle", "noise"),
anova_table=list(correction = "none", es = "none"))
# Default output
aov_ez("id", "rt", md_12.1, within = c("angle", "noise"))
# examples using obk.long (see ?obk.long), a long version of the OBrienKaiser
# dataset (car package). Data is a split-plot or mixed design: contains both
# within- and between-subjects factors.
data(obk.long, package = "afex")
# estimate mixed ANOVA on the full design:
aov_car(value ~ treatment * gender + Error(id/(phase*hour)),
data = obk.long, observed = "gender")
aov_4(value ~ treatment * gender + (phase*hour|id),
data = obk.long, observed = "gender")
aov_ez("id", "value", obk.long, between = c("treatment", "gender"),
within = c("phase", "hour"), observed = "gender")
# the three calls return the same ANOVA table:
# Anova Table (Type 3 tests)
#
# Response: value
# Effect df MSE F ges p.value
# 1 treatment 2, 10 22.81 3.94 + .198 .055
# 2 gender 1, 10 22.81 3.66 + .115 .085
# 3 treatment:gender 2, 10 22.81 2.86 .179 .104
# 4 phase 1.60, 15.99 5.02 16.13 *** .151 <.001
# 5 treatment:phase 3.20, 15.99 5.02 4.85 * .097 .013
# 6 gender:phase 1.60, 15.99 5.02 0.28 .003 .709
# 7 treatment:gender:phase 3.20, 15.99 5.02 0.64 .014 .612
# 8 hour 1.84, 18.41 3.39 16.69 *** .125 <.001
# 9 treatment:hour 3.68, 18.41 3.39 0.09 .002 .979
# 10 gender:hour 1.84, 18.41 3.39 0.45 .004 .628
# 11 treatment:gender:hour 3.68, 18.41 3.39 0.62 .011 .641
# 12 phase:hour 3.60, 35.96 2.67 1.18 .015 .335
# 13 treatment:phase:hour 7.19, 35.96 2.67 0.35 .009 .930
# 14 gender:phase:hour 3.60, 35.96 2.67 0.93 .012 .449
# 15 treatment:gender:phase:hour 7.19, 35.96 2.67 0.74 .019 .646
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1
#
# Sphericity correction method: GG
# "numeric" variables are per default converted to factors (as long as factorize
# = TRUE):
obk.long$hour2 <- as.numeric(as.character(obk.long$hour))
# gives same results as calls before
aov_car(value ~ treatment * gender + Error(id/phase*hour2),
data = obk.long, observed = c("gender"))
# ANCOVA: adding a covariate (necessary to set factorize = FALSE)
aov_car(value ~ treatment * gender + age + Error(id/(phase*hour)),
data = obk.long, observed = c("gender", "age"), factorize = FALSE)
aov_4(value ~ treatment * gender + age + (phase*hour|id),
data = obk.long, observed = c("gender", "age"), factorize = FALSE)
aov_ez("id", "value", obk.long, between = c("treatment", "gender"),
within = c("phase", "hour"), covariate = "age",
observed = c("gender", "age"), factorize = FALSE)
# aggregating over one within-subjects factor (phase), with warning:
aov_car(value ~ treatment * gender + Error(id/hour), data = obk.long,
observed = "gender")
aov_ez("id", "value", obk.long, c("treatment", "gender"), "hour",
observed = "gender")
# aggregating over both within-subjects factors (again with warning),
# only between-subjects factors:
aov_car(value ~ treatment * gender + Error(id), data = obk.long,
observed = c("gender"))
aov_4(value ~ treatment * gender + (1|id), data = obk.long,
observed = c("gender"))
aov_ez("id", "value", obk.long, between = c("treatment", "gender"),
observed = "gender")
# only within-subject factors (ignoring between-subjects factors)
aov_car(value ~ Error(id/(phase*hour)), data = obk.long)
aov_4(value ~ (phase*hour|id), data = obk.long)
aov_ez("id", "value", obk.long, within = c("phase", "hour"))
### changing defaults of ANOVA table:
# no df-correction & partial eta-squared:
aov_car(value ~ treatment * gender + Error(id/(phase*hour)),
data = obk.long, anova_table = list(correction = "none", es = "pes"))
# no df-correction and no MSE
aov_car(value ~ treatment * gender + Error(id/(phase*hour)),
data = obk.long,observed = "gender",
anova_table = list(correction = "none", MSE = FALSE))
# add p-value adjustment for all effects (see Cramer et al., 2015, PB&R)
aov_ez("id", "value", obk.long, between = "treatment",
within = c("phase", "hour"),
anova_table = list(p_adjust_method = "holm"))
###########################
## 2: Follow-up Analysis ##
###########################
# use data as above
data(obk.long, package = "afex")
# 1. obtain afex_aov object:
a1 <- aov_ez("id", "value", obk.long, between = c("treatment", "gender"),
within = c("phase", "hour"), observed = "gender")
if (requireNamespace("ggplot2") & requireNamespace("emmeans")) {
# 1b. plot data using afex_plot function, for more see:
## vignette("afex_plot_introduction", package = "afex")
## default plot uses multivariate model-based CIs
afex_plot(a1, "hour", "gender", c("treatment", "phase"))
a1b <- aov_ez("id", "value", obk.long, between = c("treatment", "gender"),
within = c("phase", "hour"), observed = "gender",
include_aov = TRUE)
## you can use a univariate model and CIs if you refit the model with the aov
## slot
afex_plot(a1b, "hour", "gender", c("treatment", "phase"),
emmeans_arg = list(model = "univariate"))
## in a mixed between-within designs, no error-bars might be preferrable:
afex_plot(a1, "hour", "gender", c("treatment", "phase"), error = "none")
}
if (requireNamespace("emmeans")) {
library("emmeans") # package emmeans needs to be attached for follow-up tests.
# 2. obtain reference grid object (default uses multivariate model):
r1 <- emmeans(a1, ~treatment +phase)
r1
# 3. create list of contrasts on the reference grid:
c1 <- list(
A_B_pre = c(rep(0, 6), 0, -1, 1), # A versus B for pretest
A_B_comb = c(-0.5, 0.5, 0, -0.5, 0.5, 0, 0, 0, 0), # A vs. B for post and follow-up combined
effect_post = c(0, 0, 0, -1, 0.5, 0.5, 0, 0, 0), # control versus A&B post
effect_fup = c(-1, 0.5, 0.5, 0, 0, 0, 0, 0, 0), # control versus A&B follow-up
effect_comb = c(-0.5, 0.25, 0.25, -0.5, 0.25, 0.25, 0, 0, 0) # control versus A&B combined
)
# 4. test contrasts on reference grid:
contrast(r1, c1)
# same as before, but using Bonferroni-Holm correction for multiple testing:
contrast(r1, c1, adjust = "holm")
# 2. (alternative): all pairwise comparisons of treatment:
emmeans(a1, "treatment", contr = "pairwise")
}
#######################
## 3: Other examples ##
#######################
data(obk.long, package = "afex")
# replicating ?Anova using aov_car:
obk_anova <- aov_car(value ~ treatment * gender + Error(id/(phase*hour)),
data = obk.long, type = 2)
# in contrast to aov you do not need the within-subject factors outside Error()
str(obk_anova, 1, give.attr = FALSE)
# List of 5
# $ anova_table:Classes ‘anova’ and 'data.frame': 15 obs. of 6 variables:
# $ aov : NULL
# $ Anova :List of 14
# $ lm :List of 13
# $ data :List of 3
obk_anova$Anova
# Type II Repeated Measures MANOVA Tests: Pillai test statistic
# Df test stat approx F num Df den Df Pr(>F)
# (Intercept) 1 0.96954 318.34 1 10 6.532e-09 ***
# treatment 2 0.48092 4.63 2 10 0.0376868 *
# gender 1 0.20356 2.56 1 10 0.1409735
# treatment:gender 2 0.36350 2.86 2 10 0.1044692
# phase 1 0.85052 25.61 2 9 0.0001930 ***
# treatment:phase 2 0.68518 2.61 4 20 0.0667354 .
# gender:phase 1 0.04314 0.20 2 9 0.8199968
# treatment:gender:phase 2 0.31060 0.92 4 20 0.4721498
# hour 1 0.93468 25.04 4 7 0.0003043 ***
# treatment:hour 2 0.30144 0.35 8 16 0.9295212
# gender:hour 1 0.29274 0.72 4 7 0.6023742
# treatment:gender:hour 2 0.57022 0.80 8 16 0.6131884
# phase:hour 1 0.54958 0.46 8 3 0.8324517
# treatment:phase:hour 2 0.66367 0.25 16 8 0.9914415
# gender:phase:hour 1 0.69505 0.85 8 3 0.6202076
# treatment:gender:phase:hour 2 0.79277 0.33 16 8 0.9723693
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1