Linear Regression

The previous course measured every step of the features in twelve lessons and left two accounts open at its close. The first: throughout all those measurements the model was always chosen from the same two or three families — a decision tree of depth 6, a gradient-based linear model, a distance-based procedure — and what the family itself contributes was never asked. The second: why the depth is six and why the number of neighbors is fifteen, and which set those numbers should be read from, was never measured. The same closing counted that the test set had been looked at at least twenty-four times.

This course takes up both questions together and changes the measure for that purpose. The features are held constant; the only thing that changes is the model. In every lesson three numbers stand side by side: the baseline, the number of the family on the held-out set, and the number of candidates tried to find that number. The third is new and is the rule of the course: a model’s number cannot be read without writing how many candidates were tried to find it. A number whose tuning budget is not written is an unmeasured number. This lesson establishes the triple and starts with the simplest family, least squares.

• LM1. The set is the fictional subscriber table of the Introduction to Machine Learning and Data Preparation and Feature Engineering courses; it is rebuilt with the same seed (20260218). The split is 756 train, 252 validation and 252 test rows and does not change throughout the course.
• LM2. In the classification lessons the target is the subscriber carrying an illicit-use or fault suspicion, and the baseline is saying the most frequent class: 0.7579 on the test set.
• LM3. In the regression lessons the target is the subscriber’s mean consumption over up to three periods of readings (m³) and the baseline is the mean of the train set: saying the same number to every subscriber.
• LM4. The regression feature is not derived from any reading that enters the computation of the target. Volatility, zero readings and per-person consumption come out of the same reading series, so they are not used; what remains is region, household size, satisfaction and period count.
• LM5. The sum of the five region indicators equals one in every row, that is, the constant term. One of them (north) is left out and becomes the reference.
• LM6. The measure is the root mean squared error (m³); the mean absolute error is printed next to it. Where the two diverge is left to evaluation-and-tuning/05.
• LM7. The coefficient of determination is defined relative to the train mean: zero means making the same error as the baseline. It can also come out negative on a held-out set.
• LM8. The causal reading of a coefficient is not the subject of this course. The sign flipping under a confounding variable was measured in the Applied Statistics topic of the Statistical Analysis and Exploratory Data Analysis course; it is not repeated here.
• LM9. The test set is looked at in each lesson only after the choice is finished, and how many times it was looked at is written.

Three Numbers

The setup produces the set, splits it and prints the two baselines.

# regression.py — MODEL: the FICTIONAL subscriber table of the Introduction to Machine Learning course is rebuilt with the same seed;
# the target is the subscriber's three-period mean consumption (m3).
import math

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


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, DATA = [], []
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])
    rating = 1 + pick(r(), [0.08, 0.14, 0.27, 0.34, 0.17])
    if r() >= 0.046:
        SUBSCRIBERS.append({"id": 10001 + i, "region": b, "household": household,
                            "satisfaction": rating})
for k in SUBSCRIBERS:
    r, v = rng(SEED + 7000 + k["id"]), []
    for n in range(pick(r(), [0.05, 0.12, 0.21, 0.62])):
        v.append(0.0 if r() < 0.038 else math.floor(
            k["region"][2] * math.exp((r() + r() + r() - 1.5) * 0.62)
            * (1 - n * 0.05) * 100 + 0.5) / 100)
        r()
    if not v:
        continue
    mean, vol, r = sum(v) / len(v), max(v) - min(v), rng(SEED + 51000 + k["id"])
    z = (0.052 * (mean - 20) + 0.85 * (min(v) == 0.0) + 0.026 * vol
         + 0.24 * (k["household"] >= 5) - 0.20 * (k["satisfaction"] >= 4) + k["region"][3]
         + 0.9 * (k["id"] - 10001) / 1399 + (r() + r() + r() - 1.5) * 1.30)
    x = {"id": k["id"], "region_name": k["region"][0], "household": k["household"],
         "satisfaction": k["satisfaction"], "period": len(v), "target": round(mean, 2),
         "suspicious": int(z > 1.35)}
    for b in REGION:
        x["area_" + b[0]] = int(k["region"][0] == b[0])
    DATA.append(x)


def shuffle(data, seed):
    r, s = rng(seed), list(range(len(data)))
    for i in range(len(s) - 1, 0, -1):
        j = int(r() * (i + 1))
        s[i], s[j] = s[j], s[i]
    return [data[i] for i in s]


K = shuffle(DATA, SEED + 90000)
TRAIN, VAL, TEST = K[:756], K[756:1008], K[1008:]
MEAN = sum(x["target"] for x in TRAIN) / len(TRAIN)


def rmse(predict, s):
    return math.sqrt(sum((predict(x) - x["target"]) ** 2 for x in s) / len(s))


def mae(predict, s):
    return sum(abs(predict(x) - x["target"]) for x in s) / len(s)


print(f"subscribers {len(DATA)}, train {len(TRAIN)}, validation {len(VAL)}, test {len(TEST)}")
print(f"classification baseline (most frequent class, on test) "
      f"{sum(x['suspicious'] == 0 for x in TEST) / len(TEST):.4f}")
print(f"regression baseline (train mean) {MEAN:.4f} m3\n")
print(f"{'split':<12} {'rmse':>9} {'mean abs err':>12}")
for name, s in (("train", TRAIN), ("validation", VAL), ("test", TEST)):
    print(f"{name:<12} {rmse(lambda x: MEAN, s):>9.4f} {mae(lambda x: MEAN, s):>12.4f}")

subscribers 1260, train 756, validation 252, test 252
classification baseline (most frequent class, on test) 0.7579
regression baseline (train mean) 19.7665 m3

split             rmse mean abs err
train           6.1995       4.9017
validation      6.4193       5.1148
test            6.5324       5.0806

Both baselines are non-learning procedures and both stay constant throughout this course. In classification the number 0.7579 is carried over as it is from the previous two courses. In regression the baseline says 19.7665 m³ to every subscriber and gives a root error of 6.5324 m³ on the test set: the zero point against which all the remaining numbers will be measured.

Least Squares Is a Prediction Procedure

Least squares solves for the coefficients that minimize the sum of squared errors. The solution is in closed form: a system of equations built from the products of the columns with one another is solved by elimination.

FIELDS = ["area_south", "area_east", "area_west", "area_central", "household",
          "satisfaction", "period"]


def solve(A, b):                            # elimination with partial pivoting
    n = len(b)
    for i in range(n):
        p = max(range(i, n), key=lambda k: abs(A[k][i]))
        A[i], A[p], b[i], b[p] = A[p], A[i], b[p], b[i]
        for k in range(i + 1, n):
            f = A[k][i] / A[i][i]
            for j in range(i, n):
                A[k][j] -= f * A[i][j]
            b[k] -= f * b[i]
    w = [0.0] * n
    for i in range(n - 1, -1, -1):
        w[i] = (b[i] - sum(A[i][j] * w[j] for j in range(i + 1, n))) / A[i][i]
    return w


def row(x, fields):
    return [1.0] + [float(x[a]) for a in fields]


def least_squares(s, fields):               # MODEL: linear regression
    X = [row(x, fields) for x in s]
    y = [x["target"] for x in s]
    p = len(X[0])
    A = [[sum(X[i][a] * X[i][b] for i in range(len(X))) for b in range(p)]
         for a in range(p)]
    return solve(A, [sum(X[i][a] * y[i] for i in range(len(X))) for a in range(p)])


def predictor(w, fields):
    return lambda x: sum(a * b for a, b in zip(w, row(x, fields)))


def determination(predict, s):              # coefficient of determination
    return 1 - (sum((predict(x) - x["target"]) ** 2 for x in s)
                / sum((MEAN - x["target"]) ** 2 for x in s))


W = least_squares(TRAIN, FIELDS)
f = predictor(W, FIELDS)
print(f"intercept   {W[0]:>8.3f}   (north is the reference)")
for i, a in enumerate(FIELDS):
    print(f"{a:<12} {W[i + 1]:>8.3f}")
print(f"\n{'split':<12} {'baseline rmse':>13} {'model rmse':>10} {'difference':>10} "
      f"{'mae':>8} {'determination':>13}")
for name, s in (("train", TRAIN), ("validation", VAL), ("test", TEST)):
    t = rmse(lambda x: MEAN, s)
    print(f"{name:<12} {t:>13.4f} {rmse(f, s):>10.4f} {t - rmse(f, s):>+10.4f} "
          f"{mae(f, s):>8.4f} {determination(f, s):>13.4f}")
print(f"\ncandidates tried 1")

intercept     21.820   (north is the reference)
area_south     -3.383
area_east       4.598
area_west      -5.851
area_central    2.340
household       0.062
satisfaction   -0.148
period         -0.550

split        baseline rmse model rmse difference      mae determination
train               6.1995     5.0757    +1.1238   3.7461        0.3297
validation          6.4193     5.3225    +1.0969   4.0005        0.3125
test                6.5324     5.3096    +1.2228   4.0020        0.3393

candidates tried 1

The triple is complete: baseline 6.5324, the family’s number 5.3096, candidates tried 1. The difference is 1.2228 m³, that is, 18.7% of the baseline. That the number of candidates is one is the privilege of this lesson — the column set was written in advance, there is no hyperparameter, and the solution came in closed form in a single pass. In the remaining lessons this number will never be one again.

The same coefficient list can be put to two different jobs, and the two are not the same thing. The region coefficients give the order of the base consumption levels in the setup correctly: east is the highest, west the lowest. But the measure of this lesson is not the correctness of a coefficient, it is the error made on a subscriber not seen before. The coefficient of determination on the test set is 0.3393 and is defined relative to the baseline: zero means making the same error as saying the train mean to every subscriber.

Assumptions Are Checked with Numbers

Linear regression has assumptions: the mean of the residuals is zero, their spread is independent of the prediction, the tails are symmetric. These are checked not as theory but as numbers.

def residuals(predict, s):
    return [predict(x) - x["target"] for x in s]


def resid_mean(v):                          # a residual at machine precision is zero
    m = sum(v) / len(v)
    return m if abs(m) > 5e-5 else 0.0


print(f"{'region':<9} {'train mean':>13} {'test n':>11} {'test mean':>13}")
for b in REGION:
    e = [x for x in TRAIN if x["region_name"] == b[0]]
    v = [x for x in TEST if x["region_name"] == b[0]]
    print(f"{b[0]:<9} {resid_mean(residuals(f, e)):>+13.4f} {len(v):>11} "
          f"{resid_mean(residuals(f, v)):>+13.4f}")
ordered = sorted(TEST, key=f)
n = len(ordered) // 4
print(f"\n{'slice':<7} {'mean prediction':>16} {'mean residual':>14} {'residual spread':>16}")
for k in range(4):
    d = ordered[k * n:(k + 1) * n] if k < 3 else ordered[3 * n:]
    v = residuals(f, d)
    m = sum(v) / len(v)
    sd = math.sqrt(sum((q - m) ** 2 for q in v) / len(v))
    print(f"slice {k + 1:<1} {sum(f(x) for x in d) / len(d):>16.2f} {m:>+14.3f} "
          f"{sd:>16.3f}")
v = residuals(f, TEST)
m = sum(v) / len(v)
sd = math.sqrt(sum((q - m) ** 2 for q in v) / len(v))
print(f"\nresiduals beyond two spreads on test "
      f"{sum(1 for q in v if abs(q - m) > 2 * sd)}, "
      f"negative {sum(1 for q in v if q - m < -2 * sd)}, "
      f"positive {sum(1 for q in v if q - m > 2 * sd)}")

region       train mean      test n     test mean
north           +0.0000          68       -1.0370
south           +0.0000          58       -0.0391
east            +0.0000          36       -1.0769
west            +0.0000          32       +1.3685
central         +0.0000          58       +0.5708

slice    mean prediction  mean residual  residual spread
slice 1            15.33         +1.192            3.566
slice 2            18.67         -1.007            4.746
slice 3            21.29         -0.237            4.883
slice 4            23.85         -0.498            7.137

residuals beyond two spreads on test 17, negative 11, positive 6

The first column is not a check but an identity. On the train set the mean residual of every region is exactly zero, and the reason is not that the model is good: the region indicators are inside the model as columns, and the least squares solution makes every column orthogonal to the residuals. This number can never come out different from zero, so it says nothing. The column that speaks is the third: the same means wander between -1.0769 and +1.3685 on the test set. On the thirty-two subscribers in the west region the model says 1.37 m³ too much on average.

The second table shows a sharper breakdown. The spread of the residuals doubles as the prediction grows, from 3.566 to 7.137. The model is far more accurate on low-consumption subscribers and far coarser on high-consumption ones. The reason the equal-spread assumption fails lies in the setup itself: consumption is generated from a multiplicative distribution, so the scatter of large values is also large. The last line gives the same asymmetry from the other direction — of the seventeen residuals beyond two spreads, eleven lie in one direction and six in the other.

Broken assumptions would invalidate the margin of play of the coefficients in an explanation study; the Applied Statistics topic of the Statistical Analysis and Exploratory Data Analysis course did that computation. For a prediction target the measure is single and has already been taken: 5.3096. The assumption check here does not invalidate the number, it says where the number is bad.

A Number Whose Tuning Budget Is Not Written Is an Unmeasured Number

In this lesson a single candidate was tried. What happens if all the non-empty subsets of the seven columns are tried?

import itertools

best_val, best_test, candidates = None, None, 0
for r in range(1, len(FIELDS) + 1):
    for c in itertools.combinations(FIELDS, r):
        candidates += 1
        g = predictor(least_squares(TRAIN, list(c)), list(c))
        d, s = rmse(g, VAL), rmse(g, TEST)
        if best_val is None or d < best_val[0]:
            best_val = (d, s, c)
        if best_test is None or s < best_test[1]:
            best_test = (d, s, c)
print(f"candidates tried {candidates}")
print(f"full model            validation {rmse(f, VAL):.4f}  test {rmse(f, TEST):.4f}")
print(f"chosen on validation  validation {best_val[0]:.4f}  "
      f"test {best_val[1]:.4f}  {len(best_val[2])} columns")
print(f"chosen on test        validation {best_test[0]:.4f}  "
      f"test {best_test[1]:.4f}  {len(best_test[2])} columns")
print(f"number reported if chosen on test {best_test[1]:.4f}, "
      f"honest {best_val[1]:.4f}, difference {best_val[1] - best_test[1]:.4f} m3")
print(f"excluded column: "
      f"{[a for a in FIELDS if a not in best_val[2]]}")

candidates tried 127
full model            validation 5.3225  test 5.3096
chosen on validation  validation 5.3073  test 5.3157  6 columns
chosen on test        validation 5.4150  test 5.2931  6 columns
number reported if chosen on test 5.2931, honest 5.3157, difference 0.0226 m3
excluded column: ['satisfaction']

Three rows give three different numbers, and all three come from the same model family, the same features, the same split. The only differences among them are how many candidates were tried and which set the choice was read from.

The choice read from the validation set drops the satisfaction column and beats the full model’s 5.3225 on validation with 5.3073. The same choice gives 5.3157 on test — that is, behind the full model. Choosing among a hundred and twenty-seven candidates by looking at validation produced a worse result than not choosing at all; searching 127 times for the smallest value on a set of 252 rows also selects that set’s own scatter.

The last two lines turn the rule of the course into a number. Had the choice been made on the test set, the reported number would have been 5.2931; the honest number is 5.3157 and the difference is 0.0226 m³. The sign of the difference is fixed — the number on the set where the choice was made can never come out worse. This is why the comparison between 5.3096 and 5.2931 carries meaning only when 1 and 127 are written next to them.

Summary

• The three numbers of the course were established: the baseline, the number of the family on the held-out set, and the number of candidates tried. The baseline is 0.7579 in classification and 6.5324 m³ in regression.
• Least squares gives 5.3096 m³ on the test set with seven columns; the difference is 1.2228 m³ and the number of candidates tried is 1. This is the only lesson in the course where the number of candidates is one.
• The mean residual is exactly zero per region on the train set, and this is not a check but the identity of the least squares solution; the same mean wanders between -1.0769 and +1.3685 on the test set.
• The equal-spread assumption does not hold: the residual spread rises from 3.566 to 7.137 as the prediction grows. For a prediction target this does not invalidate the number; it says where the number is bad.
• When the 127 subsets of the seven columns are tried, the subset chosen on validation falls behind the full model with 5.3157 on test; had it been chosen on test, 5.2931 would have been reported and the difference is 0.0226 m³.

Next Step

In this lesson the model was a plane and the columns were used as given. Yet the relation between household size and consumption need not be a line; the consumption that grows in going from a two-person household to a four-person one may not grow at the same rate in going from four to six. There is a way to build a nonlinear relation with a linear model, and it consists only of adding columns: the square, the cube, the pairwise products. The next lesson measures this expansion by degree — the degree is swept from 1 to 8, the train and held-out numbers of each degree are printed side by side, the choice is read from the validation set, and which degree would have been chosen had it been read from the test set is shown in the same table.