# Comparing Multiple Means in R

## Repeated Measures ANOVA in R

The repeated-measures ANOVA is used for analyzing data where same subjects are measured more than once. This test is also referred to as a within-subjects ANOVA or ANOVA with repeated measures. The “within-subjects” term means that the same individuals are measured on the same outcome variable under different time points or conditions.

For example, you might have measured 10 individuals’ self-esteem score (the outcome or dependent variable) on three time points during a specific diet to determine whether their self-esteem improved.

This chapter describes the different types of repeated measures ANOVA, including:

• One-way repeated measures ANOVA, an extension of the paired-samples t-test for comparing the means of three or more levels of a within-subjects variable.
• two-way repeated measures ANOVA used to evaluate simultaneously the effect of two within-subject factors on a continuous outcome variable.
• three-way repeated measures ANOVA used to evaluate simultaneously the effect of three within-subject factors on a continuous outcome variable.

The main goal of two-way and three-way repeated measures ANOVA is, respectively, to evaluate if there is a statistically significant interaction effect between two and three within-subjects factors in explaining a continuous outcome variable.

You will learn how to:

• Compute and interpret the different repeated measures ANOVA in R.
• Check repeated measures ANOVA test assumptions
• Perform post-hoc tests, multiple pairwise comparisons between groups to identify which groups are different
• Visualize the data using box plots, add ANOVA and pairwise comparisons p-values to the plot

Contents:

#### Related Book

Practical Statistics in R II - Comparing Groups: Numerical Variables

## Assumptions

The repeated measures ANOVA makes the following assumptions about the data:

• No significant outliers in any cell of the design. This can be checked by visualizing the data using box plot methods and by using the function identify_outliers() [rstatix package].
• Normality: the outcome (or dependent) variable should be approximately normally distributed in each cell of the design. This can be checked using the Shapiro-Wilk normality test (shapiro_test() [rstatix]) or by visual inspection using QQ plot (ggqqplot() [ggpubr package]).
• Assumption of sphericity: the variance of the differences between groups should be equal. This can be checked using the Mauchly’s test of sphericity, which is automatically reported when using the R function anova_test() [rstatix package]. Read more in Chapter @ref(mauchly-s-test-of-sphericity-in-r).

Before computing repeated measures ANOVA test, you need to perform some preliminary tests to check if the assumptions are met.

Note that, if the above assumptions are not met there are a non-parametric alternative (Friedman test) to the one-way repeated measures ANOVA.

Unfortunately, there are no non-parametric alternatives to the two-way and the three-way repeated measures ANOVA. Thus, in the situation where the assumptions are not met, you could consider running the two-way/three-way repeated measures ANOVA on the transformed and non-transformed data to see if there are any meaningful differences.

If both tests lead you to the same conclusions, you might not choose to transform the outcome variable and carry on with the two-way/three-way repeated measures ANOVA on the original data.

It’s also possible to perform robust ANOVA test using the WRS2 R package.

No matter your choice, you should report what you did in your results.

## Prerequisites

Make sure that you have installed the following R packages:

• tidyverse for data manipulation and visualization
• ggpubr for creating easily publication ready plots
• rstatix provides pipe-friendly R functions for easy statistical analyses
• datarium: contains required data sets for this chapter

library(tidyverse)
library(ggpubr)
library(rstatix)

Key R functions:

• anova_test() [rstatix package], a wrapper around car::Anova() for making easy the computation of repeated measures ANOVA. Key arguments for performing repeated measures ANOVA:
• data: data frame
• dv: (numeric) the dependent (or outcome) variable name.
• wid: variable name specifying the case/sample identifier.
• within: within-subjects factor or grouping variable
• get_anova_table() [rstatix package]. Extracts the ANOVA table from the output of anova_test(). It returns ANOVA table that is automatically corrected for eventual deviation from the sphericity assumption. The default is to apply automatically the Greenhouse-Geisser sphericity correction to only within-subject factors violating the sphericity assumption (i.e., Mauchly’s test p-value is significant, p <= 0.05). Read more in Chapter @ref(mauchly-s-test-of-sphericity-in-r).

## One-way repeated measures ANOVA

### Data preparation

We’ll use the self-esteem score dataset measured over three time points. The data is available in the datarium package.

# Data preparation
# Wide format
data("selfesteem", package = "datarium")
head(selfesteem, 3)
## # A tibble: 3 x 4
##      id    t1    t2    t3
##   <int> <dbl> <dbl> <dbl>
## 1     1  4.01  5.18  7.11
## 2     2  2.56  6.91  6.31
## 3     3  3.24  4.44  9.78
# Gather columns t1, t2 and t3 into long format
# Convert id and time into factor variables
selfesteem <- selfesteem %>%
gather(key = "time", value = "score", t1, t2, t3) %>%
convert_as_factor(id, time)
head(selfesteem, 3)
## # A tibble: 3 x 3
##   id    time  score
##   <fct> <fct> <dbl>
## 1 1     t1     4.01
## 2 2     t1     2.56
## 3 3     t1     3.24

The one-way repeated measures ANOVA can be used to determine whether the means self-esteem scores are significantly different between the three time points.

### Summary statistics

Compute some summary statistics of the self-esteem score by groups (time): mean and sd (standard deviation)

selfesteem %>%
group_by(time) %>%
get_summary_stats(score, type = "mean_sd")
## # A tibble: 3 x 5
##   time  variable     n  mean    sd
##   <fct> <chr>    <dbl> <dbl> <dbl>
## 1 t1    score       10  3.14 0.552
## 2 t2    score       10  4.93 0.863
## 3 t3    score       10  7.64 1.14

### Visualization

Create a box plot and add points corresponding to individual values:

bxp <- ggboxplot(selfesteem, x = "time", y = "score", add = "point")
bxp

### Check assumptions

#### Outliers

Outliers can be easily identified using box plot methods, implemented in the R function identify_outliers() [rstatix package].

selfesteem %>%
group_by(time) %>%
identify_outliers(score)
## # A tibble: 2 x 5
##   time  id    score is.outlier is.extreme
##   <fct> <fct> <dbl> <lgl>      <lgl>
## 1 t1    6      2.05 TRUE       FALSE
## 2 t2    2      6.91 TRUE       FALSE

There were no extreme outliers.

Note that, in the situation where you have extreme outliers, this can be due to: 1) data entry errors, measurement errors or unusual values.

You can include the outlier in the analysis anyway if you do not believe the result will be substantially affected. This can be evaluated by comparing the result of the ANOVA with and without the outlier.

It’s also possible to keep the outliers in the data and perform robust ANOVA test using the WRS2 package.

#### Normality assumption

The normality assumption can be checked by computing Shapiro-Wilk test for each time point. If the data is normally distributed, the p-value should be greater than 0.05.

selfesteem %>%
group_by(time) %>%
shapiro_test(score)
## # A tibble: 3 x 4
##   time  variable statistic     p
##   <fct> <chr>        <dbl> <dbl>
## 1 t1    score        0.967 0.859
## 2 t2    score        0.876 0.117
## 3 t3    score        0.923 0.380

The self-esteem score was normally distributed at each time point, as assessed by Shapiro-Wilk’s test (p > 0.05).

Note that, if your sample size is greater than 50, the normal QQ plot is preferred because at larger sample sizes the Shapiro-Wilk test becomes very sensitive even to a minor deviation from normality.

QQ plot draws the correlation between a given data and the normal distribution. Create QQ plots for each time point:

ggqqplot(selfesteem, "score", facet.by = "time")

From the plot above, as all the points fall approximately along the reference line, we can assume normality.

#### Assumption of sphericity

As mentioned in previous sections, the assumption of sphericity will be automatically checked during the computation of the ANOVA test using the R function anova_test() [rstatix package]. The Mauchly’s test is internally used to assess the sphericity assumption.

By using the function get_anova_table() [rstatix] to extract the ANOVA table, the Greenhouse-Geisser sphericity correction is automatically applied to factors violating the sphericity assumption.

### Computation

res.aov <- anova_test(data = selfesteem, dv = score, wid = id, within = time)
get_anova_table(res.aov)
## ANOVA Table (type III tests)
##
##   Effect DFn DFd    F        p p<.05   ges
## 1   time   2  18 55.5 2.01e-08     * 0.829

The self-esteem score was statistically significantly different at the different time points during the diet, F(2, 18) = 55.5, p < 0.0001, eta2[g] = 0.83.

where,

• F Indicates that we are comparing to an F-distribution (F-test); (2, 18) indicates the degrees of freedom in the numerator (DFn) and the denominator (DFd), respectively; 55.5 indicates the obtained F-statistic value
• p specifies the p-value
• ges is the generalized effect size (amount of variability due to the within-subjects factor)

### Post-hoc tests

You can perform multiple pairwise paired t-tests between the levels of the within-subjects factor (here time). P-values are adjusted using the Bonferroni multiple testing correction method.

# pairwise comparisons
pwc <- selfesteem %>%
pairwise_t_test(
score ~ time, paired = TRUE,
)
pwc
## # A tibble: 3 x 10
## * <chr> <chr>  <chr>  <int> <int>     <dbl> <dbl>       <dbl>    <dbl> <chr>
## 1 score t1     t2        10    10     -4.97     9 0.000772    0.002    **
## 2 score t1     t3        10    10    -13.2      9 0.000000334 0.000001 ****
## 3 score t2     t3        10    10     -4.87     9 0.000886    0.003    **

All the pairwise differences are statistically significant.

### Report

We could report the result as follow:

The self-esteem score was statistically significantly different at the different time points, F(2, 18) = 55.5, p < 0.0001, generalized eta squared = 0.82.

Post-hoc analyses with a Bonferroni adjustment revealed that all the pairwise differences, between time points, were statistically significantly different (p <= 0.05).

# Visualization: box plots with p-values
pwc <- pwc %>% add_xy_position(x = "time")
bxp +
stat_pvalue_manual(pwc) +
labs(
subtitle = get_test_label(res.aov, detailed = TRUE),
caption = get_pwc_label(pwc)
)

## Two-way repeated measures ANOVA

### Data preparation

We’ll use the selfesteem2 dataset [in datarium package] containing the self-esteem score measures of 12 individuals enrolled in 2 successive short-term trials (4 weeks): control (placebo) and special diet trials.

Each participant performed all two trials. The order of the trials was counterbalanced and sufficient time was allowed between trials to allow any effects of previous trials to have dissipated.

The self-esteem score was recorded at three time points: at the beginning (t1), midway (t2) and at the end (t3) of the trials.

The question is to investigate if this short-term diet treatment can induce a significant increase of self-esteem score over time. In other terms, we wish to know if there is significant interaction between diet and time on the self-esteem score.

The two-way repeated measures ANOVA can be performed in order to determine whether there is a significant interaction between diet and time on the self-esteem score.

Load and show one random row by treatment group:

# Wide format
set.seed(123)
data("selfesteem2", package = "datarium")
selfesteem2 %>% sample_n_by(treatment, size = 1)
## # A tibble: 2 x 5
##   id    treatment    t1    t2    t3
##   <fct> <fct>     <dbl> <dbl> <dbl>
## 1 4     ctr          92    92    89
## 2 10    Diet         90    93    95
# Gather the columns t1, t2 and t3 into long format.
# Convert id and time into factor variables
selfesteem2 <- selfesteem2 %>%
gather(key = "time", value = "score", t1, t2, t3) %>%
convert_as_factor(id, time)
# Inspect some random rows of the data by groups
set.seed(123)
selfesteem2 %>% sample_n_by(treatment, time, size = 1)
## # A tibble: 6 x 4
##   id    treatment time  score
##   <fct> <fct>     <fct> <dbl>
## 1 4     ctr       t1       92
## 2 10    ctr       t2       84
## 3 5     ctr       t3       68
## 4 11    Diet      t1       93
## 5 12    Diet      t2       80
## 6 1     Diet      t3       88

In this example, the effect of “time” on self-esteem score is our focal variable, that is our primary concern.

However, it is thought that the effect “time” will be different if treatment is performed or not. In this setting, the “treatment” variable is considered as moderator variable.

### Summary statistics

Group the data by treatment and time, and then compute some summary statistics of the score variable: mean and sd (standard deviation).

selfesteem2 %>%
group_by(treatment, time) %>%
get_summary_stats(score, type = "mean_sd")
## # A tibble: 6 x 6
##   treatment time  variable     n  mean    sd
##   <fct>     <fct> <chr>    <dbl> <dbl> <dbl>
## 1 ctr       t1    score       12  88    8.08
## 2 ctr       t2    score       12  83.8 10.2
## 3 ctr       t3    score       12  78.7 10.5
## 4 Diet      t1    score       12  87.6  7.62
## 5 Diet      t2    score       12  87.8  7.42
## 6 Diet      t3    score       12  87.7  8.14

### Visualization

Create box plots of the score colored by treatment groups:

bxp <- ggboxplot(
selfesteem2, x = "time", y = "score",
color = "treatment", palette = "jco"
)
bxp

### Chek assumptions

#### Outliers

selfesteem2 %>%
group_by(treatment, time) %>%
identify_outliers(score)
## [1] treatment  time       id         score      is.outlier is.extreme
## <0 rows> (or 0-length row.names)

There were no extreme outliers.

#### Normality assumption

Compute Shapiro-Wilk test for each combinations of factor levels:

selfesteem2 %>%
group_by(treatment, time) %>%
shapiro_test(score)
## # A tibble: 6 x 5
##   treatment time  variable statistic      p
##   <fct>     <fct> <chr>        <dbl>  <dbl>
## 1 ctr       t1    score        0.828 0.0200
## 2 ctr       t2    score        0.868 0.0618
## 3 ctr       t3    score        0.887 0.107
## 4 Diet      t1    score        0.919 0.279
## 5 Diet      t2    score        0.923 0.316
## 6 Diet      t3    score        0.886 0.104

The self-esteem score was normally distributed at each time point (p > 0.05), except for ctr treatment at t1, as assessed by Shapiro-Wilk’s test.

Create QQ plot for each cell of design:

ggqqplot(selfesteem2, "score", ggtheme = theme_bw()) +
facet_grid(time ~ treatment, labeller = "label_both")

From the plot above, as all the points fall approximately along the reference line, we can assume normality.

### Computation

res.aov <- anova_test(
data = selfesteem2, dv = score, wid = id,
within = c(treatment, time)
)
get_anova_table(res.aov)
## ANOVA Table (type III tests)
##
##           Effect  DFn  DFd    F        p p<.05   ges
## 1      treatment 1.00 11.0 15.5 2.00e-03     * 0.059
## 2           time 1.31 14.4 27.4 5.03e-05     * 0.049
## 3 treatment:time 2.00 22.0 30.4 4.63e-07     * 0.050

There is a statistically significant two-way interactions between treatment and time, F(2, 22) = 30.4, p < 0.0001.

### Post-hoc tests

A significant two-way interaction indicates that the impact that one factor (e.g., treatment) has on the outcome variable (e.g., self-esteem score) depends on the level of the other factor (e.g., time) (and vice versa). So, you can decompose a significant two-way interaction into:

• Simple main effect: run one-way model of the first variable (factor A) at each level of the second variable (factor B),
• Simple pairwise comparisons: if the simple main effect is significant, run multiple pairwise comparisons to determine which groups are different.

For a non-significant two-way interaction, you need to determine whether you have any statistically significant main effects from the ANOVA output.

#### Procedure for a significant two-way interaction

Effect of treatment. In our example, we’ll analyze the effect of treatment on self-esteem score at every time point.

Note that, the treatment factor variable has only two levels (“ctr” and “Diet”); thus, ANOVA test and paired t-test will give the same p-values.

# Effect of treatment at each time point
one.way <- selfesteem2 %>%
group_by(time) %>%
anova_test(dv = score, wid = id, within = treatment) %>%
get_anova_table() %>%
one.way
## # A tibble: 3 x 9
##   time  Effect      DFn   DFd      F       p p<.05      ges   p.adj
##   <fct> <chr>     <dbl> <dbl>  <dbl>   <dbl> <chr>      <dbl>   <dbl>
## 1 t1    treatment     1    11  0.376 0.552   ""      0.000767 1
## 2 t2    treatment     1    11  9.03  0.012   *       0.052    0.036
## 3 t3    treatment     1    11 30.9   0.00017 *       0.199    0.00051
# Pairwise comparisons between treatment groups
pwc <- selfesteem2 %>%
group_by(time) %>%
pairwise_t_test(
score ~ treatment, paired = TRUE,
)
pwc
## # A tibble: 3 x 11
##   time  .y.   group1 group2    n1    n2 statistic    df       p   p.adj p.adj.signif
## * <fct> <chr> <chr>  <chr>  <int> <int>     <dbl> <dbl>   <dbl>   <dbl> <chr>
## 1 t1    score ctr    Diet      12    12     0.613    11 0.552   0.552   ns
## 2 t2    score ctr    Diet      12    12    -3.00     11 0.012   0.012   *
## 3 t3    score ctr    Diet      12    12    -5.56     11 0.00017 0.00017 ***

Considering the Bonferroni adjusted p-value (p.adj), it can be seen that the simple main effect of treatment was not significant at the time point t1 (p = 1). It becomes significant at t2 (p = 0.036) and t3 (p = 0.00051).

Pairwise comparisons show that the mean self-esteem score was significantly different between ctr and Diet group at t2 (p = 0.12) and t3 (p = 0.00017) but not at t1 (p = 0.55).

Effect of time. Note that, it’s also possible to perform the same analysis for the time variable at each level of treatment. You don’t necessarily need to do this analysis.

The R code:

# Effect of time at each level of treatment
one.way2 <- selfesteem2 %>%
group_by(treatment) %>%
anova_test(dv = score, wid = id, within = time) %>%
get_anova_table() %>%
# Pairwise comparisons between time points
pwc2 <- selfesteem2 %>%
group_by(treatment) %>%
pairwise_t_test(
score ~ time, paired = TRUE,
)
pwc2

After executing the R code above, you can see that the effect of time is significant only for the control trial, F(2, 22) = 39.7, p < 0.0001. Pairwise comparisons show that all comparisons between time points were statistically significant for control trial.

#### Procedure for non-significant two-way interaction

If the interaction is not significant, you need to interpret the main effects for each of the two variables: treatment and time. A significant main effect can be followed up with pairwise comparisons.

In our example (see ANOVA table in res.aov), there was a statistically significant main effects of treatment (F(1, 11) = 15.5, p = 0.002) and time (F(2, 22) = 27.4, p < 0.0001) on the self-esteem score.

Pairwise paired t-test comparisons:

# comparisons for treatment variable
selfesteem2 %>%
pairwise_t_test(
score ~ treatment, paired = TRUE,
)
# comparisons for time variable
selfesteem2 %>%
pairwise_t_test(
score ~ time, paired = TRUE,
)

All pairwise comparisons are significant.

### Report

We could report the result as follow:

A two-way repeated measures ANOVA was performed to evaluate the effect of different diet treatments over time on self-esteem score.

There was a statistically significant interaction between treatment and time on self-esteem score, F(2, 22) = 30.4, p < 0.0001. Therefore, the effect of treatment variable was analyzed at each time point. P-values were adjusted using the Bonferroni multiple testing correction method. The effect of treatment was significant at t2 (p = 0.036) and t3 (p = 0.00051) but not at the time point t1 (p = 1).

Pairwise comparisons, using paired t-test, show that the mean self-esteem score was significantly different between ctr and Diet trial at time points t2 (p = 0.012) and t3 (p = 0.00017) but not at t1 (p = 0.55).

# Visualization: box plots with p-values
pwc <- pwc %>% add_xy_position(x = "time")
bxp +
stat_pvalue_manual(pwc, tip.length = 0, hide.ns = TRUE) +
labs(
subtitle = get_test_label(res.aov, detailed = TRUE),
caption = get_pwc_label(pwc)
)

## Three-way repeated measures ANOVA

### Data preparation

We’ll use the weightloss dataset [datarium package]. In this study, a researcher wanted to assess the effects of Diet and Exercises on weight loss in 10 sedentary individuals.

The participants were enrolled in four trials: (1) no diet and no exercises; (2) diet only; (3) exercises only; and (4) diet and exercises combined.

Each participant performed all four trials. The order of the trials was counterbalanced and sufficient time was allowed between trials to allow any effects of previous trials to have dissipated.

Each trial lasted nine weeks and the weight loss score was measured at the beginning (t1), at the midpoint (t2) and at the end (t3) of each trial.

Three-way repeated measures ANOVA can be performed in order to determine whether there is a significant interaction between diet, exercises and time on the weight loss score.

Load the data and show some random rows by groups:

# Wide format
set.seed(123)
data("weightloss", package = "datarium")
weightloss %>% sample_n_by(diet, exercises, size = 1)
## # A tibble: 4 x 6
##   id    diet  exercises    t1    t2    t3
##   <fct> <fct> <fct>     <dbl> <dbl> <dbl>
## 1 4     no    no         11.1   9.5  11.1
## 2 10    no    yes        10.2  11.8  17.4
## 3 5     yes   no         11.6  13.4  13.9
## 4 11    yes   yes        12.7  12.7  15.1
# Gather the columns t1, t2 and t3 into long format.
# Convert id and time into factor variables
weightloss <- weightloss %>%
gather(key = "time", value = "score", t1, t2, t3) %>%
convert_as_factor(id, time)
# Inspect some random rows of the data by groups
set.seed(123)
weightloss %>% sample_n_by(diet, exercises, time, size = 1)
## # A tibble: 12 x 5
##   id    diet  exercises time  score
##   <fct> <fct> <fct>     <fct> <dbl>
## 1 4     no    no        t1     11.1
## 2 10    no    no        t2     10.7
## 3 5     no    no        t3     12.3
## 4 11    no    yes       t1     10.2
## 5 12    no    yes       t2     13.2
## 6 1     no    yes       t3     15.8
## # … with 6 more rows

In this example, the effect of the “time” is our focal variable, that is our primary concern.

It is thought that the effect of “time” on the weight loss score will depend on two other factors, “diet” and “exercises”, which are called moderator variables.

### Summary statistics

Group the data by diet, exercises and time, and then compute some summary statistics of the score variable: mean and sd (standard deviation)

weightloss %>%
group_by(diet, exercises, time) %>%
get_summary_stats(score, type = "mean_sd")
## # A tibble: 12 x 7
##   diet  exercises time  variable     n  mean    sd
##   <fct> <fct>     <fct> <chr>    <dbl> <dbl> <dbl>
## 1 no    no        t1    score       12  10.9 0.868
## 2 no    no        t2    score       12  11.6 1.30
## 3 no    no        t3    score       12  11.4 0.935
## 4 no    yes       t1    score       12  10.8 1.27
## 5 no    yes       t2    score       12  13.4 1.01
## 6 no    yes       t3    score       12  16.8 1.53
## # … with 6 more rows

### Visualization

Create box plots:

bxp <- ggboxplot(
weightloss, x = "exercises", y = "score",
color = "time", palette = "jco",
facet.by = "diet", short.panel.labs = FALSE
)
bxp

### Check assumptions

#### Outliers

weightloss %>%
group_by(diet, exercises, time) %>%
identify_outliers(score)
## # A tibble: 5 x 7
##   diet  exercises time  id    score is.outlier is.extreme
##   <fct> <fct>     <fct> <fct> <dbl> <lgl>      <lgl>
## 1 no    no        t3    2      13.2 TRUE       FALSE
## 2 yes   no        t1    1      10.2 TRUE       FALSE
## 3 yes   no        t1    3      13.2 TRUE       FALSE
## 4 yes   no        t1    4      10.2 TRUE       FALSE
## 5 yes   no        t2    10     15.3 TRUE       FALSE

There were no extreme outliers.

#### Normality assumption

Compute Shapiro-Wilk test for each combinations of factor levels:

weightloss %>%
group_by(diet, exercises, time) %>%
shapiro_test(score)
## # A tibble: 12 x 6
##   diet  exercises time  variable statistic     p
##   <fct> <fct>     <fct> <chr>        <dbl> <dbl>
## 1 no    no        t1    score        0.917 0.264
## 2 no    no        t2    score        0.957 0.743
## 3 no    no        t3    score        0.965 0.851
## 4 no    yes       t1    score        0.922 0.306
## 5 no    yes       t2    score        0.912 0.229
## 6 no    yes       t3    score        0.953 0.674
## # … with 6 more rows

The weight loss score was normally distributed, as assessed by Shapiro-Wilk’s test of normality (p > .05).

Create QQ plot for each cell of design:

ggqqplot(weightloss, "score", ggtheme = theme_bw()) +
facet_grid(diet + exercises ~ time, labeller = "label_both")

From the plot above, as all the points fall approximately along the reference line, we can assume normality.

### Computation

res.aov <- anova_test(
data = weightloss, dv = score, wid = id,
within = c(diet, exercises, time)
)
get_anova_table(res.aov)
## ANOVA Table (type III tests)
##
##                Effect  DFn  DFd       F        p p<.05   ges
## 1                diet 1.00 11.0   6.021 3.20e-02     * 0.028
## 2           exercises 1.00 11.0  58.928 9.65e-06     * 0.284
## 3                time 2.00 22.0 110.942 3.22e-12     * 0.541
## 4      diet:exercises 1.00 11.0  75.356 2.98e-06     * 0.157
## 5           diet:time 1.38 15.2   0.603 5.01e-01       0.013
## 6      exercises:time 2.00 22.0  20.826 8.41e-06     * 0.274
## 7 diet:exercises:time 2.00 22.0  14.246 1.07e-04     * 0.147

From the output above, it can be seen that there is a statistically significant three-way interactions between diet, exercises and time, F(2, 22) = 14.24, p = 0.00011.

Note that, if the three-way interaction is not statistically significant, you need to consult the two-way interactions in the output.

In our example, there was a statistically significant two-way diet:exercises interaction (p < 0.0001), and two-way exercises:time (p < 0.0001). The two-way diet:time interaction was not statistically significant (p = 0.5).

### Post-hoc tests

If there is a significant three-way interaction effect, you can decompose it into:

• Simple two-way interaction: run two-way interaction at each level of third variable,
• Simple simple main effect: run one-way model at each level of second variable, and
• simple simple pairwise comparisons: run pairwise or other post-hoc comparisons if necessary.

#### Compute simple two-way interaction

You are free to decide which two variables will form the simple two-way interactions and which variable will act as the third (moderator) variable. In the following R code, we have considered the simple two-way interaction of exercises*time at each level of diet.

Group the data by diet and analyze the simple two-way interaction between exercises and time:

# Two-way ANOVA at each diet level
two.way <- weightloss %>%
group_by(diet) %>%
anova_test(dv = score, wid = id, within = c(exercises, time))
two.way
## # A tibble: 2 x 2
##   diet  anova
##   <fct> <list>
## 1 no    <anov_tst>
## 2 yes   <anov_tst>
# Extract anova table
get_anova_table(two.way)
## # A tibble: 6 x 8
##   diet  Effect           DFn   DFd     F        p p<.05   ges
##   <fct> <chr>          <dbl> <dbl> <dbl>    <dbl> <chr>   <dbl>
## 1 no    exercises          1    11 72.8  3.53e- 6 *       0.526
## 2 no    time               2    22 71.7  2.32e-10 *       0.587
## 3 no    exercises:time     2    22 28.9  6.92e- 7 *       0.504
## 4 yes   exercises          1    11 13.4  4.00e- 3 *       0.038
## 5 yes   time               2    22 20.5  9.30e- 6 *       0.49
## 6 yes   exercises:time     2    22  2.57 9.90e- 2 ""      0.06

There was a statistically significant simple two-way interaction between exercises and time for “diet no” trial, F(2, 22) = 28.9, p < 0.0001, but not for “diet yes” trial, F(2, 22) = 2.6, p = 0.099.

Note that, it’s recommended to adjust the p-value for multiple testing. One common approach is to apply a Bonferroni adjustment to downgrade the level at which you declare statistical significance.

This can be done by dividing the current level you declare statistical significance at (i.e., p < 0.05) by the number of simple two-way interaction you are computing (i.e., 2).

Thus, you only declare a two-way interaction as statistically significant when p < 0.025 (i.e., p < 0.05/2). Applying this to our current example, we would still make the same conclusions.

#### Compute simple simple main effect

A statistically significant simple two-way interaction can be followed up with simple simple main effects.

In our example, you could therefore investigate the effect of time on the weight loss score at every level of exercises or investigate the effect of exercises at every level of time.

You will only need to consider the result of the simple simple main effect analyses for the “diet no” trial as this was the only simple two-way interaction that was statistically significant (see previous section).

Group the data by diet and exercises, and analyze the simple main effect of time. The Bonferroni adjustment will be considered leading to statistical significance being accepted at the p < 0.025 level (that is 0.05 divided by the number of tests (here 2) considered for “diet:no” trial.

# Effect of time at each diet X exercises cells
time.effect <- weightloss %>%
group_by(diet, exercises) %>%
anova_test(dv = score, wid = id, within = time)
time.effect
## # A tibble: 4 x 3
##   diet  exercises anova
##   <fct> <fct>     <list>
## 1 no    no        <anov_tst>
## 2 no    yes       <anov_tst>
## 3 yes   no        <anov_tst>
## 4 yes   yes       <anov_tst>
# Extract anova table
get_anova_table(time.effect) %>%
filter(diet == "no")
## # A tibble: 2 x 9
##   diet  exercises Effect   DFn   DFd     F        p p<.05   ges
##   <fct> <fct>     <chr>  <dbl> <dbl> <dbl>    <dbl> <chr>   <dbl>
## 1 no    no        time       2    22  1.32 2.86e- 1 ""      0.075
## 2 no    yes       time       2    22 78.8  9.30e-11 *       0.801

There was a statistically significant simple simple main effect of time on weight loss score for “diet:no,exercises:yes” group (p < 0.0001), but not for when neither diet nor exercises was performed (p = 0.286).

#### Compute simple simple comparisons

A statistically significant simple simple main effect can be followed up by multiple pairwise comparisons to determine which group means are different.

Group the data by diet and exercises, and perform pairwise comparisons between time points with Bonferroni adjustment:

# Pairwise comparisons
pwc <- weightloss %>%
group_by(diet, exercises) %>%
pairwise_t_test(score ~ time, paired = TRUE, p.adjust.method = "bonferroni") %>%
select(-df, -statistic) # Remove details
# Show comparison results for "diet:no,exercises:yes" groups
pwc %>% filter(diet == "no", exercises == "yes") %>%
select(-p)     # remove p columns
## # A tibble: 3 x 9
##   <fct> <fct>     <chr> <chr>  <chr>  <int> <int>        <dbl> <chr>
## 1 no    yes       score t1     t2        12    12 0.000741     ***
## 2 no    yes       score t1     t3        12    12 0.0000000121 ****
## 3 no    yes       score t2     t3        12    12 0.000257     ***

In the pairwise comparisons table above, we are interested only in the simple simple comparisons for “diet:no,exercises:yes” groups. In our example, there are three possible combinations of group differences. We could report the pairwise comparison results as follow.

All simple simple pairwise comparisons were run between the different time points for “diet:no,exercises:yes” trial. The Bonferroni adjustment was applied. The mean weight loss score was significantly different in all time point comparisons when exercises are performed (p < 0.05).

### Report

A three-way repeated measures ANOVA was performed to evaluate the effects of diet, exercises and time on weight loss. There was a statistically significant three-way interaction between diet, exercises and time, F(2, 22) = 14.2, p = 0.00011.

For the simple two-way interactions and simple simple main effects analyses, a Bonferroni adjustment was applied leading to statistical significance being accepted at the p < 0.025 level.

There was a statistically significant simple two-way interaction between exercises and time for “diet no” trial, F(2, 22) = 28.9, p < 0.0001, but not for “diet yes”" trial, F(2, 22) = 2.6, p = 0.099.

There was a statistically significant simple simple main effect of time on weight loss score for “diet:no,exercises:yes” trial (p < 0.0001), but not for when neither diet nor exercises was performed (p = 0.286).

All simple simple pairwise comparisons were run between the different time points for “diet:no,exercises:yes” trial with a Bonferroni adjustment applied. The mean weight loss score was significantly different in all time point comparisons when exercises are performed (p < 0.05).

# Visualization: box plots with p-values
pwc <- pwc %>% add_xy_position(x = "exercises")
pwc.filtered <- pwc %>%
filter(diet == "no", exercises == "yes")
bxp +
stat_pvalue_manual(pwc.filtered, tip.length = 0, hide.ns = TRUE) +
labs(
subtitle = get_test_label(res.aov, detailed = TRUE),
caption = get_pwc_label(pwc)
)

## Summary

This chapter describes how to compute, interpret and report repeated measures ANOVA in R. We also explain the assumptions made by repeated measures ANOVA tests and provide practical examples of R codes to check whether the test assumptions are met.

Version: Français

• Hanna

What an excellent resource! Thank you for providing this website to the public. I will recommend it to my students, too.

• Kassambara

Thank you for your positive feedback, highly appreciated!

• Tom

Thank you so much, helps a lot! One thing that might be relevant for some people is when to use between subjects instead of within subjects. But, the main thing is that this is the best source for doing these types of analyses I could find on the web, so great job!

• Kassambara

Thank you for your positive feedback! I’ll take your suggestion into account in the next update

• TomO

Hi! Thanks for such a wonderfully written intro to the topic with easy to read R code – a breath of fresh air!

“It’s also possible to keep the outliers in the data and perform robust ANOVA test using the WRS2 package.”

I’ve been reading the Mair et al (2018) paper that introduces WRS2 titled “Robust Statistical Methods Using WRS2” and specifically section 6.1 on Repeated Measurement and Mixed ANOVA Designs that introduces the function I believe you are referring to (rmanova). However, the authors state that the function provides a “robust heteroscedastic repeated measurement ANOVA based on the trimmed means” which defaults to 20%. Would this trimmed mean not remove the outliers, rather than keep them?

I’d be grateful for any help you could offer, maybe I’ve got the wrong function or the wrong end of the stick?

Cheers!

• Richard

Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, …) :
0 (non-NA) cases

• Kassambara

Hi, could you please provide a reproducible example?

• Anton

I am also having this problem, though the dataset you use in the example works fine :(.

Data:

id treatment t1 t2 t3
1 Diet 3 7 7
2 Diet 5 7 7
3 Diet 2 7 8
4 Diet 4 3 4
5 Diet 5 4 4
6 Diet 3 4 4
7 Diet 5 9 9
8 Diet 3 4 4
9 Diet 3 8 8
10 Diet 5 4 4
147 ctr 5 3 3
148 ctr 4 6 6
149 ctr 7 8 8
150 ctr 3 2 2
151 ctr 7 4 4
152 ctr 3 6 6
153 ctr 5 5 5
154 ctr 4 7 7
155 ctr 5 1 1
157 ctr 6 5 5

zna % gather(key = “time”, value = “score”, t1, t2, t3) %>%
convert_as_factor(id, time)

res.aov <- anova_test(
data = zna, dv = score, wid = id,
within = c(treatment, time)
)
get_anova_table(res.aov)

• Jake

I got this error message when I tried to do two-way anova. My case was different from the example because I had different participants in all (three) groups; I think this is why I got an error.

Eventually, I analysed my data using lme function and setting a unique random effect to each participant.

• Aziz

I got the same error again and again. I was following all the steps. Later I discovered. that my data frame consisted of 2 additional variables that were not used. For some strange reasons, if I deleted the unwanted variables, the error no longer came. It is difficult to explain this situation.

• May

Hello! Thank’s a lot for this article, it helps a lot. I’m using three-way anova for my dataset and when I run function ‘anova_test’ (opy from this text here) i get an error: ‘Each row of output must be identified by a unique combination of keys.
Keys are shared for 715 rows:’. Maybe you could help and tell me what’s going on?

My function looks like:
aov <- anova_test(data=rtdata, dv=rt, wid=part, within=c(mod_sm, con, resp))
where:
rtdata – data in long format
part – factor, participant code from 1 to 12; each participant got 60 trials, so his code repeats 60 times in dataset
mod_sm – text, levels 'pos' and 'neg'
con – text, 'n' and 'y'
resp – factor, 1 and 2

Thank you!

• Marie

Hey May, did you manage to figure out how to fix this issue? I have exactly the same issue and my data structure is similar.

• Anton

I have been trying to do a two-way repeated samples and spent many hours until I worked out what was wrong – my control and treatment groups aren’t the same, so don’t have the same ids repeated. So it is not “ctr” followed by “Diet” for the same subjects but separate subject

Can you comment on how this would work? What kind of test would it be? Thanks!

• Odi Makov

Very informative!
I’ll recommend this guide everywhere

• Kassambara

Thank you for your positive feedback, highly appreciated!

• Hannah PW

Hi, I’ve been trying for days to run the 3-way repeated measures ANOVA for my dataset in a long format with 2x treatment groups (temp & region) and then a time group (T0, T1, T2, T3). Whenever I run the ANOVA, I get the same error message that has been mentioned above and in many other R forums (StackExchange, etc.), although I’ve not as yet seen a solution as my data does not contain any null data and all my groups are in factor format (I’ve also tried them in numeric format without luck):

Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, …) :
0 (non-NA) cases

Can you help or offer any suggestions that I might not have tried?

• Hans Bourgeois

Hi Kassambara,

It seems there are a conflict with the packtage “tidyverse” and “rstatix”.
This conflict may explain the errors :
– “Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, …) :
0 (non-NA) cases”
– “Error: Unknown column ID”

To remedy this i use the package “conflict_prefer”.

Moreover, the first error also appears when i have unused columns (Factor or Numeric) in my dataframe. When i remove unused columns, the anova_test work fine.

• Jason Weber

Hi great article!
Tried to produce graph report of one-way anova and error at bottom of code came up.
PS_tried to install pubr but got nothing.
Thx JW

# Repeated One-Way ANOVA
res.aov <- anova_test(data = V3, dv=Data, wid=Name, within = Rep)
get_anova_table(res.aov)

# pairwise comparisons
pwc %
pairwise_t_test(
Data ~ Rep, paired = TRUE,
)
pwc

# Visualization: box plots with p-values
bxp +
stat_pvalue_manual(pwc) +
labs(
subtitle = get_test_label(res.aov, detailed = TRUE),
caption = get_pwc_label(pwc)

Error in bxp + stat_pvalue_manual(pwc) + labs(subtitle = get_test_label(res.aov, :
non-numeric argument to binary operator

• Kassambara

Hi,

It should be possible to install pubr using devtools::install_github(“kassambara/pubr”).

In order efficiently I need a reproducible example with a demo data.

Thanks

• Elizabeth

Really great examples that are well explained.

My issue is when trying to visualize the final box plot with the p-values for any of the examples. I receive the following error each time:

Error in f(…) :
Can only handle data with groups that are plotted on the x-axis
1: Ignoring unknown parameters: hide.ns
2: Ignoring unknown aesthetics: xmin, xmax, annotations, y_position

• Kassambara

• Sophie

Thanks for the excellent tutorial! I got the same error with data visualization. Have you figured out the potential issue with this error?

• Qian Yang

Hi Kassambara,
Thanks a lot for your input, it’s really helpful. I have one question that the output of two-way repeated ANOVA analysis in R is not consistent with the results obtained in SPSS, what kind of reason do you think can explain this difference? thanks in advance.

• Kassambara

Hi!

After checking, I found that the output of SPSS and rstatix are the same when considering the rstatix raw output (res.aov).

When displaying the ANOVA table using the function get_anova_table(res.aov), the sphericity correction is automatically applied when the sphericity can’t be assumed.

The sphericity assumption can be checked using the Mauchly’s test of sphericity, which is automatically reported when using the R function rstatix::anova_test()

library(rstatix)

# Data preparation
data("selfesteem2", package = "datarium")
selfesteem2 <- selfesteem2 %>%
gather(key = "time", value = "score", t1, t2, t3) %>%
convert_as_factor(id, time)
head(selfesteem2, 3)
## # A tibble: 3 x 4
##   id    treatment time  score
##   <fct> <fct>     <fct> <dbl>
## 1 1     ctr       t1       83
## 2 2     ctr       t1       97
## 3 3     ctr       t1       93
# Anova result
res.aov <- anova_test(
data = selfesteem2, dv = score, wid = id,
within = c(treatment, time)
)

# Display raw results
res.aov
## ANOVA Table (type III tests)
##
## $ANOVA ## Effect DFn DFd F p p<.05 ges ## 1 treatment 1 11 15.541 2.00e-03 * 0.059 ## 2 time 2 22 27.369 1.08e-06 * 0.049 ## 3 treatment:time 2 22 30.424 4.63e-07 * 0.050 ## ##$Mauchly's Test for Sphericity
##           Effect     W     p p<.05
## 1           time 0.469 0.023     *
## 2 treatment:time 0.616 0.089
##
## $Sphericity Corrections ## Effect GGe DF[GG] p[GG] p[GG]<.05 HFe DF[HF] p[HF] p[HF]<.05 ## 1 time 0.653 1.31, 14.37 5.03e-05 * 0.705 1.41, 15.52 2.81e-05 * ## 2 treatment:time 0.723 1.45, 15.9 1.25e-05 * 0.803 1.61, 17.66 4.82e-06 * # Display ANOVA table # Sphericity correction is applied by default if # the the Mauchly's test is significant # Read this: https://www.datanovia.com/en/lessons/mauchlys-test-of-sphericity-in-r/ get_anova_table(res.aov) ## ANOVA Table (type III tests) ## ## Effect DFn DFd F p p<.05 ges ## 1 treatment 1.00 11.00 15.541 2.00e-03 * 0.059 ## 2 time 1.31 14.37 27.369 5.03e-05 * 0.049 ## 3 treatment:time 2.00 22.00 30.424 4.63e-07 * 0.050 # No correction get_anova_table(res.aov, correction = "none") ## ANOVA Table (type III tests) ## ## Effect DFn DFd F p p<.05 ges ## 1 treatment 1 11 15.541 2.00e-03 * 0.059 ## 2 time 2 22 27.369 1.08e-06 * 0.049 ## 3 treatment:time 2 22 30.424 4.63e-07 * 0.050 • Leon Hey, great tutorial! This is by far the easiest and most understandable approach I have found. I am working on my bachelor thesis and using R instead of SPSS. I need to calculate several 2 way repeated Anovas and rstatix seems more than handy atm. However, I cannot find any explanation of the output: What is the value of DFd? I calculated a oneway repeated Anova with 1 factor of 3 levels. DFn seems to be the degree of freedom between the levels (2) and DFd’s value is 30. I have 16 participants and I dont know how the degrees of freedom could go up this far… Though, my statistics courses were 2 years ago… Thanks in advance! • Kassambara Hi, thank you for the positive feedback! DFn is the Degrees of Freedom in the numerator (i.e. DF for the effect of the experiment); DFd Degrees of Freedom in the denominator (i.e., DF for the residual, DFd Degrees of Freedom in the denominator (i.e., DF error)). This is explained in the documentation of anova_summary(), but obviously, we should clarify it in the current blog post. If only one factor is repeated measures, the number of degrees of freedom (DFd) can be calculated as (n-1)(a-1), where n is the number of participants and a is the number of levels of the repeated measures factor. If two factors are repeated measures, the number of degrees of freedom (DFd) for the two-way interaction term can be calculated as (n-1)(a-1)(b-1) where n is the number of participants, a is the number of levels one factor, and b is the number of levels of the other factor. Note that, a sphericity correction is applied, when the sphericity assumtion is not met, resulting in a fractional degres of freedom. • Leon Thank you so much for the fast response! (16-1) (3-2) = 30 should be fine then. Due to the pandemic I was not able to test all participants and I am allowed to set the alpha level to 0.10 instead of 0.05. Is it possible to adjust it within your package, so that get_anova_table() marks p<.10 with * ? Even if not, your package worked great so far. Thanks again! • Kassambara It’s possible to do something like this: res <- get_anova_table(res.aov) res[, "p<.05"] <- ifelse(res$p <= 0.1, "*", "")
res
• Gene Lu

Hi,
Thank you very much for your website! It helps a lot with my analysis works.
I have one question, after running the rstatix::anova_test() code, it only shows the ANOVA table. I would like to check the sphericity assumption for my data but could not find it in the anova_test result. (I have tried to use your example data (selfesteem2) to do the analysis and can check the sphericity assumption in the ANOVA table. So I think the problem is from my data.) Do you have any idea what would be the possible reasons? (I run the two-wat rmANOVA) Thank you very much again!!

• Kassambara

Note that the the Mauchly’s Test of Sphericity is only reported for variables or effects with >2 levels because sphericity necessarily holds for effects with only 2 levels Read more: Mauchly’s Test of Sphericity in R. Does your within-measures variable contains only two levels?

• Gene Lu

Hi, thank you very much for your response! I found the solution from your answer in the previous question asked by Jack. My data is a mixed two-way ANOVA design (have 3 treatments). After changing the analysis method I get the Mauchly’s Test of Sphericity in my result. Thank you very much for your help!

• X32

There are non-parametric method for 2/3 way repeated ANOVA. It’s WTS (Wald-Type Statistics), permuted WTS, ATS (ANOVA-Type Statistics) and ART (Aligned-Rank Transform) ANOVA. Both methods are available in R, not to mention Conover’s ANOVA (on “classic” ranks). Please check: nparLD, GFD, rankFD, ARTool and MANOVA packages.

• Kassambara

Thank you for your input, very useful and instructive. Obviously, we’ll write on these packages!

• Sara

Hi! Thanks for this beautiful and useful explanation.
I have a question regarding treating outliers:

You say that if outliers are NOT extreme you keep in the analysis right?

If I want to remove outliers from my dataset and compute anova, which is the command to remove them?

Sara

• Kassambara

Thank you Sara, for your positive feedback, highly appreciated!

Once you identified outliers, you can filter out the corresponding samples (before computing ANOVA), by using the R function dplyr::filter().

• Viper

Excellent guide and your hard work is very much appreciated.
I am getting the same error code as several others:
Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, …) :
0 (non-NA) cases

I have tried using the Mixed ANOVA, but that will not run either. I also tried using the “conflict_prefer”as a poster suggested, but that does not help either. There are no NA’s in any of my data columns. I have three categorical IVs (speaker, complexity, format) and not sure where I am going wrong.

I will try and post a reproducible R script here later.

Thanks for the help and thanks for the excellent guides.

• Pi

Great article! Only one aspect, a non-significant p-value for the Shapiro test doesn’t guarantee that the data is normally distributed. It may happen that the data you have is not enough to reject the null. The q-q plot is always a good double check as you did!

• Kassambara

• Daniel

Hi,
Thank for build so amazing tool!
My question is: My data is not perfectly normal and I have some outliers, but results are consistent. In these cases, I am more calm if I can see model residuals. Any way to obtain anova_test residuals?

• Kassambara

The returned object by the anova_test() function has an attribute called args, which is a list holding the arguments used to fit the ANOVA model, including: data, dv, within, between, type, model, etc.

You can access to the residuals as follow: residuals(attr(res.aov, "args")$model) • Daniel Are you the best? maybe… jejeje Thanks! • Danny Hi, I’m trying to run Two-way repeated measures ANOVA but I kept getting this message: Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, …) : 0 (non-NA) cases And I don’t see any missing data for my data. This is what I did: res.aov <- anova_test( data = Tab14_4, dv = Behavior, wid = Subject, within = c(Group, Interval) ) Pls help! • Peggy Hi, I’m also trying to run Two-way repeated measures ANOVA like Danny but I kept getting this message: Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, …) : 0 (non-NA) cases And I don’t see any missing data for my data. This is what I did: anova_test(data = mydata2,dv = Score,wid = Subjects,within = c(Drinks,Gender)) May I sent my demo to you? I really need your help! • Kassambara Hi, In order to fix this issue I need a reproducible example with a data. Would you like to send your data with a reproducible R script? • Lauren Hi there, I am also experiencing the same issue as Peggy and Danny. I have spent hours trying to get this to work (after it had successfully worked many times the other day) and have had no luck. It works perfectly fine with the data presented on this site, but is no longer working with my data. I also have no missing cases present within my dataset. res.aov <- anova_test(data = newanova, dv = score, wid = id, within = c(treatment, time)) Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, …) : 0 (non-NA) cases Please help! • Kassambara Hi, In order to fix this issue I need a reproducible example with a data. Would you like to send your data with a reproducible R script? • Laurianne Have you find a solution? I am in the same situation. • Dongze chen If the data is an imbalance between groups, Does this still work? Cause in my data, the first group has 100 samples, the second group only has 10 samples I am looking forward to your reply! Thank you very much! • Laurianne Schoenberger Hi ! Thank you so much for this explanation. I found this so helpful. However, I have some issues when doing this with my data. I removed the NA from my dataset. I tried to do as your “example 1: Reproducible R script using R built-in data” but I am not sure this works… I had error mesage while running “pubr::render_r()”. My independent variable is secondhand smoke (SHS yes or no) and cognitive test score at 2 points of time. When I compute the anova test. I got this message: “Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, …) : 0 (non-NA) cases” # Computation res.aov2 <- anova_test(data = ds, dv = score, wid = ID, within = c(SHSCAT, time)) get_anova_table(res.aov2) Then when I want to see the effect of independent variable at each time point, I got an error message: "Error: Problem with mutate() input data. x 0 (non-NA) cases i Input data is map(.data$data, .f, …)"

# Effect of SHS smoke at each time point
one.way % group_by(time) %>%
anova_test(dv = score, wid = ID, within = SHSCAT) %>%
get_anova_table() %>%

I really need your help. I tried to find other codes, to watch videos, but nothing fixed my problem.
Thank you for time.

• Kassambara

Please, provide your data and the R script you used so that we can help, thanks

• Chiara

Hello,
thank you for the great article! I would like to use this procedure in my master’s thesis. Do you have some references to literature for me ?
Thank you so much!

• Chiara

I need references for the pairwise paired t-tests and the one way ANOVA with repeated measures. Thank you so much!

• Berman

Hi there, is there any way to hide ” eta square” value when adding pvalues on ggplot?

• Berman

• Stine

Hey Kassambara,

Thanks for your amazing instructions on repeated measures ANOVA. However, like many others it seem, I also receive the error message:

“Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, …) :
0 (non-NA) cases”

when trying to run the 2-way repeated ANOVA. I have removed any unused variables from the data set and tried running the: “suppressPackageStartupMessages(library(rstatix))” as well. And if I try to install your rstatix package, I get this message ‘kassambara/rstatix’ is not available (for R version 4.0.1) . I have just updated my Rstudio and believe I have the newest version?

My data looks like this:
ID Site Day Sediment

1 1 Site 1 Day 1 5.5
2 2 Site 1 Day 1 -2.75
3 3 Site 1 Day 1 -9.65
4 1 Site 1 Day 3 9.3
5 2 Site 1 Day 3 -0.2
6 3 Site 1 Day 3 -2.45

I have three Sites (site 1, Site 2, Site 3) which are measured 5 times (Day 1, Day 3, Day 7, Day 14, Day 33) and my dependent variable is sediment levels (Sediment).

I would be very grateful if you can help me solve this problem.

Kind regards

• Kassambara

Thank you for the feedback! This has been fixed now in the latest dev version (see https://github.com/kassambara/rstatix/issues/55), which can be installed using devtools::install_github(“kassambara/rstatix”)

• Anna

Hi Kassambara, thank you for your great explanations and examples! I would have one question on you. How to compare more than two treatment groups in 2-way repeated-measures ANOVA?

I am trying to do a two-way repeated-measures ANOVA, only I have 3 treatment groups (independent variable 1) over time (independent variable 2). I got significant two-way interaction and a simple main effect. You used paired t-test in your example to compare the treatment groups so I tried anova_test instead. That was significant for some of the time points. But now, I cannot figure out how to do a posthoc Tuckey to know how they differ? I can’t get TukeyHSD working with the output from anova_test (I always used aov function, but there you cannot use group_by()).

I ran my anova as follows:
comparison %
group_by(time) %>%
anova_test(dependent variable ~ treatment) %>%
get_anova_table() %>%

Kind regards

• Charles

Hi Kassambara,

Thanks for the amazing tutorial!

I have a question that whether partial counterbalance design can be analyzed with anova_test().

I have two IV: tool (N, R, S) and topic (T, F, K), each one has three levels. If it is a complete counterbalancing design, each participants should has 9 treatment. But that is impratical. So in my study, each participant experienced 3 treatment (e.g. [NT, RF, SK], [RT, NK, SF]), each tool and topic only appear once for each participant.

When I used anova_test(), I got the same error as above: Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, …) : 0 (non-NA) cases. I do not have NA in my dataset.

Should I use mixed model instead? But every IV is within subject variable, I think I should use repeated measurement ANOVA.

If it does not make sense I will provide reproducible demo data.

• Kassambara

Hi, this issue seems to be related to the issue described at https://github.com/kassambara/rstatix/issues/55, which has been fixed in the rstatix dev version available on Github. Please make sure you have installed the latest dev version of the rstatix package.

• Charles

Hi Kassambara,

I installed the package using devtools::install_github(“kassambara/rstatix”)

But still get this error.

Here is a script:

demo <- data.frame( id = c(1L,1L,1L,2L,2L, 2L,3L,3L,3L,4L,4L,4L,5L,5L,5L,6L,6L,6L, 7L,7L,7L,8L,8L,8L,9L,9L,9L,10L,10L,10L, 11L,11L,11L,12L,12L,12L,13L,13L,13L,14L,14L, 14L,15L,15L,15L,16L,16L,16L,17L,17L,17L, 18L,18L,18L,19L,19L,19L,20L,20L,20L,21L,21L, 21L,22L,22L,22L,23L,23L,23L,24L,24L,24L, 25L,25L,25L,26L,26L,26L,27L,27L,27L,28L, 28L,28L,29L,29L,29L,30L,30L,30L,31L,31L,31L, 32L,32L,32L,33L,33L,33L,34L,34L,34L,35L, 35L,35L,36L,36L,36L), result = c(4,4,4,4,4.5,3.5, 4,2,4,4,2,4,4,4,4,2,3.5,3,6,1,4,4, 4.5,3.5,3.5,2,5.5,5.5,6,6.5,3,2,5,1,6, 6.5,3.5,4.5,4.5,3,5,5.5,2,2,3.5,2,5,5,4, 6,7,5,5,5,3,5,5,4,4.5,5,4.5,6.5,6.5,5, 4.5,3,5,5,6,6,6,7,4,5.5,5,6,6,5.5,2, 4.5,5,1,7,7,1,4,4.5,3,5,6,6,2.5,7,3,3, 5,2,5,7,1.5,2,2,1.5,4.5,5,1,3,1), topic = c(1L,2L,3L,2L,3L, 1L,3L,1L,2L,3L,2L,1L,2L,1L,3L,1L,3L,2L, 2L,3L,1L,3L,1L,2L,1L,2L,3L,1L,3L,2L,3L, 2L,1L,2L,1L,3L,3L,1L,2L,1L,2L,3L,2L,3L, 1L,2L,1L,3L,1L,3L,2L,3L,2L,1L,1L,2L,3L, 2L,3L,1L,3L,1L,2L,3L,2L,1L,2L,1L,3L,1L, 3L,2L,2L,3L,1L,3L,1L,2L,1L,2L,3L,1L,3L, 2L,3L,2L,1L,2L,1L,3L,3L,1L,2L,1L,2L,3L, 2L,3L,1L,2L,1L,3L,1L,3L,2L,3L,2L,1L), tool = c(1L,2L,3L,2L,3L, 1L,3L,1L,2L,3L,2L,1L,2L,1L,3L,1L,3L,2L, 1L,2L,3L,2L,3L,1L,3L,1L,2L,3L,2L,1L,2L, 1L,3L,1L,3L,2L,1L,2L,3L,2L,3L,1L,3L,1L, 2L,3L,2L,1L,2L,1L,3L,1L,3L,2L,1L,2L,3L, 2L,3L,1L,3L,1L,2L,3L,2L,1L,2L,1L,3L,1L, 3L,2L,1L,2L,3L,2L,3L,1L,3L,1L,2L,3L,2L, 1L,2L,1L,3L,1L,3L,2L,1L,2L,3L,2L,3L,1L, 3L,1L,2L,3L,2L,1L,2L,1L,3L,1L,3L,2L), trial = c(1L,2L,3L,1L,2L, 3L,1L,2L,3L,1L,2L,3L,1L,2L,3L,1L,2L,3L, 1L,2L,3L,1L,2L,3L,1L,2L,3L,1L,2L,3L,1L, 2L,3L,1L,2L,3L,1L,2L,3L,1L,2L,3L,1L,2L, 3L,1L,2L,3L,1L,2L,3L,1L,2L,3L,1L,2L,3L, 1L,2L,3L,1L,2L,3L,1L,2L,3L,1L,2L,3L,1L, 2L,3L,1L,2L,3L,1L,2L,3L,1L,2L,3L,1L,2L, 3L,1L,2L,3L,1L,2L,3L,1L,2L,3L,1L,2L,3L, 1L,2L,3L,1L,2L,3L,1L,2L,3L,1L,2L,3L) )

 

anova_test( data = demo, dv = result, wid = id, within = c(topic, tool, trial) )

Great tutorial and I’m getting some great plots. However, I have a small issue that I cannot resolve.

When doing a pairwise comparison and trying to add the significance to the boxplot using pairwise_t_test() and add_xy_position() the groups are ordered along the x axis alphabetically.
I have 3 groups: Normal, DblEl, ON. I want to display in that order but by default they are ordered DblEL, Normal, ON. I can alter the order in ggboxplot using order = c(‘Normal’, ‘DblEl’, ‘ON’), but cannot get the paired comparisons or their associated x numbering to match. Any ideas other than relabelling groups alphabetically? Thanks!

Hi,

I would suggest to prepare your data by defining the desired order of the factor variable levels before using any plotting or stat functions. For example, I would go as follow:

1/ Data preparation: defining the order of the variable

my_data$grouping_var = factor(my_data$grouping_var, levels = c(“Normal”, “DblEl”, “ON”))

2/ now you can use ggboxplot() or any stat functions. The defined factor levels at step 1 will be taken into account by all functions.

Thanks Alboukadel. Turns out r was automatically ascribing levels when I used read.csv(), seemingly alphabetically. I reassigned as you suggested and now works perfectly. Woop! Thanks.

• Clara

Hi,

For one-way repeated measures ANOVA, is it possible to include the significance levels in the boxplot, but NOT when it is nonsignificant? In other words, is it possible to not to include it when it is “ns”.

Thanks,

Clara

• Teresa

Hi Kassambara,
I’m facing the same issue as many people here. I get the following error Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, …) :
0 (non-NA) cases
I don’t have any NAs in my dataset and I have downloaded the newest version you have suggested to other comments on here. Is there possibly a limit on data that can be managed? I have 5000+ data points. I want to compare two groups of stations within the same water body measured in 3 seasons and with a certain treatment. Hope you can help. Thanks.

Did you installed the dev version from Github? if yes and you still have the issue, please send to meby e-mail the data and the script you used so that I can help.

• SP

Hi,

I keep getting this error when I try to run two-way ANOVA: Error: Problem with mutate() input data.
x Can’t subset columns that don’t exist.

Can anyone help me with this??

• Dai

I was getting the same error message, but I found it was due to group_by(). ungroup() stopped the error. I hope this will help.

• Mehmet gazan

Hi,
Thx a lot for the highly valuable info you’ve shared. I got a question mark that might also interest you for my own research purposes:
I’ve a dataset which has 2 tasks and each task comprises of 3 items. 40 participants complete these 6 items in total respectively. Their performance is assessed by 3 raters. My question is whether this task type have any impact on the participants performance. In this case, how should my SPSS data entry look like?

• Iza

Hi, thank you so much for this tutorial!

I have a question about Mauchly’s test. I am doing a two-way repeated measures ANOVA. Both the independent variables have more than two levels (5 and 6), but I cannot find any output for the sphericity assumptions. I’m running it like this:

res.aov <- anova_test(
data = G1_noid, dv = rating, wid = ppno,
within = c(instr_pair, pitch)
)
get_anova_table(res.aov)

Is there something I'm doing wrong here?

I also had a question for the normality assumption. As you can see there are a lot of levels, so the Shaprio Wilk is applied 30 times. Six of those came out significant. Does that mean I will have to run robust testing, even if relatively few groups were non-normal?

Thank you!

• Debanjan Borthakur

The pairwise t test results don’t match in R and SPSS. In SPSS, the pairwise comparison across time seems to be effected by the levels in the factor group’s levels which is not the case with R pairwise_t test function you mentioned. Any suggestion?

• Aline

First of all, thank you very much Kassambara for this work! I am trying to run a repeated measure two-way ANOVA with my own data and I still seem to run into the error previously described by others;

Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, …) : )
0 (non-NA) cases

I have installed the latest dev version of the rstatix, however, that does not seem to solve the problem. Do you have any idea how to solve this issue? I would very much appreciate your help!

the code looks as followed:

my_data % sample_n_by(treatment, size = 1)

# Gather the columns t1, t2 and t3 into long format.
# Convert id and time into factor variables
my_data %
gather(key = “minute”, value = “score”, t1, t2) %>%
convert_as_factor(id, minute, treatment)

#Inspect some random rows of the data by groups
set.seed(123)
my_data %>% sample_n_by(treatment, minute, size = 1)

#summery of data
my_data %>%
group_by(treatment, minute) %>%
get_summary_stats(score, type = “mean_sd”)

#check for outlayers
my_data %>%
group_by(treatment, minute) %>%
identify_outliers(score)

#check for normality
my_data %>%
group_by(treatment, minute) %>%
shapiro_test(score)

#QQplot also check for normality
ggqqplot(my_data, “score”, ggtheme = theme_bw()) +
facet_grid(minute ~ treatment, labeller = “label_both”)

#if assumptions are met run repeated two way ANOVA
res.aov <- rstatix::anova_test(
data = my_data, dv = score, wid = id,
within = c(treatment, minute)
)
get_anova_table(res.aov)

data:

score id treatment minute
1 5.66666 1 ctrl t1
2 75.16670 3 ctrl t1
3 57.61120 5 ctrl t1
4 33.16670 7 ctrl t1
5 66.16660 9 ctrl t1
6 57.61110 11 ctrl t1
7 75.11110 13 ctrl t1
8 95.27780 15 ctrl t1
9 72.61110 17 ctrl t1
10 64.66670 19 ctrl t1
11 25.38890 21 ctrl t1
12 10.16670 23 ctrl t1
13 46.44440 4 Gs t1
14 3.05555 8 Gs t1
15 61.38880 12 Gs t1
16 74.61110 16 Gs t1
17 39.33340 18 Gs t1
18 81.55550 20 Gs t1
19 15.72220 22 Gs t1
20 49.72230 24 Gs t1
21 35.33330 25 Gs t1
22 37.72220 26 Gs t1
23 32.94440 1 ctrl t2
24 87.38890 3 ctrl t2
25 72.72230 5 ctrl t2
26 88.50000 7 ctrl t2
27 44.66660 9 ctrl t2
28 81.11120 11 ctrl t2
29 94.00000 13 ctrl t2
30 94.16650 15 ctrl t2
31 50.27760 17 ctrl t2
32 78.94450 19 ctrl t2
33 0.00000 21 ctrl t2
34 8.88888 23 ctrl t2
35 30.61110 4 Gs t2
36 0.00000 8 Gs t2
37 75.83330 12 Gs t2
38 69.77780 16 Gs t2
39 36.88890 18 Gs t2
40 3.05555 20 Gs t2
41 5.55555 22 Gs t2
42 29.88880 24 Gs t2
43 6.88888 25 Gs t2
44 23.27800 26 Gs t2

• Aline Mak

data.frame(
+ score = c(5.66666,75.1667,57.6112,33.1667,66.1666,
+ 57.6111,75.1111,95.2778,72.6111,64.6667,25.3889,10.1667,
+ 46.4444,3.05555,61.3888,74.6111,39.3334,81.5555,15.7222,
+ 49.7223,35.3333,37.7222,32.9444,87.3889,72.7223,88.5,
+ 44.6666,81.1112,94,94.1665,50.2776,78.9445,0,8.88888,
+ 30.6111,0,75.8333,69.7778,36.8889,3.05555,5.55555,29.8888,
+ 6.88888,23.278),
+ id = as.factor(c(“1″,”3″,”5″,”7″,”9”,
+ “11”,”13″,”15″,”17″,”19″,”21″,”23″,”4″,”8″,
+ “12”,”16″,”18″,”20″,”22″,”24″,”25″,”26″,”1″,
+ “3”,”5″,”7″,”9″,”11″,”13″,”15″,”17″,”19″,
+ “21”,”23″,”4″,”8″,”12″,”16″,”18″,”20″,”22″,
+ “24”,”25″,”26″)),
+ treatment = as.factor(c(“ctrl”,”ctrl”,”ctrl”,
+ “ctrl”,”ctrl”,”ctrl”,”ctrl”,”ctrl”,”ctrl”,
+ “ctrl”,”ctrl”,”ctrl”,”Gs”,”Gs”,”Gs”,”Gs”,”Gs”,
+ “Gs”,”Gs”,”Gs”,”Gs”,”Gs”,”ctrl”,”ctrl”,”ctrl”,
+ “ctrl”,”ctrl”,”ctrl”,”ctrl”,”ctrl”,”ctrl”,”ctrl”,
+ “ctrl”,”ctrl”,”Gs”,”Gs”,”Gs”,”Gs”,”Gs”,”Gs”,
+ “Gs”,”Gs”,”Gs”,”Gs”)),
+ minute = as.factor(c(“t1″,”t1″,”t1″,”t1”,
+ “t1″,”t1″,”t1″,”t1″,”t1″,”t1″,”t1″,”t1″,”t1”,
+ “t1″,”t1″,”t1″,”t1″,”t1″,”t1″,”t1″,”t1”,
+ “t1″,”t2″,”t2″,”t2″,”t2″,”t2″,”t2″,”t2″,”t2”,
+ “t2″,”t2″,”t2″,”t2″,”t2″,”t2″,”t2″,”t2″,”t2”,
+ “t2″,”t2″,”t2″,”t2″,”t2”))
+ )

• Aline

First of all, thank you Kassambara for this work! I am trying to run a two way repeated measurement ANOVA with my own data but I still seem to run into error described above by others;

– Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, …) :
0 (non-NA) cases

I have updated the lasted dev version by running devtools::install_github(“kassambara/rstatix”). However, it still will give me the error. Do you have any idea how to fix this? I would very much appreciate your help!

Data:

data.frame(
+ score = c(5.66666,75.1667,57.6112,33.1667,66.1666,
+ 57.6111,75.1111,95.2778,72.6111,64.6667,25.3889,10.1667,
+ 46.4444,3.05555,61.3888,74.6111,39.3334,81.5555,15.7222,
+ 49.7223,35.3333,37.7222,32.9444,87.3889,72.7223,88.5,
+ 44.6666,81.1112,94,94.1665,50.2776,78.9445,0,8.88888,
+ 30.6111,0,75.8333,69.7778,36.8889,3.05555,5.55555,29.8888,
+ 6.88888,23.278),
+ id = as.factor(c(“1″,”3″,”5″,”7″,”9”,
+ “11”,”13″,”15″,”17″,”19″,”21″,”23″,”4″,”8″,
+ “12”,”16″,”18″,”20″,”22″,”24″,”25″,”26″,”1″,
+ “3”,”5″,”7″,”9″,”11″,”13″,”15″,”17″,”19″,
+ “21”,”23″,”4″,”8″,”12″,”16″,”18″,”20″,”22″,
+ “24”,”25″,”26″)),
+ treatment = as.factor(c(“ctrl”,”ctrl”,”ctrl”,
+ “ctrl”,”ctrl”,”ctrl”,”ctrl”,”ctrl”,”ctrl”,
+ “ctrl”,”ctrl”,”ctrl”,”Gs”,”Gs”,”Gs”,”Gs”,”Gs”,
+ “Gs”,”Gs”,”Gs”,”Gs”,”Gs”,”ctrl”,”ctrl”,”ctrl”,
+ “ctrl”,”ctrl”,”ctrl”,”ctrl”,”ctrl”,”ctrl”,”ctrl”,
+ “ctrl”,”ctrl”,”Gs”,”Gs”,”Gs”,”Gs”,”Gs”,”Gs”,
+ “Gs”,”Gs”,”Gs”,”Gs”)),
+ minute = as.factor(c(“t1″,”t1″,”t1″,”t1”,
+ “t1″,”t1″,”t1″,”t1″,”t1″,”t1″,”t1″,”t1″,”t1”,
+ “t1″,”t1″,”t1″,”t1″,”t1″,”t1″,”t1″,”t1”,
+ “t1″,”t2″,”t2″,”t2″,”t2″,”t2″,”t2″,”t2″,”t2”,
+ “t2″,”t2″,”t2″,”t2″,”t2″,”t2″,”t2″,”t2″,”t2”,
+ “t2″,”t2″,”t2″,”t2″,”t2”))
+ )

code:
library(ggpubr)
devtools::install_github(“kassambara/rstatix”)
library(rstatix)

my_data % sample_n_by(treatment, size = 1)

# Gather the columns t1, t2 and t3 into long format.
# Convert id and time into factor variables
my_data %
gather(key = “minute”, value = “score”, t1, t2) %>%
convert_as_factor(id, minute, treatment)

#Inspect some random rows of the data by groups
set.seed(123)
my_data %>% sample_n_by(treatment, minute, size = 1)

#summery of data
my_data %>%
group_by(treatment, minute) %>%
get_summary_stats(score, type = “mean_sd”)

#check for outlayers
my_data %>%
group_by(treatment, minute) %>%
identify_outliers(score)

#check for normality
my_data %>%
group_by(treatment, minute) %>%
shapiro_test(score)

#QQplot also check for normality
ggqqplot(my_data, “score”, ggtheme = theme_bw()) +
facet_grid(minute ~ treatment, labeller = “label_both”)

suppressPackageStartupMessages(library(rstatix))

#if assumptions are met run repeated two way ANOVA
res.aov <- rstatix::anova_test(
data = my_data, dv = score, wid = id,
within = c(treatment, minute)
)
get_anova_table(res.aov)

• Bartek

Thank you for this useful tutorial. What about sutation when I detect outliers and want to delete them from the dataset but only from a specific condition (3×2 design)? I.e. do not want to lose whole subject if he is an outlier in only 1 of 6 conditions. If I remove him from this single condition, and get a significant effect in the 3-level factor, I cannot perform a paired t test with the method that you proposed – I get an error due to unequal number of observations (Error in complete.cases(x, y) : not all arguments have the same length). Any tips of what should I do in such a case?