Analysis of covariance the modern way — compare adjusted group means while controlling for a continuous covariate, check the ANCOVA assumptions, and run the post-hoc comparisons
Learn to run a one-way and two-way ANCOVA (analysis of covariance) in R with rstatix. Compare adjusted group means while controlling for a covariate, check every ANCOVA assumption — linearity, homogeneity of regression slopes, normality of residuals, equal variance and outliers — compute the adjusted (estimated marginal) means, run Bonferroni post-hoc comparisons, and put the p-values on the plot.
anova_test(outcome ~ covariate + group) — the covariate goes first (you remove its effect before testing the group), then read the adjusted means with get_emmeans().emmeans_test() pairwise comparisons on the adjusted means, Bonferroni-corrected — then plot the adjusted means with the p-values on them.ANCOVA — analysis of covariance — compares the means of an outcome variable across groups while controlling for one or more continuous variables, called covariates. It blends ANOVA and regression: the covariate’s effect is removed first, and the groups are compared on the adjusted means that remain. A one-way ANCOVA is a one-way ANOVA plus a covariate; a two-way ANCOVA adds a second factor and the interaction between them.
The classic scenario: a team tests whether an exercise programme reduces anxiety, but participants start at different baseline anxiety levels — and any improvement plausibly depends on where you began. ANCOVA compares the post-test anxiety of the exercise groups after adjusting for that baseline, so the group difference isn’t just an echo of who started calmer. This lesson runs the full workflow with rstatix (with the base-R aov() alongside), using the anxiety and stress datasets, and covers both the one-way and two-way cases.
ANCOVA hypotheses (one-way). H₀: after adjusting for the covariate, the adjusted group means are equal. Hₐ: at least one adjusted group mean differs. The covariate is not the thing being tested — it is removed so the group comparison is fair.
ANCOVA assumptions — check these before trusting the test. Beyond independent observations: a linear relationship between covariate and outcome in each group; homogeneity of regression slopes (the covariate × group interaction is not significant — the lines are parallel); normally distributed residuals; equal residual variance across groups (homoscedasticity); and no extreme outliers (standardized residuals within ±3). If they fail badly, use a robust ANCOVA (WRS2 package).
We use the anxiety dataset from datarium: anxiety score for people in three physical-exercise groups (grp1 low, grp2 moderate, grp3 high), measured before and after a six-month programme. We take the pre-test score as the covariate (pretest) and the post-test score as the outcome (posttest), and ask whether the groups differ at post-test once baseline is accounted for. Prepare the columns with base R:
library(rstatix)
data("anxiety", package = "datarium")
anx <- as.data.frame(anxiety[, c("id", "group", "t1", "t3")])
names(anx)[names(anx) %in% c("t1", "t3")] <- c("pretest", "posttest")
anx[14, "posttest"] <- 19 # one edited value, as in the source example
set.seed(123)
anx %>% sample_n_by(group, size = 1)# A tibble: 3 × 4
id group pretest posttest
<fct> <fct> <dbl> <dbl>
1 15 grp1 19.8 19.4
2 30 grp2 19.3 17.7
3 33 grp3 15.5 11 ANCOVA assumes the covariate and outcome are linearly related within each group. A grouped scatter plot with a regression line per group is the check — roughly straight, similarly-sloped lines are what you want:
library(rstatix)
library(ggpubr)
data("anxiety", package = "datarium")
anx <- as.data.frame(anxiety[, c("id", "group", "t1", "t3")])
names(anx)[names(anx) %in% c("t1", "t3")] <- c("pretest", "posttest")
anx[14, "posttest"] <- 19
ggscatter(
anx, x = "pretest", y = "posttest", color = "group", palette = "jco",
add = "reg.line"
) +
stat_regline_equation(aes(label = after_stat(rr.label), color = group))
The relationship is linear in every group.
The slopes of those lines must be parallel — formally, there is no covariate × group interaction. Test it by adding the interaction term and confirming it is not significant:
library(rstatix)
data("anxiety", package = "datarium")
anx <- as.data.frame(anxiety[, c("id", "group", "t1", "t3")])
names(anx)[names(anx) %in% c("t1", "t3")] <- c("pretest", "posttest")
anx[14, "posttest"] <- 19
anx %>% 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.006The group:pretest interaction is not significant, F(2, 39) = 0.13, p = 0.88 — the slopes are homogeneous, so ANCOVA is valid.
Fit the ANCOVA model with lm() (covariate first), then derive the residuals and standardized residuals with base R (residuals() and rstandard() — no extra package). Shapiro-Wilk tests normality, Levene’s test checks equal residual variance across groups, and any |std.resid| > 3 flags an outlier:
library(rstatix)
data("anxiety", package = "datarium")
anx <- as.data.frame(anxiety[, c("id", "group", "t1", "t3")])
names(anx)[names(anx) %in% c("t1", "t3")] <- c("pretest", "posttest")
anx[14, "posttest"] <- 19
model <- lm(posttest ~ pretest + group, data = anx)
model.metrics <- data.frame(
group = anx$group,
resid = residuals(model),
std.resid = rstandard(model)
)
shapiro_test(model.metrics$resid) # normality of residuals# A tibble: 1 × 3
variable statistic p.value
<chr> <dbl> <dbl>
1 model.metrics$resid 0.975 0.444model.metrics %>% levene_test(resid ~ group) # equal variance# A tibble: 1 × 4
df1 df2 statistic p
<int> <int> <dbl> <dbl>
1 2 42 2.27 0.116model.metrics[abs(model.metrics$std.resid) > 3, ] # outliers (|std.resid| > 3)[1] group resid std.resid
<0 lignes> (ou 'row.names' de longueur nulle)Shapiro-Wilk is non-significant (p = 0.44, residuals normal), Levene’s test is non-significant (p = 0.12, equal variance), and no row has a standardized residual beyond ±3 — all assumptions hold.
The order of terms matters: put the covariate first. You want to remove the covariate’s effect before testing the group, so write outcome ~ covariate + group.
library(rstatix)
data("anxiety", package = "datarium")
anx <- as.data.frame(anxiety[, c("id", "group", "t1", "t3")])
names(anx)[names(anx) %in% c("t1", "t3")] <- c("pretest", "posttest")
anx[14, "posttest"] <- 19
res.aov <- anx %>% 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.321 4.48e-26 * 0.936
2 group 2 41 218.629 1.35e-22 * 0.914After adjusting for the pre-test score, there is a statistically significant difference in post-test anxiety between the groups, F(2, 41) = 218.63, p < 0.0001.
aov() gives the same adjusted group effect — provided the covariate is entered first (base R uses sequential, Type-I sums of squares, so order matters):
data("anxiety", package = "datarium")
anx <- as.data.frame(anxiety[, c("id", "group", "t1", "t3")])
names(anx)[names(anx) %in% c("t1", "t3")] <- c("pretest", "posttest")
anx[14, "posttest"] <- 19
summary(aov(posttest ~ pretest + group, data = anx)) Df Sum Sq Mean Sq F value Pr(>F)
pretest 1 99.40 99.40 587.2 <2e-16 ***
group 2 74.02 37.01 218.6 <2e-16 ***
Residuals 41 6.94 0.17
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The group row reads F = 218.63 — identical to the rstatix result. This match is not automatic:
Put the covariate first, or you get a different answer. Base R’s aov() adjusts each term only for the terms before it. With posttest ~ pretest + group the group effect is adjusted for the covariate (F = 218.63, matching rstatix’s Type-II anova_test()). Flip the order to posttest ~ group + pretest and the group F changes to 213.05 — the group is no longer adjusted for baseline, which is not an ANCOVA. Always write outcome ~ covariate + group.
A significant group effect is followed up with pairwise comparisons of the adjusted (estimated marginal) means — emmeans_test() with the covariate argument, Bonferroni-corrected. get_emmeans() prints the adjusted means themselves:
library(rstatix)
library(emmeans)
data("anxiety", package = "datarium")
anx <- as.data.frame(anxiety[, c("id", "group", "t1", "t3")])
names(anx)[names(anx) %in% c("t1", "t3")] <- c("pretest", "posttest")
anx[14, "posttest"] <- 19
pwc <- anx %>%
emmeans_test(posttest ~ group, covariate = pretest, p.adjust.method = "bonferroni")
pwc# A tibble: 3 × 9
term .y. group1 group2 df statistic p p.adj p.adj.signif
* <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <chr>
1 pretest*gr… post… grp1 grp2 41 4.24 1.26e- 4 3.77e- 4 ***
2 pretest*gr… post… grp1 grp3 41 19.9 1.19e-22 3.58e-22 ****
3 pretest*gr… post… grp2 grp3 41 15.5 9.21e-19 2.76e-18 **** get_emmeans(pwc)# A tibble: 3 × 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 testThe adjusted mean anxiety score was significantly greater in grp1 (≈ 16.4) than in grp2 (≈ 15.8) and grp3 (≈ 13.5), and all pairwise differences are significant (p.adj < 0.001) — the higher the exercise intensity, the lower the adjusted anxiety.
The reporting figure plots the adjusted means (not the raw means) with their confidence intervals, the ANCOVA result in the subtitle, and the significant comparisons as brackets:
library(rstatix)
library(ggpubr)
library(emmeans)
data("anxiety", package = "datarium")
anx <- as.data.frame(anxiety[, c("id", "group", "t1", "t3")])
names(anx)[names(anx) %in% c("t1", "t3")] <- c("pretest", "posttest")
anx[14, "posttest"] <- 19
res.aov <- anx %>% anova_test(posttest ~ pretest + group)
pwc <- anx %>%
emmeans_test(posttest ~ group, covariate = pretest, p.adjust.method = "bonferroni") %>%
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)
)
An ANCOVA was run to determine the effect of exercise group on post-test anxiety after controlling for baseline anxiety. After adjustment for the pre-test score, there was a statistically significant difference in post-test anxiety between the groups, F(2, 41) = 218.63, p < 0.0001. Post-hoc analysis with a Bonferroni adjustment showed all pairwise differences were significant (p < 0.001); the adjusted mean anxiety score was highest in grp1 (≈ 16.4) and lowest in grp3 (≈ 13.5).
A two-way ANCOVA tests two factors and their interaction on the outcome while adjusting for a covariate. We use the stress dataset (datarium): stress-reduction score by treatment (no, yes) and exercise (low, moderate, high), adjusting for age.
library(rstatix)
data("stress", package = "datarium")
stress %>% sample_n_by(treatment, exercise, size = 1)# A tibble: 6 × 5
id score treatment exercise age
<dbl> <dbl> <fct> <fct> <dbl>
1 3 97.2 yes low 70
2 20 85.4 yes moderate 57
3 22 65.8 yes high 56
4 36 87 no low 63
5 45 88.4 no moderate 65
6 54 75.7 no high 58The assumptions are the same, checked the same way (a faceted scatter for linearity, a non-significant covariate × factor interaction for homogeneity of slopes, Shapiro/Levene on the residuals, |std.resid| within 3). The homogeneity-of-slopes check crosses the covariate with both factors:
library(rstatix)
data("stress", package = "datarium")
stress %>% anova_test(
score ~ age + treatment + exercise +
treatment:exercise + age:treatment + age:exercise + age:treatment:exercise
)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.003000None of the age:* interaction terms is significant (p > 0.05), so the slopes are homogeneous.
Covariate first, then the two factors and their interaction:
library(rstatix)
data("stress", package = "datarium")
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.110 4.00e-03 * 0.147
2 treatment 1 53 11.096 2.00e-03 * 0.173
3 exercise 2 53 20.820 2.13e-07 * 0.440
4 treatment:exercise 2 53 4.446 1.60e-02 * 0.144After adjusting for age, there is a statistically significant treatment:exercise interaction, F(2, 53) = 4.45, p = 0.016 — the effect of exercise on stress depends on treatment, so (exactly as in a two-way ANOVA) we decompose it with simple main effects rather than reading the main effects alone.
Analyse the effect of treatment within each exercise level — group by exercise, run a one-way ANCOVA of treatment adjusting for age:
library(rstatix)
data("stress", package = "datarium")
stress %>%
group_by(exercise) %>%
anova_test(score ~ age + treatment)# A tibble: 6 × 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.524At a Bonferroni-adjusted level of 0.05 / 3 = 0.0167, the effect of treatment is significant only in the high-intensity exercise group (p = 0.00045), not in the low (p = 0.52) or moderate (p = 0.53) groups. Pairwise comparisons confirm where treatment matters:
library(rstatix)
library(emmeans)
data("stress", package = "datarium")
pwc <- stress %>%
group_by(exercise) %>%
emmeans_test(score ~ treatment, covariate = age, p.adjust.method = "bonferroni")
pwc[pwc$exercise == "high", c("exercise", "group1", "group2", "p.adj", "p.adj.signif")]# A tibble: 1 × 5
exercise group1 group2 p.adj p.adj.signif
<fct> <chr> <chr> <dbl> <chr>
1 high yes no 0.0000597 **** For high-intensity exercisers, treatment yes differs significantly from no (p.adj < 0.0001) — so the treatment helps most for those already exercising hard.
The decomposition is symmetric — you can also test the effect of exercise within each treatment level. Group by treatment and run a one-way ANCOVA of exercise adjusting for age:
library(rstatix)
data("stress", package = "datarium")
stress %>%
group_by(treatment) %>%
anova_test(score ~ age + exercise)# A tibble: 4 × 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.235At a Bonferroni-adjusted level of 0.05 / 2 = 0.025, exercise has a significant effect in the treatment = yes group (F(2, 26) = 17.3, p < 0.0001) but not the no group (p = 0.031) — the two simple-main-effect views together describe where the interaction lives.
A two-way ANCOVA examined the effects of treatment and exercise on stress reduction, controlling for age. There was a statistically significant interaction between treatment and exercise, F(2, 53) = 4.45, p = 0.016. Simple main-effect analyses (Bonferroni-adjusted) showed the effect of treatment was significant in the high-intensity exercise group (p < 0.001) but not the low (p = 0.52) or moderate (p = 0.53) groups.
Plot the adjusted means by exercise, coloured by treatment, with the significant comparison marked — the reporting figure for a two-way ANCOVA:
library(rstatix)
library(ggpubr)
library(emmeans)
data("stress", package = "datarium")
res.aov <- stress %>% anova_test(score ~ age + treatment * exercise)
pwc <- stress %>%
group_by(exercise) %>%
emmeans_test(score ~ treatment, covariate = age, p.adjust.method = "bonferroni") %>%
add_xy_position(x = "exercise", fun = "mean_se", step.increase = 0.2)
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) +
stat_pvalue_manual(pwc[pwc$exercise == "high", ], hide.ns = TRUE, tip.length = 0, bracket.size = 0) +
labs(
subtitle = get_test_label(res.aov, detailed = TRUE),
caption = get_pwc_label(pwc)
)
ANCOVA fits a linear model with the covariate and the factor together, \(y_{ij} = \mu + \tau_i + \beta\,(x_{ij} - \bar{x}) + \varepsilon_{ij}\), where \(\tau_i\) is the group effect and \(\beta\) is the common within-group slope on the covariate \(x\). The group is tested after the covariate has absorbed its share of the variation — the adjusted-mean comparison. The homogeneity-of-slopes assumption is exactly the claim that a single \(\beta\) fits all groups (no group-specific slope), which is why a significant covariate × group interaction invalidates the model.
Run a one-way ANCOVA on a fresh question — does mpg differ by number of cylinders once we adjust for car weight? In mtcars, wt is the covariate, cyl the group. Edit and re-run; the sandbox boots on first Run.
Ask Prova “I want to compare group means but I have a covariate to control for — run a one-way ANCOVA on my data, check the homogeneity-of-slopes assumption, and give me the adjusted means” — it answers with rstatix code you can run on your own data, then helps you read the adjusted means and the post-hoc comparisons. The runtime is the judge. Ask Prova →
You entered the covariate after the group. ANCOVA adjusts the group effect for the covariate, which in base R’s sequential aov() only happens if the covariate is listed first (outcome ~ covariate + group). The wrong order silently gives a different, un-adjusted group F. rstatix anova_test() defaults to Type-II so it is less order-sensitive, but write the covariate first regardless for clarity.
The homogeneity-of-slopes assumption fails. If the covariate × group interaction is significant, the groups have different slopes and a single adjusted-mean comparison is misleading — the group difference depends on the covariate value. Don’t run a standard ANCOVA; model the interaction explicitly (a moderation analysis) or compare groups at specific covariate values.
You reported the raw group means. Post-hoc comparisons and the reporting plot use the adjusted (estimated marginal) means from get_emmeans(), not mean() of the raw outcome — the adjusted means are the whole point of ANCOVA.
ANCOVA (analysis of covariance) compares group means on an outcome while controlling for a continuous covariate. Use it when a nuisance or baseline variable (pre-test score, age, dose) influences the outcome and you want a fair group comparison adjusted for it — e.g. comparing post-treatment scores across groups that started at different baselines.
ANOVA compares group means directly; ANCOVA compares adjusted group means after removing the effect of one or more continuous covariates. Mechanically, ANCOVA adds the covariate to the model (outcome ~ covariate + group) and tests the group effect on the residual, covariate-adjusted variation — which usually increases power by shrinking the error term.
It requires the covariate–outcome relationship to have the same slope in every group — i.e. no covariate × group interaction. Test it by fitting outcome ~ group * covariate and confirming the interaction term is not significant. If it is significant, the groups respond to the covariate differently and a single adjusted-mean comparison is not valid.
Base R’s aov() uses sequential (Type-I) sums of squares, adjusting each term only for the terms before it. Writing outcome ~ covariate + group adjusts the group for the covariate (true ANCOVA); the reverse order does not. rstatix anova_test() defaults to Type-II, which is order-insensitive, but listing the covariate first keeps the intent unambiguous.
They are the predicted group means at a common value of the covariate (usually its overall mean) — the means you would expect if every group had the same average baseline. get_emmeans() returns them, and the post-hoc comparisons and reporting plot use them rather than the raw group means.
posttest anxiety by group, adjusting for pretest. Fill the blanks (covariate first!). Is the group effect significant after adjustment?Fill the blanks so the formula is posttest ~ pretest + group — the covariate (pretest) comes first, then + group. Look at the group row’s p column. For question 2, recall that ANCOVA assumes the covariate’s slope is the same in every group.
library(rstatix)
library(datarium)
data("anxiety", package = "datarium")
anx <- as.data.frame(anxiety[, c("id", "group", "t1", "t3")])
names(anx)[names(anx) %in% c("t1", "t3")] <- c("pretest", "posttest")
anx %>% anova_test(posttest ~ pretest + group)
#> group is significant after adjusting for pretest: F(2, 41) ≈ 218, p < 0.0001.The group effect is significant after adjustment.
For question 2: you must check homogeneity of regression slopes — fit posttest ~ group * pretest and confirm the group:pretest interaction is not significant. If it is, the groups have different slopes and a standard ANCOVA is invalid.
You can now run a one-way and two-way ANCOVA in R: put the covariate first (outcome ~ covariate + group), check the ANCOVA-specific assumptions — linearity and homogeneity of regression slopes — alongside the usual normality, equal-variance and outlier checks, then compare the adjusted (estimated marginal) means with emmeans_test() and report them on the plot. ANCOVA gives a fairer group comparison than ANOVA whenever a continuous covariate shapes the outcome — and a two-way ANCOVA decomposes a significant interaction with the same simple-main-effects workflow as a two-way ANOVA.
@online{2026,
author = {},
title = {ANCOVA in {R:} {Compare} {Group} {Means} {Adjusted} for a
{Covariate}},
date = {2026-06-23},
url = {https://www.datanovia.com/learn/biostatistics/anova/ancova-in-r},
langid = {en}
}