| SOCR ≫ | DSPA ≫ | DSPA3 Topics ≫ |
DSPA3: Data Science and Predictive Analytics (UMich HS650)
Chapter 9: Model Performance Assessment, Validation, and
Improvement
Chapter 9: Model Performance Assessment, Validation, and Improvement
SOCR/MCAIM (Ivo Dinov)
SOCR/MCAIM (Ivo Dinov)
September 04, 2026
library(ggplot2)
library(dplyr)
library(tidyr)
library(patchwork)
library(plotly) # interactive figures and ALL 3-D graphics
library(rsample) # splitting and resampling
library(yardstick) # metrics
library(caret) # confusionMatrix, train
library(pROC)How this chapter uses graphics
Every two-dimensional figure is drawn with
ggplot2and rendered statically. Immediately after each one, the equivalentplot_ly()code appears in a chunk markedeval=FALSE, echo=TRUE.Every three-dimensional figure is drawn with
plot_ly()and evaluated. Nearly every quantity in this chapter is a function of two variables, optimism of \((n, p)\), the cross-validation estimator’s variance of \((k, n)\), net benefit of (threshold, prevalence), the learning curve of (sample size, complexity). Each is a surface, and reading one from a fixed viewpoint discards the interaction that makes it interesting.
1 Learning objectives
After completing this chapter you will be able to:
- Define generalization error, conditional test error, and training error, and explain which one resampling actually estimates.
- Derive the expected optimism of the training error and connect it to \(C_p\), AIC, and effective degrees of freedom.
- Distinguish Type I/II errors from the false discovery and false omission rates by what each conditions on.
- Explain why Cohen’s \(\kappa\) can be low despite high agreement, and demonstrate both kappa paradoxes.
- Choose between ROC and precision–recall curves from the class balance, and state the chance level of each.
- Assess calibration with a reliability diagram and a proper scoring rule, and explain why AUC cannot detect miscalibration.
- Compute net benefit and read a decision curve.
- Analyze the bias–variance tradeoff of the cross-validation estimator in \(k\), and use the closed-form LOOCV shortcut for linear smoothers.
- Implement the .632 and .632+ bootstrap estimators correctly.
- Design a nested, grouped, or rolling-origin resampling scheme, and compare models with a test that accounts for correlated folds.
Estimated time: 11–14 hours including exercises. Prerequisites: Chapter 5 (the metric family and the leakage taxonomy of §5.4.1), Chapter 6 (regression metrics and temporal validation, §6.9), and Chapter 8 (the adjusted Rand index, §8.4.2). This chapter develops all three properly.
1.1 Chapter 9 Live Webapps
The DSPA3 Labs include a number of interactive webapps, which are build as portable, self-contained, and browser accessible applications. These apps demonstrate various concepts discussed in this chapter, including
2 PART I: WHAT ARE WE ESTIMATING?
Every number in this chapter estimates something. Before computing any of them, it is worth writing down precisely what.
3 Three different errors
Let \((X,Y)\sim P\) and let \(\hat f\) be a model fitted to a training sample \(\mathcal{T}=\{(x_i,y_i)\}_{i=1}^{n}\) drawn from \(P\). Given a loss \(L\), three distinct quantities are routinely conflated.
\[ \begin{aligned} \textbf{Training (apparent) error: }\quad &\overline{\mathrm{err}}=\frac{1}{n}\sum_{i=1}^{n}L\big(y_i,\hat f(x_i)\big)\\[2mm] \textbf{Conditional test error: }\quad &\mathrm{Err}_{\mathcal{T}}=\mathbb{E}_{(X,Y)\sim P}\Big[L\big(Y,\hat f(X)\big)\ \Big|\ \mathcal{T}\Big]\\[2mm] \textbf{Expected test error: }\quad &\mathrm{Err}=\mathbb{E}_{\mathcal{T}}\big[\mathrm{Err}_{\mathcal{T}}\big] \end{aligned} \]
\(\overline{\mathrm{err}}\) scores the model on the data that produced it. \(\mathrm{Err}_{\mathcal{T}}\) is what you actually want: the error of this particular fitted model on new data. \(\mathrm{Err}\) averages that over the randomness in the training sample too.
Common misconception: “cross-validation estimates the error of my model.” It does not, quite. Cross-validation estimates \(\mathrm{Err}\), the error of the fitting procedure, averaged over training sets, rather than \(\mathrm{Err}_{\mathcal{T}}\), the error of the one model you actually have (Hastie, Tibshirani & Friedman, ESL, §7.12). Each fold fits a different model on a different subset, and their errors are averaged.
The distinction matters when it is the deployed model whose performance you must certify. For that, a genuinely held-out set is the only direct estimate, and cross-validation is the tool for choosing the procedure rather than for certifying the artifact.
set.seed(11)
sim_errors <- function(n = 60, p = 8, sigma = 1, n_big = 20000, reps = 300) {
beta <- c(2, -1.5, 1, rep(0, p - 3))
out <- replicate(reps, {
X <- matrix(rnorm(n * p), n, p)
y <- as.vector(X %*% beta) + rnorm(n, sd = sigma)
fit <- lm(y ~ X)
# Training error on the fitting sample
tr <- mean(residuals(fit)^2)
# Conditional test error: THIS model on a very large fresh sample
Xb <- matrix(rnorm(n_big * p), n_big, p)
yb <- as.vector(Xb %*% beta) + rnorm(n_big, sd = sigma)
te <- mean((yb - cbind(1, Xb) %*% coef(fit))^2)
c(train = tr, cond_test = te)
})
c(mean_training_error = mean(out["train", ]),
mean_conditional_test_error = mean(out["cond_test", ]),
sd_of_conditional_test_error = sd(out["cond_test", ]),
irreducible = sigma^2)
}
round(sim_errors(), 4)#> mean_training_error mean_conditional_test_error
#> 0.8719 1.1780
#> sd_of_conditional_test_error irreducible
#> 0.0815 1.0000Three facts to read off. The training error is below the irreducible noise \(\sigma^2\), impossible for genuine predictive error, and the signature of fitting the noise. The conditional test error is above it. And that conditional error has real spread across training samples, which is what \(\mathrm{Err}\) averages away.
4 Optimism
The gap between \(\overline{\mathrm{err}}\) and the test error has a name and a closed form.
Define the in-sample error, the error at the same \(x_i\) but with fresh responses \(Y_i^{\text{new}}\) drawn from \(P(Y\mid x_i)\):
\[\mathrm{Err}_{\mathrm{in}}=\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}_{Y^{\text{new}}}\Big[L\big(Y_i^{\text{new}},\hat f(x_i)\big)\ \Big|\ \mathcal{T}\Big].\]
The optimism is \(\mathrm{op}=\mathrm{Err}_{\mathrm{in}}-\overline{\mathrm{err}}\), and its expectation has a remarkably clean form.
Expected optimism. For squared-error loss (and, more generally, for log-likelihood and 0–1 loss), \[\boxed{\;\mathbb{E}_{\mathbf y}\big[\mathrm{op}\big]=\frac{2}{n}\sum_{i=1}^{n}\operatorname{Cov}\big(\hat y_i,\,y_i\big)\;}\]
Sketch. Expand \(\mathbb{E}[(y_i^{\text{new}}-\hat y_i)^2]-\mathbb{E}[(y_i-\hat y_i)^2]\). The squared terms in \(\hat y_i\) cancel, \(y_i^{\text{new}}\) is independent of \(\hat y_i\) so its cross term vanishes, and what remains is \(2\big(\mathbb{E}[y_i\hat y_i]-\mathbb{E}[y_i]\mathbb{E}[\hat y_i]\big)=2\operatorname{Cov}(\hat y_i,y_i)\). \(\blacksquare\)
The optimism is exactly twice the covariance between each fitted value and its own response, summed and scaled. The harder a model works to match its own \(y_i\), the more optimistic its training error.
For a linear fit with \(d\) effective parameters and homoscedastic noise \(\sigma^2\), \(\sum_i\operatorname{Cov}(\hat y_i,y_i)=d\sigma^2\), giving
\[\boxed{\;\mathbb{E}\big[\mathrm{Err}_{\mathrm{in}}\big]=\mathbb{E}\big[\overline{\mathrm{err}}\big]+\frac{2d\sigma^2}{n}\;}\]
set.seed(21)
verify_optimism <- function(n, p, sigma = 1, reps = 4000) {
beta <- c(2, -1.5, 1, rep(0, max(0, p - 3)))[1:p]
res <- replicate(reps, {
X <- matrix(rnorm(n * p), n, p)
y <- as.vector(X %*% beta) + rnorm(n, sd = sigma)
H <- cbind(1, X); fit <- lm.fit(H, y)
yhat <- fit$fitted.values
ynew <- as.vector(X %*% beta) + rnorm(n, sd = sigma) # SAME x, fresh y
c(err_bar = mean((y - yhat)^2), err_in = mean((ynew - yhat)^2))
})
d <- p + 1 # intercept counts
c(n = n, p = p,
observed_optimism = mean(res["err_in", ] - res["err_bar", ]),
theoretical_2d_sigma2_over_n = 2 * d * sigma^2 / n)
}
as.data.frame(do.call(rbind, list(
verify_optimism(50, 5), verify_optimism(50, 15),
verify_optimism(200, 5), verify_optimism(200, 15)))) |>
mutate(across(where(is.numeric), \(z) round(z, 4)))#> n p observed_optimism theoretical_2d_sigma2_over_n
#> 1 50 5 0.2388 0.24
#> 2 50 15 0.6461 0.64
#> 3 200 5 0.0638 0.06
#> 4 200 15 0.1592 0.16Observed and theoretical agree to three decimals across all four configurations. Optimism grows linearly in the number of parameters and falls as \(1/n\), which is why training error is most misleading exactly where data are scarce and models are rich.
That is a surface over \((n,p)\):
n_grid <- round(seq(20, 400, length.out = 40))
p_grid <- 1:40
sigma2 <- 1
Zop <- outer(p_grid, n_grid, function(p, n) 2 * (p + 1) * sigma2 / n)
plot_ly(x = n_grid, y = p_grid, z = Zop, type = "surface",
colorscale = "Inferno", reversescale = TRUE,
colorbar = list(title = "E[optimism]")) |>
layout(title = "Expected optimism of the training error, 2(p+1)sigma^2 / n",
scene = list(xaxis = list(title = "Sample size n"),
yaxis = list(title = "Parameters p"),
zaxis = list(title = "Optimism")))Rotate toward the small-\(n\), large-\(p\) corner. The surface climbs steeply, and that corner is where most biomedical studies sit.
5 Information criteria as optimism corrections
If the optimism can be estimated, it can be added back. That single idea generates the whole family of information criteria.
\[ \begin{aligned} \textbf{Mallows' } C_p:\quad & C_p=\overline{\mathrm{err}}+\frac{2d\hat\sigma^2}{n}\\[1mm] \textbf{AIC:}\quad & \mathrm{AIC}=-\frac{2}{n}\,\mathrm{loglik}+\frac{2d}{n}\\[1mm] \textbf{BIC:}\quad & \mathrm{BIC}=-2\,\mathrm{loglik}+d\log n \end{aligned} \]
\(C_p\) is the optimism identity written as an estimator. AIC is its log-likelihood analogue. BIC replaces the penalty \(2d\) with \(d\log n\), so for \(n>e^2\approx7.4\) it penalizes complexity more heavily, BIC is consistent (it selects the true model as \(n\to\infty\) when the true model is in the candidate set), while AIC is efficient (it minimizes prediction error asymptotically when it is not).
Effective degrees of freedom generalizes \(d\) beyond parameter counting:
\[\boxed{\;\mathrm{df}(\hat f)=\frac{1}{\sigma^2}\sum_{i=1}^{n}\operatorname{Cov}(\hat y_i,y_i)\;}\]
For any linear smoother \(\hat{\mathbf y}=S\mathbf y\), this is \(\mathrm{df}=\operatorname{tr}(S)\), the trace of the hat matrix. For OLS, \(\operatorname{tr}(H)=p+1\), recovering the parameter count. For ridge regression it is \(\sum_j\frac{d_j^2}{d_j^2+\lambda}\), a continuous function of \(\lambda\) that falls smoothly from \(p\) to 0, so complexity is not a count but a dial.
set.seed(31)
n_df <- 100; p_df <- 20
Xdf <- scale(matrix(rnorm(n_df * p_df), n_df, p_df))
sv <- svd(Xdf)
lambdas <- 10^seq(-3, 4, length.out = 60)
edf <- vapply(lambdas, \(l) sum(sv$d^2 / (sv$d^2 + l)), numeric(1))
# Verify against the covariance definition by simulation at one lambda
lam0 <- 10
beta0 <- c(rep(1, 5), rep(0, p_df - 5))
S0 <- Xdf %*% solve(crossprod(Xdf) + lam0 * diag(p_df)) %*% t(Xdf)
sim_cov <- rowMeans(replicate(3000, {
y <- as.vector(Xdf %*% beta0) + rnorm(n_df)
(S0 %*% y) * y
})) - rowMeans(replicate(3000, {
y <- as.vector(Xdf %*% beta0) + rnorm(n_df); S0 %*% y
})) * as.vector(Xdf %*% beta0)
c(trace_of_smoother = round(sum(diag(S0)), 4),
formula_sum_d2_over_d2_plus_lambda = round(sum(sv$d^2 / (sv$d^2 + lam0)), 4),
simulated_covariance_sum = round(sum(sim_cov), 3))#> trace_of_smoother formula_sum_d2_over_d2_plus_lambda
#> 17.7906 17.7906
#> simulated_covariance_sum
#> 19.1780ggplot(data.frame(lambda = lambdas, edf = edf), aes(lambda, edf)) +
geom_line(linewidth = 1, color = "steelblue") +
geom_hline(yintercept = p_df, linetype = "dashed", color = "grey45") +
annotate("text", x = 1e-3, y = p_df - 0.9, hjust = 0, size = 3.2,
color = "grey35", label = "OLS: df = p") +
scale_x_log10() +
labs(title = "Effective degrees of freedom for ridge regression",
subtitle = expression(df(lambda) == sum(d[j]^2/(d[j]^2 + lambda))~", a continuous dial from p down to 0"),
x = expression(lambda~"(log scale)"), y = "Effective df") +
theme_dspa()
# --- Interactive equivalent ------------------------------------------------
plot_ly(x = ~lambdas, y = ~edf, type = "scatter", mode = "lines",
name = "Effective df") |>
add_lines(x = range(lambdas), y = c(p_df, p_df), name = "OLS df = p",
line = list(dash = "dash", color = "gray")) |>
layout(title = "Effective degrees of freedom vs. ridge penalty",
xaxis = list(title = "lambda", type = "log"),
yaxis = list(title = "Effective df"),
legend = list(orientation = "h"))The trace, the closed-form sum, and the simulated covariance all agree. This is why “number of parameters” is the wrong notion of complexity for regularized, ensemble, or early-stopped models, and why resampling, which requires no such count, is the general tool.
6 PART II: METRICS
7 The confusion matrix, and what each rate conditions on
Fix a positive class. For binary predictions:
| Reference: Positive | Reference: Negative | ||
|---|---|---|---|
| Predicted: Positive | \(TP\) | \(FP\) | \(\mathrm{PPV}=\frac{TP}{TP+FP}\) |
| Predicted: Negative | \(FN\) | \(TN\) | \(\mathrm{NPV}=\frac{TN}{TN+FN}\) |
| \(\mathrm{Sens}=\frac{TP}{TP+FN}\) | \(\mathrm{Spec}=\frac{TN}{TN+FP}\) |
caret::confusionMatrix(data, reference)takes predictions first. The printed matrix places predictions in rows and truth in columns. Reversing the arguments transposes the matrix, which exchanges sensitivity with specificity and PPV with NPV, silently, because the output still looks entirely plausible. Always passpositive =explicitly too; otherwisecaretuses the first factor level, which may not be the class you mean.
set.seed(41)
truth <- factor(rep(c("pos", "neg"), times = c(30, 70)), levels = c("neg", "pos"))
pred <- factor(c(rep("pos", 24), rep("neg", 6),
rep("pos", 14), rep("neg", 56)), levels = c("neg", "pos"))
right <- caret::confusionMatrix(data = pred, reference = truth, positive = "pos")
wrong <- caret::confusionMatrix(data = truth, reference = pred, positive = "pos")
rbind(`correct order (pred, truth)` = right$byClass[c("Sensitivity", "Specificity",
"Pos Pred Value", "Neg Pred Value")],
`reversed order (truth, pred)` = wrong$byClass[c("Sensitivity", "Specificity",
"Pos Pred Value", "Neg Pred Value")]) |>
round(4)#> Sensitivity Specificity Pos Pred Value
#> correct order (pred, truth) 0.8000 0.8000 0.6316
#> reversed order (truth, pred) 0.6316 0.9032 0.8000
#> Neg Pred Value
#> correct order (pred, truth) 0.9032
#> reversed order (truth, pred) 0.8000Sensitivity and PPV have swapped places. Neither output signals an error.
7.1 Row conditioning versus column conditioning
The single most useful organizing idea for this whole family: every rate is a conditional probability, and what distinguishes them is what they condition on.
\[ \begin{aligned} \textbf{Column rates (condition on the TRUTH):}\quad &\alpha=\text{Type I}=P(\hat Y{=}1\mid Y{=}0)=\tfrac{FP}{FP+TN}=1-\mathrm{Spec}\\ &\beta=\text{Type II}=P(\hat Y{=}0\mid Y{=}1)=\tfrac{FN}{TP+FN}=1-\mathrm{Sens}\\[3mm] \textbf{Row rates (condition on the PREDICTION):}\quad &\mathrm{FDR}=P(Y{=}0\mid \hat Y{=}1)=\tfrac{FP}{TP+FP}=1-\mathrm{PPV}\\ &\mathrm{FOR}=P(Y{=}1\mid \hat Y{=}0)=\tfrac{FN}{TN+FN}=1-\mathrm{NPV} \end{aligned} \]
Common misconception: “\(1-\text{precision}\) is the Type I error rate.” It is the false discovery rate. Type I error conditions on the case being truly negative; FDR conditions on the case being predicted positive. Bayes’ rule connects them through the prevalence \(\pi=P(Y{=}1)\):
\[\mathrm{PPV}=\frac{\mathrm{Sens}\cdot\pi}{\mathrm{Sens}\cdot\pi+(1-\mathrm{Spec})(1-\pi)} .\]
They coincide only under very particular conditions, and they diverge dramatically when the positive class is rare, which is the situation in screening, rare-disease diagnosis, and fraud detection. Sensitivity and specificity are properties of the test; PPV and NPV are properties of the test in a population.
sens_fixed <- 0.95; spec_fixed <- 0.95
prev <- 10^seq(-4, -0.3, length.out = 200)
ppv <- sens_fixed * prev / (sens_fixed * prev + (1 - spec_fixed) * (1 - prev))
data.frame(prevalence = prev, PPV = ppv,
`Type I error` = 1 - spec_fixed, check.names = FALSE) |>
pivot_longer(-prevalence, names_to = "quantity", values_to = "value") |>
ggplot(aes(prevalence, value, color = quantity)) +
geom_line(linewidth = 1) +
scale_x_log10(labels = scales::percent) +
scale_y_continuous(labels = scales::percent) +
scale_color_manual(values = c(PPV = "#3B7DD8", `Type I error` = "#D8433B")) +
labs(title = "A test with 95% sensitivity and 95% specificity",
subtitle = "Type I error is fixed at 5% by construction; PPV collapses as the disease becomes rare",
x = "Prevalence (log scale)", y = NULL, color = NULL) +
theme_dspa()
ppv_at <- function(p, se = 0.95, sp = 0.95) se * p / (se * p + (1 - sp) * (1 - p))
data.frame(prevalence = c(0.5, 0.1, 0.01, 0.001),
sensitivity = 0.95, specificity = 0.95,
type_I_error = 0.05,
PPV = round(ppv_at(c(0.5, 0.1, 0.01, 0.001)), 4),
FDR = round(1 - ppv_at(c(0.5, 0.1, 0.01, 0.001)), 4))#> prevalence sensitivity specificity type_I_error PPV FDR
#> 1 0.500 0.95 0.95 0.05 0.9500 0.0500
#> 2 0.100 0.95 0.95 0.05 0.6786 0.3214
#> 3 0.010 0.95 0.95 0.05 0.1610 0.8390
#> 4 0.001 0.95 0.95 0.05 0.0187 0.9813At a prevalence of 1 in 1,000, a test with 95% specificity, Type I error fixed at exactly 5%, produces a false discovery rate above 98%. The test has not changed; only what we condition on has.
binary_counts <- function(pred, truth, positive) {
pred <- as.character(pred); truth <- as.character(truth)
c(TP = sum(pred == positive & truth == positive),
FP = sum(pred == positive & truth != positive),
FN = sum(pred != positive & truth == positive),
TN = sum(pred != positive & truth != positive))
}
binary_metrics <- function(pred, truth, positive) {
k <- binary_counts(pred, truth, positive); n <- sum(k)
sens <- k[["TP"]] / (k[["TP"]] + k[["FN"]])
spec <- k[["TN"]] / (k[["TN"]] + k[["FP"]])
ppv <- k[["TP"]] / (k[["TP"]] + k[["FP"]])
npv <- k[["TN"]] / (k[["TN"]] + k[["FN"]])
acc <- (k[["TP"]] + k[["TN"]]) / n
p_e <- ((k[["TP"]] + k[["FP"]]) * (k[["TP"]] + k[["FN"]]) +
(k[["TN"]] + k[["FN"]]) * (k[["TN"]] + k[["FP"]])) / n^2
mcc_den <- sqrt(prod(c(k[["TP"]] + k[["FP"]], k[["TP"]] + k[["FN"]],
k[["TN"]] + k[["FP"]], k[["TN"]] + k[["FN"]])))
c(k, n = n, prevalence = (k[["TP"]] + k[["FN"]]) / n,
accuracy = acc, NIR = max(mean(truth == positive), 1 - mean(truth == positive)),
sensitivity = sens, specificity = spec, PPV = ppv, NPV = npv,
type_I = 1 - spec, type_II = 1 - sens, FDR = 1 - ppv, FOR = 1 - npv,
balanced_accuracy = (sens + spec) / 2,
F1 = 2 * ppv * sens / (ppv + sens),
kappa = (acc - p_e) / (1 - p_e),
MCC = if (mcc_den > 0) (k[["TP"]] * k[["TN"]] - k[["FP"]] * k[["FN"]]) / mcc_den else NA_real_)
}
round(binary_metrics(pred, truth, "pos"), 4)#> TP FP FN TN
#> 24.0000 14.0000 6.0000 56.0000
#> n prevalence accuracy NIR
#> 100.0000 0.3000 0.8000 0.7000
#> sensitivity specificity PPV NPV
#> 0.8000 0.8000 0.6316 0.9032
#> type_I type_II FDR FOR
#> 0.2000 0.2000 0.3684 0.0968
#> balanced_accuracy F1 kappa MCC
#> 0.8000 0.7059 0.5575 0.5665Every rate is computed from a named count vector.
Extracting counts by position from table() is fragile: a
perfectly classified fold has only one logical level, and
table(x)[1] then silently returns the wrong number.
8 Cohen’s kappa, and its two paradoxes
Accuracy does not account for agreement expected by chance. Cohen’s \(\kappa\) does:
\[\boxed{\;\kappa=\frac{p_o-p_e}{1-p_e}\;},\qquad p_o=\text{observed agreement},\quad p_e=\sum_{k}\hat P(\hat Y{=}k)\,\hat P(Y{=}k).\]
\(\kappa=1\) is perfect agreement with the ground truth; \(\kappa=0\) is chance-level agreement. Note the range:
\(\kappa\) is not bounded below by zero. It reaches \(-1\) in the symmetric case and can be negative whenever agreement is worse than chance. Its attainable minimum depends on the marginals, so a \(\kappa\) of \(-0.1\) on unbalanced data may already be near the floor.
#> Reference
#> Predicted neg pos
#> neg 56 6
#> pos 14 24n_tot <- sum(cm)
p_o <- sum(diag(cm)) / n_tot
p_e <- sum(rowSums(cm) * colSums(cm)) / n_tot^2
c(observed_agreement = round(p_o, 4),
expected_agreement = round(p_e, 4),
kappa_manual = round((p_o - p_e) / (1 - p_e), 4),
kappa_caret = round(caret::confusionMatrix(cm)$overall[["Kappa"]], 4))#> observed_agreement expected_agreement kappa_manual kappa_caret
#> 0.8000 0.5480 0.5575 0.5575Note that \(p_o\) is the accuracy (0.80 here), while \(\kappa\) is a different number entirely. Reporting one as the other is a common slip.
8.1 Paradox 1: kappa depends on prevalence
Two classifiers with identical sensitivity and specificity can have very different \(\kappa\) if the base rate differs.
kappa_from <- function(sens, spec, prev, n = 10000) {
TP <- n * prev * sens; FN <- n * prev * (1 - sens)
TN <- n * (1 - prev) * spec; FP <- n * (1 - prev) * (1 - spec)
p_o <- (TP + TN) / n
p_e <- ((TP + FP) * (TP + FN) + (TN + FN) * (TN + FP)) / n^2
(p_o - p_e) / (1 - p_e)
}
pv <- seq(0.02, 0.98, length.out = 200)
kp <- data.frame(prevalence = pv,
`sens = spec = 0.90` = kappa_from(0.90, 0.90, pv),
`sens = spec = 0.75` = kappa_from(0.75, 0.75, pv),
check.names = FALSE)
kp |> pivot_longer(-prevalence, names_to = "test", values_to = "kappa") |>
ggplot(aes(prevalence, kappa, color = test)) +
geom_line(linewidth = 1) +
scale_color_manual(values = c("#3B7DD8", "#D8433B")) +
labs(title = "Kappa depends on prevalence even when the test does not",
subtitle = "Sensitivity and specificity are held fixed along each curve; only the base rate changes",
x = "Prevalence", y = expression(kappa), color = NULL) +
theme_dspa()
data.frame(prevalence = c(0.5, 0.2, 0.05, 0.01),
sensitivity = 0.90, specificity = 0.90,
kappa = round(kappa_from(0.90, 0.90, c(0.5, 0.2, 0.05, 0.01)), 4))#> prevalence sensitivity specificity kappa
#> 1 0.50 0.9 0.9 0.8000
#> 2 0.20 0.9 0.9 0.7191
#> 3 0.05 0.9 0.9 0.4318
#> 4 0.01 0.9 0.9 0.1367The same test, 90% sensitive, 90% specific at every point, has \(\kappa=0.80\) at balanced prevalence and \(\kappa\) near 0.15 at 1% prevalence. \(\kappa\) is not a property of the classifier alone.
8.2 Paradox 2: high agreement with low kappa
paradox <- function(tab) {
n <- sum(tab); p_o <- sum(diag(tab)) / n
p_e <- sum(rowSums(tab) * colSums(tab)) / n^2
c(accuracy = round(p_o, 4), expected = round(p_e, 4),
kappa = round((p_o - p_e) / (1 - p_e), 4))
}
A <- matrix(c(45, 15, 25, 15), 2, byrow = TRUE) # balanced marginals
B <- matrix(c(86, 4, 10, 0), 2, byrow = TRUE) # very skewed marginals
rbind(`balanced marginals` = paradox(A),
`skewed marginals` = paradox(B))#> accuracy expected kappa
#> balanced marginals 0.60 0.540 0.1304
#> skewed marginals 0.86 0.868 -0.0606The second table has higher accuracy (0.86 versus 0.60) and a lower \(\kappa\), because heavily skewed marginals make chance agreement itself very high, leaving little room above it. This is the Feinstein–Cicchetti paradox.
The standard interpretive bands are Landis & Koch (1977): poor \(<0.20\), fair \(0.20\)–\(0.40\), moderate \(0.40\)–\(0.60\), good \(0.60\)–\(0.80\), very good \(0.80\)–\(1.00\). Use them with the prevalence in hand, and prefer the Matthews correlation coefficient when marginals are skewed, MCC is the Pearson correlation between the binary prediction and truth vectors and is symmetric in the two classes.
Weighted kappa is for ordinal outcomes. It assigns partial credit for near-misses, which requires an ordering on the categories. With two nominal classes there are no near-misses, so weighting has nothing to express. The weight matrix must also encode disagreement, with zeros on the diagonal; a matrix with non-zero diagonal entries does not define a valid scheme, and the resulting “kappa” can leave \([-1,1]\) for that reason rather than any property of the statistic.
9 Threshold-free evaluation
Every plug-in classifier produces a score \(\hat p(x)\) and then thresholds it. Metrics computed at one threshold describe one operating point; curves describe all of them.
9.1 ROC and AUC
The ROC curve traces \((\text{FPR},\text{TPR})\) as the threshold sweeps. Its area has an exact probabilistic meaning:
\[\boxed{\;\mathrm{AUC}=P\big(\hat p(X^{+})>\hat p(X^{-})\big)\;}\]
the probability that a randomly chosen positive scores above a randomly chosen negative. It is threshold-free and prevalence-invariant: change the class balance and the ROC curve does not move.
9.2 Precision–recall
The PR curve traces \((\text{Recall},\text{Precision})\). Precision is a row rate, so it does depend on prevalence, and that is the point.
Common misconception: “AUC of 0.9 means the model is good.” ROC AUC is prevalence-invariant, which is a virtue when comparing tests and a liability when deploying one. Under heavy imbalance, the FPR denominator \(FP+TN\) is dominated by \(TN\), so even a large number of false positives moves the ROC curve very little. The precision–recall curve, whose denominator is \(TP+FP\), registers them immediately.
The chance levels differ too: ROC AUC has chance level \(0.5\) regardless of balance, while PR-AUC has chance level equal to the positive prevalence. A PR-AUC of 0.3 is excellent at 1% prevalence and terrible at 50%. Always report the prevalence alongside.
library(PRROC)
set.seed(51)
make_scores <- function(n, prev, sep = 1.2) {
y <- rbinom(n, 1, prev)
s <- rnorm(n, mean = ifelse(y == 1, sep, 0))
list(y = y, s = s)
}
curves <- bind_rows(lapply(c(0.5, 0.1, 0.02), function(pv) {
d <- make_scores(4000, pv)
r <- PRROC::roc.curve(scores.class0 = d$s[d$y == 1],
scores.class1 = d$s[d$y == 0], curve = TRUE)
p <- PRROC::pr.curve(scores.class0 = d$s[d$y == 1],
scores.class1 = d$s[d$y == 0], curve = TRUE)
bind_rows(
data.frame(x = r$curve[, 1], y = r$curve[, 2], type = "ROC",
panel = sprintf("prevalence %.0f%%\nROC-AUC %.3f | PR-AUC %.3f",
100 * pv, r$auc, p$auc.integral)),
data.frame(x = p$curve[, 1], y = p$curve[, 2], type = "Precision-Recall",
panel = sprintf("prevalence %.0f%%\nROC-AUC %.3f | PR-AUC %.3f",
100 * pv, r$auc, p$auc.integral)))
}))
ggplot(curves, aes(x, y, color = type)) +
geom_line(linewidth = 0.9) +
facet_grid(type ~ panel, switch = "y") +
scale_color_manual(values = c(ROC = "#3B7DD8",
`Precision-Recall` = "#D8433B"), guide = "none") +
coord_fixed() +
labs(title = "The same classifier at three prevalences",
subtitle = "ROC is unchanged; the PR curve collapses as positives become rare",
x = NULL, y = NULL) +
theme_dspa(9)
The top row barely moves while the bottom row collapses. The classifier is identical in all three columns, only the base rate differs.
# --- Interactive equivalent ------------------------------------------------
d <- make_scores(4000, 0.1)
r <- PRROC::roc.curve(scores.class0 = d$s[d$y == 1],
scores.class1 = d$s[d$y == 0], curve = TRUE)
p <- PRROC::pr.curve(scores.class0 = d$s[d$y == 1],
scores.class1 = d$s[d$y == 0], curve = TRUE)
plot_ly() |>
add_lines(x = r$curve[, 1], y = r$curve[, 2],
name = sprintf("ROC (AUC = %.3f)", r$auc)) |>
add_lines(x = c(0, 1), y = c(0, 1), name = "Chance",
line = list(dash = "dash", color = "black")) |>
layout(title = "ROC curve",
xaxis = list(title = "False positive rate", scaleanchor = "y"),
yaxis = list(title = "True positive rate"),
legend = list(orientation = "h"))
plot_ly() |>
add_lines(x = p$curve[, 1], y = p$curve[, 2],
name = sprintf("PR (AUC = %.3f)", p$auc.integral)) |>
add_lines(x = c(0, 1), y = rep(mean(d$y), 2), name = "Chance = prevalence",
line = list(dash = "dash", color = "black")) |>
layout(title = "Precision-Recall curve",
xaxis = list(title = "Recall"), yaxis = list(title = "Precision"),
legend = list(orientation = "h"))
ROCR::performance(pred, "auc")computes the ROC AUC. Displaying it next to a precision–recall plot reports the wrong quantity. UsePRROC::pr.curve()oryardstick::pr_auc()for the area under the PR curve.
set.seed(53)
d <- make_scores(3000, 0.05)
r <- PRROC::roc.curve(scores.class0 = d$s[d$y == 1], scores.class1 = d$s[d$y == 0])
p <- PRROC::pr.curve(scores.class0 = d$s[d$y == 1], scores.class1 = d$s[d$y == 0])
c(ROC_AUC = round(r$auc, 4), PR_AUC = round(p$auc.integral, 4),
prevalence = round(mean(d$y), 4),
PR_chance_level = round(mean(d$y), 4), ROC_chance_level = 0.5)#> ROC_AUC PR_AUC prevalence PR_chance_level
#> 0.8080 0.2154 0.0513 0.0513
#> ROC_chance_level
#> 0.500010 Calibration
AUC measures ranking. It says nothing about whether a predicted probability of 0.8 corresponds to an event that happens 80% of the time. Those are separate properties, and a model can have one without the other.
Common misconception: “a high AUC means the predicted probabilities are trustworthy.” AUC is invariant to any strictly increasing transformation of the scores. Replace every \(\hat p\) by \(\hat p^{10}\) and the AUC is unchanged, while every probability is now wildly wrong. If the output feeds a threshold chosen by cost, a risk communication, or a downstream calculation, ranking is not enough.
A model is calibrated if, among cases assigned probability \(q\), the event occurs a fraction \(q\) of the time:
\[\boxed{\;P\big(Y=1\ \big|\ \hat p(X)=q\big)=q\qquad\text{for all }q\in[0,1]\;}\]
10.1 Reliability diagrams and the calibration line
Bin the predictions, plot observed frequency against mean predicted probability, and compare to the diagonal. The parametric summary is the calibration intercept and slope from the logistic recalibration model
\[\operatorname{logit}\big(P(Y=1)\big)=\alpha+\gamma\,\operatorname{logit}(\hat p),\]
where perfect calibration is \(\alpha=0\), \(\gamma=1\). A slope \(\gamma<1\) means the predictions are too extreme (over-confident); \(\gamma>1\) means they are too conservative; \(\alpha\ne0\) means a systematic shift in the overall risk level.
set.seed(61)
n_cal <- 4000
x_cal <- rnorm(n_cal)
p_true <- plogis(1.2 * x_cal)
y_cal <- rbinom(n_cal, 1, p_true)
# Three score vectors with IDENTICAL ranking, hence identical AUC
p_good <- p_true
p_over <- plogis(3.0 * qlogis(p_true)) # over-confident: slope < 1
p_shift <- plogis(qlogis(p_true) + 1.2) # shifted: intercept != 0
auc_of <- function(p) as.numeric(pROC::auc(pROC::roc(y_cal, p, quiet = TRUE)))
c(AUC_well_calibrated = round(auc_of(p_good), 4),
AUC_over_confident = round(auc_of(p_over), 4),
AUC_shifted = round(auc_of(p_shift), 4))#> AUC_well_calibrated AUC_over_confident AUC_shifted
#> 0.7795 0.7795 0.7795All three AUCs are identical to four decimals, because all three are monotone transformations of the same score.
calib_summary <- function(p, y, bins = 12) {
eps <- 1e-12; p <- pmin(pmax(p, eps), 1 - eps)
fit <- glm(y ~ qlogis(p), family = binomial())
c(Brier = mean((p - y)^2),
log_loss = -mean(y * log(p) + (1 - y) * log(1 - p)),
calib_intercept = unname(coef(fit)[1]),
calib_slope = unname(coef(fit)[2]),
mean_predicted = mean(p), observed_rate = mean(y))
}
rbind(`well calibrated` = calib_summary(p_good, y_cal),
`over-confident` = calib_summary(p_over, y_cal),
`shifted` = calib_summary(p_shift, y_cal)) |> round(4)#> Brier log_loss calib_intercept calib_slope mean_predicted
#> well calibrated 0.1913 0.5631 -0.0400 1.0395 0.4974
#> over-confident 0.2187 0.7444 -0.0400 0.3465 0.4979
#> shifted 0.2497 0.7079 -1.2873 1.0395 0.7168
#> observed_rate
#> well calibrated 0.4898
#> over-confident 0.4898
#> shifted 0.4898The Brier score, log-loss, and calibration slope all separate the three models that AUC cannot tell apart.
reliability <- function(p, y, label, bins = 12) {
b <- cut(p, breaks = seq(0, 1, length.out = bins + 1), include.lowest = TRUE)
data.frame(p = p, y = y, b = b) |>
summarise(predicted = mean(p), observed = mean(y), n = dplyr::n(), .by = b) |>
filter(n >= 20) |> mutate(model = label)
}
bind_rows(reliability(p_good, y_cal, "Well calibrated"),
reliability(p_over, y_cal, "Over-confident"),
reliability(p_shift, y_cal, "Shifted")) |>
ggplot(aes(predicted, observed, color = model)) +
geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "grey50") +
geom_line(linewidth = 0.9) + geom_point(aes(size = n)) +
scale_color_manual(values = c(`Well calibrated` = "#3B7DD8",
`Over-confident` = "#D8433B",
Shifted = "#7FB069")) +
scale_size_continuous(range = c(1.2, 4.5), guide = "none") +
coord_fixed(xlim = c(0, 1), ylim = c(0, 1)) +
labs(title = "Reliability diagram: three models with identical AUC",
subtitle = "Distance from the diagonal is miscalibration, which ranking metrics cannot see",
x = "Mean predicted probability", y = "Observed frequency", color = NULL) +
theme_dspa()
# --- Interactive equivalent ------------------------------------------------
rg <- reliability(p_good, y_cal, "Well calibrated")
ro <- reliability(p_over, y_cal, "Over-confident")
plot_ly() |>
add_trace(x = rg$predicted, y = rg$observed, type = "scatter",
mode = "lines+markers", name = "Well calibrated") |>
add_trace(x = ro$predicted, y = ro$observed, type = "scatter",
mode = "lines+markers", name = "Over-confident") |>
add_lines(x = c(0, 1), y = c(0, 1), name = "Perfect calibration",
line = list(dash = "dash", color = "black")) |>
layout(title = "Reliability diagram",
xaxis = list(title = "Predicted probability", scaleanchor = "y"),
yaxis = list(title = "Observed frequency"),
legend = list(orientation = "h"))Miscalibration is a two-dimensional phenomenon, its magnitude depends on both the predicted probability and the strength of the distortion, so it is worth rotating:
p_grid <- seq(0.01, 0.99, length.out = 60)
gamma_grid <- seq(0.3, 3.0, length.out = 60)
# Observed frequency when predictions are distorted with slope gamma
Zcal <- outer(gamma_grid, p_grid, function(g, q) plogis(qlogis(q) / g))
plot_ly(x = p_grid, y = gamma_grid, z = Zcal, type = "surface",
colorscale = "RdBu", reversescale = TRUE,
colorbar = list(title = "Observed\nfrequency")) |>
add_trace(x = p_grid, y = rep(1, length(p_grid)), z = p_grid,
type = "scatter3d", mode = "lines", name = "Perfect calibration",
line = list(width = 8, color = "black")) |>
layout(title = "Observed frequency as a function of predicted probability and calibration slope",
scene = list(xaxis = list(title = "Predicted probability"),
yaxis = list(title = "Calibration slope gamma"),
zaxis = list(title = "Observed frequency")))The black ridge at \(\gamma=1\) is the diagonal. Away from it the surface bends: \(\gamma<1\) pulls observed frequencies toward the middle (predictions too extreme), \(\gamma>1\) pushes them outward.
10.2 Proper scoring rules
A scoring rule is proper if it is minimized in expectation by reporting the true probability, and strictly proper if that minimizer is unique. Accuracy is not proper, it can be improved by distorting probabilities toward 0 and 1.
\[ \begin{aligned} \textbf{Brier score: }\quad &\mathrm{BS}=\frac1n\sum_{i=1}^{n}\big(\hat p_i-y_i\big)^2 &&\text{(quadratic, strictly proper)}\\[1mm] \textbf{Log-loss: }\quad &\mathrm{LL}=-\frac1n\sum_{i=1}^{n}\Big[y_i\log\hat p_i+(1-y_i)\log(1-\hat p_i)\Big] &&\text{(logarithmic, strictly proper)} \end{aligned} \]
Log-loss is unbounded, a confident wrong prediction (\(\hat p\to0\) when \(y=1\)) costs infinitely, so it punishes over-confidence far more harshly than Brier.
Murphy’s decomposition splits the Brier score into three interpretable pieces:
\[\boxed{\;\mathrm{BS}=\underbrace{\frac1n\sum_k n_k(\bar p_k-\bar y_k)^2}_{\text{reliability (calibration), lower is better}} -\underbrace{\frac1n\sum_k n_k(\bar y_k-\bar y)^2}_{\text{resolution (discrimination), higher is better}} +\underbrace{\bar y(1-\bar y)}_{\text{uncertainty (irreducible)}}\;}\]
The third term depends only on the base rate, it is the Brier score of the best constant predictor and cannot be improved by any model.
brier_decomp <- function(p, y, bins = 15) {
b <- cut(p, breaks = seq(0, 1, length.out = bins + 1), include.lowest = TRUE)
d <- data.frame(p = p, y = y, b = b) |>
summarise(n = dplyr::n(), pbar = mean(p), ybar = mean(y), .by = b)
N <- length(y); ybar_all <- mean(y)
rel <- sum(d$n * (d$pbar - d$ybar)^2) / N
res <- sum(d$n * (d$ybar - ybar_all)^2) / N
unc <- ybar_all * (1 - ybar_all)
c(reliability = rel, resolution = res, uncertainty = unc,
reconstructed = rel - res + unc, actual_brier = mean((p - y)^2))
}
rbind(`well calibrated` = brier_decomp(p_good, y_cal),
`over-confident` = brier_decomp(p_over, y_cal),
`shifted` = brier_decomp(p_shift, y_cal)) |> round(5)#> reliability resolution uncertainty reconstructed actual_brier
#> well calibrated 0.00086 0.05892 0.24989 0.19184 0.19130
#> over-confident 0.02632 0.05643 0.24989 0.21979 0.21866
#> shifted 0.05911 0.05892 0.24989 0.25008 0.24966The decomposition reconstructs the Brier score exactly. Read the columns: all three models have similar resolution (they rank equally well, as the AUCs showed) and very different reliability. That is the distinction the whole section is about, now quantified.
10.3 Recalibration
Miscalibration is fixable without touching the model, by learning a monotone map from scores to probabilities. Two standard methods:
Platt scaling fits a one-parameter logistic model \(P(Y=1\mid s)=\sigma(a+bs)\) on held-out scores. Parametric, needs little data, assumes a sigmoidal distortion.
Isotonic regression fits the best non-decreasing step function by pool-adjacent-violators. Non-parametric and more flexible, but it needs more data and can overfit in small samples.
The recalibration map must be fitted on data the model did not train on. Fitting it on the training predictions learns the model’s in-sample over-confidence, which is not the distortion present out of sample.
set.seed(63)
idx_cal <- sample(n_cal, n_cal / 2) # calibration set
idx_ev <- setdiff(seq_len(n_cal), idx_cal) # evaluation set
# Platt: logistic on the LOGIT of the raw score, fitted on the calibration half
platt <- glm(y_cal[idx_cal] ~ qlogis(p_over[idx_cal]), family = binomial())
p_platt <- plogis(predict(platt, newdata = data.frame(
`qlogis(p_over[idx_cal])` = qlogis(p_over[idx_ev]), check.names = FALSE)))
p_platt <- plogis(coef(platt)[1] + coef(platt)[2] * qlogis(p_over[idx_ev]))
# Isotonic: monotone step function, also fitted on the calibration half
iso <- isoreg(p_over[idx_cal], y_cal[idx_cal])
iso_fn <- approxfun(iso$x, iso$yf, method = "linear", rule = 2)
p_iso <- pmin(pmax(iso_fn(p_over[idx_ev]), 1e-6), 1 - 1e-6)
rbind(`uncalibrated` = calib_summary(p_over[idx_ev], y_cal[idx_ev]),
`Platt scaling` = calib_summary(p_platt, y_cal[idx_ev]),
`isotonic` = calib_summary(p_iso, y_cal[idx_ev])) |>
round(4)#> Brier log_loss calib_intercept calib_slope mean_predicted
#> uncalibrated 0.2180 0.7365 -0.0432 0.3511 0.4929
#> Platt scaling 0.1901 0.5602 -0.0057 1.0267 0.4868
#> isotonic 0.2988 0.8220 -0.0586 -0.0379 0.4942
#> observed_rate
#> uncalibrated 0.4855
#> Platt scaling 0.4855
#> isotonic 0.4855bind_rows(reliability(p_over[idx_ev], y_cal[idx_ev], "Uncalibrated"),
reliability(p_platt, y_cal[idx_ev], "Platt scaling"),
reliability(p_iso, y_cal[idx_ev], "Isotonic")) |>
ggplot(aes(predicted, observed, color = model)) +
geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "grey50") +
geom_line(linewidth = 0.9) + geom_point(size = 2) +
scale_color_manual(values = c(Uncalibrated = "#D8433B",
`Platt scaling` = "#3B7DD8",
Isotonic = "#7FB069")) +
coord_fixed(xlim = c(0, 1), ylim = c(0, 1)) +
labs(title = "Recalibration on held-out data",
subtitle = "The ranking, and hence the AUC, is unchanged -- only the probabilities move",
x = "Mean predicted probability", y = "Observed frequency", color = NULL) +
theme_dspa()