One-way multivariate analysis of variance the modern way — test a grouping factor against two or more continuous outcomes together, check the multivariate assumptions, and follow up with univariate ANOVAs
Learn to run a one-way MANOVA in R. Test whether groups differ on several continuous outcome variables simultaneously with car’s Manova() and Pillai’s trace, check the MANOVA assumptions — multivariate normality, multivariate outliers, multicollinearity, homogeneity of covariances (Box’s M) and of variances — then follow a significant result with univariate ANOVAs and pairwise comparisons, with the p-values on the plot.
lm(cbind(y1, y2) ~ group) then car::Manova(); read Pillai’s trace (the robust default, best when assumptions are shaky or the design is unbalanced).mshapiro_test()), no multivariate outliers (Mahalanobis distance), moderate correlation between outcomes (no multicollinearity), and homogeneity of covariance matrices (Box’s M test).A one-way MANOVA extends the one-way ANOVA to several outcome variables measured together. Instead of asking “do the groups differ on this score?” it asks “do the groups differ on the combination of these scores?” — testing a single grouping factor against a vector of continuous outcomes.
Use it when your outcomes measure different facets of one cohesive theme — several exam scores for overall academic performance, or several flower measurements for plant size — and you want one test of the joint difference rather than many separate ANOVAs (which inflate the error rate). The classic scenario: do iris species differ on the combination of Sepal.Length and Petal.Length? This lesson runs the full workflow with rstatix and car, using the built-in iris data.
One-way MANOVA hypotheses. H₀: the group mean vectors are equal (the groups do not differ on the combined outcomes). Hₐ: at least one group’s mean vector differs. MANOVA builds a linear combination of the outcomes and compares its group means.
We use the built-in iris dataset, keeping two outcomes — Sepal.Length and Petal.Length — and the grouping factor Species (setosa, versicolor, virginica, 50 each):
library(rstatix)
iris2 <- data.frame(id = 1:nrow(iris), iris[, c("Sepal.Length", "Petal.Length", "Species")])
head(iris2) id Sepal.Length Petal.Length Species
1 1 5.1 1.4 setosa
2 2 4.9 1.4 setosa
3 3 4.7 1.3 setosa
4 4 4.6 1.5 setosa
5 5 5.0 1.4 setosa
6 6 5.4 1.7 setosaA merged box plot shows both outcomes by species side by side:
library(ggpubr)
iris2 <- data.frame(id = 1:nrow(iris), iris[, c("Sepal.Length", "Petal.Length", "Species")])
ggboxplot(
iris2, x = "Species", y = c("Sepal.Length", "Petal.Length"),
merge = TRUE, palette = "jco"
)
Summarize each outcome by group with get_summary_stats():
library(rstatix)
iris2 <- data.frame(id = 1:nrow(iris), iris[, c("Sepal.Length", "Petal.Length", "Species")])
iris2 %>%
group_by(Species) %>%
get_summary_stats(Sepal.Length, Petal.Length, type = "mean_sd")# A tibble: 6 × 5
Species variable n mean sd
<fct> <fct> <dbl> <dbl> <dbl>
1 setosa Sepal.Length 50 5.01 0.352
2 setosa Petal.Length 50 1.46 0.174
3 versicolor Sepal.Length 50 5.94 0.516
4 versicolor Petal.Length 50 4.26 0.47
5 virginica Sepal.Length 50 6.59 0.636
6 virginica Petal.Length 50 5.55 0.552MANOVA assumes adequate sample size (n per group > number of outcomes — here 50 > 2 ✓), independent observations, no multivariate outliers, multivariate normality, moderate (not excessive) correlation between outcomes, a linear relationship between the outcomes within each group, and homogeneity of covariance matrices.
First, check each outcome for extreme values within each group with identify_outliers() (an empty result means none):
library(rstatix)
iris2 <- data.frame(id = 1:nrow(iris), iris[, c("Sepal.Length", "Petal.Length", "Species")])
iris2 %>% group_by(Species) %>% identify_outliers(Sepal.Length)# A tibble: 1 × 6
Species id Sepal.Length Petal.Length is.outlier is.extreme
<fct> <int> <dbl> <dbl> <lgl> <lgl>
1 virginica 107 4.9 4.5 TRUE FALSE iris2 %>% group_by(Species) %>% identify_outliers(Petal.Length)# A tibble: 5 × 6
Species id Sepal.Length Petal.Length is.outlier is.extreme
<fct> <int> <dbl> <dbl> <lgl> <lgl>
1 setosa 14 4.3 1.1 TRUE FALSE
2 setosa 23 4.6 1 TRUE FALSE
3 setosa 25 4.8 1.9 TRUE FALSE
4 setosa 45 5.1 1.9 TRUE FALSE
5 versicolor 99 5.1 3 TRUE FALSE There were no univariate extreme outliers in either outcome.
Multivariate outliers are unusual combinations of values. mahalanobis_distance() flags them per group:
library(rstatix)
iris2 <- data.frame(id = 1:nrow(iris), iris[, c("Sepal.Length", "Petal.Length", "Species")])
iris2 %>%
group_by(Species) %>%
mahalanobis_distance(-id) %>%
filter(is.outlier == TRUE) %>%
as.data.frame()[1] id Sepal.Length Petal.Length mahal.dist is.outlier
<0 lignes> (ou 'row.names' de longueur nulle)The empty result means no multivariate outliers (p > 0.001).
Test each outcome per group with Shapiro-Wilk, and the outcomes jointly with mshapiro_test():
library(rstatix)
iris2 <- data.frame(id = 1:nrow(iris), iris[, c("Sepal.Length", "Petal.Length", "Species")])
# univariate normality, each outcome per group
iris2 %>%
group_by(Species) %>%
shapiro_test(Sepal.Length, Petal.Length)# A tibble: 6 × 4
Species variable statistic p
<fct> <chr> <dbl> <dbl>
1 setosa Petal.Length 0.955 0.0548
2 setosa Sepal.Length 0.978 0.460
3 versicolor Petal.Length 0.966 0.158
4 versicolor Sepal.Length 0.978 0.465
5 virginica Petal.Length 0.962 0.110
6 virginica Sepal.Length 0.971 0.258 # multivariate normality
iris2[, c("Sepal.Length", "Petal.Length")] %>% mshapiro_test()# A tibble: 1 × 2
statistic p.value
<dbl> <dbl>
1 0.995 0.855Both outcomes are normal in every group (p > 0.05), and mshapiro_test() is non-significant (p = 0.86) — multivariate normality holds. Back the test with a Q-Q plot per group:
library(ggpubr)
iris2 <- data.frame(id = 1:nrow(iris), iris[, c("Sepal.Length", "Petal.Length", "Species")])
ggqqplot(iris2, "Sepal.Length", facet.by = "Species", ggtheme = theme_bw())
With n > 50 per group the Q-Q plot is the better guide — Shapiro-Wilk becomes oversensitive to tiny departures at large samples. Here both agree.
The outcomes should be moderately correlated — not above r = 0.9 (multicollinearity), not near zero (then run separate ANOVAs):
library(rstatix)
iris2 <- data.frame(id = 1:nrow(iris), iris[, c("Sepal.Length", "Petal.Length", "Species")])
iris2 %>% cor_test(Sepal.Length, Petal.Length)# A tibble: 1 × 9
var1 var2 cor statistic df p conf.low conf.high method
<chr> <chr> <dbl> <dbl> <int> <dbl> <dbl> <dbl> <chr>
1 Sepal.Length Petal.L… 0.87 21.6 148 1.04e-47 0.827 0.906 Pears…The outcomes correlate at r = 0.87 (p < 0.0001) — high but below 0.9, so acceptable.
MANOVA builds a linear combination of the outcomes, so it assumes a linear relationship between the outcome variables within each group — a curved relationship would distort the multivariate test. Check it visually with a scatterplot matrix per group, using ggpairs() from the GGally package and colouring the points by species:
library(GGally)
iris2 <- data.frame(id = 1:nrow(iris), iris[, c("Sepal.Length", "Petal.Length", "Species")])
ggpairs(
iris2, columns = c("Sepal.Length", "Petal.Length"),
mapping = aes(color = Species)
)
The off-diagonal panel plots Sepal.Length against Petal.Length for each species: within every group the points fall along a straight line with no curved or bent pattern, so the linearity assumption holds. (The upper panel prints the per-group correlations too — these vary in strength, but linearity is about the shape of each cloud, not how tight it is.)
If instead a group’s cloud were curved or U-shaped, the relationship would be non-linear. In that case, transform or remove the offending outcome variable, or run the MANOVA anyway and accept that you lose some statistical power.
Box’s M is the multivariate analogue of equal variances. It is very sensitive, so judge it at α = 0.001:
library(rstatix)
iris2 <- data.frame(id = 1:nrow(iris), iris[, c("Sepal.Length", "Petal.Length", "Species")])
box_m(iris2[, c("Sepal.Length", "Petal.Length")], iris2$Species)# A tibble: 1 × 4
statistic p.value parameter method
<dbl> <dbl> <dbl> <chr>
1 58.4 9.62e-11 6 Box's M-test for Homogeneity of Covariance Matri…Box’s M is significant (p < 0.001), so the covariance matrices are not homogeneous. With a balanced design (50 per group) this is not fatal — but it is the reason to report Pillai’s trace, the robust statistic, rather than Wilks’ lambda.
For each outcome separately, reshape to long form with base R, then Levene’s test per outcome:
library(rstatix)
iris2 <- data.frame(id = 1:nrow(iris), iris[, c("Sepal.Length", "Petal.Length", "Species")])
long <- data.frame(
Species = rep(iris2$Species, 2),
variable = rep(c("Sepal.Length", "Petal.Length"), each = nrow(iris2)),
value = c(iris2$Sepal.Length, iris2$Petal.Length)
)
long %>% group_by(variable) %>% levene_test(value ~ Species)# A tibble: 2 × 5
variable df1 df2 statistic p
<chr> <int> <int> <dbl> <dbl>
1 Petal.Length 2 147 19.5 0.0000000313
2 Sepal.Length 2 147 6.35 0.00226 Levene’s test is significant for both outcomes — variances are not equal across species, which (with Box’s M) again points to Pillai’s trace and to Games-Howell for the post-hoc step.
Fit a multivariate linear model with cbind() on the left, then Manova() with Pillai’s trace:
library(rstatix)
library(car)
iris2 <- data.frame(id = 1:nrow(iris), iris[, c("Sepal.Length", "Petal.Length", "Species")])
model <- lm(cbind(Sepal.Length, Petal.Length) ~ Species, data = iris2)
Manova(model, test.statistic = "Pillai")
Type II MANOVA Tests: Pillai test statistic
Df test stat approx F num Df den Df Pr(>F)
Species 2 0.9885 71.829 4 294 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1There is a statistically significant difference between species on the combined outcomes, F(4, 294) = 71.83, p < 0.0001 — species differ on the sepal-length/petal-length combination.
Which multivariate statistic? Manova() offers Pillai, Wilks, Hotelling-Lawley and Roy. Wilks’ lambda is the common default, but Pillai’s trace is more robust and is preferred when the design is unbalanced or Box’s M is significant (our case) — so we report Pillai.
A significant MANOVA is decomposed into a univariate one-way ANOVA per outcome (Bonferroni-adjusted to α = 0.05 / 2 = 0.025 for two outcomes), then pairwise comparisons. Because variances are unequal, use the Games-Howell test:
library(rstatix)
iris2 <- data.frame(id = 1:nrow(iris), iris[, c("Sepal.Length", "Petal.Length", "Species")])
long <- data.frame(
Species = rep(iris2$Species, 2),
variable = rep(c("Sepal.Length", "Petal.Length"), each = nrow(iris2)),
value = c(iris2$Sepal.Length, iris2$Petal.Length)
)
# univariate one-way ANOVA per outcome
long %>% group_by(variable) %>% anova_test(value ~ Species)# A tibble: 2 × 8
variable Effect DFn DFd F p `p<.05` ges
* <chr> <chr> <dbl> <dbl> <dbl> <dbl> <chr> <dbl>
1 Petal.Length Species 2 147 1180. 2.86e-91 * 0.941
2 Sepal.Length Species 2 147 119. 1.67e-31 * 0.619# pairwise comparisons per outcome (Games-Howell, unequal variances)
long %>% group_by(variable) %>% games_howell_test(value ~ Species)# A tibble: 6 × 9
variable .y. group1 group2 estimate conf.low conf.high p.adj p.adj.signif
* <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <chr>
1 Petal.L… value setosa versi… 2.80 2.63 2.97 1.85e-11 ****
2 Petal.L… value setosa virgi… 4.09 3.89 4.29 1.68e-11 ****
3 Petal.L… value versi… virgi… 1.29 1.05 1.54 4.45e-10 ****
4 Sepal.L… value setosa versi… 0.93 0.719 1.14 2.86e-10 ****
5 Sepal.L… value setosa virgi… 1.58 1.34 1.83 0 ****
6 Sepal.L… value versi… virgi… 0.652 0.376 0.928 5.58e- 7 **** Both outcomes differ significantly by species — Sepal.Length F(2, 147) = 119 and Petal.Length F(2, 147) = 1180 (both p < 0.0001) — and every pairwise species comparison is significant for each outcome.
The reporting figure overlays the pairwise comparisons on the merged box plot, with the MANOVA result in the subtitle (built with create_test_label()):
library(rstatix)
library(ggpubr)
iris2 <- data.frame(id = 1:nrow(iris), iris[, c("Sepal.Length", "Petal.Length", "Species")])
long <- data.frame(
Species = rep(iris2$Species, 2),
variables = rep(c("Sepal.Length", "Petal.Length"), each = nrow(iris2)),
value = c(iris2$Sepal.Length, iris2$Petal.Length)
)
pwc <- long %>%
group_by(variables) %>%
games_howell_test(value ~ Species) %>%
add_xy_position(x = "Species")
test.label <- create_test_label(
description = "MANOVA", statistic.text = quote(italic("F")),
statistic = 71.83, p = "<0.0001", parameter = "4,294",
type = "expression", detailed = TRUE
)
ggboxplot(
iris2, x = "Species", y = c("Sepal.Length", "Petal.Length"),
merge = TRUE, palette = "jco"
) +
stat_pvalue_manual(
pwc, hide.ns = TRUE, tip.length = 0,
step.increase = 0.1, step.group.by = "variables", color = "variables"
) +
labs(subtitle = test.label, caption = get_pwc_label(pwc, type = "expression"))
A one-way MANOVA was performed to determine the effect of iris species on sepal length and petal length. There was a statistically significant difference between species on the combined dependent variables, F(4, 294) = 71.83, p < 0.0001 (Pillai’s trace). Follow-up univariate ANOVAs (Bonferroni-adjusted α = 0.025) showed a significant difference in sepal length, F(2, 147) = 119, and petal length, F(2, 147) = 1180 (both p < 0.0001). All pairwise comparisons between species were significant for each outcome.
Where ANOVA partitions a single sum of squares, MANOVA partitions matrices: a between-groups matrix \(H\) (hypothesis) and a within-groups matrix \(E\) (error). The multivariate statistics summarise the \(H\)-vs-\(E\) relationship through the eigenvalues \(\lambda_i\) of \(E^{-1}H\). Pillai’s trace is \(V = \sum_i \frac{\lambda_i}{1 + \lambda_i}\) — a bounded, robust measure of how much of the combined variation lies between groups. Each multivariate statistic is converted to an approximate \(F\) for the p-value.
Run a MANOVA on a fresh question — do iris species differ on the petal measurements (Petal.Length and Petal.Width) together? Edit the outcomes or the test statistic and re-run; the sandbox boots on first Run.
Ask Prova “I have one grouping factor and several numeric outcomes — run a one-way MANOVA on my data, check the multivariate assumptions, and follow up a significant result with univariate ANOVAs” — it answers with rstatix + car code you can run on your own data, then helps you read Pillai’s trace and the follow-up tests. The runtime is the judge. Ask Prova →
You ran a separate ANOVA per outcome instead of a MANOVA. Testing each outcome separately inflates the family-wise error rate and ignores the correlation between outcomes. Run the MANOVA first for the joint test; only then decompose into univariate ANOVAs (with a Bonferroni-adjusted α).
Box’s M is significant, so you used Wilks’ lambda. When Box’s M fails (unequal covariance matrices) or the design is unbalanced, Wilks’ lambda is not robust. Report Pillai’s trace instead — it is far less sensitive to the violation.
Your outcomes are almost perfectly correlated (r > 0.9). That is multicollinearity — the outcomes carry nearly the same information, which destabilises the MANOVA. Drop one of the redundant outcomes, or combine them into a single composite score.
ANOVA tests a grouping factor against one continuous outcome; MANOVA tests it against two or more outcomes simultaneously, accounting for the correlation between them. MANOVA gives a single joint test (avoiding the inflated error rate of many separate ANOVAs) and can detect group differences that appear only in the combination of outcomes.
Wilks’ lambda is the traditional default, but Pillai’s trace is more robust to violations of the homogeneity-of-covariance (Box’s M) assumption and to unbalanced designs. When Box’s M is significant — as in most real data — report Pillai’s trace. car::Manova() uses Pillai by default.
Adequate sample size (n per group greater than the number of outcomes), independent observations, no multivariate outliers (Mahalanobis distance), multivariate normality (mshapiro_test()), a linear relationship between outcomes, no multicollinearity (outcome correlations below 0.9), and homogeneity of covariance matrices (Box’s M). MANOVA is fairly robust to mild departures, especially with a balanced design.
Draw a scatterplot matrix of the outcome variables for each group with ggpairs() from the GGally package (colour the points by group). The relationship between each pair of outcomes should look linear — a straight cloud with no curve — within every group. If a group shows a curved pattern, transform or remove that outcome variable, or run the MANOVA anyway and accept a loss of power.
Decompose it into a univariate one-way ANOVA for each outcome, using a Bonferroni-adjusted significance level (α divided by the number of outcomes). For each outcome that differs, run pairwise comparisons — Tukey HSD if variances are equal, Games-Howell if they are not — to see which groups differ.
Sepal.Length and Sepal.Width by Species and report Pillai’s trace. Fill the blanks. Is the multivariate effect significant?Fill the outcomes as cbind(Sepal.Length, Sepal.Width) and the statistic as "Pillai". Read the Pr(>F) column. For question 2, recall which statistic is robust when the covariance matrices are unequal.
library(rstatix); library(car)
model <- lm(cbind(Sepal.Length, Sepal.Width) ~ Species, data = iris)
Manova(model, test.statistic = "Pillai")
#> Species effect is highly significant (p < 2e-16).The multivariate effect is significant. For question 2: report Pillai’s trace — it is the robust statistic, least affected by a significant Box’s M (unequal covariance matrices).
You can now run a one-way MANOVA in R: fit lm(cbind(y1, y2) ~ group), test it with car::Manova() reporting Pillai’s trace, and check the multivariate assumptions — multivariate normality, multivariate outliers, multicollinearity and Box’s M. When the joint test is significant, decompose it into univariate ANOVAs (Bonferroni-adjusted) and pairwise comparisons to see which groups differ on which outcome, then put the comparisons on a merged box plot. MANOVA is the right tool whenever several outcomes measure one cohesive theme.
@online{2026,
author = {},
title = {MANOVA in {R:} {Compare} {Groups} on {Several} {Outcomes} at
{Once}},
date = {2026-06-23},
url = {https://www.datanovia.com/learn/biostatistics/anova/manova-in-r},
langid = {en}
}