chg grp age
1 17.05 aerobic 31
2 4.96 aerobic 23
3 10.40 aerobic 27
4 11.05 aerobic 28
5 0.26 aerobic 22
6 2.51 aerobic 24
7 -0.87 running 23
8 -10.74 running 22
9 -3.27 running 22
10 -1.97 running 25
11 7.50 running 27
12 -7.25 running 20STAT 516 Lec 11
Analysis of covariance (ANCOVA)
Karl Gregory
2026-10-26
VO\(_2\) max data from Kuehl (2000)
Change in VO\(_2\) max of 12 males (after-minus-before) randomly assigned to two exercise programs (running, step aerobics). Ages recorded.
Which exercise program will lead to a greater average change in VO\(_2\) max?
What role does age play?
VO\(_2\) max
VO\(_2\) max
VO\(_2\) max
VO\(_2\) max
Analysis of covariance
- Useful when EUs are not homogeneous.
- Measurement capturing EU inhomogeneity is called a covariate.
- Can isolate treatment effects even when EUs differ across treatment groups.
Single-slope analysis of covariance (ANCOVA) model
Assume \[ Y_{ij} = \mu + \tau_i + \beta x_{ij}+ \varepsilon_{ij} \] for \(i = 1,\dots,a\), \(j = 1,\dots,n_i\), where
- \(Y_{ij}\) is the response of EU \(j\) in treatment group \(i\).
- \(\mu\) is a baseline or overall mean.
- the \(\tau_i\) are treatment effects.
- the \(x_{ij}\) are covariate values measured on the EUs.
- \(\beta\) is a slope coefficient expressing the effect of the covariate.
- the \(\varepsilon_{ij}\) are independent \(\text{Normal}(0,\sigma_{\varepsilon}^2)\) error terms.
- \(n_1 + \dots + n_a = N\). Unbalancedness not an issue.
Set \(\mu_i = \mu + \tau_i\) for \(i =1,\dots,a\).
Goals in analysis of covariance
- Estimate the parameters \(\mu\), \(\tau_1,\dots,\tau_a\), and \(\beta\).
- Visualize the data.
- Fit full and reduced models and collect error sums of squares.
- Test whether the covariate has any effect.
- Test for a treatment effect.
- Adjust treatment group means for the inhomogeneity of the EUs.
- Compare the adjusted group means.
Parameter constraints
- There are \(a\) treatment groups and one covariate slope.
- The model has \(a + 1\) parameters, which is one too-many.
- R will set \(\tau_1 = 0\) so that it can estimate all the parameters.
VO\(_2\) max
Call:
lm(formula = chg ~ grp + age, data = vo2max)
Residuals:
Min 1Q Median 3Q Max
-5.77 -0.99 0.14 1.82 3.04
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -41.014 7.714 -5.32 0.00048 ***
grprunning -5.443 1.796 -3.03 0.01425 *
age 1.886 0.295 6.39 0.00013 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.8 on 9 degrees of freedom
Multiple R-squared: 0.902, Adjusted R-squared: 0.88
F-statistic: 41.4 on 2 and 9 DF, p-value: 2.89e-05VO\(_2\) max
VO\(_2\) max
Full and reduced models in ANCOVA
We construct test statistics by fitting full and reduced models:
- Full model: \(Y_{ij} = \mu + \tau_i + \beta x_{ij} + \varepsilon_{ij}\)
- Reduced model (no covariate effect): \(Y_{ij} = \mu + \tau_i + \varepsilon_{ij}\)
- Reduced model (no treatment effect): \(Y_{ij} = \mu + \beta x_{ij} + \varepsilon_{ij}\)
Full- and reduced-model error sums of squares
| SS | Formula |
|---|---|
| \(\operatorname{SS}_{\operatorname{Error}}(\text{Full})\) | \(\sum_{i=1}^a\sum_{j=1}^{n_{ij}}(Y_{ij} - (\hat \mu + \hat \tau_i + \hat \beta x_{ij}))^2\) |
| \(\operatorname{SS}_{\operatorname{Error}}(\text{No covariate})\) | \(\sum_{i=1}^a\sum_{j=1}^{n_{ij}}(Y_{ij} - (\hat \mu + \hat \tau_i) )^2\) |
| \(\operatorname{SS}_{\operatorname{Error}}(\text{No treatment})\) | \(\sum_{i=1}^a\sum_{j=1}^{n_{ij}}(Y_{ij} - (\hat \mu + \hat \beta x_{ij}))^2\) |
Now define these difference in error sums of squares:
- \(\operatorname{SS}_{\operatorname{Cov}}= \operatorname{SS}_{\operatorname{Error}}(\text{No covariate}) - \operatorname{SS}_{\operatorname{Error}}(\text{Full})\)
- \(\operatorname{SS}_{\operatorname{Trt}}= \operatorname{SS}_{\operatorname{Error}}(\text{No treatment}) - \operatorname{SS}_{\operatorname{Error}}(\text{Full})\)
ANCOVA table
| Source | Df | SS | MS | F value |
|---|---|---|---|---|
| Covariate | \(1\) | \(\operatorname{SS}_{\operatorname{Cov}}\) | \(\operatorname{MS}_{\operatorname{Cov}}\) | \(\operatorname{F}_{\operatorname{Cov}}= \operatorname{MS}_{\operatorname{Cov}}/ \operatorname{MS}_{\operatorname{Error}}\) |
| Treatment | \(a-1\) | \(\operatorname{SS}_{\operatorname{Trt}}\) | \(\operatorname{MS}_{\operatorname{Trt}}\) | \(\operatorname{F}_{\operatorname{Trt}}= \operatorname{MS}_{\operatorname{Trt}}/ \operatorname{MS}_{\operatorname{Error}}\) |
| Error | \(N - a - 1\) | \(\operatorname{SS}_{\operatorname{Error}}\) | \(\operatorname{MS}_{\operatorname{Error}}\) |
- Reject \(H_0\): \(\beta = 0\) at \(\alpha\) if \(\operatorname{F}_{\operatorname{Cov}}> F_{1,N-a-1,\alpha}\).
- Reject \(H_0\): \(\mu_1 = \dots = \mu_a\) at \(\alpha\) if \(\operatorname{F}_{\operatorname{Trt}}> F_{a-1,N-a-1,\alpha}\).
VO\(_2\) max
Using anova() on the lm() output gives the wrong SS (gives sequential).
Use Anova() from R package car on the lm() output.
Anova Table (Type III tests)
Response: chg
Sum Sq Df F value Pr(>F)
(Intercept) 221 1 28.27 0.00048 ***
grp 72 1 9.18 0.01425 *
age 319 1 40.78 0.00013 ***
Residuals 70 9
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1VO\(_2\) max
lm_nocov <- lm(chg ~ grp, data = vo2max) # reduced model with no covariate
lm_notrt <- lm(chg ~ age, data = vo2max) # reduced model with no treatment
lm_full <- lm(chg ~ age + grp, data = vo2max) # full model
SSE_nocov <- sum(lm_nocov$resid^2)
SSE_notrt <- sum(lm_notrt$resid^2)
SSE_full <- sum(lm_full$resid^2)
SSCov <- SSE_nocov - SSE_full
SSTrt <- SSE_notrt - SSE_full
SSE <- SSE_full
a <- 2
N <- nrow(vo2max)
MSCov <- SSCov / 1
MSTrt <- SSTrt / (a-1)
MSE <- SSE / (N - a - 1)
FCov <- MSCov / MSE
FTrt <- MSTrt / MSE
pCov <- 1 - pf(FCov,1,N - a - 1)
pTrt <- 1 - pf(FTrt,a-1,N - a - 1)Treatment means in terms of parameter estimates
We can write the treatment group means as \[ \bar Y_{i.} = \hat \mu + \hat \tau_i + \hat \beta \bar x_{i.} \quad \text{ for } i=1,\dots,a, \]
We also define the covariate-adjusted means \[ \bar Y_{i.}^{\operatorname{adj}} = \hat \mu + \hat \tau_i + \hat \beta \bar x_{..} \quad \text{ for } i=1,\dots,a, \]
We can equivalently write covariate-adjusted means as \[ \bar Y_{i.}^{\operatorname{adj}} = \bar Y_{i.} - \hat\beta(\bar x_{i.} - \bar x_{..}) \quad \text{ for } i = 1,\dots,a. \]
Variances of covariate-adjusted means and their differences
| Contrast | Variance |
|---|---|
| \(\bar Y_{i.}^{\operatorname{adj}}\) | \(\sigma^2\left[\dfrac{1}{n_i} + \dfrac{(\bar x_{i.} - \bar x_{..})^2}{E_{xx}}\right]\) |
| \(\bar Y_{i.}^{\operatorname{adj}} - \bar Y_{i'.}^{\operatorname{adj}}\) | \(\sigma^2\left[\dfrac{1}{n_i} + \dfrac{1}{n_{i'}} + \dfrac{(\bar x_{i.} - \bar x_{i'.})^2}{E_{xx}}\right]\) |
where \(E_{xx} = \sum_{i=1}^a\sum_{j=1}^{n_i} (x_{ij} - \bar x_{i.})^2\).
Some (unadjusted) CIs in the ANCOVA model
Define the true or population-level covariate-adjusted means as
\[ \mu_i^{\operatorname{adj}} = \mu + \tau_i + \beta \bar x_{..}\quad \text{ for } i = 1,\dots,a. \]
| Target | \((1-\alpha)100\%\) confidence interval |
|---|---|
| \(\mu_i^{\operatorname{adj}}\) | \(\bar Y_{i.}^{\operatorname{adj}} \pm t_{N - a -1,\alpha/2}\sqrt{\operatorname{MS}_{\operatorname{Error}}}\sqrt{\dfrac{1}{n_i} + \dfrac{(\bar x_{i.} - \bar x_{..})^2}{E_{xx}}}\) |
| \(\mu_i^{\operatorname{adj}} - \mu_{i'}^{\operatorname{adj}}\) | \(\bar Y_{i.}^{\operatorname{adj}} - \bar Y_{i'.}^{\operatorname{adj}} \pm t_{N - a -1,\alpha/2} \sqrt{\operatorname{MS}_{\operatorname{Error}}}\sqrt{\dfrac{1}{n_i} + \dfrac{1}{n_{i'}} + \dfrac{(\bar x_{i.} - \bar x_{i'.})^2}{E_{xx}}}\) |
VO\(_2\) max
y1. <- mean(vo2max$chg[vo2max$grp == "aerobic"])
y2. <- mean(vo2max$chg[vo2max$grp == "running"])
x1. <- mean(vo2max$age[vo2max$grp == "aerobic"])
x2. <- mean(vo2max$age[vo2max$grp == "running"])
x.. <- mean(vo2max$age)
bhat <- coef(lm_out)[3]
y1.adj <- y1. - bhat * (x1. - x..)
y2.adj <- y2. - bhat * (x2. - x..)VO\(_2\) max
n1 <- 6
n2 <- 6
alpha <- 0.05
tval <- qt(1-alpha/2,N-a-1)
age1 <- vo2max$age[vo2max$grp == "aerobic"]
age2 <- vo2max$age[vo2max$grp == "running"]
Ex1 <- sum((age1 - mean(age1))^2)
Ex2 <- sum((age2 - mean(age2))^2)
Exx <- Ex1 + Ex2
se1 <- sqrt(MSE) * sqrt(1/n1 + (x1. - x..)^2 / Exx)
lo1 <- y1.adj - tval * se1
up1 <- y1.adj + tval * se1
se2 <- sqrt(MSE) * sqrt(1/n2 + (x2. - x..)^2 / Exx)
lo2 <- y2.adj - tval * se2
up2 <- y2.adj + tval * se2
se12 <- sqrt(MSE) * sqrt(1/n1 + 1/n2 + (x1. - x2.)^2 / Exx)
lo12 <- y1.adj - y2.adj - tval * se12
up12 <- y1.adj - y2.adj + tval * se12VO\(_2\) max
| Unadjusted | Age-adjusted(CI) | |
|---|---|---|
| Aerobic | 7.71 | 5.19 (2.46,7.92) |
| Running | -2.77 | -0.25 (-2.98,2.48) |
| Aerobic - Running | 10.47 | 5.44 (1.38,9.51) |
Soybean data from Dr. Longnecker’s notes
Soybean plants assigned to three greenhouse conditions:
- Supplemental lighting (SL)
- Partial shading (PS)
- Control (C).
The response was seed yield. The pre-treatment height of each plant was also recorded.
Do the greenhouse conditions effect the seed yield?
What is the role of plant height (proxy for plant vigor)?
Soybean data
soybean <- data.frame(yield = c(12.2,12.1,12.3,16.6,15.6,15.0,9.5,9.8,8.6,
12.4,13.1,12.3,15.8,15.8,16.2,9.5,9.1,10.4,
11.9,12.7,12.6,16.5,15.8,16.7,9.6,10.3,9.4,
11.3,12.4,13.2,15.0,16.0,15.8,8.8,9.5,10.2,
11.8,11.4,12.3,15.4,15.8,15.9,9.5,8.5,9.3),
trt = as.factor(rep(c(rep("C",3),rep("SL",3),rep("PS",3)),5)),
height = c(45,48,48,63,45,35,52,62,41,
52,60,51,50,50,50,54,52,67,
42,61,56,63,48,62,58,67,55,
35,50,65,33,50,49,45,55,66,
40,33,51,38,49,52,57,40,56))
head(soybean,n = 12) yield trt height
1 12.2 C 45
2 12.1 C 48
3 12.3 C 48
4 16.6 SL 63
5 15.6 SL 45
6 15.0 SL 35
7 9.5 PS 52
8 9.8 PS 62
9 8.6 PS 41
10 12.4 C 52
11 13.1 C 60
12 12.3 C 51Soybean data
Soybean data
Soybean data
ANCOVA allowing different slopes in each group
Assume \[ Y_{ij} = \mu + \tau_i + (\beta + (\tau\beta)_i)x_{ij} + \varepsilon_{ij} \] for \(i = 1,\dots,a\), \(j = 1,\dots,n_i\), where
- \(Y_{ij}\) is the response of EU \(j\) in treatment group \(i\).
- \(\mu\) is a baseline or overall mean.
- the \(\tau_i\) are treatment effects.
- the \(x_{ij}\) are covariate values measured on the EUs.
- \(\beta\) is a slope coefficient expressing the effect of the covariate.
- the \((\tau\beta)_i\) allow interaction between the treatment and the covariate.
- the \(\varepsilon_{ij}\) are independent \(\text{Normal}(0,\sigma_{\varepsilon}^2)\) error terms.
- \(n_1 + \dots + n_a = N\). Unbalancedness not an issue.
Set \(\mu_i = \mu + \tau_i\) and \(\beta_i = \beta + (\tau\beta)_i\) for \(i =1,\dots,a\).
Parameter constraints in multiple slopes model
- There are \(a\) treatment groups and \(a\) covariate slopes.
- The model has \(2(a + 1)\) parameters, which is two too-many.
- R will set \(\tau_1 = 0\) and \((\tau\beta)_1 = 0\) to make all parameters estimable.
Soybean data
Call:
lm(formula = yield ~ trt + height + trt:height, data = soybean)
Residuals:
Min 1Q Median 3Q Max
-0.23473 -0.08875 -0.00395 0.05764 0.29332
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 9.50057 0.18193 52.22 < 2e-16 ***
trtPS -3.63959 0.28553 -12.75 1.7e-15 ***
trtSL 3.71305 0.25858 14.36 < 2e-16 ***
height 0.05630 0.00364 15.45 < 2e-16 ***
trtPS:height 0.00910 0.00537 1.69 0.098 .
trtSL:height -0.00244 0.00518 -0.47 0.641
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.126 on 39 degrees of freedom
Multiple R-squared: 0.998, Adjusted R-squared: 0.998
F-statistic: 4.05e+03 on 5 and 39 DF, p-value: <2e-16Soybean data
Soybean data
An F-test for equal slopes
Define the error sums of squares:
| SS | Formula |
|---|---|
| \(\operatorname{SS}_{\operatorname{Error}}(\text{Full})\) | \(\sum_{i=1}^a\sum_{j=1}^{n_{ij}}(Y_{ij} - (\hat \mu + \hat \tau_i + (\hat \beta + (\hat{\tau\beta})_i) x_{ij}))^2\) |
| \(\operatorname{SS}_{\operatorname{Error}}(\text{Equal slopes})\) | \(\sum_{i=1}^a\sum_{j=1}^{n_{ij}}(Y_{ij} - (\hat \mu + \hat \tau_i + \hat \beta x_{ij}))^2\) |
Now set
\[ \operatorname{F}_{\operatorname{T}\times\operatorname{C}}= \frac{[\operatorname{SS}_{\operatorname{Error}}(\text{Equal slopes}) - \operatorname{SS}_{\operatorname{Error}}(\text{Full})] / (a-1)}{\operatorname{SS}_{\operatorname{Error}}(\text{Full})/(N-2a)} \]
Reject \(H_0\): Equal slopes at \(\alpha\) if \(\operatorname{F}_{\operatorname{T}\times\operatorname{C}}> F_{a-1,N-2a,\alpha}\).
Soybean data
lm_eqslp <- lm(yield ~ trt + height, data = soybean) # equal slopes model
lm_full <- lm(yield ~ trt + height + trt:height, data = soybean) # full model
SSE_eqslp <- sum(lm_eqslp$resid^2)
SSE_full <- sum(lm_full$resid^2)
a <- 3
N <- nrow(soybean)
FTrtCov <- ((SSE_eqslp - SSE_full) / (a-1)) / (SSE_full / (N - 2*a))
pTrtCov <- 1 - pf(FTrtCov,a-1,N - 2*a)Soybean data
Same as interaction p value from anova() on the lm() output.
When slopes are unequal
If we reject \(H_0\): Equal slopes, then
- We can compute covariate-adjusted means with the different slopes.
- We can compare the covariate adjusted means—but the CI formulas are different from the ones above (details omitted).
References
Kuehl, R. O. 2000. Design of Experiments: Statistical Principles of Research Design and Analysis. Duxbury/Thomson Learning.