Principal component analysis

A simulated example, with the code and its output

The numbers below come from made-up data, made by the code on this page. They show the method. They are not output from this tool and not a finding about your data.

Question: If six measures move together, can they be summarised by fewer dimensions, and how many are real rather than noise?

Inputs: Simulated data, not real measurements. Built in by construction: 6 measures and 2 hidden factors. Measures m1 to m3 each equal 0.8 times factor A plus noise (sd 0.6); m4 to m6 each equal 0.8 times factor B plus noise. Seed 3 makes 200 rows for the main example. The repeat checks use seeds 1000 to 1299, 300 seeds each, with 200 or 30 rows, on two-factor data and on pure random noise with no structure at all. The random-data comparison uses 200 random datasets of the same size per row count.

Method: Principal component analysis. Standardise each measure, take the correlation matrix, and find its eigenvalues and eigenvectors. Each eigenvalue is how much variance a component carries; each loading shows how strongly a measure goes with a component. To judge how many components are more than noise, the example uses parallel analysis, one common way to compare against random data: make random data of the same size, take the 95th percentile of its eigenvalues, and keep leading components, in order, until one falls at or below that. This example does not confirm how the tool decides how many components to keep.

Result: In seed 3 the two largest components carried 0.387 and 0.380 of the variance, with eigenvalues 2.32 and 2.28 against random-data thresholds of 1.32 and 1.19. Both were kept and the third (0.42) was not. Loadings recovered the built-in groups: one component loaded 0.87 to 0.88 on m4 to m6 and at most 0.11 on m1 to m3; the other loaded 0.85 to 0.88 on m1 to m3. Across 300 seeds the rule kept exactly 2 components in 1.00 of seeds with 200 rows and 0.96 with 30 rows. On pure noise it kept 0 components in 0.95 (200 rows) and 0.93 (30 rows) of seeds and kept 1 in 0.04 and 0.07, close to the 5% level the 95th percentile threshold implies for the first component.

Limit: The factors here are simple, equal in strength and uncorrelated with each other, and the measures all have similar noise. Real data has skewed measures, outliers, missing values and factors that overlap, which this example does not show. It does not show rotation, how the tool names what drives each component, or how it reports results.

Python script (numpy and scipy only):

import numpy as np
# Simulated data. Built in: 6 measures, 2 hidden factors. m1-m3 load on factor A, m4-m6 on factor B (loading 0.8, noise sd 0.6).
def make(rng, n, two=True):
    A = rng.standard_normal(n); B = rng.standard_normal(n)
    X = np.empty((n, 6))
    for j in range(6):
        f = (A if j < 3 else B) if two else 0
        X[:, j] = 0.8 * f + 0.6 * rng.standard_normal(n)
    if not two: X = rng.standard_normal((n, 6))
    return X

def pca(X):
    Z = (X - X.mean(0)) / X.std(0, ddof=1)
    ev, V = np.linalg.eigh(np.corrcoef(Z, rowvar=False))
    o = np.argsort(-ev); return ev[o], V[:, o]

def kept(ev, thr):
    k = 0
    for a, b in zip(ev, thr):
        if a > b: k += 1
        else: break
    return k

def parallel(n, p, rng, reps=200):
    r = np.array([pca(rng.standard_normal((n, p)))[0] for _ in range(reps)])
    return np.percentile(r, 95, axis=0)       # 95th percentile of random-data eigenvalues

rng = np.random.default_rng(3)
X = make(rng, 200); ev, V = pca(X)
thr = parallel(200, 6, np.random.default_rng(99))
print("seed 3, n=200")
print("eigenvalues:", np.round(ev, 2))
print("share of variance:", np.round(ev / ev.sum(), 3))
print("random-data 95th pct:", np.round(thr, 2))
print("components kept (leading components above random):", kept(ev, thr))
for k in range(2):
    v = V[:, k] * np.sqrt(ev[k]); v = v if v[np.argmax(abs(v))] > 0 else -v
    print("PC%d loadings m1..m6:" % (k + 1), np.round(v, 2))
cache = {n: parallel(n, 6, np.random.default_rng(n)) for n in (30, 200)}
for n, two in ((200, True), (30, True), (200, False), (30, False)):
    k = [kept(pca(make(np.random.default_rng(1000 + s), n, two))[0], cache[n]) for s in range(300)]
    k = np.array(k)
    print("n=%3d %s: kept 0 in %.2f, 1 in %.2f, 2 in %.2f, 3+ in %.2f of seeds" % (n, "two factors" if two else "pure noise ", np.mean(k == 0), np.mean(k == 1), np.mean(k == 2), np.mean(k >= 3)))

Printed output from running the code above:

seed 3, n=200
eigenvalues: [2.32 2.28 0.42 0.35 0.34 0.3 ]
share of variance: [0.387 0.38  0.069 0.058 0.057 0.05 ]
random-data 95th pct: [1.32 1.19 1.08 1.   0.93 0.85]
components kept (leading components above random): 2
PC1 loadings m1..m6: [0.08 0.05 0.11 0.87 0.88 0.87]
PC2 loadings m1..m6: [ 0.88  0.87  0.85 -0.09 -0.05 -0.1 ]
n=200 two factors: kept 0 in 0.00, 1 in 0.00, 2 in 1.00, 3+ in 0.00 of seeds
n= 30 two factors: kept 0 in 0.00, 1 in 0.04, 2 in 0.96, 3+ in 0.00 of seeds
n=200 pure noise : kept 0 in 0.95, 1 in 0.04, 2 in 0.01, 3+ in 0.00 of seeds
n= 30 pure noise : kept 0 in 0.93, 1 in 0.07, 2 in 0.00, 3+ in 0.00 of seeds

What's in the report?

Reduces many measures to a few dimensions and shows which measures go with each.