Multivariate Analysis

The previous lesson examined two columns together and wrote its margin of variation next to the coefficient of a single pair. In a real analysis, though, pairs are not asked one by one; a table is taken and the coefficients of all numeric columns are computed at the same time. Six columns produce 15 pairs, twenty columns 190 pairs. Each pair is a separate question, and each question gives chance one more shot.

The unit of measure of this lesson is the spurious relationship produced by the number of comparisons. In columns generated completely independently of one another, the number of pairs exceeding the threshold is counted and it is shown how this number grows with the number of columns. The second measure is the gain of dimensionality reduction, and it is reduced to a single question: how many components carry how much of the variation.

• EA16. Multivariate analysis looks at more than two columns at once. A correlation matrix with $k$ columns carries $k(k-1)/2$ pairs; the number of pairs grows quadratically with the number of columns.
• EA17. Every cell in the matrix is a separate question. A threshold that is defensible for a single pair automatically inflates the number of pairs exceeding the threshold when it is applied to all pairs at once.
• EA18. Here the threshold comes not from a distribution table but from permutation: a column is permuted within itself 2,000 times, and the 95th percentile of the resulting coefficient distribution is taken as the threshold.
• EA19. Dimensionality reduction is covered in this lesson through a single question: how many components carry how much of the variation. Choice of method, rotation and interpretation of the components are not written.
• EA20. Component shares are read from the eigenvalues of the correlation matrix; the sum of the eigenvalues equals the number of columns, and when written as shares the total is 100. If the columns do not repeat one another, the shares are distributed close to equal.

The Matrix and the Threshold

# exploration.py — MODEL: the fictional measurement network of K03 and K04 is rebuilt
# from the same seed. The dataset is fictional; no real subscriber, meter or region appears.
import math

SEED, RAW, M32 = 20260218, 1400, 0xFFFFFFFF
PERIODS = ["2026-01", "2026-02", "2026-03"]
TARIFFS = [10, 25, 40]
REGION = [("north", 0.28, 21), ("south", 0.22, 17), ("east", 0.18, 26),
          ("west", 0.14, 14), ("central", 0.18, 23)]


def rng(t: int):
    """Linear congruential generator: seed written in the code, same output every run."""
    x = ((t ^ (t >> 16)) * 2246822507) & M32
    x = ((x ^ (x >> 13)) * 3266489909) & M32
    s = [(x ^ (x >> 16)) & M32]

    def next_value() -> float:
        s[0] = (s[0] * 1664525 + 1013904223) & M32
        return s[0] / 4294967296
    return next_value


def pick(u: float, w: list[float]) -> int:
    t = 0.0
    for i, x in enumerate(w):
        t += x
        if u < t:
            return i
    return len(w) - 1


_subscribers, _dropped = [], []
for i in range(RAW):
    r = rng(SEED + i)
    b = REGION[pick(r(), [x[1] for x in REGION])]
    k = {"subscriber_id": 10001 + i, "region": b[0], "base": b[2],
         "household_size": pick(r(), [0.06, 0.24, 0.30, 0.24, 0.11, 0.05]),
         "satisfaction": pick(r(), [0.08, 0.14, 0.27, 0.34, 0.17]) + 1}
    (_dropped if r() < 0.046 else _subscribers).append(k)
_valid = {a["subscriber_id"] for a in _subscribers}

READINGS = []
for k in _subscribers + _dropped:
    r = rng(SEED + 7000 + k["subscriber_id"])
    for d in range(pick(r(), [0.05, 0.12, 0.21, 0.62])):
        t = 0.0 if r() < 0.038 else math.floor(k["base"] * math.exp(
            (r() + r() + r() - 1.5) * 0.62) * (1 - d * 0.05) * 100 + 0.5) / 100
        if k["subscriber_id"] not in _valid:
            r()
            continue
        READINGS.append({"subscriber_id": k["subscriber_id"], "region": k["region"],
                         "period": PERIODS[d], "consumption_m3": t,
                         "household_size": k["household_size"],
                         "satisfaction": k["satisfaction"],
                         "tier": 1 if t < TARIFFS[0] else 2 if t < TARIFFS[1]
                         else 3 if t < TARIFFS[2] else 4})
        r()

def mean(v: list[float]) -> float:
    return sum(v) / len(v)


def corr(x: list[float], y: list[float]) -> float:
    mx, my = mean(x), mean(y)
    sx = math.sqrt(sum((a - mx) ** 2 for a in x))
    sy = math.sqrt(sum((b - my) ** 2 for b in y))
    return sum((a - mx) * (b - my) for a, b in zip(x, y)) / (sx * sy)


def align(w: list[int], left: int = 1):
    """The first `left` columns are left-aligned, the rest right-aligned."""
    return lambda h: "".join(str(v).ljust(w[i]) if i < left else str(v).rjust(w[i])
                             for i, v in enumerate(h))


COLUMNS = [("consumption", [o["consumption_m3"] for o in READINGS]),
           ("tier", [float(o["tier"]) for o in READINGS]),
           ("household_size", [float(o["household_size"]) for o in READINGS]),
           ("satisfaction", [float(o["satisfaction"]) for o in READINGS]),
           ("period", [float(PERIODS.index(o["period"])) for o in READINGS]),
           ("subscriber_id", [float(o["subscriber_id"]) for o in READINGS])]
K = len(COLUMNS)
m = align([15] + [11] * K)
print(m(["column"] + [a[:9] for a, _ in COLUMNS]))
for i, (name, x) in enumerate(COLUMNS):
    print(m([name] + [f"{corr(x, y):+.3f}" for _, y in COLUMNS[:i + 1]]))

RUNS = 2000


def threshold_95(x: list[float], seed: int) -> float:
    """The column is permuted against itself; the 95th percentile of the |coefficient| distribution is the threshold."""
    r, n = rng(seed), len(x)
    mx = mean(x)
    d = [a - mx for a in x]
    total_sq = sum(a * a for a in d)
    values = []
    for _ in range(RUNS):
        h = d[:]
        for i in range(n - 1, 0, -1):
            j = int(r() * (i + 1))
            h[i], h[j] = h[j], h[i]
        values.append(abs(sum(a * b for a, b in zip(d, h)) / total_sq))
    return sorted(values)[int(0.95 * RUNS)]


THRESHOLD = threshold_95(COLUMNS[0][1], SEED + 901)
pairs = [(COLUMNS[i][0], COLUMNS[j][0], corr(COLUMNS[i][1], COLUMNS[j][1]))
         for i in range(K) for j in range(i + 1, K)]
exceeded = [c for c in pairs if abs(c[2]) >= THRESHOLD]
print(f"\n{K} columns, {len(pairs)} pairs; in {RUNS} permutations the 95th percentile of the "
      f"|coefficient| distribution is {THRESHOLD:.4f}")
print(f"pairs exceeding the threshold: {len(exceeded)} / {len(pairs)}")
for a, b, v in sorted(exceeded, key=lambda c: -abs(c[2])):
    print(f"  {a} / {b}: {v:+.3f}")

column           consumpti       tier  household  satisfact     period  subscribe
consumption         +1.000
tier                +0.883     +1.000
household_size      +0.010     +0.007     +1.000
satisfaction        -0.010     -0.015     -0.004     +1.000
period              -0.111     -0.116     +0.000     -0.006     +1.000
subscriber_id       +0.031     +0.041     +0.046     -0.001     -0.005     +1.000

6 columns, 15 pairs; in 2000 permutations the 95th percentile of the |coefficient| distribution is 0.0354
pairs exceeding the threshold: 5 / 15
  consumption / tier: +0.883
  tier / period: -0.116
  consumption / period: -0.111
  household_size / subscriber_id: +0.046
  tier / subscriber_id: +0.041

Six columns produced 15 pairs and the threshold came out as 0.0354. The threshold was not read from a table: the consumption column was permuted within itself 2,000 times, the coefficient was recomputed in every run, and the 95th percentile of the resulting distribution was taken as the threshold. The meaning of this threshold is narrow and is exactly this — when there is no relationship at all, in a single pair, a coefficient of this size is seen in five of a hundred runs.

Five pairs exceed the threshold, but not all five are of the same kind. The +0.883 between consumption and tier is not a discovery but a definition; tier is computed from consumption with the limits 10, 25 and 40 m³. The −0.116 between tier and period is the same as the −0.111 between consumption and period, because tier is a rewriting of consumption — two lines count the same fact twice. What remains are the two pairs that include the subscriber number: the subscriber number is a record number and can be the cause of nothing.

Five of the 15 pairs exceeded the threshold and the number of independent findings that will enter the report is at most one. The matrix itself does not produce a finding; it enlarges the number of questions to be asked, and the larger it grows, the easier it becomes to find a cell that exceeds the threshold.

The Number of Pairs Produces Findings

# multiplicity.py — FICTIONAL: columns with no relation to one another are generated and the
# same threshold is applied to all pairs. What is sought is the number of pairs that exceed it.
def independent(k: int, n: int, seed: int) -> list[list[float]]:
    r = rng(seed)
    return [[r() + r() + r() for _ in range(n)] for _ in range(k)]


INDEPENDENT = independent(20, len(READINGS), SEED + 902)
_m = [mean(s) for s in INDEPENDENT]
_d = [[a - t for a in s] for s, t in zip(INDEPENDENT, _m)]
_n = [math.sqrt(sum(a * a for a in s)) for s in _d]
pairwise = lambda i, j: sum(a * b for a, b in zip(_d[i], _d[j])) / (_n[i] * _n[j])

g = align([8, 8, 12, 12, 12])
print(g(["columns", "pairs", "exceeded", "expected", "largest"]) + "   largest pair")
for k in (4, 6, 8, 12, 20):
    combos = [(i, j) for i in range(k) for j in range(i + 1, k)]
    values = {(i, j): abs(pairwise(i, j)) for i, j in combos}
    best = max(values, key=values.get)
    print(g([k, len(combos), sum(1 for v in values.values() if v >= THRESHOLD),
             f"{0.05 * len(combos):.1f}", f"{values[best]:.3f}"])
          + f"   column{best[0] + 1} / column{best[1] + 1}")
print(f"\nthreshold {THRESHOLD:.4f}; none of the columns is derived from another")

columns    pairs    exceeded    expected     largest   largest pair
4              6           0         0.3       0.021   column1 / column2
6             15           1         0.8       0.037   column2 / column6
8             28           2         1.4       0.048   column1 / column8
12            66           5         3.3       0.049   column6 / column12
20           190          12         9.5       0.049   column6 / column12

threshold 0.0354; none of the columns is derived from another

None of the columns in this table is derived from another; all of them come from separate draws and there is no relationship among them. Even so, 0 pairs in four columns, 1 in six columns, 2 in eight columns and 12 in twenty columns exceed the threshold.

Where the number comes from is written in the definition of the threshold. The threshold is the value exceeded in five of a hundred runs for a single pair; when the same threshold is applied to 190 pairs, about 190 × 0.05 = 9.5 pairs are expected to exceed it. The observed 12 is near the expected 9.5. So the pairs exceeding the threshold are not a defect but a direct consequence of the definition of the threshold itself.

The danger lies in pulling a single row out of this table. The sentence “A relationship of 0.049 was found between the sixth and the twelfth column, the threshold is 0.0354” contains correct numbers and is entirely made up: that pair was selected because it is the largest among 190 pairs. Next to a coefficient, not only its margin but also how many pairs it was selected from must be written. “One of a single pair” and “the largest of 190 pairs” are not the same sentence.

One way to reduce the number of comparisons is not to raise the threshold but to lower the number of columns. This is the subject of the next section.

How Many Components

# component.py — the single question of dimensionality reduction: how many components carry how
# much of the variation. Shares are read from the eigenvalues of the correlation matrix.
def eigenvalues(A: list[list[float]]) -> list[float]:
    """Eigenvalues of a symmetric matrix: rotated until the largest off-diagonal entry drops to zero."""
    n, M = len(A), [s[:] for s in A]
    for _ in range(200):
        p, q, largest = 0, 1, 0.0
        for i in range(n):
            for j in range(i + 1, n):
                if abs(M[i][j]) > largest:
                    p, q, largest = i, j, abs(M[i][j])
        if largest < 1e-12:  # the rotation stops when the off-diagonal entries fall below this threshold
            break
        th = (M[q][q] - M[p][p]) / (2 * M[p][q])
        t = (1.0 if th >= 0 else -1.0) / (abs(th) + math.sqrt(th * th + 1))
        c = 1 / math.sqrt(t * t + 1)
        s = t * c
        for k in range(n):
            kp, kq = M[k][p], M[k][q]
            M[k][p], M[k][q] = c * kp - s * kq, s * kp + c * kq
        for k in range(n):
            pk, qk = M[p][k], M[q][k]
            M[p][k], M[q][k] = c * pk - s * qk, s * pk + c * qk
    return sorted((M[i][i] for i in range(n)), reverse=True)


def shared_driver(seed: int) -> list[list[float]]:
    """FICTIONAL: six columns are generated from two hidden drivers and noise."""
    r = rng(seed)
    f1 = [r() + r() + r() - 1.5 for _ in range(len(READINGS))]
    f2 = [r() + r() + r() - 1.5 for _ in range(len(READINGS))]
    WEIGHTS = [(1.0, 0.0), (0.9, 0.2), (0.8, 0.3), (0.1, 1.0), (0.2, 0.9), (0.3, 0.8)]
    return [[a * x + b * y + 0.55 * (r() + r() + r() - 1.5) for x, y in zip(f1, f2)]
            for a, b in WEIGHTS]


DATASETS = [("measurement network", [x for _, x in COLUMNS]),
            ("shared driver", shared_driver(SEED + 903)),
            ("independent", independent(6, len(READINGS), SEED + 904))]
b = align([20] + [9] * 6)
print(b(["dataset"] + [f"c{i + 1}" for i in range(6)]) + "   first two components")
for name, S in DATASETS:
    eig = eigenvalues([[corr(a, c) for c in S] for a in S])
    total = sum(eig)
    print(b([name] + [f"{100 * v / total:.1f}" for v in eig])
          + f"   {100 * (eig[0] + eig[1]) / total:.1f} percent")
print(f"\nshares sum to 100 in every row; the sum of eigenvalues is the number of columns ({len(COLUMNS)})")

dataset                    c1       c2       c3       c4       c5       c6   first two components
measurement network      31.9     17.4     16.7     16.2     15.9      1.9   49.3 percent
shared driver            55.2     26.9      4.9      4.7      4.3      4.0   82.1 percent
independent              17.8     17.2     16.7     16.4     16.3     15.5   35.0 percent

shares sum to 100 in every row; the sum of eigenvalues is the number of columns (6)

The three rows show when dimensionality reduction does work. In the dataset generated from a shared driver, six columns come from two hidden variables and the first two components carry 82.1% of the variation; writing two numbers instead of six columns means leaving 17.9% of the variation behind. In the independent columns the shares are almost equal, from 17.8 to 15.5, and the first two components carry only 35.0%: there is nothing to reduce, because the columns never repeat one another.

The measurement network row lies between the two. The first component carries 31.9%, while the next four are squeezed between 17.4 and 15.9. The real information is in the last column: the share of the sixth component has dropped to 1.9%. This drop is a sign of repetition, and what the repetition is can be known from the first table — consumption and tier write the same thing twice. The share table does not say which columns repeat one another, but it gives whether there is repetition at a glance.

The two sections meet here. Six columns produce 15 pairs; if the columns are reduced to two components the number of pairs falls to 1. This is the measurable gain of dimensionality reduction: it lowers the number of cells to look at and thus also lowers the chance that coincidence exceeds the threshold. Its cost is that the components no longer have a name like “consumption” or “household size”. Choice of method and interpretation of the components are outside this lesson.

Summary

• A matrix of six columns carries 15 pairs; the threshold comes not from a distribution table but from permutation, and as the 95th percentile of the 2,000-run distribution it comes out as 0.0354.
• Of the five pairs exceeding the threshold, three exceed it by definition (tier is computed from consumption) and two come from the subscriber number column; the number of independent findings is at most one.
• In columns with no relation to one another, when the number of pairs rises from 6 to 190 the number of pairs exceeding the threshold rises from 0 to 12, and the expected value is 9.5.
• Next to a coefficient, its margin is written together with how many pairs it was selected from; “one of a single pair” and “the largest of 190 pairs” are not the same claim.
• Component shares tell the gain of dimensionality reduction in advance: the first two components carry 82.1% in the shared-driver dataset, 49.3% in the measurement network and 35.0% in the independent columns.

Next Step

The four lessons of this topic approached the dataset always through computation: a mean was taken, a distribution was counted, a coefficient was computed, a matrix was built. Their common assumption was that the values in the dataset are measured values. Yet during exploration the first things caught are often not found by computation but seen while looking at the table: the repetition of the same value with suspicious frequency, records piling up at an upper limit, field combinations that should never occur together, one period being markedly sparser than the others. The next lesson takes up this first inspection and counts which summary each caught problem moves and by how much.