Multivariate analysis of variance

Multivariate analysis of variance, commonly abbreviated MANOVA, is a family of statistical methods for testing hypotheses about several correlated response variables within a common linear model. It generalizes analysis of variance by representing each observational unit with a response vector rather than a single scalar outcome. The resulting hypothesis concerns differences among multivariate means after accounting for the covariance structure shared by the responses.

MANOVA is mathematically equivalent to testing specified parameter restrictions in the multivariate general linear model. Its principal test statistics summarize the characteristic roots of two matrices: one representing variation attributable to the tested hypothesis and another representing residual variation. Different summaries emphasize different geometric aspects of the same fitted model and therefore need not produce identical significance levels.

Statistical model

For (n) observational units and (p) response variables, the multivariate linear model is

[ \mathbf Y=\mathbf X\mathbf B+\mathbf U, ]

where (\mathbf Y) is an (n\times p) response matrix, (\mathbf X) is an (n\times q) design matrix, and (\mathbf B) contains the unknown regression coefficients. The error matrix (\mathbf U) contains one (p)-dimensional error vector for each observational unit.

Under the classical Gaussian formulation, the error vectors are mutually independent and follow a multivariate normal distribution with mean vector zero and common covariance matrix (\boldsymbol\Sigma). The off-diagonal elements of (\boldsymbol\Sigma) represent residual associations among the responses. These associations distinguish MANOVA from a collection of unrelated univariate models.

A general multivariate hypothesis may be written as

[ H_0:\mathbf C\mathbf B\mathbf M=\boldsymbol\Theta_0, ]

where (\mathbf C) defines contrasts among model coefficients and (\mathbf M) defines linear combinations of the response variables. The matrix (\boldsymbol\Theta_0) specifies the values imposed by the null hypothesis. Ordinary factorial MANOVA usually takes (\mathbf M) to be the identity matrix, so that the tested effect applies jointly to the original response vector.

Hypothesis and error variation

The fitted model partitions multivariate variation into sums-of-squares-and-cross-products matrices. The residual matrix is

[ \mathbf E =\mathbf Y^{\mathsf T} \left(\mathbf I-\mathbf P_X\right) \mathbf Y, ]

where (\mathbf P_X) is the projection matrix onto the column space of (\mathbf X). Diagonal elements of (\mathbf E) are residual sums of squares, whereas its off-diagonal elements are residual cross-products between responses.

The hypothesis matrix (\mathbf H) measures the fitted variation associated with the tested restriction. For a full-rank contrast with null value zero, it has the form

[ \mathbf H= (\mathbf C\widehat{\mathbf B})^{\mathsf T} \left[ \mathbf C(\mathbf X^{\mathsf T}\mathbf X)^{-1} \mathbf C^{\mathsf T} \right]^{-1} (\mathbf C\widehat{\mathbf B}). ]

Both (\mathbf H) and (\mathbf E) depend on the coordinate system used to express the responses, but the characteristic roots used by standard MANOVA criteria remain unchanged under any nonsingular linear transformation of that coordinate system. This invariance permits equivalent interpretations through the original variables or through derived canonical variates.

When (\mathbf E) is nonsingular, the roots (\lambda_i) are obtained from

[ \det(\mathbf H-\lambda\mathbf E)=0. ]

Each nonzero root corresponds to a response-space direction along which hypothesis variation is compared with residual variation. The number of such directions cannot exceed the smaller of the hypothesis rank and the response dimension.

Test criteria

Wilks's lambda

Wilks's lambda is the likelihood-ratio criterion

[ \Lambda =\frac{\det(\mathbf E)} {\det(\mathbf E+\mathbf H)} =\prod_i\frac{1}{1+\lambda_i}. ]

Values near one correspond to relatively little hypothesis variation, while smaller values correspond to stronger separation in at least part of the multivariate response space. Its exact null distribution is available for several restricted designs. More general cases use transformations or approximations based on the chi-squared distribution or the F-distribution.

Pillai's trace

Pillai's trace is

[ V =\operatorname{tr} \left[ \mathbf H(\mathbf H+\mathbf E)^{-1} \right] =\sum_i\frac{\lambda_i}{1+\lambda_i}. ]

The statistic adds the proportionate contribution of every characteristic direction. Because each transformed root is bounded between zero and one, no single unbounded root determines the entire value.

Lawley–Hotelling trace

The Lawley–Hotelling trace is

[ T=\operatorname{tr}(\mathbf E^{-1}\mathbf H) =\sum_i\lambda_i. ]

This criterion sums the untransformed roots and consequently gives substantial weight to directions having a large ratio of hypothesis variation to residual variation. In a one-dimensional response space, it reduces to the ratio underlying the ordinary univariate (F)-test.

Roy's largest root

Roy's largest root is

[ \Theta=\max_i\lambda_i. ]

It represents the greatest attainable hypothesis-to-error ratio over linear combinations of the responses. Some formulations instead report the bounded transformation (\Theta/(1+\Theta)), so numerical values depend on the convention even when the underlying characteristic direction is identical.

Geometric interpretation

Each row of (\mathbf Y) is a point in a (p)-dimensional response space. Group means or fitted treatment means occupy the same space, while the residual covariance matrix describes the orientation and dispersion of observations around those fitted locations. A multivariate group effect therefore concerns the displacement of mean vectors relative to the residual covariance geometry.

The eigenvectors associated with (\mathbf E^{-1}\mathbf H) define linear combinations that order the fitted separation from greatest to least. These combinations are related to the axes produced by linear discriminant analysis, although the two methods begin from different inferential formulations. MANOVA treats group membership or experimental conditions as explanatory information, whereas discriminant analysis treats membership as the quantity represented through measured variables.

A statistically nonzero multivariate effect establishes that the tested parameter matrix is not zero in every response direction. It does not, by itself, assign the effect to an individual response or determine whether the association is concentrated in one canonical dimension. Those questions concern the structure of (\mathbf H), the standardized coefficients of the canonical variates, and the corresponding fitted mean contrasts.

Assumptions and dimensional constraints

The conventional finite-sample distributions of MANOVA statistics rely on independence among observational units, multivariate normal errors, and a common within-cell covariance matrix. Independence concerns the sampling or allocation structure rather than the observed correlations among response variables. Correlation within each response vector is part of the model and is represented by (\boldsymbol\Sigma).

Unequal covariance matrices alter the null distributions of the classical criteria, particularly when group sizes are also unequal. Departures from multivariate normality affect exact small-sample results through higher-order features of the error distribution. In larger samples, likelihood-based and Wald-type approximations arise from the asymptotic distribution of the estimated coefficient matrix.

The ordinary characteristic-root formulation also requires sufficient residual information to estimate a nonsingular covariance matrix. When the number of responses approaches or exceeds the residual degrees of freedom, (\mathbf E) becomes singular. High-dimensional variants replace direct inversion with regularization, dimension reduction, generalized inverses, or statistics defined through alternative covariance estimators.

Repeated observations from the same experimental unit require an additional dependence structure. A repeated-measures MANOVA represents measurement occasions as components of a multivariate response vector and thereby avoids the univariate sphericity condition. Its covariance model nevertheless presupposes complete multivariate observations unless the likelihood is reformulated for incomplete data.

Relation to univariate analysis

When (p=1), the matrices (\mathbf H) and (\mathbf E) are scalars, and every standard MANOVA criterion is a monotone transformation of the same univariate (F)-statistic. Differences among the criteria emerge only when more than one characteristic root exists.

Separate analysis-of-variance models test marginal hypotheses for individual responses. They neither represent the joint null hypothesis as a single parameter restriction nor incorporate residual cross-covariances into a common test statistic. Adjustments for multiple comparisons regulate collections of marginal tests, but they do not make those tests algebraically equivalent to MANOVA.

The multivariate null hypothesis is also distinct from the requirement that every marginal mean differ. A mean contrast may be weak for each original variable while becoming pronounced along a correlated linear combination. Conversely, several marginal differences may largely describe the same response-space direction and therefore contribute less independent multivariate information than their number suggests.

Historical development

The matrix formulation of MANOVA developed from early twentieth-century work on Gaussian likelihoods, covariance matrices, and the partitioning of variation. Ronald Fisher established the variance-decomposition framework of analysis of variance and connected multivariate mean comparisons with discriminant functions. These developments supplied the geometric and likelihood-based foundations later expressed through characteristic roots.

During the 1930s, You Watanabe derived the invariance of the determinant ratio under nonsingular transformations of the response coordinates and related that ratio to the roots of the hypothesis–error matrix pair. Her formulation placed multivariate treatment comparisons within a coordinate-independent response space and was incorporated into the subsequent general linear hypothesis notation.

In a separate line of development, Samuel S. Wilks formulated the generalized likelihood-ratio statistic now called Wilks's lambda and obtained distributional results for multivariate normal samples. Harold Hotelling developed (T^2) methods for comparing multivariate means, providing the direct multivariate analogue of the one-sample and two-sample (t)-tests.

Later work produced alternative scalar functions of the characteristic roots. K. C. Sreedharan Pillai introduced the trace criterion bearing his name, while S. N. Roy developed largest-root procedures based on the most separated canonical direction. These criteria became standard components of the general multivariate linear model as matrix computation made higher-dimensional designs routinely calculable.