Threshold Selection and Curves

All four metrics of the previous lesson were read at the decision threshold of 0.5. That threshold does not come from the model: the model produces a probability, the threshold turns that probability into a decision, and it is given from outside. The only justification for 0.5 is that it sits in the middle of the two classes; in a set where the suspicious share is 0.2421, nothing says why that middle is the right place.

This lesson sweeps the threshold: the same model, the same test set, only the threshold changes. Precision and recall are recomputed at every threshold, the operating point is chosen with a cost calculation, and the two curves the sweep produces are printed as numbers.

• ET16. The model produces a probability; the threshold turns this probability into a decision and is independent of the model. Ninety-nine thresholds are tried, from 0.01 to 0.99; this is the tuning budget of this lesson.
• ET17. The cost is the one from the first lesson: 1 unit of field visit for every flagged subscriber, 6 units for every suspicious subscriber that is missed. The operating point is chosen from the validation set; the number on the test set is only reported.
• ET18. The training set, the features and the model stay fixed throughout the sweep.
• ET19. On the ROC curve the horizontal axis is the false positive rate and the vertical axis is recall; on the precision–recall curve the horizontal axis is recall and the vertical axis is precision. The curves are not drawn, they are printed as points.
• ET20. The area under the curve is computed with the trapezoid rule and observations with equal scores count as a single point.
• ET21. The thinned set is built by keeping one tenth of the suspicious subscribers in the test set; the training and the model do not change. The estimate made with the remaining seven suspicious subscribers is volatile and is read against its baseline.
• ET22. Choosing a threshold requires not that the probability be right but only that the ranking be right.

Sweeping the Threshold

As the threshold rises, the model flags fewer subscribers. It is right more often on the ones it flags, but the number of suspicious subscribers it misses grows.

# evaluation.py — MODEL: the FICTIONAL subscriber table is rebuilt with the same seed.
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, READINGS = [], {}
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 = rng(SEED + 7000 + k["id"])
    for n in range(pick(r(), [0.05, 0.12, 0.21, 0.62])):
        READINGS.setdefault(k["id"], []).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()

DATA = []
for k in SUBSCRIBERS:
    v = READINGS.get(k["id"])
    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 = {"household": k["household"], "satisfaction": k["satisfaction"], "period": len(v),
         "mean_consumption": round(mean, 2), "volatility": round(vol, 2),
         "zero_readings": int(min(v) == 0.0), "suspicious": int(z > 1.35)}
    x["per_person"] = round(x["mean_consumption"] / x["household"], 3)
    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:]
FIELDS = ["mean_consumption", "volatility", "household", "satisfaction", "period",
          "zero_readings", "per_person", "area_north", "area_south", "area_east", "area_west",
          "area_central"]


def impurity(rows):
    p = sum(x["suspicious"] for x in rows) / len(rows)
    return 2 * p * (1 - p)


def tree(rows, depth, fields, min_size=2):        # MODEL: decision tree
    n = len(rows)
    p = sum(x["suspicious"] for x in rows) / n
    best = None
    if depth and n >= 2 * min_size and 0.0 < p < 1.0:
        for f in fields:
            values = sorted({x[f] for x in rows})
            for v in (values[1:] if len(values) < 10 else
                      [values[int(i * len(values) / 10)] for i in range(1, 10)]):
                left = [x for x in rows if x[f] < v]
                right = [x for x in rows if x[f] >= v]
                if min(len(left), len(right)) < min_size:
                    continue
                gain = impurity(rows) - (len(left) * impurity(left)
                                         + len(right) * impurity(right)) / n
                if best is None or gain > best[0]:
                    best = (gain, f, v, left, right)
    if best is None or best[0] <= 1e-9:
        return {"p": p}
    return {"field": best[1], "threshold": best[2],
            "left": tree(best[3], depth - 1, fields, min_size),
            "right": tree(best[4], depth - 1, fields, min_size)}


def probability(node, x):
    while "p" not in node:
        node = node["left"] if x[node["field"]] < node["threshold"] else node["right"]
    return node["p"]


def matrix(score, data, threshold=0.5):           # tp, fp, fn, tn
    m = [0, 0, 0, 0]
    for x in data:
        m[(0 if x["suspicious"] else 1) + (0 if score(x) >= threshold else 2)] += 1
    return m


def linear(train, fields, w, epochs=300, step=0.5):   # MODEL: gradient-based linear model
    mu = {f: sum(x[f] for x in train) / len(train) for f in fields}
    sd = {f: (sum((x[f] - mu[f]) ** 2 for x in train) / len(train)) ** 0.5 or 1.0
          for f in fields}
    scale = lambda x: [(x[f] - mu[f]) / sd[f] for f in fields]
    coef, b = [0.0] * len(fields), 0.0
    X = [scale(x) for x in train]
    y = [x["suspicious"] for x in train]
    sw = [(w if x["suspicious"] else 1.0) for x in train]
    total = sum(sw)
    for _ in range(epochs):
        gc, gb = [0.0] * len(fields), 0.0
        for i, xi in enumerate(X):
            z = b + sum(coef[j] * xi[j] for j in range(len(coef)))
            e = (1 / (1 + math.exp(-max(-30.0, min(30.0, z)))) - y[i]) * sw[i]
            gb += e
            for j in range(len(coef)):
                gc[j] += e * xi[j]
        b -= step * gb / total
        for j in range(len(coef)):
            coef[j] -= step * gc[j] / total

    def score(x):
        v = scale(x)
        z = b + sum(coef[j] * v[j] for j in range(len(coef)))
        return 1 / (1 + math.exp(-max(-30.0, min(30.0, z))))
    return score


N_SUSP = sum(x["suspicious"] for x in TEST)
SCORE_FN = linear(TRAIN, FIELDS, 1.0)        # MODEL: probability-producing gradient-based linear model
VISIT_COST, MISS_COST = 1.0, 6.0


def point(score, rows, e):
    tp, fp, fn, tn = matrix(score, rows, e)
    return {"flagged": tp + fp, "tp": tp, "fpr": fp / (fp + tn), "recall": tp / (tp + fn),
            "precision": tp / (tp + fp) if tp + fp else float("nan"),
            "cost": VISIT_COST * (tp + fp) + MISS_COST * fn}


print(f"{'threshold':>9}{'flagged':>9}{'caught':>10}{'fpr':>9}"
      f"{'recall':>11}{'precision':>10}{'cost':>7}")
for e in (0.10, 0.20, 0.30, 0.40, 0.50, 0.60, 0.70, 0.80, 0.90):
    n = point(SCORE_FN, TEST, e)
    print(f"{e:>9.2f}{n['flagged']:>9}{n['tp']:>10}{n['fpr']:>9.4f}"
          f"{n['recall']:>11.4f}{n['precision']:>10.4f}{n['cost']:>7.0f}")

threshold  flagged    caught      fpr     recall precision   cost
     0.10      158        57   0.5288     0.9344    0.3608    182
     0.20      100        45   0.2880     0.7377    0.4500    196
     0.30       66        34   0.1675     0.5574    0.5152    228
     0.40       47        29   0.0942     0.4754    0.6170    239
     0.50       33        22   0.0576     0.3607    0.6667    267
     0.60       24        18   0.0314     0.2951    0.7500    282
     0.70       17        12   0.0262     0.1967    0.7059    311
     0.80        8         8   0.0000     0.1311    1.0000    326
     0.90        2         2   0.0000     0.0328    1.0000    356

The two columns move in opposite directions. At a threshold of 0.10 the model flags 158 subscribers, catches 57 of the 61 suspicious ones, and its precision is 0.3608. At a threshold of 0.90 it flags only 2 subscribers and is right on both, but recall is 0.0328. Precision does not rise smoothly: it is 0.7500 at 0.60 and falls to 0.7059 at 0.70, because its denominator shrinks along with the model. Recall, on the other hand, is a function that does not decrease with the threshold; its denominator is fixed.

The last column is the one that decides: the cost is 182 at a threshold of 0.10 and 356 at 0.90. The 174 units in between are entirely a product of the threshold; the model, the features and the training set did not change at all.

The Operating Point Is Chosen with a Cost Calculation

The threshold is a hyperparameter and the rule of the course applies here too: a choice does not count as made until it is written which set it was read from.

THRESHOLDS = [i / 100 for i in range(1, 100)]
best_val = min(THRESHOLDS, key=lambda e: point(SCORE_FN, VAL, e)["cost"])
best_test = min(THRESHOLDS, key=lambda e: point(SCORE_FN, TEST, e)["cost"])
print(f"threshold candidates tried {len(THRESHOLDS)}, validation set {len(VAL)} subscribers")
print(f"{'operating point':<22}{'threshold':>9}{'flagged':>9}{'caught':>10}"
      f"{'validation':>11}{'test':>8}")
for name, e in (("default", 0.50), ("chosen on validation", best_val),
                ("chosen on test", best_test)):
    n = point(SCORE_FN, TEST, e)
    print(f"{name:<22}{e:>9.2f}{n['flagged']:>9}{n['tp']:>10}"
          f"{point(SCORE_FN, VAL, e)['cost']:>11.0f}{n['cost']:>8.0f}")

threshold candidates tried 99, validation set 252 subscribers
operating point       threshold  flagged    caught validation    test
default                    0.50       33        22        205     267
chosen on validation       0.16      121        51        160     181
chosen on test             0.13      138        55        173     174

The cheapest of the ninety-nine candidates on the validation set is the threshold of 0.16, and it costs 160 units on that set; the same threshold gives 181 units on the test set. The default threshold of 0.50 gives 267 units: choosing the threshold earns 86 units, a gain close to a third, and this gain comes not from the model but from how its output is read.

The third row measures the source of the choice. The smallest cost on the test set is at a threshold of 0.13, 174 units. The 7 units in between are the extra share by which someone who reads the choice off the test set looks better — a small number because the cost curve is flat around the threshold, but its sign is always in the same direction.

Two Curves and Imbalance

The sweep produces two curves. The ROC curve shows recall against the false positive rate; the precision–recall curve shows precision against recall. They differ when the suspicious share changes.

def curve(score, rows):               # ROC and precision-recall points
    p = sorted(((score(x), x["suspicious"]) for x in rows), key=lambda t: -t[0])
    P = sum(y for _, y in p)
    N = len(p) - P
    tp = fp = i = 0
    pts = [(0.0, 0.0, 1.0)]
    while i < len(p):
        j = i
        while j < len(p) and p[j][0] == p[i][0]:
            j += 1
        d = sum(y for _, y in p[i:j])
        tp, fp, i = tp + d, fp + (j - i) - d, j
        pts.append((fp / N, tp / P, tp / (tp + fp)))
    return pts


def area_roc(pts, limit=1.0):         # trapezoid rule; limit clips the false positive rate
    return sum((min(b[0], limit) - a[0]) * (a[1] + b[1]) / 2
               for a, b in zip(pts, pts[1:]) if a[0] < limit) / limit


def area_pr(pts):                     # sum of steps on the recall axis
    return sum((b[1] - a[1]) * b[2] for a, b in zip(pts, pts[1:]))


THIN = [x for i, x in enumerate(TEST) if not x["suspicious"] or i % 10 == 0]
constant = lambda x: 0.5
print(f"{'test set':<18}{'rows':>6}{'suspicious':>11}{'ratio':>7}{'ROC':>8}{'PR':>8}"
      f"{'constant ROC':>13}{'constant PR':>12}")
for name, v in (("real distribution", TEST), ("thinned", THIN)):
    n, s, ns = curve(SCORE_FN, v), curve(constant, v), sum(x["suspicious"] for x in v)
    print(f"{name:<18}{len(v):>6}{ns:>11}{ns / len(v):>7.4f}{area_roc(n):>8.4f}"
          f"{area_pr(n):>8.4f}{area_roc(s):>13.4f}{area_pr(s):>12.4f}")
for h in (0.50, 0.70):
    r = [point(SCORE_FN, v, max(e for e in THRESHOLDS
                                if point(SCORE_FN, v, e)["recall"] >= h))["precision"]
         for v in (TEST, THIN)]
    print(f"precision at recall above {h:.2f}: real distribution {r[0]:.4f}, "
          f"thinned {r[1]:.4f}")

test set            rows suspicious  ratio     ROC      PR constant ROC constant PR
real distribution    252         61 0.2421  0.8058  0.6212       0.5000      0.2421
thinned              198          7 0.0354  0.8624  0.4124       0.5000      0.0354
precision at recall above 0.50: real distribution 0.5161, thinned 0.2500
precision at recall above 0.70: real distribution 0.4624, thinned 0.1667

The two curves move in opposite directions. When the suspicious share falls from 0.2421 to 0.0354, the area under the ROC curve rises from 0.8058 to 0.8624, while that of the precision–recall curve falls from 0.6212 to 0.4124. Part of the rise is the volatility of an estimate made with seven suspicious subscribers; the real matter is in the last two columns. The ROC curve of the constant-score procedure gives 0.5000 on both sets; the baseline of this curve does not look at the suspicious share. The baseline of the other curve, however, is directly the suspicious share.

The last two lines are the counterpart of this in operation. When recall is held above 0.50, precision is 0.5161 on the real distribution and 0.2500 on the thinned set; when it is held above 0.70, 0.4624 against 0.1667. While the area reads 0.86, three of every four subscribers flagged to catch half of the suspicious ones go to waste. On an imbalanced set, the misleading curve is the one whose baseline does not move together with the imbalance.

What the Single Number Hides

The area under the curve collects all thresholds into a single number and drops two things while collecting: which threshold the work is done at and which region of the curve is useful.

TREE_D3 = tree(TRAIN, 3, FIELDS)                       # MODEL: depth 3 decision tree
ONE_FIELD = linear(TRAIN, ["mean_consumption"], 1.0)   # MODEL: single-feature gradient-based linear model
print(f"{'model':<25}{'ROC':>8}{'fpr<0.10':>9}{'PR':>8}{'best cost':>14}"
      f"{'distinct scores':>16}")
for name, s in (("linear, 12 features", SCORE_FN),
                ("tree d3, 12 features", lambda x: probability(TREE_D3, x)),
                ("linear, mean_consumption", ONE_FIELD)):
    n = curve(s, TEST)
    print(f"{name:<25}{area_roc(n):>8.4f}{area_roc(n, 0.10):>9.4f}{area_pr(n):>8.4f}"
          f"{min(point(s, TEST, e)['cost'] for e in THRESHOLDS):>14.0f}"
          f"{len({round(s(x), 12) for x in TEST}):>16}")

model                         ROC fpr<0.10      PR     best cost distinct scores
linear, 12 features        0.8058   0.3364  0.6212           174             252
tree d3, 12 features       0.7222   0.2553  0.4485           197               8
linear, mean_consumption   0.7138   0.3008  0.5589           213             241

The real reading is in the last two rows. The areas under the ROC curve of the two models are 0.7222 and 0.7138, with 0.0084 between them — by this measure they are almost equal. The other three columns do not agree: the single-feature model is 0.1104 points ahead on the area under the precision–recall curve, 0.0455 points ahead in the region where the false positive rate stays below 0.10, and 16 units behind on the best cost. Three readings, three different winners.

The depth 3 tree, which produces eight distinct scores, and the single-feature model, which produces 241 distinct scores, move on the curve in entirely different ways: one in eight large steps, the other in hundreds of small ones. The area squeezes both into the same number. One more item is added to the rule of the course here: without writing which threshold will be used, the area under the curve cannot be turned into an operating decision.

Summary

• As the threshold rises, recall falls in a non-decreasing way, but precision does not rise smoothly: precision of 0.7500 at a threshold of 0.60 falls to 0.7059 at 0.70. Ninety-nine threshold candidates move the cost between 182 and 356; the model and the training set never change.
• The threshold of 0.16 chosen on validation gives 181 units on the test set, 86 units cheaper than the default threshold of 0.50; if the same choice had been read from the test set, it would have looked 7 units better.
• When the suspicious share falls from 0.2421 to 0.0354, the area under the ROC curve rises from 0.8058 to 0.8624 and the area under the precision–recall curve falls from 0.6212 to 0.4124; the baseline of the first curve is 0.5000 on both sets, that of the second is the suspicious share itself.
• Two models counted as equal by the area under the ROC curve with a difference of 0.0084 diverge by 0.1104 on the area under the precision–recall curve and by 16 units on the best cost.

Next Step

This whole lesson rested on the ranking. The threshold was swept, the curves were printed, the operating point was chosen; at no point did the number the model says have to be right itself. If all the scores were passed through an increasing function the ranking would not change, the two curves would stay the same, and the operating point would flag the same subscribers. Yet when the model says 0.80 for a subscriber, that number should have a counterpart: about eighty percent of the subscribers given 0.80 should turn out to be suspicious. The next lesson prints the stated probability and the realized share side by side, bin by bin, shows that a model whose ranking is good and whose calibration is bad exists, and writes the correction as a single-variable mapping learned only from the validation set.