7.1 Principal components and latent variables

Many measurements often move together because a few underlying directions drive them. Principal components find the directions in which the data vary most; a handful of them can carry most of the information in dozens of measurements.

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from matplotlib.axes import Axes
from sklearn import datasets
from sklearn.datasets import make_swiss_roll
from sklearn.decomposition import PCA, FactorAnalysis
from sklearn.manifold import MDS, TSNE, Isomap
from sklearn.preprocessing import StandardScaler

from statly import style

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

wine = datasets.load_wine(as_frame=True)  # 178 Italian wines, 13 chemical measurements
measurements = wine.data
NAMES = [
    c.replace("_", " ").replace("od280/od315 of diluted wines", "od280")
    for c in measurements.columns
]
cultivar = wine.target.to_numpy()  # three grape cultivars, used only to check results
CULTIVAR_COLORS = np.array([style.BLUE, style.GREEN, style.CORAL])
scaled = StandardScaler().fit_transform(measurements)
pca = PCA().fit(scaled)
scores = pca.transform(scaled)
digits = datasets.load_digits()  # 1,797 handwritten digits, 8 × 8 pixels each

7.1.1 The idea

The 178 wines from Chapter 3.5 have 13 chemical measurements each. Thirteen numbers per wine cannot be drawn, but many of them move together (Chapter 1.2 showed their correlation matrix). Principal component analysis (PCA) finds the direction in which the wines differ most, then the direction at right angles to it that captures most of what is left, and so on. With two measurements it is easy to see:

pair = StandardScaler().fit_transform(measurements[["total_phenols", "flavanoids"]])
pair_pca = PCA().fit(pair)


def plot_pair(ax: Axes) -> None:
    """The wines, the two principal directions, and projections onto the first."""
    first = pair_pca.components_[0]
    projected = np.outer(pair @ first, first)
    some = rng.choice(len(pair), 40, replace=False)
    for a, b in zip(pair[some], projected[some]):
        ax.plot([a[0], b[0]], [a[1], b[1]], color=style.GREY, lw=0.7)
    ax.scatter(pair[:, 0], pair[:, 1], s=12, color=style.BLUE, alpha=0.6)
    for direction, length, color in zip(
        pair_pca.components_,
        np.sqrt(pair_pca.explained_variance_) * 2,
        (style.CORAL, style.GOLD),
    ):
        ax.annotate(
            "",
            xy=direction * length,
            xytext=(0, 0),
            arrowprops=dict(arrowstyle="->", color=color, lw=2.2),
        )
    share = pair_pca.explained_variance_ratio_[0]
    style.note(ax, f"first component:\n{share:.0%} of the variation", x=0.03, color=style.CORAL)
    ax.set(
        xlabel="total phenols (standardised)",
        ylabel="flavanoids (standardised)",
        aspect="equal",
    )


fig, ax = plt.subplots(figsize=(5.4, 4.4))
plot_pair(ax)
plt.show()

Figure 7.1.1: Two measurements of the wines that move together (both standardised). The long arrow is the first principal component, the direction of largest spread; the short one is the second. Projecting each wine onto the first arrow (grey lines) keeps most of the information in one number, its score.

7.1.2 In symbols

With p standardised measurements $x_1, \dots, x_p$, the first principal component is the weighted sum

$$z_1 = \phi_{11} x_1 + \phi_{21} x_2 + \dots + \phi_{p1} x_p, \qquad \text{with } \sum_j \phi_{j1}^2 = 1$$

whose weights make the variance of $z_1$ across the cases as large as possible. The weights $\phi_{j1}$ are the loadings: what the component is made of. The values $z_{i1}$ for each case are its scores: where the case sits along the component. The second component has the largest variance among all directions at right angles to the first, and so on.

There is a neat shortcut: the loadings are the eigenvectors of the correlation matrix (the covariance matrix of the standardised measurements, Chapter 1.2), and the variance of each component is its eigenvalue $\lambda_m$. The share of the total variation that component m carries, the proportion of variance explained, is

$$\text{PVE}_m = \frac{\lambda_m}{\sum_{k=1}^{p} \lambda_k}$$

These are the arrows of the ellipses in Chapter 2.4: principal components are the axes of the data cloud.

7.1.3 Reading the components

def plot_scree(ax: Axes) -> None:
    """Share of variation per component, and the running total."""
    shares = pca.explained_variance_ratio_
    ax.bar(range(1, 14), shares, color=style.BLUE, alpha=0.6)
    ax.plot(range(1, 14), np.cumsum(shares), "o-", ms=3, color=style.CORAL, lw=1.2)
    ax.yaxis.set_major_formatter(style.PERCENT)
    ax.set(
        title="Variation",
        xlabel="component",
        ylabel="share of variation",
        xticks=range(1, 14, 2),
    )


def plot_biplot(ax: Axes, shown: int = 9) -> None:
    """Scores of the wines and loading arrows of the strongest measurements."""
    ax.scatter(scores[:, 0], scores[:, 1], c=CULTIVAR_COLORS[cultivar], s=10, alpha=0.5)
    loadings = pca.components_[:2].T * 4  # stretched to the scale of the scores
    strongest = np.argsort(-np.hypot(loadings[:, 0], loadings[:, 1]))[:shown]
    for i in strongest:
        ax.annotate(
            "",
            xy=loadings[i],
            xytext=(0, 0),
            arrowprops=dict(arrowstyle="->", color=style.INK, lw=1),
        )
        ax.text(*(loadings[i] * 1.12), NAMES[i], fontsize=7.5, ha="center", color=style.INK)
    ax.set(
        title="Biplot",
        xlabel=f"component 1 ({pca.explained_variance_ratio_[0]:.0%})",
        ylabel=f"component 2 ({pca.explained_variance_ratio_[1]:.0%})",
    )


fig = plt.figure(figsize=(10.5, 4))
grid = fig.add_gridspec(1, 2, width_ratios=[1, 1.5])
plot_scree(fig.add_subplot(grid[0]))
plot_biplot(fig.add_subplot(grid[1]))
fig.tight_layout()
plt.show()

Figure 7.1.2: Principal components of all 13 standardised measurements. Left: the share of variation each component carries (bars) and the running total (line); the first two carry over half. Right: a biplot. Every wine is a dot at its scores on the first two components, coloured by cultivar, which the method never saw: the cultivars separate. The arrows show the loadings of the measurements: an arrow pointing towards a group of wines means those wines are high on that measurement.

loading_table = pd.DataFrame(
    pca.components_[:2].T, index=NAMES, columns=["component 1", "component 2"]
).round(2)
print(loading_table.sort_values("component 1").to_markdown(colalign=["left", "right", "right"]))

Table 7.1.1: Loadings of the first two components: the weight of each standardised measurement. Sorted by the first component.

	component 1	component 2
nonflavanoid phenols	-0.3	0.03
malic acid	-0.25	0.22
alcalinity of ash	-0.24	-0.01
color intensity	-0.09	0.53
ash	-0	0.32
alcohol	0.14	0.48
magnesium	0.14	0.3
proline	0.29	0.36
hue	0.3	-0.28
proanthocyanins	0.31	0.04
od280	0.38	-0.16
total phenols	0.39	0.07
flavanoids	0.42	-0

To read a component, look at its largest loadings and their signs. The first component sets flavanoids, total phenols, od280 (how much light the diluted wine absorbs), proanthocyanins and hue against nonflavanoid phenols, alcalinity of ash and malic acid: roughly “phenolic richness”. The second is dominated by colour intensity, alcohol, proline and ash: roughly “body”. The signs themselves are arbitrary: flipping all loadings of a component flips its scores and changes nothing. Two numbers per wine instead of 13, and the three cultivars fall into three clouds: the typical use of PCA is a map of high-dimensional data, to look at before modelling anything.

7.1.4 How many components?

There is no single answer. Common guides: keep components until the running total passes a level such as 80%; look for an elbow in the scree plot; if the components feed a regression (PCR, Chapter 6.4), choose by cross-validation. A useful check is parallel analysis: shuffle each measurement separately, which destroys all correlation, and see which eigenvalues of the real data beat those of the shuffled data.

shuffled_eigenvalues = []
for _ in range(50):
    shuffled = np.column_stack([rng.permutation(scaled[:, j]) for j in range(scaled.shape[1])])
    shuffled_eigenvalues.append(PCA().fit(shuffled).explained_variance_)
shuffled_eigenvalues = np.array(shuffled_eigenvalues)

fig, ax = plt.subplots(figsize=(8, 3.1))
components = np.arange(1, 14)
ax.fill_between(
    components,
    shuffled_eigenvalues.min(axis=0),
    shuffled_eigenvalues.max(axis=0),
    color=style.GREY,
    alpha=0.3,
    lw=0,
    label="shuffled data",
)
ax.plot(components, pca.explained_variance_, "o-", color=style.BLUE, label="real data")
ax.legend(fontsize=8)
ax.set(xlabel="component", ylabel="eigenvalue (variance)", xticks=components)
plt.show()

Figure 7.1.3: Eigenvalues of the 13 wine components (blue) against the eigenvalues of 50 data sets made by shuffling each measurement on its own (grey band: their range). The first three components carry more variation than correlation-free data would; the rest do not.

7.1.5 Scale first

Principal components chase variation, so a measurement in large units dominates. Proline is measured in the hundreds and thousands, most others below 10:

raw_scores = PCA(2).fit_transform(measurements)

fig, axes = style.subplots(2, height=3.2, width=9.5)
for ax, values, title in zip(axes, (raw_scores, scores), ("Raw measurements", "Standardised")):
    ax.scatter(values[:, 0], values[:, 1], c=CULTIVAR_COLORS[cultivar], s=12, alpha=0.7)
    ax.set(title=title, xlabel="component 1", ylabel="component 2")
fig.tight_layout()
plt.show()

Figure 7.1.4: The first two components of the wines without standardising (left) and with (right). Unscaled, the first component is nearly just proline, the measurement with the largest numbers, and two of the cultivars mix.

Unless the measurements share a unit and their spread matters, standardise them first, as for the penalties of Chapter 6.4 and the distances of k nearest neighbours.

7.1.6 Components as compression

PCA is also the best way to approximate the data with fewer numbers: among all ways to describe each case by M numbers that combine linearly back into the measurements, the first M principal components lose the least (in squared error). The digits make it visible: each 8 × 8 image is 64 numbers, and a few components rebuild it surprisingly well.

digit_pca = PCA().fit(digits.data)
COUNTS = [1, 3, 10, 30]
examples = digits.data[[0, 7, 22, 42]]


def rebuild(images: np.ndarray, m: int) -> np.ndarray:
    """Images rebuilt from their first m principal components."""
    z = (images - digit_pca.mean_) @ digit_pca.components_[:m].T
    return digit_pca.mean_ + z @ digit_pca.components_[:m]


fig = plt.figure(figsize=(9, 5.4))
grid = fig.add_gridspec(5, 5, height_ratios=[1, 1, 1, 1, 2.6])
for row, image in enumerate(examples):
    versions = [rebuild(image[None], m)[0] for m in COUNTS] + [image]
    for col, version in enumerate(versions):
        ax = fig.add_subplot(grid[row, col])
        ax.imshow(version.reshape(8, 8), cmap="gray_r", vmin=0, vmax=16)
        ax.set_axis_off()
        if row == 0:
            ax.set_title(f"{COUNTS[col]} comp." if col < 4 else "original", fontsize=9)
curve = fig.add_subplot(grid[4, :])
curve.plot(range(1, 65), np.cumsum(digit_pca.explained_variance_ratio_), color=style.CORAL)
curve.yaxis.set_major_formatter(style.PERCENT)
curve.set(xlabel="components kept (M)", ylabel="variation kept", ylim=(0, 1.02))
fig.tight_layout()
plt.show()

Figure 7.1.5: Top: four digits rebuilt from their first 1, 3, 10 and 30 principal components, and the originals (right). Bottom: the share of the variation in all 1,797 images kept by the first M components; 10 components keep about three quarters of 64 pixels’ worth.

7.1.7 Filling in missing values

The compression view leads to a practical trick. When some entries of a table are missing, the patterns among the measurements can fill them in: a missing value is predicted from the components, which are estimated from what is there. The simple version iterates: fill the gaps with column means, fit M components, replace the gaps with the components’ reconstruction, refit, and repeat until the filled values settle. This is matrix completion, the idea behind recommendation systems that guess a viewer’s rating of films they have not seen.

missing = rng.uniform(size=scaled.shape) < 0.2
holey = np.where(missing, np.nan, scaled)


def complete(table: np.ndarray, m: int = 3, rounds: int = 100) -> np.ndarray:
    """Fill the gaps of a table by repeatedly reconstructing it from m principal components."""
    filled = np.where(np.isnan(table), np.nanmean(table, axis=0), table)
    for _ in range(rounds):
        model = PCA(m).fit(filled)
        rebuilt = model.inverse_transform(model.transform(filled))
        filled = np.where(np.isnan(table), rebuilt, table)
    return filled


completed = complete(holey)
truth, guess = scaled[missing], completed[missing]
mean_guess = np.nanmean(holey, axis=0)[np.where(missing)[1]]


def plot_guesses(ax: Axes, guesses: np.ndarray, title: str, color: str) -> None:
    """Filled-in values against the deleted true values."""
    ax.scatter(truth, guesses, s=6, color=color, alpha=0.4)
    ax.plot([-3, 4], [-3, 4], color=style.INK, ls="--", lw=1)
    r = np.corrcoef(truth, guesses)[0, 1] if guesses.std() > 0 else 0.0
    error = np.sqrt(np.mean((truth - guesses) ** 2))
    style.note(ax, f"error {error:.2f}\ncorrelation {r:.2f}", x=0.03, y=0.98, color=style.INK)
    ax.set(
        title=title,
        xlabel="true value (standardised)",
        ylabel="filled-in value",
        xlim=(-3, 4),
        ylim=(-3, 4),
        aspect="equal",
    )


fig, axes = style.subplots(2, height=3.6, width=8)
plot_guesses(axes[0], mean_guess, "Column means", style.GREY)
plot_guesses(axes[1], guess, "Iterative PCA", style.BLUE)
fig.tight_layout()
plt.show()

Figure 7.1.6: One in five values of the standardised wine table deleted at random, then filled in. Left: filling each gap with its column’s mean (all at 0) ignores everything else known about the wine. Right: iterative PCA with 3 components uses the other measurements, and its guesses follow the true values much better, though far from perfectly.

7.1.8 Latent variables: factor analysis

Chapter 2.4 met latent variables: hidden quantities that cause several measured ones to move together. Factor analysis is the model built on that idea. Each measurement is a weighted sum of a few hidden factors plus its own noise:

$$x_j = \lambda_{j1} f_1 + \dots + \lambda_{jk} f_k + \varepsilon_j$$

PCA and factor analysis look similar but answer different questions. PCA describes the variation, all of it, with no model behind it. Factor analysis explains the correlations by a few causes and leaves each measurement its own private noise. Its factors are only defined up to a rotation, so they are usually rotated (varimax is the common choice) to make each measurement load mainly on one factor, which makes them easier to name.

factors = FactorAnalysis(3, rotation="varimax", random_state=0).fit(scaled)


def plot_loading_map(ax: Axes, loadings: np.ndarray, title: str, labels: list[str]) -> None:
    """A heat map of loadings, measurements by row."""
    limit = abs(loadings).max()
    ax.imshow(loadings, cmap="RdBu_r", vmin=-limit, vmax=limit, aspect="auto")
    for i in range(loadings.shape[0]):
        for j in range(loadings.shape[1]):
            ink = "white" if abs(loadings[i, j]) > 0.6 * limit else "black"
            ax.text(
                j,
                i,
                f"{loadings[i, j]:.1f}".replace("-0.0", "0.0"),
                ha="center",
                va="center",
                fontsize=7.5,
                color=ink,
            )
    ax.set_xticks(range(loadings.shape[1]), labels)
    ax.set_yticks(range(len(NAMES)), NAMES, fontsize=8)
    ax.grid(False)
    ax.set_title(title)


fig, axes = style.subplots(2, height=4.6, width=8.5, sharey=True)
plot_loading_map(
    axes[0], pca.components_[:3].T, "Principal components", ["PC 1", "PC 2", "PC 3"]
)
plot_loading_map(axes[1], factors.components_.T, "Rotated factors", ["F 1", "F 2", "F 3"])
fig.tight_layout()
plt.show()

Figure 7.1.7: Loadings of the 13 wine measurements on three principal components (left) and on three factors after varimax rotation (right). Red is a positive weight, blue negative. The components mix many measurements; the rotated factors are cleaner, each led by a few measurements: the phenols, colour intensity and alcohol (with proline), and ash with alcalinity of ash.

Naming factors is interpretation, not discovery: the data say which measurements move together, and it is up to us to decide what, if anything, that common cause is.

7.1.9 Nonlinear maps: MDS, Isomap, t-SNE and UMAP

Principal components can only rotate and flatten. When the data lie on a curved surface, or form groups that no flat view separates, nonlinear maps can do better. They all place each case as a point in two dimensions, but they try to keep different things:

• Multidimensional scaling (MDS) keeps the distances between all pairs as well as possible.
• Isomap keeps distances measured along the data, by hopping from neighbour to neighbour, so a rolled-up sheet can be unrolled.
• t-SNE keeps only neighbourhoods: points that are close in the data should be close on the map. Distances between far-apart points are deliberately ignored.
• UMAP has the same aim as t-SNE, from a different theory; it is faster and keeps a little more of the large-scale layout. It needs a separate package, so it is not drawn here, but its maps look and must be read like t-SNE’s.

roll, position = make_swiss_roll(1500, noise=0.3, random_state=0)
MAPS = {
    "PCA": PCA(2),
    "MDS": MDS(2, n_init=1, random_state=0, normalized_stress="auto"),
    "Isomap": Isomap(n_neighbors=12),
    "t-SNE": TSNE(perplexity=30, random_state=0),
}

fig = plt.figure(figsize=(11, 2.6))
ax = fig.add_subplot(1, 5, 1, projection="3d")
ax.scatter(roll[:, 0], roll[:, 1], roll[:, 2], c=position, cmap="viridis", s=3)
ax.set(title="the data", xticks=[], yticks=[], zticks=[])
for k, (name, method) in enumerate(MAPS.items(), start=2):
    ax = fig.add_subplot(1, 5, k)
    xy = method.fit_transform(roll)
    ax.scatter(xy[:, 0], xy[:, 1], c=position, cmap="viridis", s=3)
    ax.set(title=name, xticks=[], yticks=[])
fig.tight_layout()
plt.show()

Figure 7.1.8: A ‘Swiss roll’: 1,500 points on a rolled-up sheet in three dimensions (left), coloured by their position along the sheet. PCA and MDS, keeping straight-line distances, show the roll end-on: the sheet stays rolled up, and points on neighbouring turns sit close although they are far apart along the sheet. Isomap, measuring along the sheet, unrolls it. t-SNE keeps neighbours together, so the colours run smoothly, but it tears the sheet into pieces and bends them: local structure kept, global shape lost.

t-SNE turns distances into neighbour probabilities. In the data, point j is a neighbour of point i with probability proportional to a normal curve around i; on the map, with a heavy-tailed curve (a t distribution with one degree of freedom, hence the name). The map is moved by gradient descent until the two sets of probabilities agree, measured by the Kullback-Leibler divergence:

$$p_{j \mid i} \propto e^{-\|x_i - x_j\|^2 / 2\sigma_i^2}, \qquad q_{ij} \propto \big(1 + \|y_i - y_j\|^2\big)^{-1}, \qquad \text{minimise } \sum_{i, j} p_{ij} \log \frac{p_{ij}}{q_{ij}}$$

The width $\sigma_i$ is set for each point so that it has about the same effective number of neighbours, the perplexity, typically 5 to 50. On the digits, where PCA overlaps several digits, t-SNE separates nearly all of them, without ever seeing the labels:

pixels = digits.data / 16
flat = PCA(2).fit_transform(pixels)
tsne = TSNE(perplexity=30, random_state=0).fit_transform(pixels)


def plot_map(ax: Axes, xy: np.ndarray, title: str) -> None:
    """A 2D map of the digits, each digit drawn as its own coloured numeral at the group centre."""
    for d in range(10):
        group = xy[digits.target == d]
        ax.scatter(group[:, 0], group[:, 1], s=3, color=plt.cm.tab10(d), alpha=0.6)
        centre = np.median(group, axis=0)
        ax.text(*centre, str(d), fontsize=12, fontweight="bold", ha="center", va="center")
    ax.set(title=title, xticks=[], yticks=[])


fig, axes = style.subplots(2, height=4, width=10)
plot_map(axes[0], flat, "PCA")
plot_map(axes[1], tsne, "t-SNE")
fig.tight_layout()
plt.show()

Figure 7.1.9: The 1,797 digits (64 pixels each) mapped to two dimensions, coloured by the digit, which neither method saw. Left: the first two principal components; several digits overlap. Right: t-SNE with perplexity 30; nearly every digit forms its own island (the 1s, written two ways, form two), with a few digits stranded on the wrong island.

t-SNE maps are easy to over-read. Three things on them carry no meaning: the size of a cluster, the distance between clusters, and the exact shape. And with a very small perplexity, t-SNE finds clumps even in data without any:

local = np.random.default_rng(0)
blobs = np.r_[
    local.normal(0, 1, (200, 10)),
    local.normal(0, 10, (200, 10)) + np.r_[40, np.zeros(9)],
    local.normal(0, 1, (200, 10)) + np.r_[200, np.zeros(9)],
]
group_of = np.repeat([0, 1, 2], 200)
noise = local.normal(0, 1, (600, 10))
GROUP_COLORS = np.array([style.BLUE, style.CORAL, style.GREEN])


def plot_points(ax: Axes, xy: np.ndarray, title: str, colors: np.ndarray | str) -> None:
    """A 2D scatter without axis values."""
    ax.scatter(xy[:, 0], xy[:, 1], s=3, c=colors)
    ax.set(title=title, xticks=[], yticks=[])


fig, axes = plt.subplots(2, 4, figsize=(11, 5.2))
plot_points(
    axes[0, 0], PCA(2).fit_transform(blobs), "three groups (PCA)", GROUP_COLORS[group_of]
)
for ax, perplexity in zip(axes[0, 1:], (5, 30, 100)):
    xy = TSNE(perplexity=perplexity, random_state=0).fit_transform(blobs)
    plot_points(ax, xy, f"t-SNE, perplexity {perplexity}", GROUP_COLORS[group_of])
plot_points(axes[1, 0], PCA(2).fit_transform(noise), "pure noise (PCA)", style.GREY)
for ax, perplexity in zip(axes[1, 1:], (2, 5, 30)):
    xy = TSNE(perplexity=perplexity, random_state=0).fit_transform(noise)
    plot_points(ax, xy, f"t-SNE, perplexity {perplexity}", style.GREY)
fig.tight_layout()
plt.show()

Figure 7.1.10: Top: three groups of 200 points in 10 dimensions, shown on their first two principal components (left): one tight, one ten times wider, and one far away. t-SNE (right three, perplexities 5, 30 and 100) draws the tight and the wide group at about the same size, and at perplexities 5 and 30 the two gaps look alike; only at 100 does some of the difference show. Bottom: 600 points of pure noise, without any groups (left); at a very small perplexity, 2, t-SNE breaks them into small clumps that look like structure.

So use t-SNE and UMAP to look, not to measure: cluster on the data (or on principal components), not on the map, try several perplexities and seeds, and check that groups seen on the map also show up in the original measurements.

7.1.10 In one picture

• Principal components are weighted sums $z = \sum_j \phi_j x_j$ with the largest variance; the loadings say what they are made of, the scores where each case sits.
• They are the eigenvectors of the correlation matrix; PVE = eigenvalue / total.
• Choose how many with the scree plot, a total like 80%, parallel analysis or cross-validation.
• Standardise first, or the measurement with the largest numbers takes over.
• PCA is the best low-dimensional compression, and can fill in missing values.
• Factor analysis models hidden causes of the correlations; rotate the factors to read them.
• Nonlinear maps: MDS keeps distances, Isomap distances along the data, t-SNE and UMAP neighbourhoods. On a t-SNE map, cluster sizes and gaps mean nothing; vary the perplexity.

Next: finding groups.