Modeling Relationships with Regression

The previous lesson counted in how many combinations a difference was measured, but what it counted was always the difference between two groups. How consumption changes together with a continuous variable was never asked. When a line is asked about instead of two groups, the number that comes out is not a difference but a coefficient: by how many units consumption changes on average when the variable grows by one unit.

The number of this lesson is the margin of variation of the coefficient. A second number stands next to it and it is the real subject of the lesson: the same coefficient changes sign when a single variable is added to the dataset. Prediction accuracy measures are not computed here; mean absolute error, bias and the band were measured in the Estimation and Planning lesson of the Process, Team and Delivery course. Here the subject is the interpretation and the limit of the coefficient.

• AS31. The dataset is inherited from K03 and is fictional. It is reduced to the subscriber level: the value of each subscriber is the mean of that subscriber’s periods. Of 1,329 subscribers, 1,260 have at least one reading.
• AS32. The meter age field is fictional too and carries two components together: a component tied to the region and a component tied to the within-region consumption deviation. The aim is to collect two directions in a single variable.
• AS33. Least squares is written inside the lesson: the normal equations are built and solved by elimination with partial pivoting. No ready-made solver is used.
• AS34. The margin of variation of the coefficient is produced by drawing with replacement: the subscriber list is redrawn 1,000 times and the model is rebuilt at every draw.
• AS35. The confounder is the base level of the region: the base consumption of the region determines both the meter age and the consumption.

A Coefficient and Its Margin

Least squares looks for the coefficients that make the sum of the squares of the residuals smallest. This condition gives the normal equations, and the system of equations is solved by elimination. The block below builds the dataset, writes the solver and computes the univariate model.

# regression.py — MODEL. The fictional dataset of K03 is produced with the same seed and is
# reduced to the subscriber level. The meter age field is FICTIONAL too: it carries a component
# tied to the region and a component tied to the within-region consumption deviation together.
import math

SEED, RAW, M32 = 20260218, 1400, 0xFFFFFFFF
REGION = [("north", 0.28, 21), ("south", 0.22, 17), ("east", 0.18, 26),
          ("west", 0.14, 14), ("central", 0.18, 23)]


def rng(t):
    x = ((t ^ (t >> 16)) * 2246822507) & M32
    x = ((x ^ (x >> 13)) * 3266489909) & M32
    s = [(x ^ (x >> 16)) & M32]

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


def pick(u, w):
    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])]
    household = 1 + pick(r(), [0.06, 0.24, 0.30, 0.24, 0.11, 0.05])
    satisfied = 1 + pick(r(), [0.08, 0.14, 0.27, 0.34, 0.17])
    (DROPPED if r() < 0.046 else SUBSCRIBERS).append(
        {"subscriber_id": 10001 + i, "region": b[0], "base": b[2], "household": household,
         "satisfaction": satisfied})
VALID = {k["subscriber_id"] for k 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"] in VALID:
            READINGS.setdefault(k["subscriber_id"], []).append(t)
        r()

DATA = [dict(k, consumption=sum(READINGS[k["subscriber_id"]]) / len(READINGS[k["subscriber_id"]]))
        for k in SUBSCRIBERS if k["subscriber_id"] in READINGS]
REGION_MEAN = {}
for x in DATA:
    REGION_MEAN.setdefault(x["region"], []).append(x["consumption"])
REGION_MEAN = {b: sum(v) / len(v) for b, v in REGION_MEAN.items()}
for x in DATA:
    r = rng(SEED + 91000 + x["subscriber_id"])
    y = (14.0 - 0.75 * (x["base"] - 20.5)
         + 0.15 * (x["consumption"] - REGION_MEAN[x["region"]]) + (r() + r() + r() - 1.5) * 2.4)
    x["age"] = max(1, min(30, int(y + 0.5)))


def solve(A, b):
    """Solves the normal equations by elimination with partial pivoting."""
    n = len(b)
    M = [A[i][:] + [b[i]] for i in range(n)]
    for i in range(n):
        p = max(range(i, n), key=lambda k: abs(M[k][i]))
        M[i], M[p] = M[p], M[i]
        for k in range(i + 1, n):
            f = M[k][i] / M[i][i]
            for j in range(i, n + 1):
                M[k][j] -= f * M[i][j]
    x = [0.0] * n
    for i in range(n - 1, -1, -1):
        x[i] = (M[i][n] - sum(M[i][j] * x[j] for j in range(i + 1, n))) / M[i][i]
    return x


def ols(rows, fields):
    """Least squares: the coefficients and the coefficient of determination."""
    p = len(fields) + 1
    A, b = [[0.0] * p for _ in range(p)], [0.0] * p
    for s in rows:
        v = [1.0] + [float(s[a]) for a in fields]
        for i in range(p):
            for j in range(p):
                A[i][j] += v[i] * v[j]
            b[i] += v[i] * s["consumption"]
    coefs = solve(A, b)
    avg = sum(s["consumption"] for s in rows) / len(rows)
    sst = sum((s["consumption"] - avg) ** 2 for s in rows)
    sse = sum((s["consumption"] - sum(k * x for k, x in
                                      zip(coefs, [1.0] + [float(s[a]) for a in fields]))) ** 2
              for s in rows)
    return coefs, 1 - sse / sst


RESAMPLES = 1000
rs = rng(SEED + 77000)


def coefficient_interval(rows, fields, index):
    """Drawing with replacement: the subscriber list is redrawn 1,000 times."""
    n, drawn = len(rows), []
    for _ in range(RESAMPLES):
        sample = [rows[int(rs() * n)] for _ in range(n)]
        drawn.append(ols(sample, fields)[0][index])
    drawn.sort()
    return drawn[24], drawn[974]


print(f"fictional dataset (seed {SEED}): {len(SUBSCRIBERS)} subscribers, "
      f"{sum(len(v) for v in READINGS.values())} readings")
print(f"with at least one reading: {len(DATA)} subscribers")
coefs, r2 = ols(DATA, ["age"])
low, high = coefficient_interval(DATA, ["age"], 1)
print(f"\nconsumption = {coefs[0]:.3f} + ({coefs[1]:.3f}) x meter_age")
print(f"  meter age coefficient : {coefs[1]:.3f} m3/year")
print(f"  margin of variation   : [{low:.3f}, {high:.3f}]  ({RESAMPLES} draws with replacement)")
print(f"  explained variation   : {r2:.4f}   ({len(DATA)} observations)")

fictional dataset (seed 20260218): 1329 subscribers, 3199 readings
with at least one reading: 1260 subscribers

consumption = 28.352 + (-0.605) x meter_age
  meter age coefficient : -0.605 m3/year
  margin of variation   : [-0.714, -0.503]  (1000 draws with replacement)
  explained variation   : 0.0984   (1260 observations)

This is the line to be reported: when the meter age grows by one year, mean consumption falls by 0.605 m³, the margin of variation is [-0.714, -0.503]. The interval is far from zero and narrow; it rests on 1,260 observations. This number looks as if it satisfies the rule of the course.

The third number next to the coefficient is the coefficient of determination, 0.098. The meter age carries about a tenth of the variation in consumption. This low value is often read as “the relationship is weak” and turns into a one-line warning in the report. That reading is incomplete. The coefficient of determination says something about the strength of the relationship; it says nothing about the direction of the relationship.

A Confounder Turns the Coefficient Around

The base consumption level of the region determines two things at once: the higher the consumption in a region, the newer the meters, and the higher the consumption. This definition is a confounder. When it enters the model, what the coefficient is can be measured.

# Confounder: the base level of the region determines both the meter age and the consumption.
for x in DATA:
    r = rng(SEED + 99000 + x["subscriber_id"])
    x["random"] = r() * 10

MODELS = [("age", ["age"]), ("age + base", ["age", "base"]),
          ("age + base + satisfaction", ["age", "base", "satisfaction"]),
          ("age + base + random", ["age", "base", "random"])]
print("model                        age coefficient   explained variation")
for name, fields in MODELS:
    k, r2 = ols(DATA, fields)
    print(f"{name:<28} {k[1]:>15.3f}   {r2:>19.6f}")

low2, high2 = coefficient_interval(DATA, ["age", "base"], 1)
k1, _ = ols(DATA, ["age"])
k2, _ = ols(DATA, ["age", "base"])
print(f"\ncoefficient          observed   margin of variation")
print(f"univariate           {k1[1]:>8.3f}   [{low:.3f}, {high:.3f}]")
print(f"confounder added     {k2[1]:>8.3f}   [{low2:.3f}, {high2:.3f}]")

print("\nregion    n   mean meter age   mean consumption   within-group slope")
for b in sorted(REGION_MEAN, key=lambda b: -REGION_MEAN[b]):
    sub = [x for x in DATA if x["region"] == b]
    e = ols(sub, ["age"])[0][1]
    print(f"{b:<7} {len(sub):>3}   {sum(x['age'] for x in sub)/len(sub):>14.2f}"
          f"   {REGION_MEAN[b]:>16.3f}   {e:>19.3f}")

model                        age coefficient   explained variation
age                                   -0.605              0.098357
age + base                             1.825              0.505363
age + base + satisfaction              1.827              0.505764
age + base + random                    1.825              0.505377

coefficient          observed   margin of variation
univariate             -0.605   [-0.714, -0.503]
confounder added        1.825   [1.629, 2.003]

region    n   mean meter age   mean consumption   within-group slope
east    219             9.95             24.814                 2.416
central 231            12.13             22.483                 2.154
north   363            13.60             20.351                 1.777
south   258            16.69             16.728                 1.550
west    189            18.84             13.935                 0.771

The coefficient turned from -0.605 to 1.825. The two margins of variation, [-0.714, -0.503] and [1.629, 2.003], do not overlap at all. Both numbers are correct and both come from the same dataset; the only thing that differs is which variables are in the model.

The third table shows the mechanism in the open. Between regions the direction is clear: while the mean age of the meters rises from 9.95 to 18.84, mean consumption falls from 24.814 to 13.935. Within the regions, the slope is positive in all five of the five regions; the lowest is 0.771, the highest 2.416. The univariate model collects these two directions in a single number, and because the direction between regions weighs more, the result comes out negative. When the confounder enters the model, the direction between regions is separated out and the direction within regions remains.

The first table gives two separate ways in which the coefficient of determination does not speak. The first: while the coefficient of determination rose from 0.098 to 0.505, the coefficient changed sign. The rise of the coefficient of determination does not say which sign is right; it says only that the residuals got smaller. The second: when a completely randomly generated column is added to the model, the coefficient of determination rises from 0.505363 to 0.505377. Because least squares can give every new column a coefficient of zero at worst, this coefficient never falls. A coefficient of determination whose variable count is not written cannot be compared.

A coefficient is read as “while the other variables in the model are held constant”, and the whole of this reading depends on what is inside the model. What the coefficient says is a relationship; it is not causation, and it carries no information about what a variable that is not in the model does.

Residuals and the Width of the Margin

The assumption check is done here. If the model is built correctly, the residuals should gather around zero in every breakdown and their dispersion should not wander together with the prediction level.

# Assumption check: the mean of the residuals in a breakdown and their dispersion against the prediction.
def residuals(rows, fields):
    coefs = ols(rows, fields)[0]
    result = []
    for s in rows:
        v = [1.0] + [float(s[a]) for a in fields]
        prediction = sum(k * x for k, x in zip(coefs, v))
        result.append((s, prediction, s["consumption"] - prediction))
    return result


R1, R2 = residuals(DATA, ["age"]), residuals(DATA, ["age", "base"])
print("region    n   residual mean (age)   residual mean (age + base)")
for b in sorted(REGION_MEAN, key=lambda b: -REGION_MEAN[b]):
    d1 = [a for s, _, a in R1 if s["region"] == b]
    d2 = [a for s, _, a in R2 if s["region"] == b]
    print(f"{b:<7} {len(d1):>3}   {sum(d1)/len(d1):>19.3f}   {sum(d2)/len(d2):>26.3f}")


def std_dev(d):
    o = sum(d) / len(d)
    return math.sqrt(sum((x - o) ** 2 for x in d) / len(d))


print("\nprediction band    n      residual mean   residual dispersion")
ordered = sorted(R2, key=lambda t: t[1])
step = len(ordered) // 4
for i in range(4):
    chunk = ordered[i * step:(i + 1) * step if i < 3 else len(ordered)]
    d = [a for _, _, a in chunk]
    print(f"{chunk[0][1]:>6.1f} - {chunk[-1][1]:>5.1f} {len(chunk):>5}   {sum(d)/len(d):>16.3f}"
          f"   {std_dev(d):>19.3f}")

# The same coefficient, with a quarter of the observations
r4 = rng(SEED + 88000)
QUARTER = [x for x in DATA if r4() < 0.25]
low3, high3 = coefficient_interval(QUARTER, ["age", "base"], 1)
k3 = ols(QUARTER, ["age", "base"])[0][1]
print(f"\n observations   coefficient   margin of variation      width")
print(f"{len(DATA):>13}   {k2[1]:>11.3f}   [{low2:.3f}, {high2:.3f}]        {high2 - low2:>8.3f}")
print(f"{len(QUARTER):>13}   {k3:>11.3f}   [{low3:.3f}, {high3:.3f}]        {high3 - low3:>8.3f}")

region    n   residual mean (age)   residual mean (age + base)
east    219                 2.481                       -0.229
central 231                 1.468                        0.295
north   363                 0.227                        0.031
south   258                -1.526                       -0.118
west    189                -3.020                        0.005

prediction band    n      residual mean   residual dispersion
   8.7 -  16.5   315              0.293                 3.891
  16.5 -  19.7   315             -0.168                 3.988
  19.7 -  23.3   315             -0.346                 4.852
  23.3 -  34.3   315              0.221                 4.893

 observations   coefficient   margin of variation      width
         1260         1.825   [1.629, 2.003]           0.374
          313         1.665   [1.341, 1.989]           0.648

The first table catches the flaw of the univariate model without looking at the coefficient. The regional mean of the residuals wanders between -3.020 and 2.481: the model systematically underestimates the east region and systematically overestimates the west region. A pattern that remains in a set of residuals by breakdown is the signature of a variable the model left out. When the confounder is added, the same numbers drop to between -0.229 and 0.295. The examination of residuals says which model is missing something before the result is seen.

The second table shows that an assumption is not met. The dispersion of the residuals grows together with the prediction level: 3.891 in the lowest band, 4.893 in the highest. The coefficient can still be computed in this case, but a single margin of variation is not valid for the whole prediction range; among high-consumption subscribers the margin is wider. Drawing with replacement carries this irregularity by itself, because it redraws the dataset as it is.

The third table gives what the margin depends on. When the same model is built with a quarter of the observations, the coefficient slips from 1.825 to 1.665 and the width of the margin of variation rises from 0.374 to 0.648. The observation count fell to a quarter, the width roughly doubled. Reporting a coefficient computed in a small breakdown without writing its margin is showing the two rows of this table as a single row.

Summary

• The univariate coefficient of meter age is -0.605 m³/year, the margin of variation [-0.714, -0.503] and the coefficient of determination 0.098; the number rests on 1,260 observations.
• When the base level of the region is added as a confounder, the same coefficient becomes 1.825 and its margin comes out as [1.629, 2.003]; the two intervals never overlap.
• The mechanism is visible in the breakdown: between regions the mean age rises while mean consumption falls, whereas in all five of the five regions the within-group slope is positive (between 0.771 and 2.416).
• While the coefficient of determination rose from 0.098 to 0.505, the coefficient changed sign; the coefficient of determination does not announce the direction, and when a completely random column is added it rises again from 0.505363 to 0.505377.
• The regional mean of the residuals wanders between -3.020 and 2.481 in the univariate model and drops to between -0.229 and 0.295 when the confounder is added; the dispersion of the residuals grows with the prediction from 3.891 to 4.893.
• With a quarter of the observations, the margin of the same coefficient grows in width from 0.374 to 0.648.

Next Step

This lesson showed that a coefficient changes direction by breakdown, but the breakdown was always the region and time never entered the work. Subscribers do not join the network at the same time and are not observed for the same length of time. The next lesson splits the dataset by entry period: the retention rate is computed, then by changing the entry window and the observation window it is counted how many different rates come out of the same dataset. The same change of direction appears here once more, this time not in a coefficient but in a rate: the aggregated rate shows one direction, the rate of every cohort shows the opposite.