Bivariate Analysis

The previous lesson looked at a single column and wrote its shape with five numbers and bin counts. Most questions ask for two columns, though: does consumption rise together with household size, does it fall as the periods go on, is it related to the satisfaction level. When two columns are examined together a new summary comes out — the correlation coefficient — and it too compresses a pile of numbers into a single number.

The unit of measure of this lesson is the pattern the coefficient cannot see. Three separate relationships that give the same coefficient are built, all three are reported with the same number, and it is shown that the difference between them never shows in the coefficient.

• EA11. Bivariate analysis considers two columns together. The correlation coefficient is the ratio of the joint variation to the product of the two separate variations and measures only linear joint variation.
• EA12. The margin is written next to the coefficient: one of the columns is randomly rearranged and it is counted in how many runs a coefficient of the same size appears.
• EA13. The coefficient is not a substitute for the scatter plot. The scatter is not drawn here; it is printed with slice means: the driver column is divided into five equal parts and the mean of the second column is written for each part.
• EA14. Cross tabulation is the joint count of two categorical columns. On a nominal scale the correlation coefficient is not defined; the relationship is read from cell counts.
• EA15. Rank correlation looks only at the ranking. The numbers given to levels on an ordinal scale are a coding decision — the scale definitions were set up in the Introduction to Data Analytics course and are not repeated here. The coefficient is affected by this decision; rank correlation is not.

The Coefficient and Its Margin

# 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:
    """Correlation coefficient: the ratio of the joint variation to the product of the two separate variations."""
    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 ranks(v: list[float]) -> list[float]:
    """Equal values take the average rank; this is what an ordinal scale requires."""
    positions: dict[float, list[int]] = {}
    for k, i in enumerate(sorted(range(len(v)), key=lambda j: v[j])):
        positions.setdefault(v[i], []).append(k + 1)
    return [mean(positions[x]) for x in v]


def rank_corr(x: list[float], y: list[float]) -> float:
    return corr(ranks(x), ranks(y))


RUNS = 2000
TOLERANCE = 1e-9  # the same coefficient computed two ways shifts in the last digit


def corr_margin(x: list[float], y: list[float], seed: int) -> tuple[int, float]:
    """The y column is randomly permuted; the coefficient is re-measured in every run."""
    r, observed = rng(seed), corr(x, y)
    mx, my = mean(x), mean(y)
    divisor = math.sqrt(sum((a - mx) ** 2 for a in x) * sum((b - my) ** 2 for b in y))
    dx, dy = [a - mx for a in x], [b - my for b in y]
    exceeded, extreme = 0, 0.0
    for _ in range(RUNS):
        h = dy[:]
        for i in range(len(h) - 1, 0, -1):
            j = int(r() * (i + 1))
            h[i], h[j] = h[j], h[i]
        c = sum(a * b for a, b in zip(dx, h)) / divisor
        exceeded += abs(c) >= abs(observed) - TOLERANCE
        extreme = max(extreme, abs(c))
    return exceeded, extreme


CONSUMPTION = [o["consumption_m3"] for o in READINGS]
PAIRS = [("household_size", [o["household_size"] for o in READINGS]),
         ("satisfaction", [o["satisfaction"] for o in READINGS]),
         ("period order", [PERIODS.index(o["period"]) for o in READINGS]),
         ("tier", [o["tier"] for o in READINGS])]
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))


y = align([30, 7, 11, 12, 10, 10])
print(y(["pair", "n", "coef.", "rank coef.", "extreme", "exceeded"]))
for i, (name, v) in enumerate(PAIRS):
    exceeded, extreme = corr_margin(CONSUMPTION, v, SEED + 710 + i)
    print(y([f"consumption / {name}", len(v), f"{corr(CONSUMPTION, v):+.3f}",
             f"{rank_corr(CONSUMPTION, v):+.3f}", f"{extreme:.3f}", exceeded]))
print(f"\n{RUNS} runs; extreme = the largest coefficient seen in the permuted sets")

pair                                n      coef.  rank coef.   extreme  exceeded
consumption / household_size     3199     +0.010      +0.014     0.069      1189
consumption / satisfaction       3199     -0.010      -0.010     0.058      1154
consumption / period order       3199     -0.111      -0.113     0.061         0
consumption / tier               3199     +0.883      +0.813     0.065         0

2000 runs; extreme = the largest coefficient seen in the permuted sets

If only the coefficient column is looked at, the first two rows can be read as “a weak but existing relationship”: +0.010 and −0.010 are not zero. The margin column closes this. The permuted sets produce a coefficient of this size in 1,189 and in 1,154 of 2,000 runs. Both rows are noise.

The third row goes the other way. The −0.111 between period order and consumption is small too, but the largest coefficient seen in the permuted sets is 0.061 and none of the 2,000 runs reaches the observed one. What speaks is not the size of the coefficient but where it stands against the permutation distribution.

The fourth row shows that the margin alone is not enough either. The coefficient between tariff tier and consumption is +0.883 and its margin is zero, but the tier is computed from consumption: the limits of 10, 25 and 40 m³ determine it. The coefficient here measures a definition, not a discovery.

Same Coefficient, Three Different Relationships

# same_coefficient.py — FICTIONAL: three different relationships that give the same coefficient.
# The free parameter (noise weight, distance of the extreme point) is fitted to the target by bisection.
N, TARGET = 557, 0.600


def driver(seed: int):
    r = rng(seed)
    return ([round(4 + 16 * r(), 2) for _ in range(N)],
            [r() + r() + r() - 1.5 for _ in range(N)])


def tune(f, upper: float) -> float:
    """f(p) gives the coefficient and varies monotonically with p; the target is approached by bisection."""
    lower, increasing = 0.0, f(upper) > f(0.0)
    for _ in range(50):
        mid = (lower + upper) / 2
        if (f(mid) < TARGET) == increasing:
            lower = mid
        else:
            upper = mid
    return (lower + upper) / 2


def mix(x, base, noise):
    w = tune(lambda p: corr(x, [t + p * g for t, g in zip(base, noise)]), 200.0)
    return [round(t + w * g, 2) for t, g in zip(base, noise)]


x1, g1 = driver(SEED + 801)
y1 = mix(x1, [1.4 * a for a in x1], g1)
x2, g2 = driver(SEED + 802)
y2 = mix(x2, [-0.16 * (a - 14.0) ** 2 + 40 for a in x2], g2)  # peak at 14 m3

r3 = rng(SEED + 803)
x3 = [round(4 + 16 * r3(), 2) for _ in range(N - 1)]
y3 = [round(20 + 9 * (r3() + r3() + r3() - 1.5), 2) for _ in range(N - 1)]
far = tune(lambda k: corr(x3 + [20 + k], y3 + [20 + 3 * k]), 400.0)
x3, y3 = x3 + [round(20 + far, 2)], y3 + [round(20 + 3 * far, 2)]


def slice_mean(x: list[float], y: list[float], d: int = 5) -> list[float]:
    s = sorted(zip(x, y))
    return [mean([b for _, b in s[i * len(s) // d:(i + 1) * len(s) // d]])
            for i in range(d)]


rel = align([13, 11, 12])
print(rel(["relationship", "coef.", "rank coef."]) + "   slice means (y)")
for name, x, y in [("linear", x1, y1), ("curvilinear", x2, y2), ("leverage", x3, y3)]:
    print(rel([name, f"{corr(x, y):+.3f}", f"{rank_corr(x, y):+.3f}"]) + "   "
          + " ".join(f"{v:5.1f}" for v in slice_mean(x, y)))
print(f"\nin the leverage relationship the extreme point is ({x3[-1]}, {y3[-1]}); the remaining {N - 1} points have an "
      f"x range of 4-20")
print(f"when that single point is removed the coefficient becomes {corr(x3[:-1], y3[:-1]):+.3f}")

relationship       coef.  rank coef.   slice means (y)
linear            +0.600      +0.594     6.6  12.6  15.4  21.4  24.9
curvilinear       +0.600      +0.559    27.8  34.9  38.8  39.8  37.1
leverage          +0.600      -0.052    20.1  20.0  19.9  19.4  21.7

in the leverage relationship the extreme point is (104.4, 273.2); the remaining 556 points have an x range of 4-20
when that single point is removed the coefficient becomes -0.057

The coefficients of the three relationships are the same to the third decimal, and the sentence “correlation coefficient +0.600” does not tell the three apart. What does is the last column: when the driver column is divided into five slices and the mean of the second column is written for each slice, three different patterns come out.

In the linear relationship the slice means climb smoothly from 6.6 to 24.9. In the curvilinear one they rise to 39.8 from 27.8 and return to 37.1: the second column increases up to a point, then decreases. The coefficient cannot see this turn, because it measures linear joint variation and the two sides of the turn partly cancel each other. Rank correlation cannot see it either; the difference between 0.559 and 0.594 does not count as a warning on its own.

The third row is the harshest. The slice means hold along the line 20.1 — 20.0 — 19.9 — 19.4 and jump to 21.7 in the last slice: not a trend, but the trace of a single point. That point carries the coefficient — while 556 points sit between 4 and 20, one is at (104.4, 273.2) — and when that point is removed the coefficient falls to −0.057. Rank correlation says this from the start, because ranking does not care about the distance of a point.

The summary of the measure is this: the coefficient reduces a two-dimensional shape to a single number and does not say what it reduced; five slice means tell all three relationships apart.

Categorical and Ordinal Columns

# categorical.py — a cross tabulation counts two categorical columns together; region is
# a nominal scale and the correlation coefficient is not defined for this pair.
from collections import Counter

cell = Counter((o["region"], o["tier"]) for o in READINGS)
REGIONS = sorted({o["region"] for o in READINGS})
TIERS = [1, 2, 3, 4]
c = align([9, 6] + [6] * 4 + [7] * 4)
print(c(["region", "n"] + [f"t{k}" for k in TIERS] + [f"%t{k}" for k in TIERS]))
for b in REGIONS + ["dataset"]:
    count = [sum(cell[(x, k)] for x in (REGIONS if b == "dataset" else [b]))
             for k in TIERS]
    print(c([b, sum(count)] + count + [f"{100 * s / sum(count):.1f}" for s in count]))

# On an ordinal scale the numbers given to levels are a decision: when the coding changes
# the coefficient changes; rank correlation does not, since it looks only at the ranking.
COLUMNS = [("tier", [o["tier"] for o in READINGS],
            [("equal interval", {1: 1, 2: 2, 3: 3, 4: 4}),
             ("doubling", {1: 1, 2: 2, 3: 4, 4: 8}),
             ("tier midpoint", {1: 5.0, 2: 17.5, 3: 32.5, 4: 45.0})]),
           ("satisfaction", [o["satisfaction"] for o in READINGS],
            [("equal interval", {d: d for d in range(1, 6)}),
             ("doubling", {1: 1, 2: 2, 3: 4, 4: 8, 5: 16})])]
o_ = align([13, 23, 12, 12], 2)
print("\n" + o_(["column", "coding", "coef.", "rank coef."]))
for col_name, col, codings in COLUMNS:
    coefs = [corr(CONSUMPTION, [m[d] for d in col]) for _, m in codings]
    for (name, m), value in zip(codings, coefs):
        print(o_([col_name, name, f"{value:+.4f}",
                  f"{rank_corr(CONSUMPTION, [m[d] for d in col]):+.4f}"]))
    shift = 100 * (max(map(abs, coefs)) / min(map(abs, coefs)) - 1)
    print(f"{col_name}: coefficient moved {shift:.0f} percent in absolute value, "
          f"rank correlation did not move")

region        n    t1    t2    t3    t4    %t1    %t2    %t3    %t4
central     600    30   352   198    20    5.0   58.7   33.0    3.3
east        557    21   281   225    30    3.8   50.4   40.4    5.4
north       921    42   661   210     8    4.6   71.8   22.8    0.9
south       645    59   538    48     0    9.1   83.4    7.4    0.0
west        476    88   377    11     0   18.5   79.2    2.3    0.0
dataset    3199   240  2209   692    58    7.5   69.1   21.6    1.8

column       coding                        coef.  rank coef.
tier         equal interval              +0.8832     +0.8125
tier         doubling                    +0.8417     +0.8125
tier         tier midpoint               +0.8800     +0.8125
tier: coefficient moved 5 percent in absolute value, rank correlation did not move
satisfaction equal interval              -0.0103     -0.0102
satisfaction doubling                    -0.0184     -0.0102
satisfaction: coefficient moved 79 percent in absolute value, rank correlation did not move

Region is a nominal scale: there is no order between north and east, and the correlation coefficient cannot be computed. The relationship is read from cell counts. In the west 79.2% of the 476 rows are in the second tier, and there are no rows in the fourth tier; in the east 5.4% of the 557 rows are in the fourth tier. In the whole dataset this share is 1.8%, so the share of the east is close to three times as large.

The table’s own margin shows in the cell counts. There are 8 rows in the fourth tier of the north and 11 in the third of the west; in a cell of 8 rows a single row moves the ratio by about one eighth. Unless the percentage is written together with how many observations it rests on, the table cannot be read.

The lower table measures the decision of the ordinal scale. The tariff tier can be numbered 1, 2, 3, 4; a doubling sequence and the midpoints of the tiers give the same ranking too. The coefficient nevertheless changes as 0.8832 — 0.8417 — 0.8800. On satisfaction the effect is larger: 79% in absolute value. Rank correlation stays in place on both columns.

The reason is in the definition of the scale. On an ordinal scale the distance between levels is not measured, only the order is known. Since the coefficient multiplies the level numbers, it counts the range between the first and second tier as equal to the range between the third and fourth; the scale does not support this. Rank correlation is independent of the coding.

Summary

• The coefficient is read with its margin: the coefficients +0.010 and −0.010 are exceeded in 1,189 and 1,154 of 2,000 permutations, and −0.111 in none; size alone does not tell them apart.
• A coefficient whose margin is zero may not be a finding either: +0.883 measures a definition, since the tier is computed from consumption.
• Three relationships that give the same +0.600 coefficient are told apart by slice means: a rise from 6.6 to 24.9, a curve that rises to 39.8 and returns to 37.1, and a flat line around 20.
• In the leverage relationship a single point carries the coefficient; when that point is removed the coefficient falls to −0.057.
• On an ordinal scale the coding moves the coefficient — by 79% in absolute value on satisfaction — and does not move rank correlation.

Next Step

This lesson looked at a single pair and wrote its margin next to the coefficient. In a real dataset the number of pairs does not stay at two: four numeric columns produce six pairs, six columns produce fifteen pairs. Looking at each pair separately is asking a new question each time, and as questions multiply, the chance that coincidence produces a high coefficient in at least one of them grows. The next lesson builds the correlation matrix, counts how many pairs give a high coefficient in columns produced entirely independently, and computes how much of the variation how many components carry.