MANOVA in R: manova(), summary() and Reading the Output

MANOVA in R is run with manova() from the base stats package, followed by summary() to get the multivariate test and summary.aov() to get one ANOVA table per outcome. You pass the outcomes together as a matrix on the left of the formula, for example manova(cbind(y1, y2, y3) ~ group). The part most tutorials skip is choosing the test statistic: R’s default is Pillai, and the R documentation spells out when the four options give identical answers.

What MANOVA is used for in research

MANOVA tests whether group means differ on several related outcome variables considered together, instead of running a separate ANOVA on each. Use it when you have two or more continuous outcomes measured on the same units and one or more categorical factors. If you only have one outcome, see the guide to analysis of variance; if you are unsure which test fits your design, start with how to choose a statistical test. The distinction between several predictors and several outcomes is covered in multivariable vs multivariate analysis.

Install or access (official source only)

MANOVA needs no extra package: manova() and summary.manova() are part of the stats package that ships with R. Download R from CRAN, the Comprehensive R Archive Network. As checked on 2026-10-29, the CRAN Windows base page offers R-4.6.1 (a 64-bit installer, R-4.6.1-win.exe) at cran.r-project.org/bin/windows/base; builds for other systems are linked from CRAN. The assumption check shown below uses the rstatix package (install.packages("rstatix")).

Step-by-step: fitting a MANOVA with manova()

The worked example is the plastic-film data from Krzanowski (1998, p. 381), used in R’s own documentation for summary.manova: three outcomes (tear, gloss, opacity) measured on 20 films under two factors, rate and additive.

1. Enter the data and bind the outcomes into a matrix.

   tear <- c(6.5, 6.2, 5.8, 6.5, 6.5, 6.9, 7.2, 6.9, 6.1, 6.3,
             6.7, 6.6, 7.2, 7.1, 6.8, 7.1, 7.0, 7.2, 7.5, 7.6)
   gloss <- c(9.5, 9.9, 9.6, 9.6, 9.2, 9.1, 10.0, 9.9, 9.5, 9.4,
              9.1, 9.3, 8.3, 8.4, 8.5, 9.2, 8.8, 9.7, 10.1, 9.2)
   opacity <- c(4.4, 6.4, 3.0, 4.1, 0.8, 5.7, 2.0, 3.9, 1.9, 5.7,
                2.8, 4.1, 3.8, 1.6, 3.4, 8.4, 5.2, 6.9, 2.7, 1.9)
   Y <- cbind(tear, gloss, opacity)
   rate     <- gl(2, 10, labels = c("Low", "High"))
   additive <- gl(2, 5, length = 20, labels = c("Low", "High"))

2. Fit the model. manova() passes its arguments to aov() and then adds the manova class; the two differ mainly in which summary method is used.

   fit <- manova(Y ~ rate * additive)

3. Get the multivariate tests.

   summary(fit)                  # Pillai (default)
   summary(fit, test = "Wilks")  # Wilks' lambda

4. Get the follow-up univariate tables.

   summary.aov(fit)   # one ANOVA table per outcome

(Source: R documentation, summary.manova, example section, checked 2026-10-29.)

Choosing the test statistic

The usage line is summary(object, test = c("Pillai", "Wilks", "Hotelling-Lawley", "Roy"), intercept = FALSE, tol = 1e-7, ...). The first option, "Pillai", is the default, and partial matching of the name is allowed. According to the R documentation:

test =	What the R documentation says
“Pillai”	Pillai-Bartlett trace; the default, recommended by Hand and Taylor (1987)
“Wilks”	Wilks’ lambda; the most widely used in the literature but not the default
“Hotelling-Lawley”	Hotelling-Lawley trace; one of the four available statistics
“Roy”	Roy’s greatest root; gives an upper bound unless the term has 1 df

When a term has 1 degree of freedom, all four approximate F statistics are identical, so the choice only matters for terms with more than 1 df. In the plastic-film model, rate and additive each have two levels (1 df), which is why the documentation notes that the default summary gives the same F statistics as the Wilks version for single-df terms.

Output and how to read it

summary(fit) returns an object of class summary.manova. Its printed table has one row per model term with the test statistic, the approximate F, the numerator and denominator degrees of freedom, and the p-value. Stored in the object are the sums-of-squares-and-products matrices (SS), the Eigenvalues, and the stats matrix holding those columns, so you can extract values programmatically, for example summary(fit)$stats.

Read the multivariate row first: a small p-value for a term means the vector of means differs across that factor’s levels on at least one outcome or combination of outcomes. Then use summary.aov(fit) to see which outcomes carry the effect. Because that runs one univariate test per outcome, keep the number of follow-up tests in mind when you interpret them. For pairwise group comparisons on a single outcome, see Tukey HSD test. By default intercept = FALSE, so the intercept is not tested.

Check the assumptions in R

MANOVA assumes, among other things, that the covariance matrices of the outcomes are equal across groups. Box’s M test is a common check. The rstatix package provides box_m(data, group), where data is a numeric data frame or matrix with n rows and p columns (n greater than p) and group is a factor of length n. It returns the test statistic (an approximate chi-square), its degrees of freedom, the p-value, and the method name “Box’s M-test for Homogeneity of Covariance Matrices”. The package’s own example is box_m(iris[, -5], iris[, 5]). In general, a small p-value indicates the covariance matrices differ between groups. With the plastic-film data and one grouping factor: library(rstatix); box_m(as.data.frame(Y), rate).

For outliers across several variables, see Mahalanobis distance, and for a normality check on each outcome, the Shapiro-Wilk test in R.

Troubleshooting

• You get a single ANOVA instead of a multivariate test. The left side of the formula must be a matrix of outcomes (cbind(...)), and the object passed to summary() must be a manova object, or an aov object with multiple responses. Per the documentation, summary.manova accepts either.
• You want a different statistic. Pass test = "Wilks", "Hotelling-Lawley" or "Roy"; partial matching of the name is allowed.
• The intercept row is missing. That is the default (intercept = FALSE); set intercept = TRUE to include it.
• Highly correlated outcomes. The tol argument (default 1e-7) is used in deciding whether the residuals are rank-deficient. The R documentation says the default guards against very highly correlated responses, that lowering it allows rather inaccurate results, and that it is normally better to transform the responses to remove the high correlation.
• The four statistics disagree. This can only happen for terms with more than 1 df, as above; report which statistic you chose.

Other software

R is the route covered here. If your department works in SPSS, the equivalent procedure is covered in MANOVA in SPSS.

FAQ

What is the default test in manova() in R?

Pillai’s trace. summary.manova uses test = "Pillai" unless you set another value.

Do I need a package to run MANOVA in R?

No. manova() is in the base stats package. Only the Box’s M check in this guide needs rstatix.

How do I get the univariate follow-up tests?

Call summary.aov() on the fitted object; it returns one ANOVA table per outcome.

Sources

• R Core Team. summary.manova, R documentation (stats package). Checked 2026-10-29.
• R Core Team. manova, R documentation (stats package). Checked 2026-10-29.
• rstatix documentation. box_m. Checked 2026-10-29.
• The Comprehensive R Archive Network. R for Windows. Checked 2026-10-29.