ANCOVA in R: lm(), car::Anova() and emmeans

Q: How do I run an ANCOVA in R?

Fit lm(outcome ~ covariate + group, data = dat) , set options(contrasts = c("contr.sum", "contr.poly")) , refit, and run car::Anova(m, type = "III") . Then use emmeans(m, "group") for adjusted means.

Q: Which sums of squares should I use for ANCOVA in R?

Use Type II, or Type III with sum-to-zero contrasts, via car::Anova() . The sequential table from anova() is valid only when the covariate comes first in the formula.

ANCOVA in R is a linear model with one categorical factor and one continuous covariate: fit it with lm(), test the factor with car::Anova(), and report group means adjusted for the covariate with emmeans(). Three details decide whether the result is right: the order of terms in the formula, the contrast coding when you use Type III tests, and whether the covariate slope is the same in every group. This guide shows each one in code, including a check you can run yourself to see the order effect.

What ANCOVA in R is used for in research

Analysis of covariance compares group means on an outcome after accounting for a continuous covariate, typically a baseline measurement such as a pre-test score. In R there is no dedicated ANCOVA function: it is an ordinary linear model, lm(outcome ~ covariate + group), whose terms you then test and summarise. If you are new to the underlying test, start with ANOVA (analysis of variance); for the R mechanics of plain ANOVA see running ANOVA in R.

The question ANCOVA answers is not “do the raw group means differ?” but “do the groups differ at the same value of the covariate?” That second quantity is the adjusted mean, and the section on emmeans() below shows exactly how R computes it.

Install or access (official source only)

Base R supplies lm() and anova() through the stats package. Two add-on packages are used here, both from CRAN, the Comprehensive R Archive Network at cran.r-project.org:

install.packages(c("car", "emmeans"))
library(car)
library(emmeans)

The car package provides Anova(), and emmeans provides estimated marginal means. Nothing else is needed for the workflow below.

Step-by-step: an ANCOVA with a covariate, three groups and adjusted means

The examples use synthetic data generated in the script itself, so you can paste the whole block and run it. A pre-treatment score (pre) is the covariate, a post-treatment score (post) is the outcome, and group has three levels.

Create the data.

set.seed(42)
dat = data.frame(
  group = factor(rep(c("control", "drugA", "drugB"), each = 20)),
  pre   = rnorm(60, mean = 50, sd = 10)
)
dat$post = 10 + 0.8 * dat$pre + c(0, 5, 8)[dat$group] + rnorm(60, sd = 5)

Check that the covariate slope is the same in every group. Fit the additive model and the model with a group-by-covariate interaction, then compare them. The base R documentation for anova.lm states that comparing models only makes sense when they are nested, and these two are.
```
m_main = lm(post ~ pre + group, data = dat)
m_int  = lm(post ~ pre * group, data = dat)
anova(m_main, m_int)
```
The row for m_int tests the pre:group interaction. A small p-value means the slopes differ between groups; in that case an adjusted mean at one covariate value does not summarise the effect, and you should report the interaction instead. The full reasoning is in ANCOVA assumptions: homogeneity of regression slopes.
Fit the ANCOVA model with the covariate first. lm(formula, data, ...) takes a formula of the form response ~ terms. Put the covariate before the factor:
```
m = lm(post ~ pre + group, data = dat)
```
Test the terms with Type III tests and sum-to-zero contrasts. The car documentation says Type III tests are normally sensible only with contrasts that are orthogonal across terms, such as contr.sum, contr.poly or contr.helmert, and warns against R’s default contr.treatment. Set the contrasts, refit, and then call Anova():
```
options(contrasts = c("contr.sum", "contr.poly"))
m3 = lm(post ~ pre + group, data = dat)
Anova(m3, type = "III")
```
The group row is the ANCOVA test of the factor after adjusting for pre. The pre row tests the covariate itself.
Get the adjusted means with emmeans().
```
emm = emmeans(m3, "group")
emm
```
By default emmeans evaluates the model with each covariate at its average, so every group is compared at the same covariate value. To evaluate at specific covariate values instead, use the at argument:
```
emmeans(m3, "group", at = list(pre = c(40, 50, 60)))
```
Compare the groups pairwise.
```
pairs(emm)
pairs(emm, adjust = "bonferroni")
```
These are comparisons of adjusted means, not of raw means. For effect sizes to report alongside them, see effect size in research.

Why term order matters: the covariate-first rule

This is the point most tutorials pass over. According to the R documentation for anova.lm, a single-model table reports “the reductions in the residual sum of squares as each term of the formula is added in turn” The sums of squares are sequential, so the order of terms in the formula changes the numbers. You can see it directly by fitting the same model in both orders:

anova(lm(post ~ pre + group, data = dat))   # covariate first
anova(lm(post ~ group + pre, data = dat))   # group first

With the covariate first, the group row is the reduction in residual sum of squares after pre is already in the model, which is the ANCOVA test. With the group first, the group row is the reduction before the covariate is added, so it does not adjust for pre at all. Both calls run without a warning.

The decision rule that follows from the documentation:

With plain anova(lm(...)), always write the covariate before the factor.
With car::Anova(), order does not decide the test. The car documentation describes Type II tests as respecting marginality and Type III as testing each term after all others; Type II results do not depend on contrast coding, while Type III needs the orthogonal contrasts from step 4.
With drop1(), the anova.lm documentation describes a Type II-style table in which each term is dropped while respecting hierarchy.

Output and how to read it

The anova() table lists one row per term with the columns Df, Sum Sq, Mean Sq, F value and Pr(>F), plus a Residuals row. The Anova() table from car reports Sum Sq, Df, F value and Pr(>F) per term. Read it in this order:

The slopes comparison from step 2. Only if the interaction is not supported do adjusted means make sense.
The group row of the Type III table. This is the test of whether the adjusted means differ.
The emmeans table. It lists each group’s adjusted mean, its standard error and confidence limits. The covariate value used is printed with the output.
The pairs() table. It gives each difference between adjusted means with its adjusted p-value.

To see what adjustment did, put raw and adjusted means side by side. aggregate(post ~ group, data = dat, FUN = mean) returns the raw group means; summary(emm) returns the adjusted ones. When group means on the covariate differ, the two sets differ, and that gap is what the covariate removed.

Troubleshooting

The Type III table looks wrong or changes when I recode the factor

You are most likely using the default contr.treatment. The car documentation warns against Type III tests with it. Set options(contrasts = c("contr.sum", "contr.poly")) before fitting, as in step 4, or use Anova(m, type = "II"), which does not depend on the contrast coding.

The group effect changes when I swap the formula order

You are reading a sequential anova() table. Put the covariate first, or use car::Anova(). See the covariate-first rule above.

The interaction is significant

The slopes are not parallel. Do not report a single set of adjusted means as if the group difference were constant. Report the model with the interaction and use emmeans() with at to show group differences at chosen covariate values.

Adjusted means are at the wrong covariate value

The default is the covariate’s average. Set the value yourself with at = list(pre = ...).

Alternative routes in R

The same ANCOVA model can be tested through four different calls. All four are documented in the sources below.

Function	What it reports	Takes order or contrasts into account?
`anova(lm(...))` (stats)	Sequential reductions in residual sum of squares as each term is added in turn	Order matters; put the covariate first
`drop1(lm(...), test = "F")` (stats)	Each term dropped in turn while respecting hierarchy (Type II-style)	Order does not matter
`car::Anova(m, type = "II")`	Type II tests, which respect marginality	Invariant to contrast coding
`car::Anova(m, type = "III")`	Each term tested after all others	Needs orthogonal contrasts such as `contr.sum`

Comparing the same data across statistical packages? The SPSS route is covered in ANCOVA in SPSS.

FAQ

How do I run an ANCOVA in R?

Fit lm(outcome ~ covariate + group, data = dat), set options(contrasts = c("contr.sum", "contr.poly")), refit, and run car::Anova(m, type = "III"). Then use emmeans(m, "group") for adjusted means.

Which sums of squares should I use for ANCOVA in R?

Use Type II, or Type III with sum-to-zero contrasts, via car::Anova(). The sequential table from anova() is valid only when the covariate comes first in the formula.

What does emmeans hold the covariate at?

By default, at the average of the covariate. The emmeans basics vignette says the default reference grid uses the average of each covariate. Change it with at = list(covariate = c(...)).