--- title: 'Modeling Relationships with Regression' source: 'https://academia.sh/en/courses/statistical-analysis-and-exploratory-data-analysis/modeling-relationships-with-regression' course: 'Statistical Analysis and Exploratory Data Analysis' language: en updated: '2026-10-07T04:31:13+00:00' license: 'CC BY-SA 4.0' objectives: - 'Builds least squares from the normal equations and writes the resampling margin of variation next to a coefficient by drawing with replacement' - 'Shows that the coefficient changes sign when a confounder is added: the coefficient -0.605 becomes 1.825 and the two margins of variation do not overlap' - 'Shows with numbers that the coefficient of determination carries no information about the direction of the coefficient and rises with every added variable' - 'Reads the mean of the residuals in a breakdown and their dispersion against the prediction as an assumption check' --- # Modeling Relationships with Regression Least squares is written inside the lesson, and what a coefficient says is separated from what it does not say: the univariate coefficient of meter age comes out as -0.605 m³/year, the margin of variation [-0.714, -0.503] and the coefficient of determination 0.098. When the base level of the region is added as a confounder, the same coefficient turns to 1.825, its margin becomes [1.629, 2.003] and the two intervals never overlap. In each of the five regions the within-group slope is positive (between 0.771 and 2.416). The coefficient of determination rises from 0.098 to 0.505 but does not announce the sign; 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 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. ```python # 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. ```python # 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. ```python # 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.