Question answered: when the model says 10%, does the event happen about 10% of the time?
The ROC page ended with a warning: AUC only measures ranking. A model can rank perfectly and still predict probabilities that are far off. Whenever the probabilities themselves are used, for a premium, an expected cost, a capital estimate or the cost-based threshold, they must also be calibrated.
A model is calibrated if, among all cases that receive prediction \(\hat p\), the event happens with frequency \(\hat p\):
\[ P(Y = 1 \mid \hat p = p) = p \quad \text{for all } p. \]
Calibration and discrimination are separate properties:
A good model needs both.
The simplest check: is the average prediction equal to the observed rate?
pd.DataFrame({
"Average prediction": {m: b[m].mean() for m in ("Null", "GLM", "GBM")},
"Observed rate": {m: rate for m in ("Null", "GLM", "GBM")},
}).round(4)| Average prediction | Observed rate | |
|---|---|---|
| Null | 0.0368 | 0.0368 |
| GLM | 0.0366 | 0.0368 |
| GBM | 0.0367 | 0.0368 |
All three are close to the observed 3.68%. This check is necessary but far from sufficient: a model can be right on average while too high for low risks and too low for high risks, and these errors cancel out.
A reliability diagram plots the observed rate against the predicted probability. A calibrated model follows the diagonal. Points above the diagonal mean the model predicts too low there, points below mean too high.
The usual recipe sorts the cases by prediction, splits them into bins (here 20 bins with equal numbers of policies), and plots each bin’s average prediction against its observed claim rate:
fig, axes = plt.subplots(1, 2, figsize=(9.5, 4.5))
for ax, m in zip(axes, ("GLM", "GBM")):
plot_reliability(y, {m: b[m]}, ax=ax, method="bins", title=m)
fig.tight_layout()
Both models stay close to the diagonal over most of the range. Most policies have a predicted probability below 10%, so the bins are packed there and sparse above. The few high-risk bins carry wide confidence intervals.
The picture depends on how the bins are chosen, and there’s no right choice. The same GBM with different binnings:
fig, axes = plt.subplots(1, 3, figsize=(12, 4.2))
for ax, (n, strategy) in zip(axes, [(10, "uniform"), (10, "quantile"), (100, "quantile")]):
plot_reliability(y, {"GBM": b["GBM"]}, ax=ax, method="bins", n_bins=n, strategy=strategy,
max_p=0.35, title=f"{n} {'equal-width' if strategy == 'uniform' else 'equal-count'} bins")
ax.get_legend().remove()
fig.tight_layout()
Equal-width bins put nearly all policies in the first two bins and leave the rest nearly empty. Few equal-count bins hide any detail at the high end. Many bins turn the diagram into noise. Summary statistics built on bins inherit the problem. The expected calibration error (ECE), the share-weighted average gap between bin prediction and bin rate, changes with the binning:
settings = [(10, "uniform"), (20, "uniform"), (5, "quantile"), (10, "quantile"), (20, "quantile"), (100, "quantile")]
pd.DataFrame({
f"{n} {'equal-width' if s == 'uniform' else 'equal-count'}": {
m: expected_calibration_error(y, b[m], n, s) for m in ("GLM", "GBM")
}
for n, s in settings
}).T.round(5)| GLM | GBM | |
|---|---|---|
| 10 equal-width | 0.00172 | 0.00089 |
| 20 equal-width | 0.00262 | 0.00101 |
| 5 equal-count | 0.00192 | 0.00051 |
| 10 equal-count | 0.00306 | 0.00187 |
| 20 equal-count | 0.00313 | 0.00215 |
| 100 equal-count | 0.00437 | 0.00342 |
With 100 bins, the GBM’s ECE is 7 times larger than with 5 bins. That’s not because the model got worse: smaller bins simply contain more sampling noise. ECE values are only comparable at the same binning, and even then mostly measure noise when bins are small.
The CORP approach (consistent, optimally binned, reproducible; Dimitriadis, Gneiting and Jordan, 2021) replaces hand-chosen bins with isotonic regression. It finds the best non-decreasing step function mapping predictions to observed outcomes, using the pool-adjacent-violators (PAV) algorithm. In effect, the data chooses the bins: adjacent predictions are merged exactly where their outcomes don’t support a higher rate.
fig, axes = plt.subplots(1, 2, figsize=(9.5, 4.5))
for ax, m in zip(axes, ("GLM", "GBM")):
plot_reliability(y, {m: b[m]}, ax=ax, method="corp", band=True, title=m)
fig.tight_layout()
Both curves follow the diagonal closely where most policies are, below about 6%. Above that, both tend to run above the diagonal, and the GLM’s curve leaves its band several times. Both models under-predict their highest risks, the GLM more clearly. The segment analysis below shows where this comes from. In practice, the CORP diagram should be the default, and binned diagrams kept for audiences used to them.
Two numbers summarize the most common patterns of miscalibration. Fit a logistic regression of the outcome on the model’s logit,
\[ \operatorname{logit} P(Y = 1) = a + b \cdot \operatorname{logit}(\hat p), \]
and read off:
| Statistic | Ideal | Meaning when off |
|---|---|---|
| Calibration intercept \(a\) (with \(b\) fixed at 1) | 0 | > 0: predictions too low on average; < 0: too high |
| Calibration slope \(b\) | 1 | < 1: predictions too extreme (too spread out, typical of overfitting); > 1: too timid (not spread out enough, typical of over-regularization) |
To see these patterns, we distort the GBM’s predictions in three ways:
center = logit(rate)
variants = {
"GBM": b["GBM"],
"too low (× 0.6)": 0.6 * b["GBM"],
"too extreme": expit(center + 1.6 * (logit(b["GBM"]) - center)),
"too timid": expit(center + 0.6 * (logit(b["GBM"]) - center)),
}fig, axes = plt.subplots(1, 3, figsize=(12, 4.2))
for ax, name in zip(axes, list(variants)[1:]):
plot_reliability(y, {name: variants[name]}, ax=ax, max_p=0.25, title=name,
colors={name: "#eb6834"})
ax.get_legend().remove()
fig.tight_layout()
rows = {m: b[m] for m in ("GLM", "GBM")} | {k: v for k, v in variants.items() if k != "GBM"}
pd.DataFrame({
k: pd.concat([calibration_intercept_slope(y, v), pd.Series({"AUC": auc_rank(y, v)})])
for k, v in rows.items()
}).T.round(3)| calibration intercept | calibration slope | AUC | |
|---|---|---|---|
| GLM | 0.004 | 0.987 | 0.694 |
| GBM | 0.004 | 1.025 | 0.713 |
| too low (× 0.6) | 0.544 | 1.063 | 0.713 |
| too extreme | -0.242 | 0.641 | 0.713 |
| too timid | 0.073 | 1.711 | 0.713 |
Calibration plots diagnose. To compare models with a single number that rewards both calibration and discrimination, use a proper scoring rule. It’s a score whose expected value is best when the predicted probability equals the true probability, so no forecaster can improve their expected score by reporting anything other than what they believe. The two standard ones, lower is better for both:
\[ \text{Brier score} = \frac{1}{n}\sum_{i} (\hat p_i - y_i)^2, \qquad \text{log loss} = -\frac{1}{n}\sum_{i} \big[ y_i \log \hat p_i + (1 - y_i)\log(1 - \hat p_i) \big]. \]
Accuracy, F1 and AUC are not proper scoring rules. Accuracy and F1 only see which side of a threshold a prediction falls on, and AUC only sees ranking. A model can improve on them while its probabilities get worse.
A raw score is hard to interpret on its own. The CORP decomposition splits any proper score exactly into three parts:
\[ \text{score} = \underbrace{\text{MCB}}_{\text{miscalibration}} - \underbrace{\text{DSC}}_{\text{discrimination}} + \underbrace{\text{UNC}}_{\text{uncertainty}}. \]
cases = {m: b[m] for m in ("Null", "GLM", "GBM")} | {"GBM / 10": b["GBM"] / 10, "GBM ^ 0.2": b["GBM"] ** 0.2}
murphy = pd.DataFrame({k: murphy_decomposition(y, v, "log loss") for k, v in cases.items()}).T
murphy.round(5)| log loss | MCB | DSC | UNC | |
|---|---|---|---|---|
| Null | 0.15760 | 0.00000 | 0.00000 | 0.1576 |
| GLM | 0.14822 | 0.00044 | 0.00983 | 0.1576 |
| GBM | 0.14502 | 0.00033 | 0.01291 | 0.1576 |
| GBM / 10 | 0.19791 | 0.05321 | 0.01291 | 0.1576 |
| GBM ^ 0.2 | 0.68023 | 0.53554 | 0.01291 | 0.1576 |
The table answers why one model scores better than another:
The same table for the Brier score (try murphy_decomposition(y, p, "Brier")) tells the same story. Use the decomposition whenever two models are compared on a scoring rule: it separates ‘ranks better’ from ‘is better calibrated’.
Overall calibration can hide miscalibration in subgroups, and subgroups are often what matters, such as tariff cells, customer segments or regions. Insurers check this with the actual-to-expected (A/E) ratio: observed claims divided by predicted claims, per segment. A/E above 1 means the model predicts too low for that segment.
test = train_test().test
segment = pd.cut(test["BonusMalus"], [0, 50, 60, 80, 100, 250],
labels=["50 (best)", "51–60", "61–80", "81–100", "> 100 (malus)"])
df = b.assign(segment=segment.to_numpy())
grouped = df.groupby("segment", observed=True)
ae = pd.DataFrame({
"Policies": grouped.size(),
"Observed rate": grouped["y"].mean(),
"A/E GLM": grouped["y"].sum() / grouped["GLM"].sum(),
"A/E GBM": grouped["y"].sum() / grouped["GBM"].sum(),
})
ae.round(3)| Policies | Observed rate | A/E GLM | A/E GBM | |
|---|---|---|---|---|
| segment | ||||
| 50 (best) | 95915 | 0.029 | 0.960 | 1.014 |
| 51–60 | 23713 | 0.036 | 1.189 | 0.982 |
| 61–80 | 29623 | 0.047 | 1.104 | 1.013 |
| 81–100 | 18305 | 0.051 | 0.832 | 0.970 |
| > 100 (malus) | 1948 | 0.143 | 1.305 | 1.034 |
The GLM is calibrated overall but not across bonus-malus levels: its A/E ratios range from 0.83 to 1.31. For policyholders with a malus, it predicts 31% fewer claims than actually occur, while in the 81–100 band it predicts too many. The errors cancel out in the overall average. The GLM uses a single log-linear effect for bonus-malus, which is too rigid for the real relationship. The GBM stays within a few percent in every segment. A table like this points directly to the feature that needs a better specification.
If ranking is good but probabilities are off, a model can be recalibrated by mapping its predictions through a fitted function:
Fit the recalibration on data that wasn’t used to train the model, and evaluate on yet another set (scikit-learn: CalibratedClassifierCV). A monotone recalibration leaves the ranking, and therefore AUC and DSC, essentially unchanged. It only removes MCB.