| 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()
c(AUC_before = round(auc_of2 <- as.numeric(pROC::auc(pROC::roc(
y_cal[idx_ev], p_over[idx_ev], quiet = TRUE))), 4),
AUC_after_platt = round(as.numeric(pROC::auc(pROC::roc(
y_cal[idx_ev], p_platt, quiet = TRUE))), 4))#> AUC_before AUC_after_platt
#> 0.7814 0.7814The AUC is unchanged, both recalibration maps are monotone, so the ranking is preserved exactly, while the Brier score and calibration slope improve substantially. Recalibration is free improvement in everything except ranking.
11 Decision-curve analysis
Metrics so far answer “how good is this model?” A different and often more useful question is “is using this model better than not using it?”
Suppose acting on a positive prediction has benefit \(B\) per true positive and cost \(C\) per false positive. A rational decision-maker acts when the predicted probability exceeds the threshold probability \(p_t=C/(B+C)\), so \(p_t\) encodes the cost ratio. The net benefit at threshold \(p_t\) is
\[\boxed{\;\mathrm{NB}(p_t)=\frac{TP}{n}-\frac{FP}{n}\cdot\frac{p_t}{1-p_t}\;}\]
measured in units of true positives per patient, with false positives discounted by the odds of the threshold (Vickers & Elkin, 2006).
Two reference strategies bracket any model:
\[\mathrm{NB}_{\text{treat all}}(p_t)=\pi-(1-\pi)\frac{p_t}{1-p_t}, \qquad \mathrm{NB}_{\text{treat none}}=0 .\]
A model is worth using only where its curve lies above both.
net_benefit <- function(p, y, thresholds = seq(0.01, 0.60, by = 0.005)) {
n <- length(y)
vapply(thresholds, function(t) {
pos <- p >= t
TP <- sum(pos & y == 1); FP <- sum(pos & y == 0)
TP / n - (FP / n) * (t / (1 - t))
}, numeric(1))
}
set.seed(65)
n_dc <- 3000
x_dc <- rnorm(n_dc)
p_dc_true <- plogis(-2.2 + 1.4 * x_dc) # ~12% prevalence
y_dc <- rbinom(n_dc, 1, p_dc_true)
p_model <- p_dc_true
p_weak <- plogis(-2.2 + 0.4 * x_dc) # poorly discriminating
th <- seq(0.01, 0.60, by = 0.005)
pi_dc <- mean(y_dc)
dc <- bind_rows(
data.frame(t = th, nb = net_benefit(p_model, y_dc, th), strategy = "Good model"),
data.frame(t = th, nb = net_benefit(p_weak, y_dc, th), strategy = "Weak model"),
data.frame(t = th, nb = pi_dc - (1 - pi_dc) * th / (1 - th), strategy = "Treat all"),
data.frame(t = th, nb = 0, strategy = "Treat none"))
ggplot(dc, aes(t, nb, color = strategy, linetype = strategy)) +
geom_line(linewidth = 0.9) +
coord_cartesian(ylim = c(-0.05, max(dc$nb) * 1.05)) +
scale_color_manual(values = c(`Good model` = "#3B7DD8", `Weak model` = "#D8433B",
`Treat all` = "grey40", `Treat none` = "black")) +
scale_linetype_manual(values = c(`Good model` = "solid", `Weak model` = "solid",
`Treat all` = "dashed", `Treat none` = "dotted")) +
labs(title = "Decision curve: is the model worth using?",
subtitle = sprintf("Prevalence %.1f%%. A model earns its place only where it rises above BOTH references",
100 * pi_dc),
x = expression("Threshold probability "*p[t]), y = "Net benefit",
color = NULL, linetype = NULL) +
theme_dspa()
# --- Interactive equivalent ------------------------------------------------
plot_ly() |>
add_lines(x = th, y = net_benefit(p_model, y_dc, th), name = "Good model") |>
add_lines(x = th, y = net_benefit(p_weak, y_dc, th), name = "Weak model") |>
add_lines(x = th, y = pi_dc - (1 - pi_dc) * th / (1 - th), name = "Treat all",
line = list(dash = "dash", color = "gray")) |>
add_lines(x = th, y = rep(0, length(th)), name = "Treat none",
line = list(dash = "dot", color = "black")) |>
layout(title = "Decision curve analysis",
xaxis = list(title = "Threshold probability"),
yaxis = list(title = "Net benefit", range = c(-0.05, 0.15)),
legend = list(orientation = "h"))Net benefit depends on both the threshold and the prevalence, which makes it a surface:
prev_grid <- seq(0.02, 0.45, length.out = 40)
th_grid <- seq(0.02, 0.55, length.out = 40)
# Net benefit of a fixed-quality model (sens 0.85, spec 0.80 at each threshold's
# operating point, approximated by a bi-normal ROC) minus treat-all
Znb <- outer(prev_grid, th_grid, function(pv, t) {
sens <- 0.85; spec <- 0.80
nb_model <- pv * sens - (1 - pv) * (1 - spec) * (t / (1 - t))
nb_all <- pv - (1 - pv) * (t / (1 - t))
pmax(nb_model - pmax(nb_all, 0), -0.05)
})
plot_ly(x = th_grid, y = prev_grid, z = Znb, type = "surface",
colorscale = "RdBu",
colorbar = list(title = "NB(model)\n- NB(best reference)")) |>
layout(title = "Where does a model add net benefit over treat-all / treat-none?",
scene = list(xaxis = list(title = "Threshold probability"),
yaxis = list(title = "Prevalence"),
zaxis = list(title = "Excess net benefit")))The region where the surface is positive is where the model helps. Rotate to see that it shrinks toward the low-prevalence, low-threshold corner: when the disease is rare and the cost of missing it is high, treating everyone is hard to beat.
12 Regression metrics
\[ \begin{aligned} \mathrm{RMSE}&=\sqrt{\tfrac1n\textstyle\sum_i(y_i-\hat y_i)^2} & \mathrm{MAE}&=\tfrac1n\textstyle\sum_i|y_i-\hat y_i|\\[1mm] R^2&=1-\frac{\sum_i(y_i-\hat y_i)^2}{\sum_i(y_i-\bar y)^2} & \mathrm{CCC}&=\frac{2\operatorname{Cov}(\hat y,y)}{\operatorname{Var}(\hat y)+\operatorname{Var}(y)+(\bar{\hat y}-\bar y)^2} \end{aligned} \]
Common misconception: “report the correlation between predicted and observed.” Correlation is invariant to affine transformation: \(\operatorname{cor}(a+b\hat y,\,y)=\operatorname{cor}(\hat y,y)\) for any \(b>0\). A model returning \(\hat y=100+0.01y\) correlates exactly 1 with the truth and is useless. Correlation cannot detect bias, cannot detect scale error, and is not a loss anyone minimizes.
Report RMSE and MAE in the response’s units, \(R^2\) computed on held-out data, which can be negative, informatively so, when the model is worse than the mean, and the calibration slope and intercept from regressing \(y\) on \(\hat y\), where agreement means intercept 0 and slope 1. Lin’s concordance correlation coefficient combines the last two into one number.
reg_metrics <- function(pred, obs) {
cal <- coef(lm(obs ~ pred))
c(RMSE = sqrt(mean((pred - obs)^2)),
MAE = mean(abs(pred - obs)),
R2 = 1 - sum((pred - obs)^2) / sum((obs - mean(obs))^2),
correlation = cor(pred, obs),
calib_intercept = unname(cal[1]), calib_slope = unname(cal[2]),
CCC = 2 * cov(pred, obs) /
(var(pred) + var(obs) + (mean(pred) - mean(obs))^2))
}
set.seed(67)
y_reg <- rnorm(500, 50, 10)
rbind(`good model` = reg_metrics(y_reg + rnorm(500, 0, 3), y_reg),
`biased by +20` = reg_metrics(y_reg + 20, y_reg),
`wrong scale` = reg_metrics(100 + 0.01 * y_reg, y_reg),
`mean predictor` = reg_metrics(rep(mean(y_reg), 500), y_reg)) |>
round(4)#> RMSE MAE R2 correlation calib_intercept calib_slope
#> good model 2.9527 2.3633 0.9106 0.959 4.5349 0.9094
#> biased by +20 20.0000 20.0000 -3.1009 1.000 -20.0000 1.0000
#> wrong scale 51.3557 50.4164 -26.0394 1.000 -10000.0000 100.0000
#> mean predictor 9.8762 7.8065 0.0000 NA 50.0845 NA
#> CCC
#> good model 0.9577
#> biased by +20 0.3283
#> wrong scale 0.0007
#> mean predictor 0.0000Row 3 is the counterexample: correlation 1.0, RMSE 50, \(R^2 = -24\). Row 4 shows what a negative \(R^2\) is measured against, the constant mean predictor scores exactly 0, so anything below that is worse than predicting the average.
13 PART III: RESAMPLING
Optimism (§9.2) has a closed form only for linear models with known noise. Resampling estimates it for anything.
14 The holdout, and the variance of one split
Split once, fit on one part, score on the other. Simple, unbiased for the smaller training size, and highly variable.
set.seed(71)
n_hv <- 200; p_hv <- 6
beta_hv <- c(1.5, -1, 0.8, rep(0, p_hv - 3))
X_hv <- matrix(rnorm(n_hv * p_hv), n_hv, p_hv)
y_hv <- as.vector(X_hv %*% beta_hv) + rnorm(n_hv, sd = 1.5)
d_hv <- data.frame(X_hv, y = y_hv)
holdout_rmse <- function(prop = 0.7) {
tr <- sample(n_hv, floor(prop * n_hv))
fit <- lm(y ~ ., data = d_hv[tr, ])
sqrt(mean((predict(fit, d_hv[-tr, ]) - d_hv$y[-tr])^2))
}
set.seed(73)
hv <- replicate(2000, holdout_rmse())
c(mean = round(mean(hv), 4), sd = round(sd(hv), 4),
q025 = round(quantile(hv, 0.025), 4), q975 = round(quantile(hv, 0.975), 4),
relative_spread = round(sd(hv) / mean(hv), 4))#> mean sd q025.2.5% q975.97.5% relative_spread
#> 1.6139 0.1511 1.3158 1.9063 0.0937ggplot(data.frame(rmse = hv), aes(rmse)) +
geom_histogram(bins = 45, fill = "steelblue", color = "white") +
geom_vline(xintercept = mean(hv), color = "firebrick", linewidth = 1) +
labs(title = "2,000 different 70/30 splits of the same 200 observations",
subtitle = sprintf("Same data, same model. The 95%% range spans %.2f to %.2f",
quantile(hv, 0.025), quantile(hv, 0.975)),
x = "Held-out RMSE", y = "Splits") +
theme_dspa()
Common misconception: “I held out a test set, so my estimate is reliable.” It is unbiased but not precise. A single 70/30 split of 200 observations here gives a held-out RMSE anywhere in a range spanning roughly 30% of its own value. Reporting one number from one split conveys none of that uncertainty. Resampling averages over many splits, which is what makes the estimate usable.
15 k-fold cross-validation
Partition into \(k\) folds; for each, train on \(k-1\) and assess on the held-out one; average.
\[\widehat{\mathrm{CV}}_{(k)}=\frac{1}{n}\sum_{i=1}^{n}L\Big(y_i,\ \hat f^{-\kappa(i)}(x_i)\Big),\]
where \(\kappa(i)\) is the fold containing \(i\) and \(\hat f^{-\kappa}\) is the model fitted with fold \(\kappa\) removed. Every observation is assessed exactly once, by a model that never saw it.
15.1 The bias–variance tradeoff in \(k\)
Bias. Each model trains on \(n(1-1/k)\) observations rather than \(n\). Since learning curves decrease in \(n\), the CV estimate is pessimistically biased, and the bias grows as \(k\) shrinks. At \(k=2\) each model sees half the data.
Variance. The \(k\) fold estimates are positively correlated, the training sets share \(n(1-2/k)\) observations. Larger \(k\) means more overlap, more correlation, and (by the ensemble variance identity of Chapter 5, §5.22.1) a variance floor that averaging cannot reduce. LOOCV has the most overlap of all.
set.seed(81)
cv_rmse <- function(dat, k, seed) {
set.seed(seed)
fold <- sample(rep(seq_len(k), length.out = nrow(dat)))
err <- vapply(seq_len(k), function(f) {
fit <- lm(y ~ ., data = dat[fold != f, ])
mean((predict(fit, dat[fold == f, ]) - dat$y[fold == f])^2)
}, numeric(1))
sqrt(mean(err))
}
# Truth: RMSE of a model fitted on ALL n, evaluated on a very large fresh sample
truth_rmse <- local({
fit <- lm(y ~ ., data = d_hv)
Xb <- matrix(rnorm(50000 * p_hv), 50000, p_hv)
yb <- as.vector(Xb %*% beta_hv) + rnorm(50000, sd = 1.5)
sqrt(mean((yb - cbind(1, Xb) %*% coef(fit))^2))
})
ks <- c(2, 3, 5, 10, 20, 50, nrow(d_hv))
bv_tab <- do.call(rbind, lapply(ks, function(k) {
est <- vapply(1:150, \(r) cv_rmse(d_hv, k, seed = 1000 + r), numeric(1))
data.frame(k = k, mean = mean(est), sd = sd(est),
bias = mean(est) - truth_rmse)
}))
bv_tab$k_label <- ifelse(bv_tab$k == nrow(d_hv), "LOOCV", as.character(bv_tab$k))
bv_tab |> mutate(across(c(mean, sd, bias), \(z) round(z, 4)))#> k mean sd bias k_label
#> 1 2 1.6343 0.0450 0.1218 2
#> 2 3 1.6168 0.0283 0.1043 3
#> 3 5 1.6116 0.0186 0.0991 5
#> 4 10 1.6079 0.0123 0.0953 10
#> 5 20 1.6064 0.0076 0.0939 20
#> 6 50 1.6049 0.0040 0.0924 50
#> 7 200 1.6048 0.0000 0.0922 LOOCV#> truth
#> 1.5125bv_tab |>
select(k, Bias = bias, `Standard deviation` = sd) |>
pivot_longer(-k, names_to = "component", values_to = "value") |>
ggplot(aes(k, value, color = component)) +
geom_hline(yintercept = 0, color = "grey65") +
geom_line(linewidth = 1) + geom_point(size = 2.2) +
scale_x_log10(breaks = ks, labels = bv_tab$k_label) +
scale_color_manual(values = c(Bias = "#D8433B",
`Standard deviation` = "#3B7DD8")) +
labs(title = "Bias and variance of the cross-validation estimator itself",
subtitle = "Small k trains on less data (pessimistic bias); large k gives overlapping training sets (variance floor)",
x = "Number of folds k (log scale)", y = NULL, color = NULL) +
theme_dspa()
Bias falls monotonically toward zero as \(k\) grows; the standard deviation does not fall correspondingly, because the fold estimates become increasingly correlated. \(k=5\) or \(k=10\) is the standard compromise for exactly this reason.
The tradeoff is a surface over \((k,n)\):
n_grid2 <- c(40, 60, 100, 160, 260, 400)
k_grid2 <- c(2, 3, 5, 10, 20)
cv_bias_at <- function(n, k, reps = 40) {
set.seed(n * 100 + k)
vapply(seq_len(reps), function(r) {
X <- matrix(rnorm(n * p_hv), n, p_hv)
y <- as.vector(X %*% beta_hv) + rnorm(n, sd = 1.5)
dd <- data.frame(X, y = y)
fold <- sample(rep(seq_len(k), length.out = n))
e <- vapply(seq_len(k), function(f) {
fit <- lm(y ~ ., data = dd[fold != f, ])
mean((predict(fit, dd[fold == f, ]) - dd$y[fold == f])^2)
}, numeric(1))
sqrt(mean(e))
}, numeric(1)) |> mean()
}
# Zcv <- outer(k_grid2, n_grid2, Vectorize(cv_bias_at)) - truth_rmse
Zcv <-
outer(n_grid2, k_grid2, Vectorize(function(n, k) cv_bias_at(n = n, k = k))) -
truth_rmse
plot_ly(x = n_grid2, y = k_grid2, z = Zcv, type = "surface",
colorscale = "Viridis", colorbar = list(title = "CV bias")) |>
layout(title = "Bias of the k-fold CV estimate, as a function of k and n",
scene = list(xaxis = list(title = "Sample size n"),
yaxis = list(title = "Folds k"),
zaxis = list(title = "Estimated - true RMSE")))The bias is largest in the small-\(n\), small-\(k\) corner and flattens as either grows. With \(n\) in the thousands the choice of \(k\) scarcely matters; with \(n\) in the dozens it matters a great deal.
Stratification preserves the class proportions in every fold. For classification, especially under imbalance, it is not optional: an unstratified fold can contain almost no minority cases, making its metric undefined or wildly variable.
set.seed(83)
y_imb <- factor(rep(c("pos", "neg"), times = c(25, 475)))
strat_check <- function(strata = TRUE) {
f <- if (strata) rsample::vfold_cv(data.frame(y = y_imb), v = 10, strata = y)
else rsample::vfold_cv(data.frame(y = y_imb), v = 10)
vapply(f$splits, \(s) sum(rsample::assessment(s)$y == "pos"), numeric(1))
}
rbind(unstratified = range(strat_check(FALSE)),
stratified = range(strat_check(TRUE))) |>
`colnames<-`(c("min positives per fold", "max positives per fold"))#> min positives per fold max positives per fold
#> unstratified 0 4
#> stratified 1 416 Leave-one-out, and a shortcut
LOOCV sets \(k=n\): each observation is held out once. It is nearly unbiased for \(\mathrm{Err}\) but costs \(n\) model fits, except that for linear smoothers it costs one.
The LOOCV shortcut. If \(\hat{\mathbf y}=S\mathbf y\) for a smoother matrix \(S\) not depending on \(\mathbf y\), then \[\boxed{\;\mathrm{CV}_{(n)}=\frac{1}{n}\sum_{i=1}^{n}\left(\frac{y_i-\hat y_i}{1-h_{ii}}\right)^{2}\;}\] where \(h_{ii}=S_{ii}\) is the \(i\)-th leverage. One fit replaces \(n\).
Why. Removing observation \(i\) changes the fit by an amount that can be written in closed form via the Sherman–Morrison update; the leave-one-out residual turns out to be the ordinary residual inflated by \(1/(1-h_{ii})\).
Generalized cross-validation replaces each \(h_{ii}\) by their average \(\operatorname{tr}(S)/n\):
\[\mathrm{GCV}=\frac{1}{n}\sum_{i=1}^{n}\left(\frac{y_i-\hat y_i}{1-\operatorname{tr}(S)/n}\right)^{2},\]
which is rotation-invariant, cheaper still, and the standard criterion for smoothing-parameter selection.
set.seed(85)
n_lo <- 120
X_lo <- matrix(rnorm(n_lo * 5), n_lo, 5)
y_lo <- as.vector(X_lo %*% c(1, -1, 0.5, 0, 0)) + rnorm(n_lo)
d_lo <- data.frame(X_lo, y = y_lo)
# Brute force: n model fits
t_brute <- system.time({
loo_brute <- mean(vapply(seq_len(n_lo), function(i) {
fit <- lm(y ~ ., data = d_lo[-i, ])
(d_lo$y[i] - predict(fit, d_lo[i, ]))^2
}, numeric(1)))
})[["elapsed"]]
# Shortcut: ONE model fit
t_fast <- system.time({
fit_all <- lm(y ~ ., data = d_lo)
h <- hatvalues(fit_all)
loo_fast <- mean((residuals(fit_all) / (1 - h))^2)
gcv <- mean((residuals(fit_all) / (1 - sum(h) / n_lo))^2)
})[["elapsed"]]
c(loocv_brute_force = round(loo_brute, 6),
loocv_closed_form = round(loo_fast, 6),
difference = signif(abs(loo_brute - loo_fast), 3),
GCV = round(gcv, 6),
brute_seconds = round(t_brute, 3), shortcut_seconds = round(t_fast, 4),
fits_required = paste(n_lo, "vs 1"))#> loocv_brute_force loocv_closed_form difference GCV
#> "1.081128" "1.081128" "2.22e-16" "1.083113"
#> brute_seconds shortcut_seconds fits_required
#> "0.09" "0" "120 vs 1"Identical to six decimals, from one fit instead of 120. This is why LOOCV is routine for ridge regression, smoothing splines, and kernel smoothers, and impractical for random forests and neural networks, the shortcut exists only when the fit is linear in \(\mathbf y\).
17 Repeated cross-validation
A single \(k\)-fold partition is one random assignment. Repeating with different partitions and averaging reduces the variance due to the partition, though not the variance due to the data.
set.seed(87)
rep_cv <- function(reps) {
mean(vapply(seq_len(reps), \(r) cv_rmse(d_hv, 10, seed = 5000 + r), numeric(1)))
}
rc <- data.frame(repeats = c(1, 2, 3, 5, 10, 20))
rc$sd_of_estimate <- vapply(rc$repeats, function(R) {
sd(vapply(1:60, function(b) {
mean(vapply(seq_len(R), \(r) cv_rmse(d_hv, 10, seed = b * 1000 + r), numeric(1)))
}, numeric(1)))
}, numeric(1))
rc |> mutate(sd_of_estimate = round(sd_of_estimate, 5))#> repeats sd_of_estimate
#> 1 1 0.01123
#> 2 2 0.00825
#> 3 3 0.00733
#> 4 5 0.00602
#> 5 10 0.00349
#> 6 20 0.00226ggplot(rc, aes(repeats, sd_of_estimate)) +
geom_line(linewidth = 1, color = "steelblue") + geom_point(size = 2.4) +
labs(title = "Repeating 10-fold CV reduces partition variance",
subtitle = "It cannot reduce the variance arising from the single dataset itself",
x = "Number of repeats", y = "SD of the CV estimate") +
theme_dspa()
The curve flattens: repetition removes the partition noise and then stops helping, because the remaining variability comes from having only one sample of \(n\) observations.
No unbiased estimator of the variance of \(k\)-fold CV exists. The fold errors are correlated in a way that depends on the unknown data distribution (Bengio & Grandvalet, 2004). This is not a gap in the literature, it is a proof of impossibility, and it is why comparing models on CV scores needs the corrected test of §9.19.
18 The bootstrap
Resample \(n\) observations with replacement, fit, and predict the omitted ones.
Each bootstrap sample omits a fraction
\[\left(1-\frac1n\right)^{n}\ \xrightarrow[n\to\infty]{}\ e^{-1}\approx0.368,\]
so about 63.2% of observations appear at least once, giving the estimator its name.
set.seed(91)
prop_unique <- function(n, B = 4000)
mean(replicate(B, length(unique(sample(n, n, replace = TRUE))) / n))
data.frame(n = c(10, 50, 200, 1000),
observed_unique = round(vapply(c(10, 50, 200, 1000), prop_unique, numeric(1)), 4),
limit = round(1 - exp(-1), 4))#> n observed_unique limit
#> 1 10 0.6518 0.6321
#> 2 50 0.6346 0.6321
#> 3 200 0.6326 0.6321
#> 4 1000 0.6323 0.632118.1 The .632 and .632+ estimators
The out-of-bag error \(\widehat{\mathrm{Err}}^{(1)}\) is pessimistic, because each model effectively trains on only 63.2% of the data. The apparent error \(\overline{\mathrm{err}}\) is optimistic. The .632 estimator blends them:
\[\boxed{\;\widehat{\mathrm{Err}}^{.632}=0.368\,\overline{\mathrm{err}}+0.632\,\widehat{\mathrm{Err}}^{(1)}\;}\]
The first term is the apparent error, the model’s error on the data it was fitted to. Substituting a held-out error there defeats the whole construction, because both terms then describe the same out-of-sample quantity and there is nothing optimistic left to balance against.
The out-of-bag set must also be built as
setdiff(seq_len(n), boot_indices). Indexing only into1:length(boot_indices)silently excludes the tail of the dataset from ever being out-of-bag.
For learners that overfit badly, .632 is itself optimistic. The .632+ estimator (Efron & Tibshirani, 1997) corrects it using the no-information error rate \(\hat\gamma\), the error when predictions and responses are independent:
\[\hat R=\frac{\widehat{\mathrm{Err}}^{(1)}-\overline{\mathrm{err}}}{\hat\gamma-\overline{\mathrm{err}}}, \qquad \hat w=\frac{0.632}{1-0.368\hat R}, \qquad \widehat{\mathrm{Err}}^{.632+}=(1-\hat w)\,\overline{\mathrm{err}}+\hat w\,\widehat{\mathrm{Err}}^{(1)} .\]
When the learner does not overfit, \(\hat R\approx0\) and \(\hat w\approx0.632\), recovering .632. When it interpolates (\(\overline{\mathrm{err}}=0\)), \(\hat w\to1\) and the estimator becomes the pure out-of-bag error.
boot_632 <- function(dat, formula, fitter, predictor, loss, B = 100, seed = 1) {
set.seed(seed); n <- nrow(dat)
# Apparent error: the model fitted on ALL the data, scored on ALL the data
full_fit <- fitter(formula, dat)
err_bar <- loss(dat, predictor(full_fit, dat))
# Leave-one-out bootstrap: for each observation, average the loss over the
# bootstrap replicates in which it was OUT of bag.
loss_i <- matrix(NA_real_, n, B)
for (b in seq_len(B)) {
idx <- sample(n, n, replace = TRUE)
oob <- setdiff(seq_len(n), idx) # the CORRECT out-of-bag set
if (!length(oob)) next
fb <- fitter(formula, dat[idx, , drop = FALSE])
loss_i[oob, b] <- loss(dat[oob, , drop = FALSE],
predictor(fb, dat[oob, , drop = FALSE]), per_case = TRUE)
}
err_oob <- mean(rowMeans(loss_i, na.rm = TRUE), na.rm = TRUE)
# No-information rate: predictions paired with SHUFFLED responses
yhat_full <- predictor(full_fit, dat)
gamma <- mean(replicate(50, {
d2 <- dat; d2[[all.vars(formula)[1]]] <- sample(d2[[all.vars(formula)[1]]])
loss(d2, yhat_full)
}))
R <- if (gamma > err_bar) (err_oob - err_bar) / (gamma - err_bar) else 0
R <- min(max(R, 0), 1)
w <- 0.632 / (1 - 0.368 * R)
c(apparent = err_bar, oob = err_oob, no_information = gamma,
R = R, weight = w,
err_632 = 0.368 * err_bar + 0.632 * err_oob,
err_632p = (1 - w) * err_bar + w * err_oob)
}library(rpart)
loss_01 <- function(dat, pred, per_case = FALSE) {
z <- as.character(pred) != as.character(dat$y)
if (per_case) as.numeric(z) else mean(z)
}
set.seed(93)
n_bs <- 300
d_bs <- data.frame(matrix(rnorm(n_bs * 6), n_bs, 6))
d_bs$y <- factor(ifelse(d_bs$X1 + d_bs$X2 + rnorm(n_bs) > 0, "a", "b"))
# A heavily overfitting learner (unpruned tree) and a stable one (logistic)
fit_tree <- function(f, d) rpart(f, d, control = rpart.control(cp = 0, minsplit = 2))
pred_tree <- function(m, d) predict(m, d, type = "class")
fit_glm <- function(f, d) glm(f, d, family = binomial())
pred_glm <- function(m, d) factor(ifelse(predict(m, d, type = "response") > 0.5,
"b", "a"), levels = c("a", "b"))
rbind(`unpruned tree` = boot_632(d_bs, y ~ ., fit_tree, pred_tree, loss_01, B = 80),
`logistic` = boot_632(d_bs, y ~ ., fit_glm, pred_glm, loss_01, B = 80)) |>
round(4)#> apparent oob no_information R weight err_632 err_632p
#> unpruned tree 0.0000 0.2933 0.4919 0.5963 0.8097 0.1854 0.2375
#> logistic 0.2033 0.2126 0.4945 0.0319 0.6395 0.2092 0.2093Read the first row. The unpruned tree has apparent error 0, it interpolates its training data, so \(\hat R\) is large, \(\hat w\) moves toward 1, and .632+ lands well above .632. For the logistic model, which does not overfit, the two estimators nearly coincide. That divergence is exactly the correction .632+ was built to supply.
19 Nested cross-validation
Common misconception: “I cross-validated while tuning, so my CV score is an honest estimate.” It is not. Selecting the hyperparameter that maximizes a CV score and then reporting that maximum is selection on the same data, the maximum of a set of noisy estimates is biased upward. The more configurations you search, the larger the bias.
Nested CV separates the two jobs: an inner loop selects hyperparameters within each training set, and an outer loop scores the whole selection-plus-fitting procedure on data neither loop has seen.
set.seed(101)
n_nc <- 250
X_nc <- matrix(rnorm(n_nc * 30), n_nc, 30)
y_nc <- factor(ifelse(X_nc[, 1] + X_nc[, 2] + rnorm(n_nc, sd = 1.2) > 0, "a", "b"))
d_nc <- data.frame(X_nc, y = y_nc)
k_grid <- c(1, 3, 5, 9, 15, 25, 41, 61)
knn_acc <- function(tr, te, k) {
p <- class::knn(tr[, -ncol(tr)], te[, -ncol(te)], tr$y, k = k)
mean(p == te$y)
}
# FLAT: tune and report on the same folds
flat_cv <- function(seed) {
set.seed(seed)
f <- sample(rep(1:5, length.out = n_nc))
acc <- vapply(k_grid, function(k)
mean(vapply(1:5, \(i) knn_acc(d_nc[f != i, ], d_nc[f == i, ], k), numeric(1))),
numeric(1))
max(acc) # the reported number
}
# NESTED: inner loop tunes, outer loop scores untouched data
nested_cv <- function(seed) {
set.seed(seed)
outer_f <- sample(rep(1:5, length.out = n_nc))
mean(vapply(1:5, function(i) {
tr <- d_nc[outer_f != i, ]; te <- d_nc[outer_f == i, ]
inner_f <- sample(rep(1:5, length.out = nrow(tr)))
inner <- vapply(k_grid, function(k)
mean(vapply(1:5, \(j) knn_acc(tr[inner_f != j, ], tr[inner_f == j, ], k),
numeric(1))), numeric(1))
knn_acc(tr, te, k_grid[which.max(inner)])
}, numeric(1)))
}
set.seed(103)
flat <- vapply(1:20, flat_cv, numeric(1))
nested <- vapply(1:20, nested_cv, numeric(1))
c(flat_CV_reported = round(mean(flat), 4),
nested_CV_estimate = round(mean(nested), 4),
optimism = round(mean(flat) - mean(nested), 4))#> flat_CV_reported nested_CV_estimate optimism
#> 0.674 0.650 0.024The flat estimate is optimistic by a measurable margin. The gap grows with the size of the search, a point worth seeing as a surface:
grid_sizes <- c(2, 4, 8, 16)
n_sizes <- c(80, 150, 250, 400)
optimism_at <- function(n, g, reps = 6) {
set.seed(n + g)
kg <- round(seq(1, 61, length.out = g))
X <- matrix(rnorm(n * 30), n, 30)
y <- factor(ifelse(X[, 1] + X[, 2] + rnorm(n, sd = 1.2) > 0, "a", "b"))
d <- data.frame(X, y = y)
vapply(seq_len(reps), function(r) {
f <- sample(rep(1:5, length.out = n))
acc <- vapply(kg, function(k)
mean(vapply(1:5, function(i) {
p <- class::knn(d[f != i, -ncol(d)], d[f == i, -ncol(d)], d$y[f != i], k = k)
mean(p == d$y[f == i])
}, numeric(1))), numeric(1))
max(acc) - mean(acc) # optimism from taking the max
}, numeric(1)) |> mean()
}
Zno <- outer(grid_sizes, n_sizes, Vectorize(optimism_at))
plot_ly(x = n_sizes, y = grid_sizes, z = Zno, type = "surface",
colorscale = "Inferno", reversescale = TRUE,
colorbar = list(title = "Selection\noptimism")) |>
layout(title = "Optimism from reporting the best of several CV scores",
scene = list(xaxis = list(title = "Sample size n"),
yaxis = list(title = "Hyperparameter grid size"),
zaxis = list(title = "max(CV) - mean(CV)")))The surface rises with the grid size and falls with \(n\). Searching a large space on a small dataset is precisely where flat CV misleads most.
20 When random folds are wrong
\(k\)-fold CV assumes observations are exchangeable. Two common structures break that assumption, and each has its own remedy.
20.1 Grouped data
Repeated measures on the same subject, multiple samples from one site, or multiple lesions in one patient are not independent. Random folds place a subject’s observations on both sides of the split, so the model can memorize the subject rather than learn the signal.
set.seed(111)
n_sub <- 40; per_sub <- 5
subj <- rep(seq_len(n_sub), each = per_sub)
sub_effect <- rnorm(n_sub, sd = 2) # strong subject effect
y_g <- factor(ifelse(sub_effect[subj] + rnorm(n_sub * per_sub, sd = 0.5) > 0,
"a", "b"))
Xg2 <- cbind(sub_effect[subj] + rnorm(n_sub * per_sub, sd = 0.5),
matrix(rnorm(n_sub * per_sub * 3), ncol = 3))
d_g <- data.frame(Xg2, subject = subj, y = y_g)
acc_folds <- function(folds) {
mean(vapply(folds, function(idx_te) {
tr <- d_g[-idx_te, ]; te <- d_g[idx_te, ]
p <- class::knn(tr[, 1:4], te[, 1:4], tr$y, k = 5)
mean(p == te$y)
}, numeric(1)))
}
set.seed(113)
f_rand <- split(sample(nrow(d_g)), rep(1:5, length.out = nrow(d_g)))
grp <- split(seq_len(n_sub), rep(1:5, length.out = n_sub))
f_grp <- lapply(grp, \(s) which(d_g$subject %in% s))
c(random_folds = round(acc_folds(f_rand), 4),
grouped_folds = round(acc_folds(f_grp), 4),
optimism = round(acc_folds(f_rand) - acc_folds(f_grp), 4))#> random_folds grouped_folds optimism
#> 0.835 0.825 0.010Random folds report substantially higher accuracy on data where the
true signal lives entirely at the subject level. Split by
group, using rsample::group_vfold_cv().
20.2 Time series
Randomly splitting a time series puts day \(t\) in training and day \(t+1\) in assessment. On autocorrelated data those are nearly duplicates, so the estimate measures interpolation rather than forecasting. The remedy is forward chaining: train on \([1,\tau]\), assess on \((\tau,\tau+h]\), and roll \(\tau\) forward.
set.seed(121)
T_len <- 400
ar <- as.numeric(arima.sim(list(ar = 0.95), T_len))
noise_pred <- as.numeric(arima.sim(list(ar = 0.95), T_len)) # UNRELATED to ar
d_t <- data.frame(x = noise_pred, y = ar)
# Random 5-fold
set.seed(123)
f_t <- sample(rep(1:5, length.out = T_len))
r2_random <- mean(vapply(1:5, function(i) {
fit <- lm(y ~ x, d_t[f_t != i, ]); te <- d_t[f_t == i, ]
1 - sum((predict(fit, te) - te$y)^2) / sum((te$y - mean(te$y))^2)
}, numeric(1)))
# Rolling origin
origins <- round(seq(0.5 * T_len, T_len - 40, length.out = 5))
r2_rolling <- mean(vapply(origins, function(tau) {
fit <- lm(y ~ x, d_t[1:tau, ]); te <- d_t[(tau + 1):(tau + 40), ]
1 - sum((predict(fit, te) - te$y)^2) / sum((te$y - mean(te$y))^2)
}, numeric(1)))
c(random_folds_R2 = round(r2_random, 4),
rolling_origin_R2 = round(r2_rolling, 4),
truth = "0 -- x and y are independent AR(1) series")#> random_folds_R2
#> "-0.0118"
#> rolling_origin_R2
#> "-0.2156"
#> truth
#> "0 -- x and y are independent AR(1) series"There is no relationship between x and
y, both are independent AR(1) processes. Random folds
report a positive \(R^2\); rolling
origin does not.
ggplot(data.frame(t = 1:T_len, y = ar), aes(t, y)) +
geom_line(linewidth = 0.4, color = "grey40") +
annotate("rect", xmin = origins[3] + 1, xmax = origins[3] + 40,
ymin = -Inf, ymax = Inf, alpha = 0.2, fill = "#D8433B") +
annotate("rect", xmin = 1, xmax = origins[3],
ymin = -Inf, ymax = Inf, alpha = 0.12, fill = "#3B7DD8") +
labs(title = "Rolling-origin evaluation of a time series",
subtitle = "Blue: training window. Red: the forecast horizon that follows it. Nothing is shuffled",
x = "Time", y = "y") +
theme_dspa()
21 Computational cost of each scheme
Let \(C(m)\) be the cost of one model fit on \(m\) observations.
| Scheme | Fits | Cost | Note |
|---|---|---|---|
| Holdout | \(1\) | \(C(0.7n)\) | Cheapest; highest variance |
| \(k\)-fold CV | \(k\) | \(k\,C\!\left(n\tfrac{k-1}{k}\right)\) | \(k=5\) or \(10\) standard |
| Repeated \(k\)-fold | \(rk\) | \(rk\,C(\cdot)\) | Reduces partition variance only |
| LOOCV (general) | \(n\) | \(n\,C(n-1)\) | Impractical for most learners |
| LOOCV (linear smoother) | \(\mathbf{1}\) | \(C(n)+O(n)\) | Hat-value shortcut |
| GCV | \(1\) | \(C(n)+O(1)\) | Uses \(\operatorname{tr}(S)/n\) |
| Bootstrap OOB / .632 | \(B\) | \(B\,C(n)\) | \(B=50\)–\(200\) typical |
| Nested CV | \(k_{\text{out}}\!\left(1+k_{\text{in}}|\Lambda|\right)\) | — | The honest cost of honest tuning |
| Grouped \(k\)-fold | \(k\) | as \(k\)-fold | Folds sized by group, not row |
| Rolling origin | \(\#\)origins | \(\sum_\tau C(\tau)\) | Training sets grow |
cost_table <- function(k_out = 5, k_in = 5, grid = 10, repeats = 3, n = 1000) {
data.frame(
scheme = c("Holdout", "5-fold", "10-fold", "3x10-fold repeated",
"LOOCV (general)", "LOOCV (linear smoother)",
"Bootstrap B=100", "Nested 5x5, grid 10"),
model_fits = c(1, 5, 10, repeats * 10, n, 1, 100,
k_out * (1 + k_in * grid)))
}
cost_table()#> scheme model_fits
#> 1 Holdout 1
#> 2 5-fold 5
#> 3 10-fold 10
#> 4 3x10-fold repeated 30
#> 5 LOOCV (general) 1000
#> 6 LOOCV (linear smoother) 1
#> 7 Bootstrap B=100 100
#> 8 Nested 5x5, grid 10 255Nested CV with a modest grid costs 255 fits where flat CV costs 5. That factor is the price of an unbiased estimate, and it is why the honest protocol is often reserved for the final reported number while flat CV guides exploration.
22 PART IV: COMPARISON AND IMPROVEMENT
23 Comparing models
Two models produce CV scores differing by two percentage points. Is that a difference?
23.2 A worked comparison
qol <- dspa_read("https://umich.instructure.com/files/481332/download?download_frd=1",
"Case06_QoL_Symptom_ChronicIllness.csv")
qol <- qol |> filter(CHRONICDISEASESCORE != -9, CHARLSONSCORE != -9)
# The threshold is COMPUTED from the data. A median split gives balanced
# classes, so accuracy is interpretable against a ~50% no-information rate.
cut_cds <- median(qol$CHRONICDISEASESCORE)
qol$cd <- factor(qol$CHRONICDISEASESCORE > cut_cds,
levels = c(FALSE, TRUE),
labels = c("minor_disease", "severe_disease"))
c(cut_point = round(cut_cds, 4))#> cut_point
#> 1.39#>
#> minor_disease severe_disease
#> 0.5009 0.4991# The feature set is defined ONCE, by name, and reused by every model.
# CHRONICDISEASESCORE determines `cd` exactly -- leaving it in would make
# every accuracy 1. ID and INTERVIEWDATE are administrative.
drop_cols <- c("ID", "CHRONICDISEASESCORE", "INTERVIEWDATE")
qol_model_df <- qol |> select(-any_of(drop_cols))
c(dropped = paste(intersect(drop_cols, names(qol)), collapse = ", "),
predictors = ncol(qol_model_df) - 1,
leakage_check_correlation_with_outcome = "CHRONICDISEASESCORE excluded")#> dropped
#> "ID, CHRONICDISEASESCORE, INTERVIEWDATE"
#> predictors
#> "38"
#> leakage_check_correlation_with_outcome
#> "CHRONICDISEASESCORE excluded"set.seed(1234)
sp_qol <- rsample::initial_split(qol_model_df, prop = 0.8, strata = cd)
qol_train <- rsample::training(sp_qol); qol_test <- rsample::testing(sp_qol)
NIR_qol <- max(prop.table(table(qol_test$cd)))
c(train = nrow(qol_train), test = nrow(qol_test),
no_information_rate = round(NIR_qol, 4))#> train test no_information_rate
#> 1751.0000 439.0000 0.5011# library(doParallel)
# n_cores <- max(1, parallel::detectCores() - 1)
# cl <- makePSOCKcluster(n_cores)
# registerDoParallel(cl)
# on.exit({ stopCluster(cl); registerDoSEQ() }, add = TRUE) # never orphan workers
#
# set.seed(1234)
# ctrl <- trainControl(method = "repeatedcv", number = 10, repeats = 3,
# classProbs = TRUE, summaryFunction = twoClassSummary,
# savePredictions = "final")
#
# fits <- list(
# RF = train(cd ~ ., data = qol_train, method = "ranger", metric = "ROC",
# trControl = ctrl, tuneLength = 3, num.trees = 300, num.threads = 1),
# kNN = train(cd ~ ., data = qol_train, method = "knn", metric = "ROC",
# trControl = ctrl, tuneLength = 6,
# preProcess = c("center", "scale")),
# SVM = train(cd ~ ., data = qol_train, method = "svmRadial", metric = "ROC",
# trControl = ctrl, tuneLength = 3,
# preProcess = c("center", "scale")),
# GBM = train(cd ~ ., data = qol_train, method = "gbm", metric = "ROC",
# trControl = ctrl, tuneLength = 3, verbose = FALSE),
# GLM = train(cd ~ ., data = qol_train, method = "glm", family = "binomial",
# metric = "ROC", trControl = ctrl))
# stopCluster(cl); registerDoSEQ()
library(doParallel)
n_cores <- max(1, parallel::detectCores() - 1)
cl <- makePSOCKcluster(n_cores)
registerDoParallel(cl)
set.seed(1234)
ctrl <- trainControl(method = "repeatedcv", number = 10, repeats = 3,
classProbs = TRUE, summaryFunction = twoClassSummary,
savePredictions = "final")
fits <- list(
RF = train(cd ~ ., data = qol_train, method = "ranger", metric = "ROC",
trControl = ctrl, tuneLength = 3, num.trees = 300, num.threads = 1),
kNN = train(cd ~ ., data = qol_train, method = "knn", metric = "ROC",
trControl = ctrl, tuneLength = 6, preProcess = c("center","scale")),
SVM = train(cd ~ ., data = qol_train, method = "svmRadial", metric = "ROC",
trControl = ctrl, tuneLength = 3, preProcess = c("center","scale")),
GBM = train(cd ~ ., data = qol_train, method = "gbm", metric = "ROC",
trControl = ctrl, tuneLength = 3, verbose = FALSE),
GLM = train(cd ~ ., data = qol_train, method = "glm", family = "binomial",
metric = "ROC", trControl = ctrl)
)
# Robust cleanup — workers may already be gone
tryCatch(stopCluster(cl), error = function(e) NULL)
tryCatch(registerDoSEQ(), error = function(e) NULL)The preProcess argument does real work:
caret recomputes the centering and scaling inside
every resample, using only that resample’s analysis set.
Standardizing the whole frame beforehand would leak (Chapter 5, §5.4.1).
#> Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
#> RF 0.6257 0.6604 0.6799 0.6837 0.7019 0.7717 0
#> kNN 0.5114 0.5282 0.5588 0.5601 0.5850 0.6342 0
#> SVM 0.5270 0.6139 0.6274 0.6322 0.6604 0.7159 0
#> GBM 0.6123 0.6676 0.6925 0.6935 0.7236 0.7729 0
#> GLM 0.5272 0.6128 0.6232 0.6265 0.6509 0.6981 0res_long <- res$values |>
select(Resample, ends_with("~ROC")) |>
pivot_longer(-Resample, names_to = "model", values_to = "ROC") |>
mutate(model = sub("~ROC$", "", model))
ggplot(res_long, aes(reorder(model, ROC, median), ROC)) +
geom_boxplot(fill = "#9EC5E8", width = 0.55, outlier.alpha = 0.4) +
geom_jitter(width = 0.08, alpha = 0.25, size = 0.8) +
coord_flip() +
labs(title = "Cross-validated AUC across 30 common resamples",
subtitle = "Identical folds for every model, so the comparison is paired",
x = NULL, y = "AUC") +
theme_dspa()
# --- Interactive equivalents, plus caret's lattice displays ----------------
plot_ly(res_long, x = ~model, y = ~ROC, color = ~model, type = "box") |>
layout(title = "Cross-validated AUC by model", showlegend = FALSE)
dens <- with(res_long, tapply(ROC, INDEX = model, density))
df_d <- data.frame(
x = unlist(lapply(dens, "[[", "x")),
y = unlist(lapply(dens, "[[", "y")),
model = rep(names(dens), each = length(dens[[1]]$x)))
plot_ly(df_d, x = ~x, y = ~y, color = ~model) |> add_lines() |>
layout(title = "Performance density plots (AUC)",
xaxis = list(title = "AUC"), yaxis = list(title = "Density"),
legend = list(orientation = "h"))
res_wide <- res_long |> pivot_wider(names_from = model, values_from = ROC)
dims <- purrr::map2(select(res_wide, -Resample),
names(select(res_wide, -Resample)),
~list(values = .x, label = .y))
plot_ly(type = "splom", dimensions = setNames(dims, NULL),
showupperhalf = FALSE, diagonal = list(visible = FALSE)) |>
layout(title = "Pairwise resample performance (AUC)")wide <- res_long |> pivot_wider(names_from = model, values_from = ROC)
model_names <- setdiff(names(wide), "Resample")
pairs_df <- do.call(rbind, combn(model_names, 2, simplify = FALSE, FUN = function(pr) {
d <- wide[[pr[1]]] - wide[[pr[2]]]
naive <- t.test(d)
corr <- corrected_t_test(d, k = 10)
data.frame(comparison = paste(pr, collapse = " - "),
mean_diff = round(mean(d), 4),
naive_p = signif(naive$p.value, 3),
corrected_p = signif(corr[["p_value"]], 3))
}))
pairs_df$corrected_q <- signif(p.adjust(pairs_df$corrected_p, "holm"), 3)
pairs_df#> comparison mean_diff naive_p corrected_p corrected_q
#> 1 RF - kNN 0.1237 1.31e-14 1.67e-07 1.50e-06
#> 2 RF - SVM 0.0515 4.38e-05 2.84e-02 8.52e-02
#> 3 RF - GBM -0.0098 3.21e-01 6.31e-01 1.00e+00
#> 4 RF - GLM 0.0572 5.69e-07 4.69e-03 2.34e-02
#> 5 kNN - SVM -0.0722 8.92e-08 2.00e-03 1.40e-02
#> 6 kNN - GBM -0.1335 1.43e-15 3.31e-08 3.31e-07
#> 7 kNN - GLM -0.0664 3.56e-07 3.80e-03 2.28e-02
#> 8 SVM - GBM -0.0613 4.88e-06 1.18e-02 4.72e-02
#> 9 SVM - GLM 0.0057 5.93e-01 7.97e-01 1.00e+00
#> 10 GBM - GLM 0.0670 5.16e-08 1.54e-03 1.23e-02Three corrections stack here, and each matters. The corrected \(t\) widens the interval for fold overlap. Holm’s adjustment accounts for testing all \(\binom{5}{2}=10\) pairs. And the naive \(p\)-values are systematically smaller than the corrected ones, which is the anti-conservatism measured in the null simulation above, now visible on real data.
holdout_tab <- data.frame(
model = c("Majority class", names(fits)),
accuracy = round(c(NIR_qol, vapply(fits,
\(m) mean(predict(m, qol_test) == qol_test$cd), numeric(1))), 4),
AUC = round(c(0.5, vapply(fits, \(m) as.numeric(pROC::auc(pROC::roc(
qol_test$cd, predict(m, qol_test, type = "prob")[, "severe_disease"],
quiet = TRUE))), numeric(1))), 4),
Brier = round(c(NA, vapply(fits, \(m) mean((
predict(m, qol_test, type = "prob")[, "severe_disease"] -
as.integer(qol_test$cd == "severe_disease"))^2), numeric(1))), 4))
holdout_tab$calib_slope <- round(c(NA, vapply(fits, function(m) {
p <- pmin(pmax(predict(m, qol_test, type = "prob")[, "severe_disease"],
1e-6), 1 - 1e-6)
unname(coef(glm(as.integer(qol_test$cd == "severe_disease") ~ qlogis(p),
family = binomial()))[2])
}, numeric(1))), 3)
holdout_tab$acc_CI <- c(NA, vapply(fits, function(m) {
ci <- binom.test(sum(predict(m, qol_test) == qol_test$cd),
nrow(qol_test))$conf.int
sprintf("[%.3f, %.3f]", ci[1], ci[2])
}, character(1)))
holdout_tab#> model accuracy AUC Brier calib_slope acc_CI
#> Majority class 0.5011 0.5000 NA NA <NA>
#> RF RF 0.6879 0.7236 0.2107 1.097 [0.642, 0.731]
#> kNN kNN 0.5649 0.6053 0.2427 0.680 [0.519, 0.614]
#> SVM SVM 0.6424 0.6969 0.2244 1.736 [0.596, 0.687]
#> GBM GBM 0.6879 0.7490 0.2051 1.299 [0.642, 0.731]
#> GLM GLM 0.6173 0.6753 0.2295 0.996 [0.570, 0.663]Read the last three columns together. AUC ranks the models; the Brier score and calibration slope say whether their probabilities can be believed; and the confidence interval says whether the AUC differences are resolvable at this sample size. A model can lead on AUC and trail badly on calibration, which is a reason to prefer a different one if the output feeds a decision rule.
24 Learning curves
A disappointing CV score has two very different causes, and the remedy differs completely. A learning curve, training and validation error against training-set size, tells you which you have.
| Pattern | Diagnosis | Remedy |
|---|---|---|
| Both curves high, converged, small gap | High bias (underfitting) | More features, more capacity, less regularization |
| Training low, validation high, wide gap still closing | High variance (overfitting) | More data, more regularization, fewer features |
| Both low, converged | Model is doing what it can | Look at the irreducible noise |
set.seed(141)
n_lc <- 800; p_lc <- 12
X_lc <- matrix(rnorm(n_lc * p_lc), n_lc, p_lc)
# A nonlinear truth, so a linear model is bias-limited
y_lc <- X_lc[, 1]^2 + 2 * X_lc[, 2] + 0.5 * X_lc[, 3] * X_lc[, 4] +
rnorm(n_lc, sd = 1)
d_lc <- data.frame(X_lc, y = y_lc)
te_idx <- 601:800; tr_pool <- 1:600
learning_curve <- function(fitter, label, sizes = c(20, 40, 80, 150, 300, 600)) {
do.call(rbind, lapply(sizes, function(m) {
reps <- vapply(1:25, function(r) {
set.seed(r * 100 + m)
idx <- sample(tr_pool, m)
fit <- fitter(d_lc[idx, ])
c(train = sqrt(mean((predict(fit, d_lc[idx, ]) - d_lc$y[idx])^2)),
valid = sqrt(mean((predict(fit, d_lc[te_idx, ]) - d_lc$y[te_idx])^2)))
}, numeric(2))
data.frame(m = m, Training = mean(reps["train", ]),
Validation = mean(reps["valid", ]), model = label)
}))
}
lc <- bind_rows(
learning_curve(\(d) lm(y ~ ., data = d), "Linear (high bias)"),
learning_curve(\(d) rpart::rpart(y ~ ., d,
control = rpart::rpart.control(cp = 0, minsplit = 2)), "Deep tree (high variance)"))
lc |> pivot_longer(c(Training, Validation), names_to = "set", values_to = "rmse") |>
ggplot(aes(m, rmse, color = set)) +
geom_line(linewidth = 1) + geom_point(size = 2) +
facet_wrap(~ model) +
scale_x_log10() +
scale_color_manual(values = c(Training = "#7FB069", Validation = "#3B7DD8")) +
labs(title = "Learning curves diagnose bias versus variance",
subtitle = "Left: curves converge high and close -- more data will not help. Right: wide gap still closing -- it will",
x = "Training set size (log scale)", y = "RMSE", color = NULL) +
theme_dspa(10)
# --- Interactive equivalent ------------------------------------------------
d <- filter(lc, model == "Deep tree (high variance)")
plot_ly(d, x = ~m, y = ~Training, type = "scatter", mode = "lines+markers",
name = "Training") |>
add_trace(y = ~Validation, name = "Validation") |>
layout(title = "Learning curve",
xaxis = list(title = "Training set size", type = "log"),
yaxis = list(title = "RMSE"), legend = list(orientation = "h"))The two panels prescribe opposite actions. On the left, training and validation error have converged at a high level, the model class cannot represent the truth, and collecting more data changes nothing. On the right the gap is wide and still narrowing, more data will help, and so will regularization.
The full picture is a surface over (training size, model complexity):
sizes_ls <- c(25, 50, 100, 200, 400, 600)
depths <- c(1, 2, 3, 5, 8, 15)
Zls <- outer(depths, sizes_ls, Vectorize(function(dp, m) {
mean(vapply(1:12, function(r) {
set.seed(r * 77 + m + dp)
idx <- sample(tr_pool, m)
fit <- rpart::rpart(y ~ ., d_lc[idx, ],
control = rpart::rpart.control(maxdepth = dp, cp = 0, minsplit = 5))
sqrt(mean((predict(fit, d_lc[te_idx, ]) - d_lc$y[te_idx])^2))
}, numeric(1)))
}))
plot_ly(x = sizes_ls, y = depths, z = Zls, type = "surface",
colorscale = "Viridis", reversescale = TRUE,
colorbar = list(title = "Validation\nRMSE")) |>
layout(title = "Validation error over training size and tree depth",
scene = list(xaxis = list(title = "Training set size"),
yaxis = list(title = "Maximum tree depth"),
zaxis = list(title = "Validation RMSE")))Rotate to the small-sample edge: the optimal depth there is shallow, and it grows as data accumulate. The right complexity is a function of the sample size, which is why a hyperparameter tuned on one dataset does not transfer to a smaller one.
25 Tuning
25.1 Grid versus random search
Grid search evaluates every combination: \(\prod_j m_j\) fits for \(m_j\) values of parameter \(j\), exponential in the number of parameters. Random search draws configurations from distributions over the same space.
Random search is usually better for the same budget (Bergstra & Bengio, 2012). The reason is effective dimensionality: most hyperparameters barely matter, and a grid spends its budget re-evaluating the same few values of the ones that do. With \(m\) points on a grid over \(d\) parameters, only \(m^{1/d}\) distinct values of each parameter are tried; random search tries \(m\) distinct values of each.
set.seed(151)
# A response depending strongly on one parameter and weakly on another
f_hyper <- function(a, b) -( (a - 0.3)^2 * 10 + (b - 0.7)^2 * 0.15 )
budget <- 36
g <- expand.grid(a = seq(0, 1, length.out = 6), b = seq(0, 1, length.out = 6))
r <- data.frame(a = runif(budget), b = runif(budget))
c(grid_best = round(max(f_hyper(g$a, g$b)), 4),
random_best = round(max(f_hyper(r$a, r$b)), 4),
distinct_a_values_grid = length(unique(g$a)),
distinct_a_values_random = length(unique(r$a)))#> grid_best random_best distinct_a_values_grid
#> -0.1015 -0.0013 6.0000
#> distinct_a_values_random
#> 36.0000bind_rows(mutate(g, search = "Grid (6 x 6)"),
mutate(r, search = "Random (36 draws)")) |>
ggplot(aes(a, b)) +
geom_raster(data = expand.grid(a = seq(0, 1, 0.01), b = seq(0, 1, 0.01)) |>
mutate(z = f_hyper(a, b)), aes(fill = z), alpha = 0.75) +
geom_point(size = 1.8, color = "white") +
facet_wrap(~ search) +
scale_fill_viridis_c(guide = "none") +
coord_fixed() +
labs(title = "Grid search wastes its budget on the parameter that does not matter",
subtitle = "Objective varies steeply in a and weakly in b. Grid tries 6 values of a; random tries 36",
x = "a (important)", y = "b (unimportant)") +
theme_dspa(10)
25.2 The one-standard-error rule
Selecting the configuration with the best CV score picks the maximum of a set of noisy estimates. The 1-SE rule picks the simplest model whose score is within one standard error of the best, trading a statistically indistinguishable amount of performance for a smaller, more stable model.
set.seed(161)
ctrl_1se <- trainControl(method = "repeatedcv", number = 10, repeats = 3,
selectionFunction = "oneSE")
grid_k <- expand.grid(k = seq(1, 41, by = 4))
knn_best <- train(cd ~ ., data = qol_train, method = "knn", metric = "Kappa",
trControl = trainControl(method = "repeatedcv", number = 10,
repeats = 3),
tuneGrid = grid_k, preProcess = c("center", "scale"))
knn_1se <- train(cd ~ ., data = qol_train, method = "knn", metric = "Kappa",
trControl = ctrl_1se, tuneGrid = grid_k,
preProcess = c("center", "scale"))
c(best_k = knn_best$bestTune$k, oneSE_k = knn_1se$bestTune$k)#> best_k oneSE_k
#> 37 41kres <- knn_best$results
best_row <- kres[which.max(kres$Kappa), ]
thresh <- best_row$Kappa - best_row$KappaSD / sqrt(30)
ggplot(kres, aes(k, Kappa)) +
geom_ribbon(aes(ymin = Kappa - KappaSD / sqrt(30),
ymax = Kappa + KappaSD / sqrt(30)), fill = "grey86") +
geom_line(linewidth = 1, color = "steelblue") + geom_point(size = 2) +
geom_hline(yintercept = thresh, linetype = "dashed", color = "firebrick") +
geom_point(data = best_row, color = "firebrick", size = 4) +
geom_point(data = kres[kres$k == knn_1se$bestTune$k, ],
color = "darkgreen", size = 4, shape = 17) +
labs(title = "The one-standard-error rule",
subtitle = "Red: best score. Dashed: one SE below it. Green triangle: simplest model above the line",
x = "k", y = "Cross-validated Kappa") +
theme_dspa()
# --- Interactive equivalent ------------------------------------------------
plot_ly(kres, x = ~k, y = ~Kappa, type = "scatter", mode = "lines+markers",
error_y = ~list(array = KappaSD / sqrt(30)), name = "CV Kappa") |>
add_lines(x = range(kres$k), y = rep(thresh, 2), name = "1-SE threshold",
line = list(dash = "dash", color = "red")) |>
layout(title = "Tuning k with the 1-SE rule",
xaxis = list(title = "k"), yaxis = list(title = "Kappa"))25.3 Leakage in tuning: a worked example
boystown <- dspa_read(
"https://umich.instructure.com/files/399119/download?download_frd=1",
"CaseStudy02_Boystown_Data.csv", sep = " ")
boystown$sex <- boystown$sex - 1
boystown$dadjob <- 2 - boystown$dadjob
boystown$momjob <- 2 - boystown$momjob
# The outcome is DERIVED from gpa
boystown$grade <- factor(boystown$gpa %in% c(3, 4, 5),
levels = c(FALSE, TRUE),
labels = c("above_avg", "avg_or_below"))
c(outcome_definition = "grade = gpa %in% c(3,4,5)",
correlation_gpa_grade = round(cor(boystown$gpa,
as.integer(boystown$grade) - 1), 4))#> outcome_definition correlation_gpa_grade
#> "grade = gpa %in% c(3,4,5)" "0.8243"grade is a deterministic function of
gpa. Any model given gpa can reproduce the
outcome exactly.
norm01 <- function(x) (x - min(x)) / (max(x) - min(x))
bt_leaky <- boystown |> select(-any_of(c("id", "ID"))) |>
mutate(across(where(is.numeric), norm01))
bt_clean <- bt_leaky |> select(-gpa) # remove the outcome's source
set.seed(171)
ctrl_bt <- trainControl(method = "repeatedcv", number = 10, repeats = 3)
grid_bt <- expand.grid(k = c(1, 3, 5, 7, 9, 15, 21))
m_leaky <- train(grade ~ ., data = bt_leaky, method = "knn",
trControl = ctrl_bt, tuneGrid = grid_bt)
m_clean <- train(grade ~ ., data = bt_clean, method = "knn",
trControl = ctrl_bt, tuneGrid = grid_bt)
data.frame(
feature_set = c("with gpa (leaky)", "without gpa (honest)"),
best_k = c(m_leaky$bestTune$k, m_clean$bestTune$k),
cv_accuracy = round(c(max(m_leaky$results$Accuracy),
max(m_clean$results$Accuracy)), 4),
cv_kappa = round(c(max(m_leaky$results$Kappa),
max(m_clean$results$Kappa)), 4),
no_information_rate = round(max(prop.table(table(boystown$grade))), 4))#> feature_set best_k cv_accuracy cv_kappa no_information_rate
#> 1 with gpa (leaky) 7 0.8683 0.6702 0.67
#> 2 without gpa (honest) 21 0.6324 0.0253 0.67Common misconception: “cross-validation protects against leakage.” It protects against overfitting the model to the training rows. It cannot protect against a feature that is the outcome. Every fold sees
gpa, every fold recoversgradefrom it, and the CV estimate is high and entirely honest about the wrong question, it correctly estimates how well the procedure predictsgradegivengpa, which is a question no one has.The check is not statistical but causal: would this feature be available, and mean the same thing, at the moment a prediction is actually needed?
26 PART V: SYNTHESIS
27 Computational complexity summary
\(n\) = observations, \(d\) = features, \(k\) = folds, \(r\) = repeats, \(B\) = bootstrap replicates, \(|\Lambda|\) = hyperparameter grid size, \(C(m)\) = cost of one fit on \(m\) rows.
| Task | Fits | Extra cost | Note |
|---|---|---|---|
| Training error | \(1\) | \(O(n)\) | Optimistic by \(2d\sigma^2/n\) |
| \(C_p\) / AIC / BIC | \(1\) | \(O(1)\) | Requires an estimate of \(d\) and \(\sigma^2\) |
| Effective df (linear smoother) | \(1\) | \(O(n^2)\) for \(\operatorname{tr}(S)\) | \(O(nd)\) if \(S\) is never formed |
| Holdout | \(1\) | — | Unbiased, high variance |
| \(k\)-fold CV | \(k\) | — | \(k=5\) or \(10\) standard |
| Repeated \(k\)-fold | \(rk\) | — | Reduces partition variance only |
| LOOCV (general) | \(n\) | — | Impractical beyond small \(n\) |
| LOOCV (linear smoother) | \(\mathbf{1}\) | \(O(n)\) for hat values | Sherman–Morrison shortcut |
| GCV | \(1\) | \(O(1)\) | Uses \(\operatorname{tr}(S)/n\) |
| Bootstrap OOB / .632 / .632+ | \(B\) | \(O(Bn)\) | \(B=50\)–\(200\); plus \(\sim50\) fits for \(\hat\gamma\) |
| Nested CV | \(k_{\text{out}}(1+k_{\text{in}}\lvert\Lambda\rvert)\) | — | 255 fits at \(5\times5\times10\) |
| Grouped \(k\)-fold | \(k\) | — | Fold size set by groups |
| Rolling origin | #origins | \(\sum_\tau C(\tau)\) | Training sets grow |
| ROC / AUC | — | \(O(n\log n)\) | Sorting the scores |
| PR curve / PR-AUC | — | \(O(n\log n)\) | Same sort |
| Reliability diagram | — | \(O(n)\) | After binning |
| Platt scaling | — | \(O(n)\) per IRLS step | One logistic fit |
| Isotonic regression | — | \(O(n)\) | Pool-adjacent-violators after sorting |
| Decision curve | — | \(O(n\lvert T\rvert)\) | \(\lvert T\rvert\) thresholds |
| Corrected resampled \(t\) | — | \(O(k)\) | Variance inflation only |
Four consequences.
The LOOCV shortcut is a three-order-of-magnitude saving, and only for linear smoothers. \(\hat{\mathbf y}=S\mathbf y\) with \(S\) independent of \(\mathbf y\) is the precondition. Ridge, smoothing splines, and kernel smoothers qualify; random forests, boosting, and neural networks do not.
Nested CV costs a factor of \(k_{\text{in}}\lvert\Lambda\rvert\) over flat CV. That is the price of an unbiased estimate after tuning, and it is why flat CV is appropriate for exploration while nested CV is reserved for the number you report.
Metrics are cheap; resampling is not. Every curve and score in Part II costs \(O(n\log n)\) at worst. All of the expense is in refitting.
Recalibration is nearly free. One extra logistic fit or one isotonic pass on held-out predictions, and it improves the Brier score and calibration slope without touching AUC.
28 Common pitfalls
| # | Pitfall | Consequence | Fix |
|---|---|---|---|
| 1 | Reporting training error | Optimistic by \(2d\sigma^2/n\) | Resample, or correct with \(C_p\)/AIC |
| 2 | Counting parameters as complexity | Wrong for regularized or ensemble models | Effective df \(=\frac1{\sigma^2}\sum\operatorname{Cov}(\hat y_i,y_i)\) |
| 3 | confusionMatrix(truth, pred) |
Transposes; swaps sensitivity and specificity | confusionMatrix(data = pred, reference = truth, positive = ...) |
| 4 | Calling \(1-\)precision the Type I error | Confuses FDR with \(\alpha\) | Type I conditions on truth; FDR on the prediction |
| 5 | Writing \(1-\text{recall}=FN/(TN+FN)\) | That is the false omission rate | \(1-\text{recall}=FN/(TP+FN)\) |
| 6 | Reading PPV as a property of the test | PPV depends on prevalence | Sens/spec are the test; PPV/NPV are the test in a population |
| 7 | Interpreting \(\kappa\) without the prevalence | Same test gives very different \(\kappa\) | Report prevalence; consider MCC |
| 8 | Assuming \(\kappa\in[0,1]\) | \(\kappa\) can be negative | Range is \([-1,1]\); the floor depends on marginals |
| 9 | Weighted \(\kappa\) on nominal categories | Nothing to weight; needs an ordering | Weighted \(\kappa\) is for ordinal outcomes |
| 10 | Reporting ROC AUC beside a PR curve | Different quantities, different chance levels | PRROC::pr.curve() /
yardstick::pr_auc() |
| 11 | Using ROC AUC under heavy imbalance | \(FP\) barely moves the FPR | PR curve; report prevalence as the PR chance level |
| 12 | Treating AUC as evidence of good probabilities | AUC is invariant to monotone distortion | Reliability diagram, Brier, calibration slope |
| 13 | Optimizing accuracy to pick probabilities | Accuracy is not a proper scoring rule | Brier or log-loss |
| 14 | Fitting the recalibration map on training predictions | Learns in-sample over-confidence | Fit on held-out scores |
| 15 | Reporting correlation as a regression metric | Affine-invariant; \(\hat y=100+0.01y\) scores 1 | RMSE, MAE, held-out \(R^2\), calibration slope |
| 16 | One holdout split as “the” estimate | Unbiased but very imprecise | Resample; report an interval |
| 17 | Choosing \(k\) without the tradeoff | Small \(k\) pessimistic; large \(k\) correlated | \(k=5\) or \(10\); repeat to stabilize |
| 18 | Unstratified folds under imbalance | A fold may contain no positives | strata = in every split |
| 19 | train.err in .632 taken from held-out data |
Both terms then estimate the same thing | The first term is the apparent error |
| 20 | OOB set built from 1:length(boot_indices) |
Tail rows never out-of-bag | setdiff(seq_len(n), idx) |
| 21 | Reporting the best of several CV scores | Maximum of noisy estimates is biased up | Nested CV |
| 22 | Random folds on grouped or temporal data | The model memorizes subject or interpolates time | group_vfold_cv; rolling origin |
| 23 | Naive paired \(t\)-test on CV folds | Folds are correlated; test is anti-conservative | Corrected resampled \(t\); adjust for multiplicity |
| 24 | Trusting CV to catch leakage | It estimates the right answer to the wrong question | Ask whether the feature exists at prediction time |
29 Practice problems
29.1 Problem 1: Verify the optimism identity for a nonlinear model
The identity \(\mathbb{E}[\mathrm{op}]=\frac2n\sum_i\operatorname{Cov}(\hat y_i,y_i)\) holds for any fitting procedure, not just linear ones. Verify it for a regression tree and compare its effective degrees of freedom to its leaf count.
Solution
set.seed(201)
n_p1 <- 100
X_p1 <- matrix(rnorm(n_p1 * 3), n_p1, 3)
mu_p1 <- X_p1[, 1] + 0.5 * X_p1[, 2]^2
sigma_p1 <- 1
fit_and_predict <- function(depth) {
reps <- 3000
yh <- matrix(NA_real_, n_p1, reps); ys <- matrix(NA_real_, n_p1, reps)
eb <- ein <- numeric(reps)
for (r in seq_len(reps)) {
y <- mu_p1 + rnorm(n_p1, sd = sigma_p1)
d <- data.frame(X_p1, y = y)
f <- rpart::rpart(y ~ ., d,
control = rpart::rpart.control(maxdepth = depth, cp = 0, minsplit = 5))
p <- predict(f, d)
yn <- mu_p1 + rnorm(n_p1, sd = sigma_p1)
yh[, r] <- p; ys[, r] <- y
eb[r] <- mean((y - p)^2); ein[r] <- mean((yn - p)^2)
}
cov_sum <- sum(vapply(seq_len(n_p1), \(i) cov(yh[i, ], ys[i, ]), numeric(1)))
c(depth = depth,
observed_optimism = mean(ein - eb),
two_over_n_times_cov = 2 * cov_sum / n_p1,
effective_df = cov_sum / sigma_p1^2,
leaves = sum(rpart::rpart(y ~ ., data.frame(X_p1, y = ys[, 1]),
control = rpart::rpart.control(maxdepth = depth, cp = 0,
minsplit = 5))$frame$var == "<leaf>"))
}
as.data.frame(do.call(rbind, lapply(c(1, 2, 4), fit_and_predict))) |>
mutate(across(where(is.numeric), \(z) round(z, 3)))#> depth observed_optimism two_over_n_times_cov effective_df leaves
#> 1 1 0.218 0.222 11.105 2
#> 2 2 0.393 0.402 20.076 4
#> 3 4 0.737 0.738 36.894 1229.2 Problem 2: Show that AUC cannot see miscalibration
Construct several monotone distortions of one score vector and confirm that AUC is invariant while every proper scoring rule is not.
Solution
set.seed(203)
n2 <- 3000
p0 <- plogis(rnorm(n2))
y2 <- rbinom(n2, 1, p0)
distortions <- list(
identity = \(p) p,
`p^3` = \(p) p^3,
`sqrt(p)` = \(p) sqrt(p),
`logit x 4`= \(p) plogis(4 * qlogis(p)),
`logit + 2`= \(p) plogis(qlogis(p) + 2))
do.call(rbind, lapply(names(distortions), function(nm) {
q <- pmin(pmax(distortions[[nm]](p0), 1e-9), 1 - 1e-9)
data.frame(distortion = nm,
AUC = round(as.numeric(pROC::auc(pROC::roc(y2, q, quiet = TRUE))), 6),
Brier = round(mean((q - y2)^2), 5),
log_loss = round(-mean(y2 * log(q) + (1 - y2) * log(1 - q)), 5),
calib_slope = round(unname(coef(glm(y2 ~ qlogis(q),
family = binomial()))[2]), 4))
}))#> distortion AUC Brier log_loss calib_slope
#> 1 identity 0.731234 0.20904 0.60374 0.9912
#> 2 p^3 0.731234 0.31251 0.98022 0.5606
#> 3 sqrt(p) 0.731234 0.24696 0.69183 1.1682
#> 4 logit x 4 0.731234 0.25956 0.92902 0.2478
#> 5 logit + 2 0.731234 0.33730 0.97787 0.991229.3 Problem 3: The kappa paradoxes
Construct two confusion matrices with the same \(\kappa\) and very different accuracy, and two with the same accuracy and very different \(\kappa\).
Solution
kstats <- function(TP, FP, FN, TN) {
n <- TP + FP + FN + TN
p_o <- (TP + TN) / n
p_e <- ((TP + FP) * (TP + FN) + (TN + FN) * (TN + FP)) / n^2
c(accuracy = round(p_o, 4), kappa = round((p_o - p_e) / (1 - p_e), 4),
prevalence = round((TP + FN) / n, 4))
}
rbind(
`same kappa, different accuracy (A)` = kstats(45, 15, 15, 25),
`same kappa, different accuracy (B)` = kstats(10, 8, 8, 174),
`same accuracy, different kappa (C)` = kstats(40, 10, 10, 40),
`same accuracy, different kappa (D)` = kstats(76, 4, 16, 4))#> accuracy kappa prevalence
#> same kappa, different accuracy (A) 0.70 0.3750 0.60
#> same kappa, different accuracy (B) 0.92 0.5116 0.09
#> same accuracy, different kappa (C) 0.80 0.6000 0.50
#> same accuracy, different kappa (D) 0.80 0.1935 0.9229.4 Problem 4, Break the .632 estimator, then fix it
Reproduce the .632 estimate using a held-out error as the first term, and compare against the correct construction.
Solution
set.seed(205)
n4 <- 250
d4 <- data.frame(matrix(rnorm(n4 * 5), n4, 5))
d4$y <- factor(ifelse(d4$X1 + d4$X2 + rnorm(n4) > 0, "a", "b"))
tr4 <- 1:180; te4 <- 181:n4
fit4 <- fit_tree(y ~ ., d4[tr4, ])
apparent <- loss_01(d4[tr4, ], pred_tree(fit4, d4[tr4, ]))
heldout <- loss_01(d4[te4, ], pred_tree(fit4, d4[te4, ]))
# Correct OOB, built with setdiff
set.seed(207)
oob_losses <- matrix(NA_real_, length(tr4), 60)
for (b in 1:60) {
idx <- sample(tr4, length(tr4), replace = TRUE)
oob <- setdiff(tr4, idx)
if (!length(oob)) next
fb <- fit_tree(y ~ ., d4[idx, ])
oob_losses[match(oob, tr4), b] <- loss_01(d4[oob, ], pred_tree(fb, d4[oob, ]),
per_case = TRUE)
}
err_oob <- mean(rowMeans(oob_losses, na.rm = TRUE), na.rm = TRUE)
# WRONG OOB: indexing only into 1:length(idx)
set.seed(207)
wrong_oob <- mean(vapply(1:60, function(b) {
idx <- sample(length(tr4), floor(0.9 * length(tr4)), replace = TRUE)
i <- seq_along(idx); i <- i[is.na(match(i, idx))] # never exceeds 0.9n
if (!length(i)) return(NA_real_)
fb <- fit_tree(y ~ ., d4[tr4[idx], ])
loss_01(d4[tr4[i], ], pred_tree(fb, d4[tr4[i], ]))
}, numeric(1)), na.rm = TRUE)
c(apparent_error = round(apparent, 4),
true_heldout_error = round(heldout, 4),
correct_632 = round(0.368 * apparent + 0.632 * err_oob, 4),
using_heldout_as_first_term = round(0.368 * heldout + 0.632 * err_oob, 4),
using_wrong_oob_index = round(0.368 * apparent + 0.632 * wrong_oob, 4))#> apparent_error true_heldout_error
#> 0.0000 0.4286
#> correct_632 using_heldout_as_first_term
#> 0.2162 0.3739
#> using_wrong_oob_index
#> 0.210729.5 Problem 5, Measure the anti-conservatism of the naive test
Estimate the false-positive rate of the paired \(t\)-test on CV folds under a true null, as a function of \(k\).
Solution
set.seed(209)
null_rate_k <- function(k, reps = 300, n = 200) {
rej <- vapply(seq_len(reps), function(r) {
X <- matrix(rnorm(n * 4), n, 4)
y <- as.vector(X %*% c(1, 0.5, 0, 0)) + rnorm(n)
d <- data.frame(X, y = y)
f <- sample(rep(seq_len(k), length.out = n))
dif <- vapply(seq_len(k), function(i) {
tr <- d[f != i, ]; te <- d[f == i, ]
# Two models that differ only trivially (an irrelevant extra feature)
a <- mean((predict(lm(y ~ X1 + X2, tr), te) - te$y)^2)
b <- mean((predict(lm(y ~ X1 + X2 + X3, tr), te) - te$y)^2)
a - b
}, numeric(1))
c(naive = abs(mean(dif) / (sd(dif) / sqrt(k))) > qt(0.975, k - 1),
corrected = abs(mean(dif) / sqrt((1/k + 1/(k-1)) * var(dif))) > qt(0.975, k - 1))
}, logical(2))
c(k = k, naive = mean(rej["naive", ]), corrected = mean(rej["corrected", ]))
}
nr <- as.data.frame(do.call(rbind, lapply(c(5, 10, 20), null_rate_k)))
round(nr, 4)#> k naive corrected
#> 1 5 0.0700 0.0067
#> 2 10 0.1633 0.0267
#> 3 20 0.2400 0.0967nr |> pivot_longer(-k, names_to = "test", values_to = "rate") |>
ggplot(aes(factor(k), rate, fill = test)) +
geom_col(position = "dodge") +
geom_hline(yintercept = 0.05, linetype = "dashed", color = "black") +
scale_fill_manual(values = c(naive = "#D8433B", corrected = "#3B7DD8")) +
labs(title = "False-positive rate under a true null",
subtitle = "Dashed line: the nominal 5% level",
x = "Folds k", y = "Rejection rate", fill = NULL) +
theme_dspa()
29.6 Problem 6: Reproduce the LOOCV shortcut for ridge
Verify the hat-value formula for ridge regression and use it to select \(\lambda\) in one pass.
Solution
set.seed(211)
n6 <- 150; p6 <- 25
X6 <- scale(matrix(rnorm(n6 * p6), n6, p6))
y6 <- as.vector(X6 %*% c(rep(1.5, 5), rep(0, p6 - 5))) + rnorm(n6, sd = 2)
ridge_loocv <- function(lambda) {
S <- X6 %*% solve(crossprod(X6) + lambda * diag(p6)) %*% t(X6)
yh <- S %*% y6
h <- diag(S)
c(loocv = mean(((y6 - yh) / (1 - h))^2),
gcv = mean(((y6 - yh) / (1 - sum(h) / n6))^2),
edf = sum(h))
}
lams6 <- 10^seq(-2, 4, length.out = 40)
r6 <- as.data.frame(t(vapply(lams6, ridge_loocv, numeric(3))))
r6$lambda <- lams6
# Brute-force check at one lambda
lam_chk <- 10
brute <- mean(vapply(seq_len(n6), function(i) {
b <- solve(crossprod(X6[-i, ]) + lam_chk * diag(p6), crossprod(X6[-i, ], y6[-i]))
(y6[i] - X6[i, ] %*% b)^2
}, numeric(1)))
c(closed_form = round(ridge_loocv(lam_chk)[["loocv"]], 6),
brute_force = round(brute, 6),
lambda_min_loocv = round(lams6[which.min(r6$loocv)], 3),
lambda_min_gcv = round(lams6[which.min(r6$gcv)], 3))#> closed_form brute_force lambda_min_loocv lambda_min_gcv
#> 4.48745 4.48745 11.93800 11.93800r6 |> select(lambda, LOOCV = loocv, GCV = gcv) |>
pivot_longer(-lambda, names_to = "criterion", values_to = "value") |>
ggplot(aes(lambda, value, color = criterion)) +
geom_line(linewidth = 1) + scale_x_log10() +
scale_color_manual(values = c(LOOCV = "#3B7DD8", GCV = "#D8433B")) +
labs(title = "LOOCV and GCV for ridge regression, from one fit each",
x = expression(lambda~"(log scale)"), y = "Criterion", color = NULL) +
theme_dspa()
29.7 Problem 7, Grouped versus random folds
Quantify the optimism of random folds as a function of the intraclass correlation.
Solution
set.seed(213)
grouped_gap <- function(icc, n_sub = 40, per = 5) {
se <- rnorm(n_sub, sd = sqrt(icc))
res <- rnorm(n_sub * per, sd = sqrt(1 - icc))
subj <- rep(seq_len(n_sub), each = per)
y <- factor(ifelse(se[subj] + res > 0, "a", "b"))
Xg <- cbind(se[subj] + rnorm(n_sub * per, sd = 0.4),
matrix(rnorm(n_sub * per * 2), ncol = 2))
d <- data.frame(Xg, subject = subj, y = y)
acc <- function(folds) mean(vapply(folds, function(te) {
tr <- d[-te, ]
mean(class::knn(tr[, 1:3], d[te, 1:3], tr$y, k = 5) == d$y[te])
}, numeric(1)))
fr <- split(sample(nrow(d)), rep(1:5, length.out = nrow(d)))
gs <- split(seq_len(n_sub), rep(1:5, length.out = n_sub))
fg <- lapply(gs, \(s) which(d$subject %in% s))
c(icc = icc, random = acc(fr), grouped = acc(fg))
}
gg <- as.data.frame(do.call(rbind, lapply(c(0.05, 0.2, 0.5, 0.8), grouped_gap)))
gg$optimism <- gg$random - gg$grouped
round(gg, 4)#> icc random grouped optimism
#> 1 0.05 0.535 0.455 0.080
#> 2 0.20 0.545 0.565 -0.020
#> 3 0.50 0.635 0.630 0.005
#> 4 0.80 0.710 0.680 0.030gg |> select(icc, Random = random, Grouped = grouped) |>
pivot_longer(-icc, names_to = "folds", values_to = "accuracy") |>
ggplot(aes(icc, accuracy, color = folds)) +
geom_line(linewidth = 1) + geom_point(size = 2.4) +
scale_color_manual(values = c(Random = "#D8433B", Grouped = "#3B7DD8")) +
labs(title = "Optimism of random folds grows with the intraclass correlation",
subtitle = "Grouped folds keep every subject's observations on one side of the split",
x = "Intraclass correlation", y = "CV accuracy", color = NULL) +
theme_dspa()
29.8 Problem 8: Does recalibration change the decision?
Compare net benefit before and after recalibration, and explain why the decision curve moves when the ROC curve does not.
Solution
set.seed(215)
n8 <- 4000
x8 <- rnorm(n8)
pt8 <- plogis(-2.0 + 1.3 * x8)
y8 <- rbinom(n8, 1, pt8)
p_bad <- plogis(3 * qlogis(pt8)) # over-confident
cal8 <- sample(n8, n8 / 2); ev8 <- setdiff(seq_len(n8), cal8)
fit8 <- glm(y8[cal8] ~ qlogis(p_bad[cal8]), family = binomial())
p_fix <- plogis(coef(fit8)[1] + coef(fit8)[2] * qlogis(p_bad[ev8]))
th8 <- seq(0.02, 0.45, by = 0.005)
pi8 <- mean(y8[ev8])
c(AUC_before = round(as.numeric(pROC::auc(pROC::roc(y8[ev8], p_bad[ev8],
quiet = TRUE))), 4),
AUC_after = round(as.numeric(pROC::auc(pROC::roc(y8[ev8], p_fix,
quiet = TRUE))), 4),
Brier_before = round(mean((p_bad[ev8] - y8[ev8])^2), 4),
Brier_after = round(mean((p_fix - y8[ev8])^2), 4))#> AUC_before AUC_after Brier_before Brier_after
#> 0.7923 0.7923 0.1221 0.1094bind_rows(
data.frame(t = th8, nb = net_benefit(p_bad[ev8], y8[ev8], th8),
model = "Uncalibrated"),
data.frame(t = th8, nb = net_benefit(p_fix, y8[ev8], th8),
model = "Recalibrated"),
data.frame(t = th8, nb = pi8 - (1 - pi8) * th8 / (1 - th8), model = "Treat all"),
data.frame(t = th8, nb = 0, model = "Treat none")) |>
ggplot(aes(t, nb, color = model)) +
geom_line(linewidth = 0.9) +
coord_cartesian(ylim = c(-0.03, NA)) +
scale_color_manual(values = c(Uncalibrated = "#D8433B",
Recalibrated = "#3B7DD8",
`Treat all` = "grey45", `Treat none` = "black")) +
labs(title = "Recalibration changes decisions even though it cannot change AUC",
subtitle = "Net benefit depends on WHERE a probability falls relative to the threshold, not just on its rank",
x = expression("Threshold probability "*p[t]), y = "Net benefit",
color = NULL) +
theme_dspa()
30 Checkpoint
- Your training RMSE is 0.8 and the irreducible noise is known to be \(\sigma=1\). What has happened?
- A screening test has 99% sensitivity and 99% specificity for a disease with prevalence 1 in 10,000. A patient tests positive. What do you tell them?
- Two models have identical AUC. Under what circumstances would you still prefer one over the other?
- Your 10-fold CV accuracy is 0.88 after tuning 20 hyperparameter configurations. Is 0.88 your performance estimate?
- You compute the .632 bootstrap and get an error lower than your resubstitution error. What went wrong?
- A colleague reports that model A beats model B with a paired \(t\)-test across 10 CV folds, \(p = 0.03\). What do you check?
Answers
- You are fitting noise. Training RMSE below the irreducible noise level is impossible for genuine predictive error, the best achievable expected test RMSE is \(\sigma=1\). The gap is optimism, which for a linear fit has expectation \(2d\sigma^2/n\): it grows with the number of effective parameters and shrinks as \(1/n\). Report a resampled estimate, or correct the training error with \(C_p\)/AIC using the effective degrees of freedom rather than a parameter count.
- That the test is far more likely wrong than right. Sensitivity and specificity condition on the true state; what the patient needs is the PPV, which conditions on the result and depends on prevalence: \[\mathrm{PPV}=\frac{0.99\times10^{-4}}{0.99\times10^{-4}+0.01\times(1-10^{-4})}\approx 0.0098 .\] Roughly 1% of positives are true positives, a false discovery rate near 99%, with the Type I error still fixed at exactly 1%. Sensitivity and specificity are properties of the test. PPV and NPV are properties of the test in a population.
- Whenever the probabilities themselves are used. AUC is invariant to any strictly increasing transformation of the scores, so it cannot distinguish a calibrated model from a badly over-confident one. If the output feeds a cost-based threshold, a risk communication, or a downstream calculation, compare Brier score, log-loss, calibration slope and intercept, and the decision curve. Also consider computational cost, interpretability, and under heavy imbalance, the PR curve rather than the ROC.
- No, it is optimistically biased. You selected the maximum over 20 noisy CV estimates and are reporting that maximum. The bias grows with the size of the search and shrinks with \(n\). Nested cross-validation separates the jobs: an inner loop selects the configuration within each training set, an outer loop scores the whole select-and-fit procedure on data neither loop has seen. Expect the honest number to be lower.
- The first term is not the apparent error. The
estimator is \(0.368\,\overline{\mathrm{err}}+0.632\,\widehat{\mathrm{Err}}^{(1)}\),
where \(\overline{\mathrm{err}}\) is
the resubstitution error and \(\widehat{\mathrm{Err}}^{(1)}\) the
out-of-bag error. Since OOB is always at least as large as
resubstitution, the blend must lie between them, it can never fall below
the resubstitution error. Getting a lower value means either a held-out
error was substituted for the apparent one, or the out-of-bag set was
constructed incorrectly (it must be
setdiff(seq_len(n), boot_indices)). - Whether the test accounted for the correlation between folds. CV training sets share \(n(k-2)/k\) observations, so the fold errors are positively correlated and the ordinary paired \(t\)-test is anti-conservative, it rejects well above its nominal rate under a true null. Use the corrected resampled \(t\)-test, which inflates the variance by \(\left(\frac1J+\frac{n_{\text{test}}}{n_{\text{train}}}\right)\). Then check whether this was one comparison among many (adjust for multiplicity) and whether the difference is practically meaningful, a statistically detectable 0.3-point gain in AUC rarely changes a decision.
31 Summary
What we are estimating
- Training error, conditional test error, and expected test error are three different quantities. Resampling estimates the third, the error of the procedure, not of the one model you hold.
- Expected optimism is \(\frac{2}{n}\sum_i\operatorname{Cov}(\hat y_i,y_i)\), equal to \(2d\sigma^2/n\) for a linear fit. It generates \(C_p\), AIC, and effective degrees of freedom.
- Effective df replaces parameter counting and is continuous in the regularization strength.
Metrics
- Every rate is a conditional probability. Column rates (sensitivity, specificity, \(\alpha\), \(\beta\)) condition on the truth; row rates (PPV, NPV, FDR, FOR) condition on the prediction, and depend on prevalence.
- \(\kappa\) ranges over \([-1,1]\), depends strongly on prevalence and marginal asymmetry, and is for nominal agreement; weighted \(\kappa\) requires an ordering.
- ROC AUC is prevalence-invariant with chance level 0.5; PR-AUC has chance level equal to the prevalence and is the more informative summary under imbalance.
- AUC measures ranking only. A monotone distortion leaves it unchanged while destroying every probability. Assess calibration with a reliability diagram, the calibration slope, and a proper scoring rule.
- The Brier decomposition separates reliability from resolution: two models can discriminate identically and calibrate very differently.
- Net benefit answers whether a model beats treating everyone or no one.
- For regression, report RMSE, MAE, held-out \(R^2\), and the calibration slope, never correlation alone.
Resampling
- A single holdout is unbiased and imprecise.
- \(k\)-fold CV trades pessimistic bias (small \(k\)) against correlated folds (large \(k\)); \(k=5\) or \(10\) is the compromise. Stratify for classification.
- LOOCV has a closed form for linear smoothers, costing one fit instead of \(n\). GCV replaces the leverages with their mean.
- The .632 estimator blends the apparent error with the out-of-bag error; .632+ corrects it for learners that interpolate.
- Nested CV is required whenever tuning and evaluation share data.
- Random folds assume exchangeability. Use grouped folds for repeated measures and rolling origin for time series.
Comparison and improvement
- CV folds are correlated, so the naive paired \(t\)-test is anti-conservative. Use the corrected resampled \(t\) and adjust for multiplicity.
- Learning curves separate bias-limited from variance-limited models and prescribe opposite remedies.
- Random search beats grid search for the same budget when only a few hyperparameters matter; the 1-SE rule trades an indistinguishable amount of performance for a simpler model.
- Recalibration is nearly free, leaves AUC untouched, and can change which cases are acted on.
- Cross-validation cannot detect leakage. It correctly estimates the answer to whatever question the features encode.
Where these threads continue
| Thread | Continues in |
|---|---|
| Regularization paths, LASSO, stability selection | Feature selection |
| Rolling origin and forecast evaluation | Longitudinal analysis |
| Hyperparameter search as optimization | Function optimization |
| Validation for representation learning | Deep learning |
32 Chapter roadmap
- Chapter 1:
Foundations. R toolchain, reproducibility conventions,
dspa_read(), simulation. - Chapter 2: Data quality and exploratory visual analytics. Missingness, robust statistics, multiplicity.
- Chapter 3: Linear algebra, matrix computing, and regression. Projection, hat matrix, conditioning.
- Chapter 4: Dimensionality reduction. PCA, truncated SVD, distance concentration.
- Chapter 5: Supervised classification. Bayes error, the metric family, the leakage taxonomy.
- Chapter 6: Black-box methods. Neural networks, kernels, ensembles, temporal validation.
- Chapter 7: Text mining and association rules. TF-IDF, similarity, multiplicity control.
- Chapter 8: Unsupervised clustering. Internal and external validation, stability, the adjusted Rand index.
- Variable importance and feature selection. Ridge, LASSO, elastic net, stability selection, FDR control.
- Longitudinal and time-series analysis. Mixed models, ARIMA, forecast evaluation.
- Function optimization. Gradient descent, duality, Bayesian optimization.
- Deep learning. Convolutional and recurrent networks, representation learning, transfer learning.
33 Session information
#> R version 4.3.3 (2024-02-29 ucrt)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows 11 x64 (build 26200)
#>
#> Matrix products: default
#>
#>
#> locale:
#> [1] LC_COLLATE=English_United States.utf8
#> [2] LC_CTYPE=English_United States.utf8
#> [3] LC_MONETARY=English_United States.utf8
#> [4] LC_NUMERIC=C
#> [5] LC_TIME=English_United States.utf8
#>
#> time zone: America/New_York
#> tzcode source: internal
#>
#> attached base packages:
#> [1] parallel stats graphics grDevices utils datasets methods
#> [8] base
#>
#> other attached packages:
#> [1] doParallel_1.0.17 iterators_1.0.14 foreach_1.5.2 rpart_4.1.23
#> [5] PRROC_1.4 rlang_1.1.5 pROC_1.18.5 caret_6.0-94
#> [9] lattice_0.22-6 yardstick_1.3.1 rsample_1.2.1 plotly_4.12.1
#> [13] patchwork_1.3.0 tidyr_1.3.1 dplyr_1.1.4 ggplot2_4.0.1
#>
#> loaded via a namespace (and not attached):
#> [1] tidyselect_1.2.1 viridisLite_0.4.2 timeDate_4032.109
#> [4] farver_2.1.2 S7_0.2.1 fastmap_1.2.0
#> [7] digest_0.6.37 timechange_0.3.0 lifecycle_1.0.5
#> [10] survival_3.7-0 kernlab_0.9-32 magrittr_2.0.3
#> [13] compiler_4.3.3 sass_0.4.9 tools_4.3.3
#> [16] yaml_2.3.10 data.table_1.16.4 knitr_1.51
#> [19] labeling_0.4.3 htmlwidgets_1.6.4 plyr_1.8.9
#> [22] RColorBrewer_1.1-3 withr_3.0.2 purrr_1.0.2
#> [25] nnet_7.3-19 grid_4.3.3 stats4_4.3.3
#> [28] e1071_1.7-14 future_1.33.2 globals_0.16.3
#> [31] scales_1.4.0 MASS_7.3-60.0.1 cli_3.6.3
#> [34] rmarkdown_2.31 generics_0.1.3 otel_0.2.0
#> [37] rstudioapi_0.18.0 future.apply_1.11.2 httr_1.4.7
#> [40] reshape2_1.4.4 cachem_1.1.0 proxy_0.4-27
#> [43] stringr_1.5.1 splines_4.3.3 vctrs_0.6.5
#> [46] hardhat_1.4.3 Matrix_1.6-5 jsonlite_1.8.9
#> [49] listenv_0.9.1 crosstalk_1.2.1 gower_1.0.1
#> [52] jquerylib_0.1.4 recipes_1.4.0 glue_1.8.0
#> [55] parallelly_1.37.1 codetools_0.2-20 lubridate_1.9.3
#> [58] stringi_1.8.4 gtable_0.3.6 tibble_3.2.1
#> [61] furrr_0.3.1 pillar_1.10.1 htmltools_0.5.8.1
#> [64] ipred_0.9-14 gbm_2.3.1 lava_1.8.0
#> [67] R6_2.6.1 evaluate_1.0.3 bslib_0.9.0
#> [70] class_7.3-22 Rcpp_1.0.14 nlme_3.1-165
#> [73] prodlim_2024.06.25 ranger_0.16.0 xfun_0.52
#> [76] pkgconfig_2.0.3 ModelMetrics_1.2.2.2