--- title: 'Multivariate Analysis' source: 'https://academia.sh/en/courses/statistical-analysis-and-exploratory-data-analysis/multivariate-examination' course: 'Statistical Analysis and Exploratory Data Analysis' language: en updated: '2026-10-07T04:31:13+00:00' license: 'CC BY-SA 4.0' objectives: - 'Relates the number of pairs in a correlation matrix to the number of columns and derives the threshold from permutation rather than from a distribution table: computes where the 95th percentile of the coefficient distribution falls in 2,000 runs' - 'Counts how the number of pairs exceeding the threshold in independent columns grows with the number of comparisons and compares it with the expected value' - 'Measures the gain of dimensionality reduction with component shares and shows that there is no gain when the columns do not repeat one another' --- # Multivariate Analysis A correlation matrix of six columns carries 15 pairs, and the threshold comes not from a distribution table but from permutation: in 2,000 permutations the 95th percentile of the coefficient distribution comes out as 0.0354, and 5 of the 15 pairs exceed this threshold — but three of the five exceed it by definition and the other two because of the subscriber number column. In columns generated completely independently of 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; the highest coefficient, 0.049, is clearly above the threshold. Dimensionality reduction is reduced to a single question: in the measurement network columns the first two components carry 49.3% of the variation, in the dataset generated from two hidden drivers 82.1%, and in independent columns only 35.0%. 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 ```python # 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 ```python # 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 ```python # 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.