<!-- README.md is generated from README.Rmd. Please edit that file -->

afex: Analysis of Factorial EXperiments

<!-- badges: start -->

<!-- badges: end -->

The main functionalities provided by afex are:

1. Interfaces for estimating standard ANOVAs with any number or combination of within-subjects or between-subjects variables (the ANOVA functions are aov_car(), aov_ez(), and aov_4() which all fit the same model but differ in the way to specify the ANOVA model).
2. Function mixed() provides an interface for mixed models analysis (estimated via lme4 lmer or glmer) that automatically obtains p-values for fixed effects model terms (i.e., main effects and interactions).
3. afex_plot() visualizes results from factorial experiments combining estimated marginal means and uncertainties associated with the estimated means in the foreground with a depiction of the raw data in the background.
4. All afex model objects (i.e., ANOVA and mixed models) can be passed to emmeans for follow-up/post-hoc/planned contrast analysis.

For afex support visit: afex.singmann.science

Installation

• afex is available from CRAN so the current stable version can be installed directly via: install.packages("afex")

• To install the latest development version you will need the devtools package: devtools::install_github("singmann/afex@master")

ANOVA functionality

To calculate an ANOVA, afex requires the data to be in the long format (i.e., one row per data point/observation). An ANOVA can then be calculated via one of three functions that only differ in how the model components are specified, but not in the output. Note that in contrast to base lm or aov, afex ANOVA functions always require the specification of a subject identifier column (the id-column), because in case there are multiple observations per participant and cell of the design, these multiple observations are aggregated (i.e., averaged) per default.

• In aov_ez the columns containing id variable, dependent variable, and factors need to be specified as character vectors.
• aov_car behaves similar to standard aov and requires the ANOVA to be specified as a formula containing an Error term (at least to identify the id variable).
• aov_4 allows the ANOVA to be specified via a formula similar to lme4::lmer (with one random effects term).

A further overview is provided by the vignette.

The following code provides a simple example for an ANOVA with both between- and within-subject factors. For this we use the lexical-decision and word naming latencies reported by Freeman, Heathcote, Chalmers, and Hockley (2010), see also ?fhch2010. As is commonly done, we use the natural logarithm of the response times, log_rt, as dependent variable. As independent variable we will consider the between-subjects factor task ("naming" or "lexdec") as well as the within-subjects-factors stimulus ("word" or "nonword") and length (with 3 levels, 3, 4, or 5 letters).

library("afex")
# examples data set with both within- and between-subjects factors (see ?fhch2010)
data("fhch2010", package = "afex")
fhch <- fhch2010[ fhch2010$correct,] # remove errors
str(fhch2010) # structure of the data
#> 'data.frame':    13222 obs. of  10 variables:
#>  $ id       : Factor w/ 45 levels "N1","N12","N13",..: 1 1 1 1 1 1 1 1 1 1 ...
#>  $ task     : Factor w/ 2 levels "naming","lexdec": 1 1 1 1 1 1 1 1 1 1 ...
#>  $ stimulus : Factor w/ 2 levels "word","nonword": 1 1 1 2 2 1 2 2 1 2 ...
#>  $ density  : Factor w/ 2 levels "low","high": 2 1 1 2 1 2 1 1 1 1 ...
#>  $ frequency: Factor w/ 2 levels "low","high": 1 2 2 2 2 2 1 2 1 2 ...
#>  $ length   : Factor w/ 3 levels "4","5","6": 3 3 2 2 1 1 3 2 1 3 ...
#>  $ item     : Factor w/ 600 levels "abide","acts",..: 363 121 202 525 580 135 42 368 227 141 ...
#>  $ rt       : num  1.091 0.876 0.71 1.21 0.843 ...
#>  $ log_rt   : num  0.0871 -0.1324 -0.3425 0.1906 -0.1708 ...
#>  $ correct  : logi  TRUE TRUE TRUE TRUE TRUE TRUE ...

# estimate mixed ANOVA on the full design:
aov_ez("id", "log_rt", fhch, between = "task", within = c("stimulus", "length"))
aov_car(log_rt ~ task * stimulus * length + Error(id/(stimulus * length)), 
        data = fhch)
## equivalent: aov_car(log_rt ~ task + Error(id/(stimulus * length)), data = fhch)

aov_4(log_rt ~ task * stimulus * length + (stimulus * length|id), data = fhch)
## equivalent: aov_4(log_rt ~ task  + (stimulus * length|id), data = fhch)

# the three calls return the same ANOVA table:

#> Warning: More than one observation per design cell, aggregating data using `fun_aggregate = mean`.
#> To turn off this warning, pass `fun_aggregate = mean` explicitly.
#> Contrasts set to contr.sum for the following variables: task
#> Anova Table (Type 3 tests)
#> 
#> Response: log_rt
#>                 Effect          df  MSE          F   ges p.value
#> 1                 task       1, 43 0.23  13.38 ***  .221   <.001
#> 2             stimulus       1, 43 0.01 173.25 ***  .173   <.001
#> 3        task:stimulus       1, 43 0.01  87.56 ***  .096   <.001
#> 4               length 1.83, 78.64 0.00  18.55 ***  .008   <.001
#> 5          task:length 1.83, 78.64 0.00       1.02 <.001    .358
#> 6      stimulus:length 1.70, 72.97 0.00       1.91 <.001    .162
#> 7 task:stimulus:length 1.70, 72.97 0.00       1.21 <.001    .298
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
#> 
#> Sphericity correction method: GG

Plotting with afex_plot

ANOVA models can be used for plotting via afex_plot:

a <- aov_ez("id", "log_rt", fhch, between = "task", within = c("stimulus", "length"))
#> Warning: More than one observation per design cell, aggregating data using `fun_aggregate = mean`.
#> To turn off this warning, pass `fun_aggregate = mean` explicitly.

afex_plot(a, "task", "stimulus", "length")
#> Warning: Panel(s) show a mixed within-between-design.
#> Error bars do not allow comparisons across all means.
#> Suppress error bars with: error = "none"

<!-- -->

afex_plot returns a ggplot2 plot object which allows simple customization:

library("ggplot2")
afex_plot(a, "task", "stimulus", "length") + 
  theme_bw()
#> Warning: Panel(s) show a mixed within-between-design.
#> Error bars do not allow comparisons across all means.
#> Suppress error bars with: error = "none"

<!-- -->

Follow-up Tests with emmeans

Follow-up tests with emmeans need to be specified in two steps.

1. Decide which factors of model should be involved in tests. Use these factors to set-up reference grid of marginal means using emmeans().
2. Specify set of tests on reference grid from step 1. Either custom contrasts as a list and using contrast() or a convenience function such as pairs().

suppressPackageStartupMessages(library("emmeans"))
## set up reference grid using only length
em1 <- emmeans(a, "length")
em1
#>  length  emmean     SE df lower.CL upper.CL
#>  X4     -0.1087 0.0299 43   -0.169 -0.04834
#>  X5     -0.0929 0.0296 43   -0.153 -0.03310
#>  X6     -0.0653 0.0290 43   -0.124 -0.00679
#> 
#> Results are averaged over the levels of: task, stimulus 
#> Confidence level used: 0.95

## test all pairwise comparisons on reference grid:
pairs(em1)
#>  contrast estimate      SE df t.ratio p.value
#>  X4 - X5   -0.0159 0.00768 43  -2.065  0.1092
#>  X4 - X6   -0.0434 0.00782 43  -5.555  <.0001
#>  X5 - X6   -0.0276 0.00602 43  -4.583  0.0001
#> 
#> Results are averaged over the levels of: task, stimulus 
#> P value adjustment: tukey method for comparing a family of 3 estimates

## only test specified tests
con <- list(
  "4vs5" = c(-1, 1, 0),
  "5vs6" = c(0, -1, 1)
)
contrast(em1, con, adjust = "holm")
#>  contrast estimate      SE df t.ratio p.value
#>  4vs5       0.0159 0.00768 43   2.065  0.0449
#>  5vs6       0.0276 0.00602 43   4.583  0.0001
#> 
#> Results are averaged over the levels of: task, stimulus 
#> P value adjustment: holm method for 2 tests

Mixed Models

Function mixed() fits a mixed model with lme4::lmer (or lme4::glmer if a family argument is passed) and then calculates p-values for fixed effects model terms using a variety of methods. The formula to mixed needs to be the same as in a call to lme4::lmer. The default method for calculation of p-values is 'S' (Satterthwaite) which only works for linear mixed models (i.e., no family argument). A similar method that provides a somewhat better control of Type I errors for small data sets is 'KR' (Kenward-Roger), but it can require considerable RAM and time. Other methods are , similar to 'KR' but requires less RAM), 'PB' (parametric bootstrap), and 'LRT' (likelihood-ratio test).

More examples are provided in the vignette, here we use the same example data as above, the lexical decision and word naming latencies collected by Freeman et al. (2010). To avoid long computation times we only consider the two factors task and length (omitting stimulus is probably not a fully sensible model). Because mixed models easily allow