Multivariate statistics reference

Each section starts with the big idea in a few sentences, then gives the results you will look up most often, with the conditions they need.

Bold lowercase letters are vectors and bold capitals are matrices. A bold capital such as $\mathbf X$ can also be a random vector, and $\boldsymbol\rho$ is the correlation matrix. A prime means transpose, $\lvert \mathbf A \rvert$ is a determinant and $\operatorname{tr}(\mathbf A)$ a trace. Critical values are upper percentiles, so $F_{a,b}(\alpha)$, $t_\nu(\alpha)$ and $\chi^2_\nu(\alpha)$ each leave probability $\alpha$ to their right. Unless a section says otherwise, $p$ is the number of variables and $n$ the number of observations.

Describing multivariate data

With several variables measured on each item, one mean and one variance are no longer enough. The data are summarized by a mean vector and a covariance matrix, and most of the course is about what these two objects can tell you.

The data matrix and its summaries

Measuring $p$ variables on $n$ items gives an $n \times p$ data matrix $\mathbf X$. Row $j$ is the observation $\mathbf x_j'$, and $x_{jk}$ is the value of variable $k$ on item $j$.

The divisor $n - 1$ makes $\mathbf S$ unbiased. The generalized variance is zero exactly when some combination of the variables has no spread in the sample, which always happens when $n \le p$.

Population quantities

For a random vector $\mathbf X = (X_1, \dots, X_p)'$, the mean vector is $\boldsymbol\mu = E(\mathbf X)$ and the covariance matrix is

with $\sigma_{ik} = \operatorname{Cov}(X_i, X_k)$ off the diagonal and the variances $\sigma_{ii}$ on it. The correlation matrix $\boldsymbol\rho$ has entries $\rho_{ik} = \sigma_{ik}/\sqrt{\sigma_{ii}\sigma_{kk}}$. The covariance matrix is symmetric and positive semidefinite, and positive definite unless some nontrivial combination of the variables has zero variance.

Linear combinations

Most methods in the course work by forming combinations of the variables, so these rules come up everywhere. For constant vectors $\mathbf b$ and $\mathbf c$ and a constant matrix $\mathbf C$,

See also Linear algebra refresher for eigenvalues, quadratic forms and positive definite matrices

The multivariate normal

The multivariate normal plays the role the bell curve plays for one variable. It is described completely by $\boldsymbol\mu$ and $\boldsymbol\Sigma$, its contours are ellipsoids, and almost every exact test in the course assumes it.

Density and contours

A random vector $\mathbf X \sim N_p(\boldsymbol\mu, \boldsymbol\Sigma)$ with $\boldsymbol\Sigma$ positive definite has density

The density is constant on the ellipsoids $(\mathbf x - \boldsymbol\mu)'\boldsymbol\Sigma^{-1}(\mathbf x - \boldsymbol\mu) = c^2$. If $\boldsymbol\Sigma\mathbf e_i = \lambda_i\mathbf e_i$ with unit eigenvectors $\mathbf e_i$, the axes of the ellipsoid point along the $\mathbf e_i$ with half lengths $c\sqrt{\lambda_i}$. With $c^2 = \chi^2_p(\alpha)$ the ellipsoid holds probability $1 - \alpha$.

See also Covariance explorer to see how the ellipse changes with the covariance matrix in two dimensions

Key properties

If $\mathbf X \sim N_p(\boldsymbol\mu, \boldsymbol\Sigma)$, then

• every linear combination is normal, with $\mathbf a'\mathbf X \sim N(\mathbf a'\boldsymbol\mu,\ \mathbf a'\boldsymbol\Sigma\mathbf a)$, and more generally $\mathbf A\mathbf X + \mathbf d \sim N_q(\mathbf A\boldsymbol\mu + \mathbf d,\ \mathbf A\boldsymbol\Sigma\mathbf A')$ for a constant $q \times p$ matrix $\mathbf A$
• every subset of the variables is normal
• two blocks of variables are independent exactly when their cross covariance is $\mathbf 0$, which relies on joint normality
• when $\boldsymbol\Sigma$ is positive definite, the squared distance $(\mathbf X - \boldsymbol\mu)'\boldsymbol\Sigma^{-1}(\mathbf X - \boldsymbol\mu)$ has a $\chi^2_p$ distribution

Conditional distributions. Split $\mathbf X$ into blocks $\mathbf X^{(1)}$ and $\mathbf X^{(2)}$, with $\boldsymbol\mu$ and $\boldsymbol\Sigma$ split to match, and suppose $\boldsymbol\Sigma_{22}$ is invertible. Given $\mathbf X^{(2)} = \mathbf x^{(2)}$, the block $\mathbf X^{(1)}$ is normal with

The conditional mean is linear in $\mathbf x^{(2)}$, which is the link to regression, and the conditional covariance does not depend on $\mathbf x^{(2)}$ at all.

Mahalanobis distance and checking normality

The Mahalanobis distance measures how far $\mathbf x$ is from $\boldsymbol\mu$ in standard deviation units, taking the correlations into account,

Points on the same density contour are at the same distance. To check a sample for normality, compute $d_j^2 = (\mathbf x_j - \bar{\mathbf x})'\mathbf S^{-1}(\mathbf x_j - \bar{\mathbf x})$ for every observation and plot the ordered values against $\chi^2_p$ quantiles. For normal data the points lie close to a straight line.

Sampling distributions

For a random sample of size $n$ from $N_p(\boldsymbol\mu, \boldsymbol\Sigma)$,

Without normality, but with $\boldsymbol\Sigma$ positive definite and $n - p$ large, $\sqrt n\,(\bar{\mathbf X} - \boldsymbol\mu)$ is approximately $N_p(\mathbf 0, \boldsymbol\Sigma)$ and $n(\bar{\mathbf X} - \boldsymbol\mu)'\mathbf S^{-1}(\bar{\mathbf X} - \boldsymbol\mu)$ is approximately $\chi^2_p$. This is what makes the large sample tests below work.

Inference about mean vectors

Testing several means at once is not the same as running a separate $t$ test for each variable. Hotelling’s $T^2$ does it in one step, accounts for the correlations, and keeps the overall error rate at the level you choose.

Hotelling’s T² for one mean vector

To test $H_0$ that $\boldsymbol\mu = \boldsymbol\mu_0$, with a random sample from $N_p(\boldsymbol\mu, \boldsymbol\Sigma)$ and $n > p$, compute

and reject at level $\alpha$ when

For one variable, $T^2$ is the square of the usual $t$ statistic. The test gives the same answer whatever units the variables are measured in, and it agrees with the likelihood ratio test. The $100(1 - \alpha)\%$ confidence region for $\boldsymbol\mu$ is the ellipsoid of all $\boldsymbol\mu$ that this test would not reject.

Simultaneous confidence intervals

After a test, you usually want intervals for individual means or combinations $\mathbf a'\boldsymbol\mu$ that all hold together. Each interval below has the form $\mathbf a'\bar{\mathbf x} \pm (\text{multiplier}) \times \sqrt{\mathbf a'\mathbf S\mathbf a / n}$.

For a single mean $\mu_k$, take $\mathbf a$ with a one in place $k$ and zeros elsewhere, so the interval is centred at $\bar x_k$ with $\sqrt{s_{kk}/n}$ in place of $\sqrt{\mathbf a'\mathbf S\mathbf a/n}$.

Paired and repeated measures

When each unit receives both treatments, work with the differences $\mathbf d_j = \mathbf x_{1j} - \mathbf x_{2j}$ and run the one sample test of $\boldsymbol\delta = \mathbf 0$ on them, using $\bar{\mathbf d}$ and $\mathbf S_d$.

When each unit receives $q$ treatments, choose a $(q-1) \times q$ matrix $\mathbf C$ of linearly independent contrasts and test $H_0$ that $\mathbf C\boldsymbol\mu = \mathbf 0$ with

rejecting when

Comparing two mean vectors

Take independent samples of sizes $n_1$ and $n_2$ from $N_p(\boldsymbol\mu_1, \boldsymbol\Sigma)$ and $N_p(\boldsymbol\mu_2, \boldsymbol\Sigma)$, so both populations share one covariance matrix. Pool the two estimates,

and test $H_0$ that $\boldsymbol\mu_1 = \boldsymbol\mu_2$ with

Reject when

If the covariance matrices differ but $n_1 - p$ and $n_2 - p$ are both large, replace the bracket by $\frac{1}{n_1}\mathbf S_1 + \frac{1}{n_2}\mathbf S_2$ and compare $T^2$ with $\chi^2_p(\alpha)$.

MANOVA

MANOVA extends the two sample comparison to $g$ groups, in the same way that analysis of variance extends the $t$ test. It splits the total variation into a part between the groups and a part within them, and asks whether the first is large compared with the second.

Model and the two matrices

Group $\ell$ gives a sample of size $n_\ell$, with $n = n_1 + \cdots + n_g$. The observations are independent and normal with group means $\boldsymbol\mu_\ell$ and a common covariance matrix $\boldsymbol\Sigma$. The null hypothesis $H_0$ is that all group means are equal.

Here $\bar{\mathbf x}_\ell$ is the mean of group $\ell$ and $\bar{\mathbf x}$ the mean of all $n$ observations.

Wilks’ lambda and the other test statistics

The usual statistic is Wilks’ lambda,

which needs $n - g \ge p$ so that $\mathbf W$ is invertible. It lies between 0 and 1, and small values mean the group means are far apart compared with the spread inside the groups. For large $n$, reject $H_0$ when

Software also reports Pillai’s trace, the Lawley-Hotelling trace and Roy’s largest root. All four are built from the eigenvalues of $\mathbf W^{-1}\mathbf B$, and in large samples they usually agree. Some small cases have exact $F$ distributions under $H_0$.

Which means differ

After a rejection, Bonferroni intervals for every variable $i$ and every pair of groups $k$ and $\ell$ are

where $w_{ii}$ is the $i$th diagonal entry of $\mathbf W$. All of them hold together with probability at least $1 - \alpha$.

Principal components

Principal components replace the original variables with new uncorrelated ones, ordered from most to least variance. When a few of them carry most of the variance, they give a simpler description of the data with little lost.

The components

Let $\boldsymbol\Sigma$ have eigenvalues $\lambda_1 \ge \lambda_2 \ge \cdots \ge \lambda_p \ge 0$ with orthonormal eigenvectors $\mathbf e_1, \dots, \mathbf e_p$. The $i$th principal component is $Y_i = \mathbf e_i'\mathbf X$.

The first component has the largest variance of all combinations $\mathbf a'\mathbf X$ with $\mathbf a'\mathbf a = 1$, and each later one has the largest variance among those uncorrelated with the earlier ones. For a multivariate normal, the eigenvectors are the axes of the density ellipsoids.

Standardized variables and sample components

When the variables are in different units or have very different variances, standardize first and use the eigenpairs of the correlation matrix $\boldsymbol\rho$. The total variance is then $p$, so component $k$ explains the share $\lambda_k/p$. Components from $\boldsymbol\Sigma$ and from $\boldsymbol\rho$ differ in general.

In a sample, use the eigenpairs $(\hat\lambda_i, \hat{\mathbf e}_i)$ of $\mathbf S$, or of $\mathbf R$ for standardized data. The score of observation $j$ on component $i$ is $\hat{\mathbf e}_i'(\mathbf x_j - \bar{\mathbf x})$. Scores on component $i$ have sample variance $\hat\lambda_i$, and scores on different components are uncorrelated.

How many to keep

• Cumulative share. Keep enough components to explain a chosen share of the total variance, often around 80 to 90 percent.
• Scree plot. Plot $\hat\lambda_i$ against $i$ and stop at the elbow, where the remaining eigenvalues are small and similar.
• Eigenvalues above 1. With $\mathbf R$, a common rule of thumb keeps the components that explain more than one standardized variable, those with $\hat\lambda_i > 1$.
• Interpretation. A component you cannot interpret or use is rarely worth keeping.

See also Linear algebra refresher for eigenvalues and the spectral decomposition

Factor analysis

Factor analysis explains the correlations among many observed variables through a few unobserved common factors. Where principal components summarize variance, factor analysis models covariance.

The model

With $m < p$ common factors,

where $\mathbf L$ is the $p \times m$ matrix of loadings, $\mathbf F$ holds the common factors and $\boldsymbol\varepsilon$ the specific factors. The factors have mean $\mathbf 0$ with $\operatorname{Cov}(\mathbf F) = \mathbf I$, the specific factors have the diagonal covariance matrix $\boldsymbol\Psi$, and $\mathbf F$ and $\boldsymbol\varepsilon$ are uncorrelated. The model implies

So each variance splits as $\sigma_{ii} = h_i^2 + \psi_i$. The communality $h_i^2 = \ell_{i1}^2 + \cdots + \ell_{im}^2$ is the part explained by the common factors, and the specific variance $\psi_i$ is the rest. The loading $\ell_{ij}$ is also the covariance of $X_i$ with factor $j$.

Estimation

Principal component method. With $(\hat\lambda_i, \hat{\mathbf e}_i)$ the eigenpairs of $\mathbf S$ or $\mathbf R$, the loadings for factor $j$ are $\sqrt{\hat\lambda_j}\,\hat{\mathbf e}_j$, and each specific variance is the diagonal entry of $\mathbf S$ minus the communality. Factor $j$ explains the share $\hat\lambda_j / (s_{11} + \cdots + s_{pp})$ of the total sample variance, or $\hat\lambda_j / p$ with $\mathbf R$.

Maximum likelihood. Assuming normality, maximize the likelihood over $\mathbf L$ and $\boldsymbol\Psi$. This also gives a large sample test of whether $m$ factors are enough. Reject at level $\alpha$ when

where $\mathbf S_n = \frac{n-1}{n}\mathbf S$. The test needs $n$ and $n - p$ large and $\nu > 0$.

Rotation

Loadings are not unique. For any orthogonal $m \times m$ matrix $\mathbf T$, the rotated loadings $\mathbf L\mathbf T$ reproduce the same $\boldsymbol\Sigma$, the same communalities and the same $\boldsymbol\Psi$. Rotation uses this freedom to make the factors easier to read. Varimax, the most common choice, pushes each loading toward 0 or toward a large value, so each factor ends up with a few clearly large loadings.

Classification

Classification assigns a new item to one of several known groups from its measurements. The best rule weighs how likely the measurements are in each group, how common each group is, and how costly each kind of mistake is.

The optimal rule

For two populations $\pi_1$ and $\pi_2$ with densities $f_1$ and $f_2$, prior probabilities $p_1$ and $p_2$, and misclassification costs $c(1 \mid 2)$ and $c(2 \mid 1)$, the rule with the smallest expected cost assigns $\mathbf x_0$ to $\pi_1$ when

and to $\pi_2$ otherwise. Here $c(1 \mid 2)$ is the cost of putting an item from $\pi_2$ into $\pi_1$. With equal costs and priors, $\mathbf x_0$ simply goes to the population with the larger density there. With $g$ populations and equal costs, choose the population with the largest $p_k f_k(\mathbf x)$, which is the largest posterior probability.

Normal populations

If the populations are $N_p(\boldsymbol\mu_i, \boldsymbol\Sigma_i)$, the rule takes one of two forms. With equal costs, assign $\mathbf x$ to the population with the largest score.

In practice replace $\boldsymbol\mu_i$ by $\bar{\mathbf x}_i$, and $\boldsymbol\Sigma$ by $\mathbf S_{\text{pooled}}$ or $\boldsymbol\Sigma_i$ by $\mathbf S_i$. The quadratic rule needs each $\mathbf S_i$ to be invertible, so $n_i > p$, and it is more sensitive to departures from normality.

Fisher’s linear discriminant

Fisher looked for the combination $y = \mathbf a'\mathbf x$ that separates two samples best relative to the spread within them, without assuming normality. The best choice is

and a new item goes to $\pi_1$ when $\hat{\mathbf a}'\mathbf x_0$ is at least the midpoint of $\hat{\mathbf a}'\bar{\mathbf x}_1$ and $\hat{\mathbf a}'\bar{\mathbf x}_2$. This is the same rule as the normal rule with equal covariances, costs and priors. The separation it achieves is the squared Mahalanobis distance between the sample means,

Judging a rule

The apparent error rate, the share of the training items the rule gets wrong, is too optimistic because the rule is judged on the data that built it. Leave-one-out cross-validation is fairer. Hold out one item, build the rule from the rest, classify the held-out item, and repeat for every item. The share misclassified this way is a nearly unbiased estimate of the error rate on new data.

Canonical correlation

Canonical correlation measures the association between two sets of variables, such as test scores and later performance. It finds the combinations of the two sets that are most strongly correlated with each other, then the next best pair, and so on.

Canonical variates

Split the variables into $\mathbf X^{(1)}$ with $p$ entries and $\mathbf X^{(2)}$ with $q$ entries, $p \le q$, with covariance blocks $\boldsymbol\Sigma_{11}$, $\boldsymbol\Sigma_{22}$ positive definite and cross covariance $\boldsymbol\Sigma_{12}$. The first canonical pair $U_1 = \mathbf a'\mathbf X^{(1)}$ and $V_1 = \mathbf b'\mathbf X^{(2)}$ maximizes

with both variates scaled to variance 1. Each later pair maximizes it among combinations uncorrelated with the earlier pairs. The maximized values $\rho_1^* \ge \cdots \ge \rho_p^*$ are the canonical correlations, and their squares are the eigenvalues of

The canonical correlations do not change when either set is replaced by nonsingular linear combinations of itself, so standardizing the variables leaves them the same. The first one is at least as large as the absolute correlation between any single variable in one set and any single variable in the other.

Testing

Compute $\hat\rho_1^* \ge \cdots \ge \hat\rho_p^*$ from the sample blocks of $\mathbf S$ or $\mathbf R$. For normal data and large $n$, reject $H_0$ that all canonical correlations are zero, meaning $\boldsymbol\Sigma_{12} = \mathbf 0$, when

To test whether only the first $k$ are nonzero, start the product at $i = k + 1$ and use $\chi^2_{(p-k)(q-k)}(\alpha)$.

Cluster analysis

Cluster analysis looks for groups in data when no groups are given in advance. The answer depends on how distance is measured and on how clusters are joined, so it is worth trying more than one choice.

Distances between items

Hierarchical clustering

Agglomerative clustering starts with every item as its own cluster and repeatedly merges the two closest clusters until one remains. A dendrogram records which clusters merged and at what distance, and cutting it at a chosen height gives the clusters. The linkage decides the distance between clusters $A$ and $B$.

For Ward’s method, merging clusters $A$ and $B$ with $N_A$ and $N_B$ items raises the within cluster sum of squares by

k-means

1. Choose the number of clusters $K$ and a starting set of $K$ centroids or groups.
2. Assign each item to the cluster with the nearest centroid.
3. Recompute each centroid as the mean of its items.
4. Repeat steps 2 and 3 until no item changes cluster.

With Euclidean distance neither step can increase the within cluster sum of squares, so the algorithm always stops, but possibly at a local minimum. Run it from several starts, standardize variables in different units, and compare a few values of $K$.