3.9 Several outcomes at once

When groups are compared on several correlated outcomes, testing each outcome alone can miss a difference that shows only in the combination. Multivariate tests compare the groups’ centres in all dimensions at once, taking the correlations into account.

Code: data and style used in this chapter

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from matplotlib.axes import Axes
from matplotlib.patches import Ellipse
from scipy import stats
from sklearn import datasets
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
from statsmodels.multivariate.manova import MANOVA

from statly import style

style.use()
rng = np.random.default_rng(1)

iris = datasets.load_iris(as_frame=True).frame  # 150 flowers, 3 species, 4 measurements
SPECIES = ["setosa", "versicolor", "virginica"]
SPECIES_COLORS = [style.BLUE, style.GREEN, style.CORAL]


def draw_ellipse(ax: Axes, points: np.ndarray, color: str) -> None:
    """An ellipse covering about 95% of a two-dimensional normal cloud fitted to the points."""
    values, vectors = np.linalg.eigh(np.cov(points.T))
    angle = np.degrees(np.arctan2(vectors[1, 1], vectors[0, 1]))
    width, height = 2 * np.sqrt(5.99 * values[::-1])
    ax.add_patch(
        Ellipse(points.mean(axis=0), width, height, angle=angle, fill=False, color=color, lw=2)
    )

3.9.1 Why not one test per outcome?

Two groups of patients are compared on two blood values that are strongly correlated in everyone. Testing each value with a t-test (Chapter 3.4) asks whether the groups differ on that value alone. But a difference can lie across the correlation: one group a little higher on the first value and a little lower on the second, a combination that is rare in either group. One-at-a-time tests then see small, unconvincing differences, while the two clouds are clearly separate.

Code

cov = np.array([[1, 0.95], [0.95, 1]])
group_a = rng.multivariate_normal([0, 0], cov, 40)
group_b = rng.multivariate_normal([0.3, -0.3], cov, 40)


def hotelling(a: np.ndarray, b: np.ndarray) -> tuple[float, float]:
    """Two-sample Hotelling T² and its p-value, through the F distribution."""
    n1, n2, p = len(a), len(b), a.shape[1]
    pooled = ((n1 - 1) * np.cov(a.T) + (n2 - 1) * np.cov(b.T)) / (n1 + n2 - 2)
    gap = a.mean(axis=0) - b.mean(axis=0)
    t2 = n1 * n2 / (n1 + n2) * gap @ np.linalg.solve(pooled, gap)
    f = (n1 + n2 - p - 1) / (p * (n1 + n2 - 2)) * t2
    return t2, stats.f.sf(f, p, n1 + n2 - p - 1)


t2, p_joint = hotelling(group_a, group_b)
fig, axes = style.subplots(3, height=3.2, width=10.5)
for points, color in ((group_a, style.BLUE), (group_b, style.CORAL)):
    axes[0].scatter(*points.T, s=10, color=color, alpha=0.7)
    draw_ellipse(axes[0], points, color)
style.note(axes[0], f"Hotelling: p = {p_joint:.1g}", x=0.03, y=0.95, color=style.INK)
axes[0].set(xlabel="outcome 1", ylabel="outcome 2", aspect="equal")
for ax, k in zip(axes[1:], (0, 1)):
    bins = np.linspace(-3, 3, 20)
    ax.hist(group_a[:, k], bins=bins, color=style.BLUE, alpha=0.5)
    ax.hist(group_b[:, k], bins=bins, color=style.CORAL, alpha=0.5)
    style.note(
        ax,
        f"t-test: p = {stats.ttest_ind(group_a[:, k], group_b[:, k]).pvalue:.2f}",
        x=0.03,
        y=0.95,
        color=style.INK,
    )
    ax.set(xlabel=f"outcome {k + 1}", ylabel="patients")
fig.tight_layout()
plt.show()

Figure 3.9.1: Two simulated groups of 40, measured on two strongly correlated outcomes. Left: the two clouds, each with an ellipse holding about 95% of its points; they lie side by side, overlapping only at the edges. Middle and right: each outcome on its own; the groups overlap heavily, and the t-tests are borderline at best. Hotelling’s test, using both outcomes and their correlation, finds the difference overwhelming.

Hotelling’s T² is the multivariate t-test. It measures the distance between the two mean vectors in units that account for the spread and correlation of the outcomes, the Mahalanobis distance:

$T^2 = \frac{n_1 n_2}{n_1 + n_2} (\bar x_1 - \bar x_2)^\top S^{-1} (\bar x_1 - \bar x_2)$

with S the pooled covariance matrix. A rescaled $T^2$ follows an F distribution under the null hypothesis. It also controls the false alarm rate across the outcomes in one test, instead of one test per outcome (Chapter 3.10).

3.9.2 More than two groups: MANOVA

Multivariate analysis of variance (MANOVA) extends one-way ANOVA (Chapter 3.5) to several outcomes. Instead of splitting one sum of squares, it splits a matrix of sums of squares and cross-products into a between-group part H and a within-group part E, and summarises how large H is relative to E. Four summaries are in common use; Pillai’s trace is the most robust, Wilks’ lambda the most traditional:

$\Lambda = \frac{\det E}{\det(H + E)}$

Λ near 0 means the groups account for almost all the variation; near 1, for almost none.

Code

data = iris.copy()
data.columns = [c.replace(" (cm)", "").replace(" ", "_") for c in iris.columns]
data["species"] = pd.Categorical.from_codes(iris.target, SPECIES)
result = MANOVA.from_formula(
    "sepal_length + sepal_width + petal_length + petal_width ~ species", data
).mv_test()
table = result.results["species"]["stat"]
print("| statistic | value | F | p-value |")
print("|:--|--:|--:|--:|")
for name, row in table.iterrows():
    print(f"| {name} | {row['Value']:.3f} | {row['F Value']:.1f} | {row['Pr > F']:.2g} |")

Table 3.9.1: MANOVA of the four flower measurements by species, for the 150 iris flowers. All four statistics agree that the species differ.

statistic	value	F	p-value
Wilks’ lambda	0.023	199.1	1.4e-112
Pillai’s trace	1.192	53.5	9.7e-53
Hotelling-Lawley trace	32.477	582.2	1.1e-135
Roy’s greatest root	32.192	1167.0	3.8e-109

3.9.3 Which combination separates the groups?

A significant MANOVA says the group centres differ somewhere in the space of outcomes. Where exactly is shown by the discriminant directions: weighted combinations of the outcomes that separate the groups as much as possible relative to the spread within groups. They are the directions of linear discriminant analysis (Chapter 6.1) and the natural follow-up to a MANOVA.

Code

lda = LinearDiscriminantAnalysis(n_components=2).fit(iris.iloc[:, :4], iris.target)
scores = lda.transform(iris.iloc[:, :4])

fig, ax = plt.subplots(figsize=(8, 3.8))
for k, (name, color) in enumerate(zip(SPECIES, SPECIES_COLORS)):
    group = scores[iris.target == k]
    ax.scatter(*group.T, s=12, color=color, alpha=0.7, label=name)
    draw_ellipse(ax, group, color)
loadings = lda.scalings_[:, :2]
for name, (x, y) in zip(iris.columns[:4], loadings):
    ax.annotate("", xy=(x, y), xytext=(0, 0), arrowprops=dict(arrowstyle="->", color=style.INK))
    ax.text(x * 1.08, y * 1.08, name.replace(" (cm)", ""), fontsize=8)
ax.legend(fontsize=8)
ax.set(xlabel="discriminant direction 1", ylabel="direction 2")
plt.show()

Figure 3.9.2: The iris flowers on the two discriminant directions found from all four measurements. The first direction separates the three species almost completely, setosa far from the other two; the second adds little. The arrows show how each measurement contributes to the two directions.

3.9.4 When the assumptions are doubtful

MANOVA and Hotelling’s test assume multivariate normal outcomes with the same covariance matrix in every group, and need more subjects than outcomes. When these fail, a permutation version (Chapter 3.1) shuffles the group labels and recomputes the statistic; with distances between subjects in place of means and covariances it becomes PERMANOVA, common in ecology and microbiome studies where there are many more outcomes than subjects.

Code

pooled_points = np.r_[group_a, group_b]
shuffled = []
for _ in range(5000):
    order = rng.permutation(80)
    shuffled.append(hotelling(pooled_points[order[:40]], pooled_points[order[40:]])[0])

fig, ax = plt.subplots(figsize=(8, 3))
style.hist(ax, shuffled, bins=60, color=style.GREY)
ax.axvline(t2, color=style.CORAL, lw=2)
style.note(ax, f"observed T² = {t2:.0f}", x=0.6, y=0.9, color=style.CORAL)
ax.set(xlabel="T² with shuffled labels", ylabel="shuffles")
plt.show()

Figure 3.9.3: Hotelling’s T² for the two groups of Figure 1.3.9.1 (red line) against 5,000 values of T² after shuffling the group labels. None of the shuffles comes close: the permutation test agrees with the F-based p-value without assuming normality.

3.9.5 In one picture

Several correlated outcomes: one-at-a-time tests can miss differences across the correlation.
Hotelling’s T²: the multivariate t-test, a Mahalanobis distance between means.
MANOVA: between versus within matrices; Pillai’s trace, Wilks’ Λ = det E / det(H + E).
Follow up with discriminant directions to see which combination separates the groups.
Permutation versions and PERMANOVA relax normality and work with many outcomes.

Next: errors, power, and what happens when we run many tests.