MANCOVA: Interpretation and APA Reporting

The Basics

When to Use It

Use it when you have two or more separate groups of different people, two or more related continuous outcomes, and a continuous covariate measured before the groups diverged that you want to adjust for. One multivariate test, rather than one per outcome, holds the overall false-positive rate at α (if its assumptions hold) and can find a difference that shows only in a combination of the outcomes. Keep covariates few and reliably measured: each costs an error degree of freedom, and measurement error in a covariate weakens the adjustment. The outcomes should be related once the covariate is out; if they barely are (in the exercise study, speed and recall correlate only .10 after age), one ANCOVA per outcome with a Bonferroni or Holm correction, chosen in advance, is an equally defensible plan. Tabachnick and Fidell (2019, Chapter 7) give a textbook treatment of MANCOVA alongside MANOVA. For example:

Do stretching, strength training and aerobic walking change older adults' processing speed and word recall, adjusting for age?
Do CBT, medication and a waitlist differ on depression and anxiety scores at 12 weeks, adjusting for baseline severity?
Do three reading interventions differ on children's reading fluency and comprehension, adjusting for pretest vocabulary?

Vocabulary You Need

Definition Covariate

A continuous variable that is related to the outcome but isn't the question of interest, such as a pretest score or age. Including it in the model sets aside the part of the outcome it predicts, which shrinks the unexplained variation and, when groups are compared, can adjust for small chance differences between them. In experiments it should be measured before any treatment, so the treatment can't have changed it. Glossary →

Definition Adjusted Mean (

\bar{y}_j'

)

A group's mean outcome, moved along the within-group regression line to where it would sit if the group's covariate mean equaled the overall covariate mean: the group's mean minus the slope times how far its covariate mean is from the overall one. With a positive slope, groups that started higher on the covariate are nudged down and groups that started lower are nudged up; with a negative slope, it's the other way around. Also called the estimated marginal mean or least-squares mean. Glossary →

Definition Sums of Squares and Cross-Products Matrix (

\mathbf{SSCP}

)

The multivariate version of a sum of squares. With two outcomes it is a 2 × 2 table: the diagonal holds each outcome's sum of squared deviations, and the off-diagonal holds the sum of the products of the two outcomes' deviations, which says whether they go up and down together. Multivariate tests split the total SSCP into a hypothesis part (H, what the groups explain) and an error part (E, what they don't), just as ANOVA splits the total sum of squares. Glossary →

Definition Pillai's Trace (

V

)

The sum, across the dimensions on which the groups can differ, of the share of variation each dimension explains, where a dimension is a weighted blend of the outcomes chosen to separate the groups (the best blend, then the next best uncorrelated with it, and so on). It runs from 0 (the group means are identical on every outcome) to s, the smaller of the number of outcomes and the effect's degrees of freedom (with groups, the number of groups minus 1; for one continuous predictor, 1), so it isn't a percentage of all the variation. Of the four usual MANOVA statistics it holds its false-positive rate best when the covariance matrices differ or the data aren't normal, which is why it is a common default. Also called the Pillai-Bartlett trace. Glossary →

Assumptions

1. Independent Observations The people (or other units) in your study don't influence each other.

What It Means: Each person's data comes from a separate, unconnected source, so knowing one person's scores tells you nothing about anyone else's. In designs that measure people more than once, their own scores are meant to be related; it's different people who must be independent. Methods built for clustered data, such as multilevel models and GEE, also allow members of the same group to be related, as long as the analysis names the grouping; then it's the groups that must be independent.
Applies To: Observations
How to Check: By design, not by a test. Did each person contribute their own data, with no clusters such as roommates, couples, therapy groups or classrooms whose members influence each other? If there are clusters, does the analysis model every one of them?
How Much It Matters: A lot. If participants exercised together in classes and influenced each other, scores are more alike than the test expects, and p comes out too small.

If It's Violated

Alternative Test Multilevel Regression (HLM). A multilevel model, which accounts for the clustering.
Design Change Assign and measure each person separately, or analyze class averages.

2. Multivariate Normality The variables are roughly bell-shaped together (within each group, when there are groups), or the sample is big enough that it matters less.

What It Means: The variables jointly follow a multivariate normal distribution: each is roughly normal, and every weighted combination of them is too, so the scatterplot of any two is an oval cloud. When the design has groups, as in MANOVA, it applies within each group; in a multivariate regression, to the outcomes' residuals; in factor analysis and SEM, it applies to the observed variables that maximum likelihood is fitted to (for outcomes, after allowing for any predictors the model takes as given). P-values and standard errors rely on it; estimates are usually less affected, and robust standard errors or a scaled test statistic can take its place.
Applies To: Residuals
How to Check: Look at a histogram of each variable (or each outcome's residuals, in a model with predictors) and a scatterplot of each pair, within each group if there are groups, for strong skew, curves or stray points. Multivariate outliers show up as large squared Mahalanobis distances from the means (each group's own, if there are groups). Don't rely on a normality test as a gatekeeper: with small samples it misses real problems, and with large ones it flags trivial ones.
How Much It Matters: Moderate. With 30 per group and equal group sizes, mild skew matters little, especially for Pillai's trace (Olson, 1976). Multivariate outliers matter more: one unusual person can pull both adjusted means and the error matrix. Our residuals look fine: the largest squared Mahalanobis distance, 9.27, is well under 13.82, the p < .001 cutoff from a χ² with 2 df (one per outcome). That uses the error SSCP matrix divided by its 86 degrees of freedom; R's cov() on the residuals divides by 89 and gives 9.60, also well under.

If It's Violated

Robust Estimation Report Pillai's trace, the most robust of the four multivariate statistics.
Transformation Log-transform right-skewed outcomes such as reaction times. Results are then on the log scale.
Proceed with Caveat Rerun without any multivariate outliers and report whether the conclusion changes.

3. Equal Covariance Matrices Your groups have about the same spread on each outcome, and the same correlations between them.

What It Means: Each group's population has the same covariance matrix: the same variance on every outcome and the same covariance (so the same correlation) between every pair. It's the multivariate version of equal variances, and it's what lets MANOVA pool the groups into one within-groups SSCP matrix.
Applies To: Covariance structure
How to Check: Compare the groups' standard deviations and within-group correlations side by side, or overlay each group's data ellipse on a scatterplot. Box's M test (Box, 1949) is the usual formal check, but it is very sensitive to non-normality and outliers (Friendly & Sigal, 2020) and flags trivial differences in large samples, so don't use it as a gatekeeper; some authors use a strict cutoff such as p < .001.
How Much It Matters: Little with equal group sizes, like our 30 per group, especially for Pillai's trace (Olson, 1976). With unequal groups it matters more: if the bigger group has the bigger spread the test is too conservative, and if the smaller one does, too liberal. Here each group's SDs are similar and Box's M is far from significant, p = .853 (Box, 1949).

If It's Violated

Robust Estimation Report Pillai's trace, which holds its false-positive rate best when covariance matrices differ.
Design Change Keep group sizes equal; unequal spreads then barely matter.

4. Linear Relationship The relationship between your variables follows a straight line, not a curve.

What It Means: As one variable goes up, the other tends to change by about the same amount across its whole range. A curved relationship, one that levels off or a U shape, breaks this.
Applies To: Relationship
How to Check: Look at a scatterplot of each pair of variables the analysis relates: a roughly oval cloud tilted one way is fine; a clear bend, a ceiling or a U is not. In a regression with several predictors, plot the residuals against the predicted values instead: a flat, even band is fine; a curve means the model is missing a bend. A component-plus-residual plot for each predictor (R: car::crPlots) shows whether one of them bends once the others are held constant.
How Much It Matters: A lot for the adjustment. The covariate's effect is removed as a straight line, so a curved relationship with age is only partly removed, and the adjusted means are off. Pairs of outcomes that relate in a curve also weaken the multivariate test.

If It's Violated

Correction Add age² as a second covariate if the relationship bends.
Transformation Transform the covariate (for example, log) so its relationship with the outcomes is straight.

5. Equal Slopes The covariate relates to the outcome about equally strongly in every group.

What It Means: The line linking the covariate to the outcome has about the same slope in each group, so the groups' lines are parallel. The adjustment uses one shared slope, and a single adjusted difference between groups only makes sense if the gap between the lines is the same all along the covariate.
Applies To: Covariate
How to Check: Plot the outcome against the covariate with a separate line for each group and look for lines that cross or fan out, and compare the groups' own slopes. Fitting the model once with a group × covariate interaction is a supporting check, not a gate: in a small study a non-significant interaction isn't evidence that the slopes are equal, and in a large one a significant interaction can flag a trivial difference. Whether to keep one shared slope is a judgment from theory, the plots and your prespecified plan, not from the interaction test alone.
How Much It Matters: A lot. If age mattered more in one program than another, the program difference would depend on age, and one set of adjusted means would misrepresent it. Here the program × age interaction is small and not significant, Pillai's V = .03, F(4, 168) = 0.70, p = .590.

If It's Violated

Alternative Test Multivariate Regression. Keep the group × covariate interaction in the model and report group differences at several covariate values (the Johnson-Neyman approach).
Proceed with Caveat If the difference in slopes is small and not meaningful, drop the interaction and say so.

Limitations

This Is Not Your Test If…

You have only one outcome ANCOVA →
You have several outcomes but no covariate MANOVA →
Your outcomes are the same measure taken repeatedly Mixed Factorial ANOVA →
Your covariate was measured after the treatment began MANOVA →
Your outcomes aren't continuous (yes/no, categories) Categorical Data and Proportions →

What It Doesn't Tell You

Which programs, on which outcome. The omnibus test only says the three programs' adjusted averages aren't all the same on speed and recall together. Which programs differ, and on which test, comes from the follow-ups.
Why. Random assignment lets us say the programs made a difference, but not why aerobic walking helped: fitness itself, more time moving, or more social contact on the walks.
Whether it matters. Significance doesn't make about 5 more symbols in 90 seconds, or one more remembered word, important for daily life.
Who else. The result describes volunteers aged 65 to 82 who stuck with a 12-week program, and only these two tests of thinking. The age adjustment sharpens the comparison; it doesn't widen it to other ages.
That there's no effect on recall. Recall missed the Bonferroni cutoff we chose in advance (p = .033 against .025; Holm, had we chosen it, would have let it through), but the aerobic group still averaged about one more word than the stretching group. That's inconclusive, not evidence that exercise leaves memory alone.

Common Mistakes

Don'tFollow a significant MANCOVA only with separate uncorrected ANCOVAs.

InsteadCorrect the follow-up tests (Bonferroni or Holm), or describe the multivariate effect with discriminant analysis.
Don'tUse a covariate measured after the treatment started.

InsteadMeasure covariates before assignment, or use ones treatment can't change, like age.
Don'tUse a covariate to “equalize” intact groups that differ on it.

InsteadRandomize if you can. With intact groups, the adjustment can't say what groups would look like if equal (Miller & Chapman, 2001).
Don'tAdd or drop covariates after seeing which choice gives a significant result.

InsteadChoose covariates before collecting data, from theory or earlier studies, and report that choice.
Don'tThrow in many loosely related outcomes to make the test multivariate.

InsteadUse outcomes that belong together conceptually; each extra outcome costs degrees of freedom.
Don'tAdd up the partial η² values as if they split the variation between them.

InsteadRemember each has its own denominator, so they can sum past 1. Report η² if you want shares of the total.
Don'tReport partial η² under the label η².

InsteadSay which one you report. SPSS prints partial η², which is never smaller than η² and is usually larger, especially with other effects or repeated measures.
Don'tRead a p-value as the chance that the null hypothesis is true.

InsteadRead it as: if the null were true, how often data at least this extreme would turn up.

Result Decoder

Significance and Effect Size Guidelines

Benchmarks Significance (p Value)

p < .001
Very strong evidence against the null hypothesis.
.001 ≤ p < .01
Strong evidence against the null hypothesis.
.01 ≤ p < .05
Statistically significant at α = .05. Moderate evidence against the null hypothesis.
.05 ≤ p < .10
Not significant. Weak evidence; sometimes called a “trend.” Treat it as inconclusive.
p ≥ .10
Not significant. Little evidence against the null, which is not the same as evidence of no effect.

These evidence labels are one common convention for describing p in graded language rather than as significant or not; the boundaries are rough guides, not rules, and other authors draw them differently. Muff et al. (2022) argue for this kind of language of evidence over a single cutoff. α is chosen before the analysis, and p never says how big or important an effect is.

More about statistical significance ↓

Benchmarks Effect Size (Multivariate Partial

\eta^2

)

η²ₚ < .01 Below Small
Smaller than the conventional small benchmark. Whether an effect this size matters depends on the setting.
.01 ≤ η²ₚ < .06 Small
A sliver of the remaining variation lines up with the effect. Hard to see without measuring it.
.06 ≤ η²ₚ < .14 Medium
A noticeable share of the remaining variation lines up with the effect.
η²ₚ ≥ .14 Large
A large share of the remaining variation lines up with the effect.

Cohen's (1988) benchmarks were written for a single outcome. For a multivariate effect they are an even rougher guide; look at the follow-up effect sizes for each outcome too.

More about effect size ↓

How It's Reported

Equation Form

Example

Pillai's V = .15, F(4, 172) = 3.44, p = .010, η²_p = .07, 95% CI [.01, .14]

Table 1. Descriptive Statistics and Age-Adjusted Means by Exercise Program

Example

	$n$	Age $M$	Age $SD$	Speed $M$	Speed $SD$	Speed $M_{\text{adj}}$	Recall $M$	Recall $SD$	Recall $M_{\text{adj}}$
Stretching	30	70.67	4.17	50.00	8.12	49.73	8.00	2.03	7.91
Strength Training	30	71.03	4.37	52.00	7.73	52.06	8.40	1.98	8.42
Aerobic Walking	30	71.20	4.92	55.00	8.37	55.21	9.00	2.44	9.07

Note. Speed is symbols coded in 90 seconds; recall is words recalled out of 15. Adjusted means are evaluated at the overall mean age, 70.97 years.

Table 2. Multivariate Tests for Speed and Recall

Example

	Value	$F$	$df_1$	$df_2$	$p$	$\eta^2_p$
Program: Pillai's trace	.148	3.44	4	172	.010	.07
Program: Wilks' Λ	.852	3.55	4	170	.008	.08
Age: Pillai's trace	.478	38.99	2	85	< .001	.48

Note. η²_p = multivariate partial eta squared. The program rows are adjusted for age. Pillai's trace is the statistic reported; Wilks' Λ agrees. The age row shows the covariate's own effect, adjusted for program.

Table 3. Follow-Up ANCOVAs for Program, Adjusted for Age

Example

	$F$	$df$	$p$	$\eta^2_p$
Processing speed	4.65	2, 86	.012	.10
Word recall	3.55	2, 86	.033	.08

Note. Each outcome is tested at a Bonferroni-corrected α of .05 ÷ 2 = .025.

Full Text (APA Style)

Example

Adjusted for age, participants who did aerobic walking scored highest on both tests (processing speed M_adj = 55.21 symbols; recall M_adj = 9.07 words), followed by strength training (52.06; 8.42) and stretching (49.73; 7.91). A one-way MANCOVA with age as the covariate showed the difference was statistically significant, Pillai's V = .15, F(4, 172) = 3.44, p = .010 (Wilks' Λ = .85, p = .008). The multivariate effect, averaged over the two dimensions on which the programs could differ, was medium, η²_p = .07, 95% CI [.01, .14] (noncentral F interval). Follow-up ANCOVAs at a Bonferroni-corrected α of .025, chosen before the analysis, showed that the programs differed in processing speed, F(2, 86) = 4.65, p = .012, η²_p = .10, but not reliably in recall, F(2, 86) = 3.55, p = .033, η²_p = .08. On speed, aerobic walkers coded 5.49 more symbols than the stretching group, Bonferroni-adjusted 95% CI [1.08, 9.90], p = .009.

Name the covariate, and give the raw M and SD for every group and outcome in a table, with the adjusted means beside them.
Say which multivariate statistic you used. Pillai's trace is the safest default; Wilks' Λ is the most familiar.
Report both degrees of freedom of the approximate F, and say how you followed up (which tests, which correction). Choose the correction before seeing the results and use it throughout.
Give intervals: a CI for η²_p (here a two-sided interval from the noncentral F; in R, effectsize::F_to_eta2(3.44, 4, 172, alternative = "two.sided"), since its default is one-sided and fixes the upper limit at 1) and for the key adjusted difference.
Italicize Latin symbols (M, F, p, V) but not Greek ones (η², Λ), and report exact p to three decimals.

Reading Each Piece

The equation form above, piece by piece, in this order: how big, then how surprising.

η²_p = .07	The multivariate effect size, averaged over the two dimensions after age. The one combination of speed and recall that separates the programs lines up with about 15% of its variation, the other with almost none. That's a medium effect by the benchmarks above; whether it matters is a question about aging and cognition, not statistics. Benchmarks ↑
Pillai's V = .15	Program explains a combined .15 of the age-adjusted variation, summed over the two dimensions the groups could differ along (so it could run up to 2). Nearly all of it lies on one combination of speed and recall.
Wilks' Λ = .85	Program leaves 85% of the combined, age-adjusted variation unexplained. Its F version, F(4, 170) = 3.55, p = .008, agrees with Pillai's.
F(4, 172) = 3.44	Pillai's trace, rescaled to an F. The 4 and 172 are its degrees of freedom: 4 = 2 outcomes × 2 group contrasts, and 172 = 2 dimensions × 86 leftover degrees of freedom (90 people − 3 groups − 1 covariate), a shortcut that works because there are two outcomes.
p = .010	If the three programs worked equally well, adjusted averages at least this far apart would turn up about 1% of the time. That's below .05, so the result is statistically significant. Benchmarks ↑

Calculations

The Running Example

Example

The exercise study: 90 older adults, 30 per program (stretching, strength training, aerobic walking), each with two cognitive scores (processing speed and word recall) and their age as the covariate.

Null Hypothesis:: $\boldsymbol{\mu}^{\text{adj}}_{\text{stretch}} = \boldsymbol{\mu}^{\text{adj}}_{\text{strength}} = \boldsymbol{\mu}^{\text{adj}}_{\text{aerobic}}$
Alternative:: not all three adjusted mean pairs are equal

Each bold $\boldsymbol{\mu}^{\text{adj}}$ is a pair of population means, speed and recall, for people in that program, adjusted to the same age. The null says all three pairs match; the alternative doesn't say which programs differ, or on which outcome.

What Your Data Should Look Like

One row per person, with one column naming their group, one numeric column for each outcome, and one numeric column for each covariate. Each person appears once, in exactly one group, with a score on every column.

Example

participant	program	age	speed	recall
1	stretching	66	40	10
2	aerobic	80	50	5
3	stretching	69	55	9
4	aerobic	68	60	10
5	strength	70	53	8
6	strength	65	65	10
… 84 more rows

	n	Age M (SD)	Speed M (SD)	Speed Adj. M	Recall M (SD)	Recall Adj. M
Stretching	30	70.67 (4.17)	50.00 (8.12)	49.73	8.00 (2.03)	7.91
Strength Training	30	71.03 (4.37)	52.00 (7.73)	52.06	8.40 (1.98)	8.42
Aerobic Walking	30	71.20 (4.92)	55.00 (8.37)	55.21	9.00 (2.44)	9.07
All	90	70.97 (4.45)	52.33 (8.25)		8.47 (2.17)

Adjusted means are each group's mean at the overall age mean, 70.97, using the pooled within-group slopes (b = −0.91 for speed, b = −0.30 for recall).

Within-group correlation between speed and recall: r = .39; after age is taken out, r = .10, the one the test uses.

Simulated data.

Download all 90 rows (data.csv)

Napkin Math

The whole test, one idea at a time. Under each heading, a line says why the step matters, and the green note reads the equation in plain English with the exercise study's numbers. Collapse all to read just the headings: they tell the whole story on their own.

1. Put Every Group on the Same Age before Comparing

Both scores drop with age, and the programs' average ages differ a little by chance. Adjusting slides each group's averages to what they'd be at the overall average age, 70.97. Because scores drop with age, a group older than average is moved up, and a younger one down.

$\bar{y}_j' = \bar{y}_j - b_w(\bar{x}_j - \bar{x})$ = $55.00 - (-0.91)(71.20 - 70.97)$

= 55.21 symbols

$\text{Strength, speed}$ = 52.06
$\text{Stretching, speed}$ = 49.73
$\text{Aerobic, recall}$ = 9.07
$\text{Strength, recall}$ = 8.42
$\text{Stretching, recall}$ = 7.91

$\bar{y}_j$ is the aerobic group's raw speed average, 55.00 symbols, and $b_w$ = −0.91 is how much speed drops per year of age within a program. The group is 0.23 years older than average, so it gets back 0.91 × 0.23 = 0.21 symbols, for 55.21.
2. Measure the Noise Left Over after Age and Program

Fit both outcomes with program and age, and keep what's left over for each person. Squaring and adding those leftovers gives each outcome's noise; multiplying the two leftovers gives how the noise on one outcome moves with the other.

$\mathbf{E} = \begin{bmatrix} \sum e_{1i}^2 & \sum e_{1i}e_{2i} \\ \sum e_{1i}e_{2i} & \sum e_{2i}^2 \end{bmatrix}$ = $\begin{bmatrix} 4199.58 & 98.92 \\ 98.92 & 245.40 \end{bmatrix}$

$E_{11}$ = 4199.58
$E_{22}$ = 245.40
$E_{12}$ = 98.92
$r_{\text{after age}}$ = .10

$e_{1i}$ and $e_{2i}$ are one person's leftover speed and recall. Squared and added over all 90 people they give 4199.58 and 245.40; their products add to 98.92. As a correlation, 98.92 ÷ √(4199.58 × 245.40) = .10. The raw scores' within-program correlation in the data box is much higher, because age pushes both scores down together; once age is out, the two outcomes barely move together.
3. Measure How Much Program Explains beyond Age

Refit with age alone and the leftovers grow, because program differences are now lumped in with the noise. The growth is what program explains, on each outcome and in how the two move together.

$\mathbf{H} = \mathbf{E}_{\text{without groups}} - \mathbf{E}$ = $\begin{bmatrix} 4653.60 & 194.82 \\ 194.82 & 265.65 \end{bmatrix} - \begin{bmatrix} 4199.58 & 98.92 \\ 98.92 & 245.40 \end{bmatrix}$

$H_{11}$ = 454.02
$H_{22}$ = 20.26
$H_{12}$ = 95.89

$\mathbf{E}_{\text{without groups}}$ is the leftover matrix with age only. Subtracting $\mathbf{E}$ gives $\mathbf{H}$ = 454.02 for speed, 20.26 for recall and 95.89 for the two together.
4. Compare What Program Explains with the Total, on Both Outcomes at Once

This is ANOVA's “explained over total,” done on both outcomes at once. The test looks for the weighted sum of the outcomes on which the programs differ most (here about 0.3 × speed + 1 × recall), then the next-best sum unrelated to the first. Each such sum is a dimension; there are as many as the smaller of the number of outcomes and the number of groups minus 1, here 2. Pillai's trace adds up the share of variation the programs explain along each. H and E are the age-adjusted matrices from the two steps above, not MANOVA's plain between- and within-groups ones.

$V = \operatorname{tr}\!\left[\mathbf{H}(\mathbf{H} + \mathbf{E})^{-1}\right]$ = $0.085 + 0.063$

= 0.148

$\mathbf{H}(\mathbf{H} + \mathbf{E})^{-1}$ is explained over total in matrix form; multiplying by the inverse, $(\mathbf{H} + \mathbf{E})^{-1}$ , is how matrices divide. Its diagonal holds 0.085 and 0.063, and the trace (tr) adds them: 0.085 + 0.063 = 0.148. The two dimensions split the same total differently, as 0.148 and almost 0: the programs differ along one combination of speed and recall and barely at all along the other.
5. Check the Same Question with Wilks' Lambda

Wilks' Λ asks it the other way round, how much is left unexplained, using determinants, which boil each matrix down to one number for its total spread.

$\Lambda = \dfrac{|\mathbf{E}|}{|\mathbf{H} + \mathbf{E}|}$ = $\dfrac{1{,}020{,}774}{1{,}198{,}295}$

= .85

$|\mathbf{E}|$ is the leftover spread with program in the model, and $|\mathbf{H} + \mathbf{E}|$ is the spread with age alone. 1,020,774 ÷ 1,198,295 = .85, so program leaves 85% unexplained.
6. Turn Pillai's Trace into an F

Pillai's trace has no ready-made table, but a simple rescaling of it follows an F distribution closely enough to give a p-value.

$F = \dfrac{V}{s - V} \times \dfrac{df_2}{df_1}$ = $\dfrac{0.148}{2 - 0.148} \times \dfrac{172}{4}$

= 3.44

$s$ = 2 dimensions (step 4). $df_1$ = 2 outcomes × (3 groups − 1) = 4. $df_2 = s(df_{\text{error}} - p + s) = 2(86 - 2 + 2) = 172$ , where $df_{\text{error}}$ = 86 = 90 people − 3 groups − 1 covariate and $p$ = 2 outcomes. With two outcomes that's simply 2 × 86; with three outcomes it would be 2 × 85 = 170. 0.080 × 43 = 3.44.
7. Measure How Big the Adjusted Program Effect Is

F grows with sample size; this doesn't. It averages Pillai's trace over the two dimensions, once age is set aside, so a dimension with nothing on it halves the number.

$\eta^2_p = \dfrac{V}{s}$ = $\dfrac{0.148}{2}$

= .07

$95\%\ \text{CI}$ = [.01, .14]

$V$ is our 0.148 and $s$ is 2. 0.148 ÷ 2 = .07. One dimension separates the programs (0.148) and the other doesn't (almost 0), so .07 is their average, not a share of the variation in either score.
8. Ask How Surprising This Would Be If the Programs Worked Equally Well

How often would chance alone spread three groups' adjusted averages at least this far apart, on the two outcomes together, if the programs made no difference?

$p = P(F_{df_1,\,df_2} \ge F)$ = $P(F_{4,\,172} \ge 3.44)$

= .010

$F_{4,\,172}$ is an F distribution with 4 and 172 df. The chance it lands at or above 3.44 is .010, about 1 in 100.
9. Find Which Outcome the Programs Moved

A significant MANCOVA says the programs differ on the outcomes together, not on which one. The diagonal entries of H and E are exactly each outcome's ANCOVA, tested here at .05 ÷ 2 = .025.

$F = \dfrac{H_{jj} / (k - 1)}{E_{jj} / (N - k - c)}$ = $\dfrac{454.02 / 2}{4199.58 / 86}, \quad \dfrac{20.26 / 2}{245.40 / 86}$

$F_{\text{speed}}$ = 4.65
$p_{\text{speed}}$ = .012
$F_{\text{recall}}$ = 3.55
$p_{\text{recall}}$ = .033
$\text{Speed, aerobic} - \text{stretching}$ = 5.49 symbols
$95\%\ \text{CI (Bonferroni)}$ = [1.08, 9.90]
$p_{\text{Bonferroni}}$ = .009

$k$ = 3 groups, $N$ = 90 people and $c$ = 1 covariate, so the error df is 90 − 3 − 1 = 86. For speed, $H_{11}$ = 454.02 and $E_{11}$ = 4199.58, so 227.01 ÷ 48.83 = 4.65, p = .012, under .025. For recall, 10.13 ÷ 2.85 = 3.55, p = .033, which misses .025. Comparing adjusted speed means pair by pair (Bonferroni over the three pairs), aerobic walking beat stretching by 5.49 symbols.

Run It in Software

The code below reads the running example's data from data.csv; save it in the same folder you run the code from.

SPSS (version 30, as documented by IBM)

Recode program into a numeric variable: Transform → Recode into Different Variables, into prog with 1 = stretching, 2 = strength, 3 = aerobic, and give it those value labels.
Analyze → General Linear Model → Multivariate.
Move speed and recall into Dependent Variables, prog into Fixed Factor(s) and age into Covariate(s).
EM Means: move prog into Display Means for, check Compare main effects and choose Bonferroni.
Options: check Descriptive statistics, Estimates of effect size and Homogeneity tests, and click OK.

Or paste this syntax:

RECODE program ('stretching'=1) ('strength'=2) ('aerobic'=3) INTO prog.
VALUE LABELS prog 1 'stretching' 2 'strength' 3 'aerobic'.
EXECUTE.

GLM speed recall BY prog WITH age
  /METHOD=SSTYPE(3)
  /EMMEANS=TABLES(prog) WITH(age=MEAN) COMPARE ADJ(BONFERRONI)
  /PRINT=DESCRIPTIVE ETASQ HOMOGENEITY
  /CRITERIA=ALPHA(.05)
  /DESIGN=age prog.

Reading the Output

Read the tables in this order: the adjusted means in Estimated Marginal Means, then Partial Eta Squared and Pillai's Trace in the prog rows of Multivariate Tests, then Sig., then Tests of Between-Subjects Effects for each outcome.

Multivariate Tests (prog and age rows)

Effect		Value	F	Hypothesis df	Error df	Sig.	Partial Eta Squared
age	Pillai's Trace	.478	38.990	2	85	<.001	.478
prog	Pillai's Trace	.148	3.440	4	172	.010	.074
	Wilks' Lambda	.852	3.547	4	170	.008	.077
	Hotelling's Trace	.174	3.652	4	168	.007	.080
	Roy's Largest Root	.174	7.477	2	86	.001	.148

Roy's Largest Root gives only an upper bound on F, so its p is too small; don't report it alone. Box's Test of Equality of Covariance Matrices above it gives Box's M = 2.723, F(6, 188643) = 0.439, Sig. = .853.

Tests of Between-Subjects Effects (abridged)

Source	Dependent Variable	Type III Sum of Squares	df	Mean Square	F	Sig.	Partial Eta Squared
age	speed	1472.422	1	1472.422	30.153	<.001	.260
	recall	159.804	1	159.804	56.004	<.001	.394
prog	speed	454.019	2	227.009	4.649	.012	.098
	recall	20.258	2	10.129	3.550	.033	.076
Error	speed	4199.578	86	48.832
	recall	245.396	86	2.853

These are the univariate ANCOVAs. SPSS doesn't correct them, so compare each Sig. with .05 ÷ 2 = .025 yourself. The Estimated Marginal Means table gives the adjusted means (speed 49.726, 52.061, 55.213; recall 7.910, 8.420, 9.070) evaluated at age = 70.97, and Pairwise Comparisons gives aerobic − stretching on speed = 5.488, Bonferroni Sig. = .009.

Output laid out from this page's computed results, not captured from a live run; labels and layout may differ slightly in your version.

R

exercise <- read.csv("data.csv")
exercise$program <- factor(exercise$program, levels = c("stretching", "strength", "aerobic"))

# Put the covariate first: R's tests are sequential, so program is then adjusted for age
fit <- manova(cbind(speed, recall) ~ age + program, data = exercise)
summary(fit, test = "Pillai")
summary(fit, test = "Wilks")

# Follow-up ANCOVAs, one per outcome
summary.aov(fit)

          Df  Pillai approx F num Df den Df   Pr(>F)    
age        1 0.46884   37.513      2     85  2.1e-12 ***
program    2 0.14815    3.440      4    172 0.009842 ** 
Residuals 86                                            
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
          Df   Wilks approx F num Df den Df   Pr(>F)    
age        1 0.53116   37.513      2     85  2.1e-12 ***
program    2 0.85186    3.547      4    170 0.008287 ** 
Residuals 86                                            
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 Response speed :
            Df Sum Sq Mean Sq F value    Pr(>F)    
age          1 1398.4 1398.40 28.6368 7.135e-07 ***
program      2  454.0  227.01  4.6488    0.0121 *  
Residuals   86 4199.6   48.83                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 Response recall :
            Df  Sum Sq Mean Sq F value    Pr(>F)    
age          1 154.746 154.746 54.2312 1.003e-10 ***
program      2  20.258  10.129  3.5498   0.03301 *  
Residuals   86 245.396   2.853                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Read the program rows: they match SPSS because program is entered after age. The age rows differ from SPSS's because here age isn't adjusted for program; that doesn't matter for the question. For adjusted means and Bonferroni pairwise comparisons, use emmeans::emmeans(lm(speed ~ age + program, data = exercise), pairwise ~ program, adjust = "bonferroni").

Output from a real run: R 4.6.1 (stats), run 2026-09-29.

Python

import pandas as pd
from statsmodels.multivariate.manova import MANOVA

exercise = pd.read_csv("data.csv")
exercise["program"] = pd.Categorical(exercise["program"], categories=["stretching", "strength", "aerobic"])

fit = MANOVA.from_formula("speed + recall ~ C(program) + age", data=exercise)
print(fit.mv_test())

                  Multivariate linear model
=============================================================

-------------------------------------------------------------
        Intercept        Value  Num DF  Den DF F Value Pr > F
-------------------------------------------------------------
           Wilks' lambda 0.3222 2.0000 85.0000 89.4096 0.0000
          Pillai's trace 0.6778 2.0000 85.0000 89.4096 0.0000
  Hotelling-Lawley trace 2.1038 2.0000 85.0000 89.4096 0.0000
     Roy's greatest root 2.1038 2.0000 85.0000 89.4096 0.0000
-------------------------------------------------------------

-------------------------------------------------------------
       C(program)       Value  Num DF  Den DF  F Value Pr > F
-------------------------------------------------------------
          Wilks' lambda 0.8519 4.0000 170.0000  3.5475 0.0083
         Pillai's trace 0.1481 4.0000 172.0000  3.4400 0.0098
 Hotelling-Lawley trace 0.1739 4.0000 100.9691  3.6815 0.0077
    Roy's greatest root 0.1739 2.0000  86.0000  7.4773 0.0010
-------------------------------------------------------------

-------------------------------------------------------------
           age           Value  Num DF  Den DF F Value Pr > F
-------------------------------------------------------------
           Wilks' lambda 0.5215 2.0000 85.0000 38.9898 0.0000
          Pillai's trace 0.4785 2.0000 85.0000 38.9898 0.0000
  Hotelling-Lawley trace 0.9174 2.0000 85.0000 38.9898 0.0000
     Roy's greatest root 0.9174 2.0000 85.0000 38.9898 0.0000
=============================================================

Skip the Intercept block, which only tests whether the scores are zero. Read the C(program) block: Pillai's trace and Wilks' Λ match SPSS; statsmodels uses a different approximation for the Hotelling-Lawley F, so its denominator df differs. It doesn't print η²_p; divide Pillai's trace by s = 2 (0.074). For the follow-up ANCOVAs, fit smf.ols("speed ~ C(program) + age", data=exercise).fit() (with import statsmodels.formula.api as smf) and pass it to statsmodels.stats.anova.anova_lm(fit, typ=2).

Output from a real run: Python 3.9.6 (statsmodels 0.14.6, pandas 2.3.3), run 2026-09-29.

Impact

Effect Size

Multivariate partial eta squared (η²_p) averages, over the s dimensions the groups can differ along, the share of variation the groups explain on each, after the covariate's share is set aside. With Pillai's trace it's just the trace divided by s. A dimension with nothing on it still counts in the average, so read η²_p as an average, not as a share of the variation in any one outcome.

\eta^2_p = \dfrac{V}{s} = \dfrac{0.148}{2} = .07

Once age is accounted for, the programs separate along one combination of speed and recall, which lines up with about 15% of its variation, and barely at all along the second. Averaged over the two, η²_p = .07, a medium effect by convention. The 95% CI, .01 to .14, is wide: anywhere from below the small benchmark to large.

It's "partial" because age's share is removed from the denominator, so it runs larger than a plain η² would, and it isn't comparable with an η² from a study without the covariate. Each statistic also has its own version: SPSS prints .077 next to Wilks' Λ (1 − Λ^1/s), so say which one you report.

A single multivariate number is hard to picture, so also give each outcome's effect in its own units. For speed, the follow-up ANCOVA gives η²_p = .10, and aerobic walkers coded 5.49 more symbols in 90 seconds than the stretching group, at the same age. For recall, η²_p = .08, about one extra word out of 15. Benchmarks are conventions. Whether five symbols or one word is worth 12 weeks of walking is a question about aging, not statistics.

Benchmarks Multivariate Partial

\eta^2

at a Glance

Below Small

η²ₚ < .01

Small

.01 ≤ η²ₚ < .06

Medium

.06 ≤ η²ₚ < .14

Large

η²ₚ ≥ .14

Full benchmarks and caveats ↑

Statistical Significance

p answers one narrow question: if every group had the same adjusted population means on every outcome, how often would the samples' adjusted means spread at least this far apart, on the outcomes taken together? It comes from the upper tail of the approximate F, so it responds to differences in any direction and on any outcome. Below the conventional α = .05, the result is called statistically significant.

If the three programs worked equally well, age-adjusted averages at least as far apart as ours would turn up about 1.0% of the time, roughly 1 in 100. That's below .05, so it's significant. Without the covariate, the same data give p = .171: age differences between people were drowning out the program differences.

Besides Pillai's trace and Wilks' Λ, software prints Hotelling's trace and Roy's largest root; Roy's uses only the best dimension. The four can disagree when the result is near the cutoff; here all four are .010 or smaller. Pick one in advance (Pillai's trace is the safest) rather than reporting whichever is smallest.

p is not the probability that the null hypothesis is true, and it says nothing about how large the differences are, which programs differ, or on which outcome. That's what η²_p and the follow-up tests are for.

Benchmarks p at a Glance

Very Strong

p < .001

Strong

.001 ≤ p < .01

Significant

.01 ≤ p < .05

Weak

.05 ≤ p < .10

Little

p ≥ .10

Strength of evidence against the null hypothesis.

Full benchmarks and caveats ↑

Power and Sample Size

Power is the chance that a study detects an effect that's really there. A real effect doesn't guarantee a significant result: with too few people, it can easily come out non-significant. For a MANCOVA, power depends on how big the adjusted group differences are on each outcome, how strongly the outcomes correlate (after the covariate), how much the covariate explains, how many people are in each group, and α.

That's too many ingredients for a simple calculator, so there isn't one on this page. Dedicated tools such as G*Power ask for the multivariate effect directly, but the most transparent approach is to simulate: write down the means, spreads, correlations and covariate slopes you expect, generate many fake studies, and count how often the MANCOVA comes out significant.

We simulated the exercise study under these assumptions: program differences of 2 and 5 symbols and 0.4 and 1 word (stretching as the baseline); scores falling 0.9 symbols and 0.2 words per year of age; ages centered on 71 (SD 4.5, kept within 65 to 82); leftover spreads, after age, of 8 symbols and 2 words, correlated .30; and 4,000 simulated studies at each sample size, each tested with Pillai's trace at α = .05. With 30 people per program, as here, the MANCOVA came out significant only about 53% of the time. It took about 50 per program, 150 in all, to reach 80%.

Those are the population values the exercise study's data were simulated from. Our sample happened to come out a little less noisy (about 7.0 symbols and 1.7 words left over after age), which is part of why it reached significance with only 30 per program. Plan with spreads from earlier studies, not from your own data; smaller spreads would mean fewer people.

A good covariate is one of the cheapest ways to buy power: age soaked up a large share of the noise, which is why the adjusted test found an effect the unadjusted one (p = .171) missed. Each extra outcome that doesn't differ between groups costs a little power.

Researchers usually aim for at least 80% power. Plan it before collecting data, using the smallest effect that would matter.

Other Info

Follow-Up and Post Hoc Tests

Post Hoc Tests

Bonferroni Correction: The most common follow-up: one ANCOVA per outcome, each at α ÷ the number of outcomes. It's step 9 of the Napkin Math: at α = .025, speed differed, F(2, 86) = 4.65, p = .012, but recall didn't quite, F(2, 86) = 3.55, p = .033. Then compare adjusted means pairwise, Bonferroni-corrected for the three pairs: on speed, aerobic walking beat stretching by 5.49 symbols, p = .009; the other two pairs didn't differ, p > .25.
Holm Correction: Gives the same protection as Bonferroni and is uniformly more powerful. Had we chosen it in advance, Holm would test the smaller p (.012) at .025 and then the larger (.033) at .05, and recall would also have counted. We chose Bonferroni before seeing the data, so the write-up uses Bonferroni. Choose the correction in advance and stick with it; switching after seeing which one gives the result you want is its own kind of p-hacking.
Descriptive Discriminant Analysis: Use when the question is about the outcomes together: which combination of speed and recall best separates the programs. Separate ANCOVAs ignore how the outcomes relate, which is what the multivariate test was about (Huberty & Morris, 1989).
Roy–Bargmann Stepdown Analysis: Use only if the outcomes had an order of importance before the data came in, say speed first. Speed would get its age-adjusted ANCOVA, as in step 9; recall would then be compared with both age and speed as covariates, asking whether the programs differ in recall among people of the same age and speed. Like discriminant analysis, it respects the link between the outcomes, but in Finch's (2007) simulations the tests after the first tended to reject true nulls too often, so we didn't use it here.

Related Methods

General CaseMultivariate Regression: Regress both outcomes at once on k − 1 group dummy variables and the covariate: the multivariate test of the dummies is this MANCOVA. MANCOVA is multivariate regression with a categorical predictor and a continuous one.
Special CaseANCOVA: MANCOVA with a single outcome. With one outcome, all four multivariate statistics reduce to the ANCOVA's F.
Special CaseMANOVA: MANCOVA with no covariates. Drop the covariate and the groups are compared on their unadjusted outcome averages.
AlternativeMultilevel Regression (HLM): An alternative if “Independent Observations” doesn't hold. A multilevel model, which accounts for the clustering.

References

Wilks, S. S. (1932). Certain generalizations in the analysis of variance. Biometrika, 24(3–4), 471–494. https://doi.org/10.1093/biomet/24.3-4.471
Pillai, K. C. S. (1955). Some new test criteria in multivariate analysis. The Annals of Mathematical Statistics, 26(1), 117–121. https://doi.org/10.1214/aoms/1177728599
“Three new test criteria are proposed for overall tests of hypotheses in multivariate analysis.” (Abstract)
Olson, C. L. (1976). On choosing a test statistic in multivariate analysis of variance. Psychological Bulletin, 83(4), 579–586. https://doi.org/10.1037/0033-2909.83.4.579
Box, G. E. P. (1949). A general distribution theory for a class of likelihood criteria. Biometrika, 36(3–4), 317–346. https://doi.org/10.1093/biomet/36.3-4.317
Huberty, C. J., & Morris, J. D. (1989). Multivariate analysis versus multiple univariate analyses. Psychological Bulletin, 105(2), 302–308. https://doi.org/10.1037/0033-2909.105.2.302
Miller, G. A., & Chapman, J. P. (2001). Misunderstanding analysis of covariance. Journal of Abnormal Psychology, 110(1), 40–48. https://doi.org/10.1037/0021-843X.110.1.40
Tabachnick, B. G., & Fidell, L. S. (2019). Using multivariate statistics (7th ed.). Pearson.
Cohen, J. (1988). Statistical power analysis for the behavioral sciences (2nd ed.). Lawrence Erlbaum Associates.
Olson, C. L. (1974). Comparative robustness of six tests in multivariate analysis of variance. Journal of the American Statistical Association, 69(348), 894–908. https://doi.org/10.1080/01621459.1974.10480224
“the Pillai-Bartlett trace test may be recommended as the most robust of the MANOVA tests, with adequate power to detect true differences in a variety of situations” (Abstract, p. 894)
Cohen, J. (1973). Eta-squared and partial eta-squared in fixed factor ANOVA designs. Educational and Psychological Measurement, 33(1), 107–112. https://doi.org/10.1177/001316447303300111
Pierce, C. A., Block, R. A., & Aguinis, H. (2004). Cautionary note on reporting eta-squared values from multifactor ANOVA designs. Educational and Psychological Measurement, 64(6), 916–924. https://doi.org/10.1177/0013164404264848
Lakens, D. (2013). Calculating and reporting effect sizes to facilitate cumulative science: A practical primer for t-tests and ANOVAs. Frontiers in Psychology, 4, 863. https://doi.org/10.3389/fpsyg.2013.00863
“using these benchmarks when interpreting the η2p effect size in designs that include covariates or repeated measures is not consistent with the considerations upon which the benchmarks were based” (Eta-squared in between and within-subjects comparisons)
Muff, S., Nilsen, E. B., O'Hara, R. B., & Nater, C. R. (2022). Rewriting results sections in the language of evidence. Trends in Ecology & Evolution, 37(3), 203–210. https://doi.org/10.1016/j.tree.2021.10.009
“we suggest language of evidence that allows for a more nuanced approach to communicate scientific findings as a simple and intuitive alternative to statistical significance testing” (Abstract)
Friendly, M., & Sigal, M. (2020). Visualizing tests for equality of covariance matrices. The American Statistician, 74(2), 144–155. https://doi.org/10.1080/00031305.2018.1497537
“Box's test is well-known to be highly sensitive to violation of (multivariate) normality and the presence of outliers.” (Section 1.3)
Warne, R. T. (2014). A primer on multivariate analysis of variance (MANOVA) for behavioral scientists. Practical Assessment, Research & Evaluation, 19(17), 1–10. https://doi.org/10.7275/sm63-7h70
Finch, W. H. (2007). Performance of the Roy-Bargmann stepdown procedure as a follow up to a significant MANOVA. Multiple Linear Regression Viewpoints, 33(1), 12–22. https://ojs.lib.ua.edu/glmj/article/view/218
“the current study does not support the use of this technique in trying to identify specific variables for which the groups might differ, other than the first one” (Discussion, p. 20)