Comparing Multiple Means in R

ANCOVA in R

The Analysis of Covariance (ANCOVA) is used to compare means of an outcome variable between two or more groups taking into account (or to correct for) variability of other variables, called covariates. In other words, ANCOVA allows to compare the adjusted means of two or more independent groups.

For example, you might want to compare “test score” by “level of education” taking into account the “number of hours spent studying”. In this example: 1) test score is our outcome (dependent) variable; 2) level of education (high school, college degree or graduate degree) is our grouping variable; 3) sudying time is our covariate.

The one-way ANCOVA can be seen as an extension of the one-way ANOVA that incorporate a covariate variable. The two-way ANCOVA is used to evaluate simultaneously the effect of two independent grouping variables (A and B) on an outcome variable, after adjusting for one or more continuous variables, called covariates.

In this article, you will learn how to:

  • Compute and interpret the one-way and the two-way ANCOVA in R
  • Check ANCOVA assumptions
  • Perform post-hoc tests, multiple pairwise comparisons between groups to identify which groups are different
  • Visualize the data using box plots, add ANCOVA and pairwise comparisons p-values to the plot

Contents:

Related Book

Practical Statistics in R II - Comparing Groups: Numerical Variables

Assumptions

ANCOVA makes several assumptions about the data, such as:

  • Linearity between the covariate and the outcome variable at each level of the grouping variable. This can be checked by creating a grouped scatter plot of the covariate and the outcome variable.
  • Homogeneity of regression slopes. The slopes of the regression lines, formed by the covariate and the outcome variable, should be the same for each group. This assumption evaluates that there is no interaction between the outcome and the covariate. The plotted regression lines by groups should be parallel.
  • The outcome variable should be approximately normally distributed. This can be checked using the Shapiro-Wilk test of normality on the model residuals.
  • Homoscedasticity or homogeneity of residuals variance for all groups. The residuals are assumed to have a constant variance (homoscedasticity)
  • No significant outliers in the groups

Many of these assumptions and potential problems can be checked by analyzing the residual errors. In the situation, where the ANCOVA assumption is not met you can perform robust ANCOVA test using the WRS2 package.

Prerequisites

Make sure you have installed the following R packages:

  • tidyverse for data manipulation and visualization
  • ggpubr for creating easily publication ready plots
  • rstatix for easy pipe-friendly statistical analyses
  • broom for printing a nice summary of statistical tests as data frames
  • datarium: contains required data sets for this chapter

Start by loading the following required packages:

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

One-way ANCOVA

Data preparation

We’ll prepare our demo data from the anxiety dataset available in the datarium package.

Researchers investigated the effect of exercises in reducing the level of anxiety. Therefore, they conducted an experiment, where they measured the anxiety score of three groups of individuals practicing physical exercises at different levels (grp1: low, grp2: moderate and grp3: high).

The anxiety score was measured pre- and 6-months post-exercise training programs. It is expected that any reduction in the anxiety by the exercises programs would also depend on the participant’s basal level of anxiety score.

In this analysis we use the pretest anxiety score as the covariate and are interested in possible differences between group with respect to the post-test anxiety scores.

# Load and prepare the data
data("anxiety", package = "datarium")
anxiety <- anxiety %>%
  select(id, group, t1, t3) %>%
  rename(pretest = t1, posttest = t3)
anxiety[14, "posttest"] <- 19
# Inspect the data by showing one random row by groups
set.seed(123)
anxiety %>% sample_n_by(group, size = 1)
## # A tibble: 3 x 4
##   id    group pretest posttest
##   <fct> <fct>   <dbl>    <dbl>
## 1 5     grp1     16.5     15.7
## 2 27    grp2     17.8     16.9
## 3 37    grp3     17.1     14.3

Check assumptions

Linearity assumption

  • Create a scatter plot between the covariate (i.e., pretest) and the outcome variable (i.e., posttest)
  • Add regression lines, show the corresponding equations and the R2 by groups
ggscatter(
  anxiety, x = "pretest", y = "posttest",
  color = "group", add = "reg.line"
  )+
  stat_regline_equation(
    aes(label =  paste(..eq.label.., ..rr.label.., sep = "~~~~"), color = group)
    )

There was a linear relationship between pre-test and post-test anxiety score for each training group, as assessed by visual inspection of a scatter plot.

Homogeneity of regression slopes

This assumption checks that there is no significant interaction between the covariate and the grouping variable. This can be evaluated as follow:

anxiety %>% anova_test(posttest ~ group*pretest)
## ANOVA Table (type II tests)
## 
##          Effect DFn DFd       F        p p<.05   ges
## 1         group   2  39 209.314 1.40e-21     * 0.915
## 2       pretest   1  39 572.828 6.36e-25     * 0.936
## 3 group:pretest   2  39   0.127 8.81e-01       0.006

There was homogeneity of regression slopes as the interaction term was not statistically significant, F(2, 39) = 0.13, p = 0.88.

Normality of residuals

You first need to compute the model using lm(). In R, you can easily augment your data to add fitted values and residuals by using the function augment(model) [broom package]. Let’s call the output model.metrics because it contains several metrics useful for regression diagnostics.

# Fit the model, the covariate goes first
model <- lm(posttest ~ pretest + group, data = anxiety)
# Inspect the model diagnostic metrics
model.metrics <- augment(model) %>%
  select(-.hat, -.sigma, -.fitted, -.se.fit) # Remove details
head(model.metrics, 3)
## # A tibble: 3 x 6
##   posttest pretest group .resid .cooksd .std.resid
##      <dbl>   <dbl> <fct>  <dbl>   <dbl>      <dbl>
## 1     14.1    14.1 grp1   0.550  0.101       1.46 
## 2     14.3    14.5 grp1   0.338  0.0310      0.885
## 3     14.9    15.7 grp1  -0.295  0.0133     -0.750
# Assess normality of residuals using shapiro wilk test
shapiro_test(model.metrics$.resid)
## # A tibble: 1 x 3
##   variable             statistic p.value
##   <chr>                    <dbl>   <dbl>
## 1 model.metrics$.resid     0.975   0.444

The Shapiro Wilk test was not significant (p > 0.05), so we can assume normality of residuals

Homogeneity of variances

ANCOVA assumes that the variance of the residuals is equal for all groups. This can be checked using the Levene’s test:

model.metrics %>% levene_test(.resid ~ group)
## # A tibble: 1 x 4
##     df1   df2 statistic     p
##   <int> <int>     <dbl> <dbl>
## 1     2    42      2.27 0.116

The Levene’s test was not significant (p > 0.05), so we can assume homogeneity of the residual variances for all groups.

Outliers

An outlier is a point that has an extreme outcome variable value. The presence of outliers may affect the interpretation of the model.

Outliers can be identified by examining the standardized residual (or studentized residual), which is the residual divided by its estimated standard error. Standardized residuals can be interpreted as the number of standard errors away from the regression line.

Observations whose standardized residuals are greater than 3 in absolute value are possible outliers.

model.metrics %>% 
  filter(abs(.std.resid) > 3) %>%
  as.data.frame()
## [1] posttest   pretest    group      .resid     .cooksd    .std.resid
## <0 rows> (or 0-length row.names)

There were no outliers in the data, as assessed by no cases with standardized residuals greater than 3 in absolute value.

Computation

The orders of variables matters when computing ANCOVA. You want to remove the effect of the covariate first - that is, you want to control for it - prior to entering your main variable or interest.

The covariate goes first (and there is no interaction)! If you do not do this in order, you will get different results.

res.aov <- anxiety %>% anova_test(posttest ~ pretest + group)
get_anova_table(res.aov)
## ANOVA Table (type II tests)
## 
##    Effect DFn DFd   F        p p<.05   ges
## 1 pretest   1  41 598 4.48e-26     * 0.936
## 2   group   2  41 219 1.35e-22     * 0.914

After adjustment for pre-test anxiety score, there was a statistically significant difference in post-test anxiety score between the groups, F(2, 41) = 218.63, p < 0.0001.

Post-hoc test

Pairwise comparisons can be performed to identify which groups are different. The Bonferroni multiple testing correction is applied. This can be easily done using the function emmeans_test() [rstatix package], a wrapper around the emmeans package, which needs to be installed. Emmeans stands for estimated marginal means (aka least square means or adjusted means).

# Pairwise comparisons
library(emmeans)
pwc <- anxiety %>% 
  emmeans_test(
    posttest ~ group, covariate = pretest,
    p.adjust.method = "bonferroni"
    )
pwc
## # A tibble: 3 x 8
##   .y.      group1 group2    df statistic        p    p.adj p.adj.signif
## * <chr>    <chr>  <chr>  <dbl>     <dbl>    <dbl>    <dbl> <chr>       
## 1 posttest grp1   grp2      41      4.24 1.26e- 4 3.77e- 4 ***         
## 2 posttest grp1   grp3      41     19.9  1.19e-22 3.58e-22 ****        
## 3 posttest grp2   grp3      41     15.5  9.21e-19 2.76e-18 ****
# Display the adjusted means of each group
# Also called as the estimated marginal means (emmeans)
get_emmeans(pwc)
## # A tibble: 3 x 8
##   pretest group emmean    se    df conf.low conf.high method      
##     <dbl> <fct>  <dbl> <dbl> <dbl>    <dbl>     <dbl> <chr>       
## 1    16.9 grp1    16.4 0.106    41     16.2      16.7 Emmeans test
## 2    16.9 grp2    15.8 0.107    41     15.6      16.0 Emmeans test
## 3    16.9 grp3    13.5 0.106    41     13.2      13.7 Emmeans test

Data are adjusted mean +/- standard error. The mean anxiety score was statistically significantly greater in grp1 (16.4 +/- 0.15) compared to the grp2 (15.8 +/- 0.12) and grp3 (13.5 +/_ 0.11), p < 0.001.

Report

An ANCOVA was run to determine the effect of exercises on the anxiety score after controlling for basal anxiety score of participants.

After adjustment for pre-test anxiety score, there was a statistically significant difference in post-test anxiety score between the groups, F(2, 41) = 218.63, p < 0.0001.

Post hoc analysis was performed with a Bonferroni adjustment. The mean anxiety score was statistically significantly greater in grp1 (16.4 +/- 0.15) compared to the grp2 (15.8 +/- 0.12) and grp3 (13.5 +/_ 0.11), p < 0.001.

# Visualization: line plots with p-values
pwc <- pwc %>% add_xy_position(x = "group", fun = "mean_se")
ggline(get_emmeans(pwc), x = "group", y = "emmean") +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2) + 
  stat_pvalue_manual(pwc, hide.ns = TRUE, tip.length = FALSE) +
  labs(
    subtitle = get_test_label(res.aov, detailed = TRUE),
    caption = get_pwc_label(pwc)
  )

Two-way ANCOVA

Data preparation

We’ll use the stress dataset available in the datarium package. In this study, a researcher wants to evaluate the effect of treatment and exercise on stress reduction score after adjusting for age.

In this example: 1) stress score is our outcome (dependent) variable; 2) treatment (levels: no and yes) and exercise (levels: low, moderate and high intensity training) are our grouping variable; 3) age is our covariate.

Load the data and show some random rows by groups:

data("stress", package = "datarium")
stress %>% sample_n_by(treatment, exercise)
## # A tibble: 6 x 5
##      id score treatment exercise   age
##   <int> <dbl> <fct>     <fct>    <dbl>
## 1     8  83.8 yes       low         61
## 2    15  86.9 yes       moderate    55
## 3    29  71.5 yes       high        55
## 4    40  92.4 no        low         67
## 5    41 100   no        moderate    75
## 6    56  82.4 no        high        53

Check assumptions

Linearity assumption

  • Create a scatter plot between the covariate (i.e., age) and the outcome variable (i.e., score) for each combination of the groups of the two grouping variables:
  • Add smoothed loess lines, which helps to decide if the relationship is linear or not
ggscatter(
  stress, x = "age", y = "score",
  facet.by  = c("exercise", "treatment"), 
  short.panel.labs = FALSE
  )+
  stat_smooth(method = "loess", span = 0.9)

There was a linear relationship between the covariate (age variable) and the outcome variable (score) for each group, as assessed by visual inspection of a scatter plot.

Homogeneity of regression slopes

This assumption checks that there is no significant interaction between the covariate and the grouping variables. This can be evaluated as follow:

stress %>%
  anova_test(
    score ~ age + treatment + exercise + 
     treatment*exercise + age*treatment +
     age*exercise + age*exercise*treatment
  )
## ANOVA Table (type II tests)
## 
##                   Effect DFn DFd      F        p p<.05      ges
## 1                    age   1  48  8.359 6.00e-03     * 0.148000
## 2              treatment   1  48  9.907 3.00e-03     * 0.171000
## 3               exercise   2  48 18.197 1.31e-06     * 0.431000
## 4     treatment:exercise   2  48  3.303 4.50e-02     * 0.121000
## 5          age:treatment   1  48  0.009 9.25e-01       0.000189
## 6           age:exercise   2  48  0.235 7.91e-01       0.010000
## 7 age:treatment:exercise   2  48  0.073 9.30e-01       0.003000

Another simple alternative is to create a new grouping variable, say group, based on the combinations of the existing variables, and then compute ANOVA model:

stress %>%
  unite(col = "group", treatment, exercise) %>%
  anova_test(score ~ group*age)
## ANOVA Table (type II tests)
## 
##      Effect DFn DFd      F        p p<.05   ges
## 1     group   5  48 10.912 4.76e-07     * 0.532
## 2       age   1  48  8.359 6.00e-03     * 0.148
## 3 group:age   5  48  0.126 9.86e-01       0.013

There was homogeneity of regression slopes as the interaction terms, between the covariate (age) and grouping variables (treatment and exercise), was not statistically significant, p > 0.05.

Normality of residuals

# Fit the model, the covariate goes first
model <- lm(score ~ age + treatment*exercise, data = stress)
# Inspect the model diagnostic metrics
model.metrics <- augment(model) %>%
  select(-.hat, -.sigma, -.fitted, -.se.fit) # Remove details
head(model.metrics, 3)
## # A tibble: 3 x 7
##   score   age treatment exercise .resid .cooksd .std.resid
##   <dbl> <dbl> <fct>     <fct>     <dbl>   <dbl>      <dbl>
## 1  95.6    59 yes       low        9.10  0.0647       1.93
## 2  82.2    65 yes       low       -7.32  0.0439      -1.56
## 3  97.2    70 yes       low        5.16  0.0401       1.14
# Assess normality of residuals using shapiro wilk test
shapiro_test(model.metrics$.resid)
## # A tibble: 1 x 3
##   variable             statistic p.value
##   <chr>                    <dbl>   <dbl>
## 1 model.metrics$.resid     0.982   0.531

The Shapiro Wilk test was not significant (p > 0.05), so we can assume normality of residuals

Homogeneity of variances

ANCOVA assumes that the variance of the residuals is equal for all groups. This can be checked using the Levene’s test:

levene_test(.resid ~ treatment*exercise, data = model.metrics)

The Levene’s test was not significant (p > 0.05), so we can assume homogeneity of the residual variances for all groups.

Outliers

Observations whose standardized residuals are greater than 3 in absolute value are possible outliers.

model.metrics %>% 
  filter(abs(.std.resid) > 3) %>%
  as.data.frame()
## [1] score      age        treatment  exercise   .resid     .cooksd    .std.resid
## <0 rows> (or 0-length row.names)

There were no outliers in the data, as assessed by no cases with standardized residuals greater than 3 in absolute value.

Computation

res.aov <- stress %>% 
  anova_test(score ~ age + treatment*exercise)
get_anova_table(res.aov)
## ANOVA Table (type II tests)
## 
##               Effect DFn DFd     F        p p<.05   ges
## 1                age   1  53  9.11 4.00e-03     * 0.147
## 2          treatment   1  53 11.10 2.00e-03     * 0.173
## 3           exercise   2  53 20.82 2.13e-07     * 0.440
## 4 treatment:exercise   2  53  4.45 1.60e-02     * 0.144

After adjustment for age, there was a statistically significant interaction between treatment and exercise on the stress score, F(2, 53) = 4.45, p = 0.016. This indicates that the effect of exercise on score depends on the level of exercise, and vice-versa.

Post-hoc test

A statistically significant two-way interactions can be followed up by simple main effect analyses, that is evaluating the effect of one variable at each level of the second variable, and vice-versa.

In the situation, where the interaction is not significant, you can report the main effect of each grouping variable.

A significant two-way interaction indicates that the impact that one factor has on the outcome variable depends on the level of the other factor (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 ANCOVA output.

In this section we’ll describe the procedure for a significant three-way interaction

Simple main effect analyses for treatment

Analyze the simple main effect of treatment at each level of exercise. Group the data by exercise and perform one-way ANCOVA for treatment controlling for age:

# Effect of treatment at each level of exercise
stress %>%
  group_by(exercise) %>%
  anova_test(score ~ age + treatment)
## # A tibble: 6 x 8
##   exercise Effect      DFn   DFd      F        p `p<.05`   ges
##   <fct>    <chr>     <dbl> <dbl>  <dbl>    <dbl> <chr>   <dbl>
## 1 low      age           1    17  2.25  0.152    ""      0.117
## 2 low      treatment     1    17  0.437 0.517    ""      0.025
## 3 moderate age           1    17  6.65  0.02     *       0.281
## 4 moderate treatment     1    17  0.419 0.526    ""      0.024
## 5 high     age           1    17  0.794 0.385    ""      0.045
## 6 high     treatment     1    17 18.7   0.000455 *       0.524

Note that, we need to apply Bonferroni adjustment for multiple testing corrections. One common approach is lowering the level at which you declare significance by dividing the alpha value (0.05) by the number of tests performed. In our example, that is 0.05/3 = 0.016667.

Statistical significance was accepted at the Bonferroni-adjusted alpha level of 0.01667, that is 0.05/3. The effect of treatment was statistically significant in the high-intensity exercise group (p = 0.00045), but not in the low-intensity exercise group (p = 0.517) and in the moderate-intensity exercise group (p = 0.526).

Compute pairwise comparisons between treatment groups at each level of exercise. The Bonferroni multiple testing correction is applied.

# Pairwise comparisons
pwc <- stress %>% 
  group_by(exercise) %>%
  emmeans_test(
    score ~ treatment, covariate = age,
    p.adjust.method = "bonferroni"
    )
pwc %>% filter(exercise == "high")
## # A tibble: 1 x 9
##   exercise .y.   group1 group2    df statistic         p     p.adj p.adj.signif
##   <fct>    <chr> <chr>  <chr>  <dbl>     <dbl>     <dbl>     <dbl> <chr>       
## 1 high     score yes    no        53     -4.36 0.0000597 0.0000597 ****

In the pairwise comparison table, you will only need the result for “exercises:high” group, as this was the only condition where the simple main effect of treatment was statistically significant.

The pairwise comparisons between treatment:no and treatment:yes group was statistically significant in participant undertaking high-intensity exercise (p < 0.0001).

Simple main effect for exercise

You can do the same post-hoc analyses for the exercise variable at each level of treatment variable.

# Effect of exercise at each level of treatment
stress %>%
  group_by(treatment) %>%
  anova_test(score ~ age + exercise)
## # A tibble: 4 x 8
##   treatment Effect     DFn   DFd     F         p `p<.05`   ges
##   <fct>     <chr>    <dbl> <dbl> <dbl>     <dbl> <chr>   <dbl>
## 1 yes       age          1    26  2.37 0.136     ""      0.083
## 2 yes       exercise     2    26 17.3  0.0000164 *       0.572
## 3 no        age          1    26  7.26 0.012     *       0.218
## 4 no        exercise     2    26  3.99 0.031     *       0.235

Statistical significance was accepted at the Bonferroni-adjusted alpha level of 0.025, that is 0.05/2 (the number of tests). The effect of exercise was statistically significant in the treatment=yes group (p < 0.0001), but not in the treatment=no group (p = 0.031).

Perform multiple pairwise comparisons between exercise groups at each level of treatment. You don’t need to interpret the results for the “no treatment” group, because the effect of exercise was not significant for this group.

pwc2 <- stress %>% 
  group_by(treatment) %>%
  emmeans_test(
    score ~ exercise, covariate = age,
    p.adjust.method = "bonferroni"
    ) %>%
  select(-df, -statistic, -p) # Remove details
pwc2 %>% filter(treatment == "yes")
## # A tibble: 3 x 6
##   treatment .y.   group1   group2         p.adj p.adj.signif
##   <fct>     <chr> <chr>    <chr>          <dbl> <chr>       
## 1 yes       score low      moderate 1           ns          
## 2 yes       score low      high     0.00000113  ****        
## 3 yes       score moderate high     0.000000466 ****

There was a statistically significant difference between the adjusted mean of low and high exercise group (p < 0.0001) and, between moderate and high group (p < 0.0001). The difference between the adjusted means of low and moderate was not significant.

Report

A two-way ANCOVA was performed to examine the effects of treatment and exercise on stress reduction, after controlling for age.

There was a statistically significant two-way interaction between treatment and exercise on score concentration, whilst controlling for age, F(2, 53) = 4.45, p = 0.016.

Therefore, an analysis of simple main effects for exercise and treatment was performed with statistical significance receiving a Bonferroni adjustment and being accepted at the p < 0.025 level for exercise and p < 0.0167 for treatment.

The simple main effect of treatment was statistically significant in the high-intensity exercise group (p = 0.00046), but not in the low-intensity exercise group (p = 0.52) and the moderate-intensity exercise group (p = 0.53).

The effect of exercise was statistically significant in the treatment=yes group (p < 0.0001), but not in the treatment=no group (p = 0.031).

All pairwise comparisons were computed for statistically significant simple main effects with reported p-values Bonferroni adjusted. For the treatment=yes group, there was a statistically significant difference between the adjusted mean of low and high exercise group (p < 0.0001) and, between moderate and high group (p < 0.0001). The difference between the adjusted means of low and moderate exercise groups was not significant.

  • Create a line plot:
# Line plot
lp <- ggline(
  get_emmeans(pwc), x = "exercise", y = "emmean", 
  color = "treatment", palette = "jco"
  ) +
  geom_errorbar(
    aes(ymin = conf.low, ymax = conf.high, color = treatment), 
    width = 0.1
    )
  • Add p-values
# Comparisons between treatment group at each exercise level
pwc <- pwc %>% add_xy_position(x = "exercise", fun = "mean_se", step.increase = 0.2)
pwc.filtered <- pwc %>% filter(exercise == "high")
lp + 
stat_pvalue_manual(
  pwc.filtered, hide.ns = TRUE, tip.length = 0,
  bracket.size = 0
  ) +
labs(
  subtitle = get_test_label(res.aov,  detailed = TRUE),
  caption = get_pwc_label(pwc)
)

# Comparisons between exercises group at each treatment level
pwc2 <- pwc2 %>% add_xy_position(x = "exercise", fun = "mean_se")
pwc2.filtered <- pwc2 %>% filter(treatment == "yes")
lp + 
stat_pvalue_manual(
  pwc2.filtered, hide.ns = TRUE, tip.length = 0,
  step.group.by = "treatment", color = "treatment"
  ) +
labs(
  subtitle = get_test_label(res.aov,  detailed = TRUE),
  caption = get_pwc_label(pwc2)
)

Summary

This article describes how to compute and interpret one-way and two-way ANCOVA in R. We also explain the assumptions made by ANCOVA tests and provide practical examples of R codes to check whether the test assumptions are met or not.

Mixed ANOVA in R (Prev Lesson)
(Next Lesson) One-Way MANOVA in R
Back to Comparing Multiple Means in R

No Comments

Give a comment

Teacher
Alboukadel Kassambara
Role : Founder of Datanovia
Read More