| SOCR ≫ | DSPA ≫ | DSPA3 Topics ≫ |
DSPA3: Data Science and Predictive Analytics (UMich HS650)
Chapter 5: Supervised Classification
Chapter 5: Supervised Classification
SOCR/MCAIM (Ivo Dinov)
SOCR/MCAIM (Ivo Dinov)
August 26, 2026
library(ggplot2)
library(dplyr)
library(tidyr)
library(patchwork)
library(plotly) # interactive figures and ALL 3-D graphics
library(class) # knn
library(caret) # confusionMatrix, train, resampling
library(pROC) # ROC / AUCHow 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, printed in these notes, ready to paste into a live session.Every three-dimensional figure is drawn with
plot_ly()and evaluated. That matters unusually much here: a classifier is a posterior probability surface over feature space, and the difference between \(k=1\) and \(k=100\), or between a single tree and a forest, is a difference in the roughness of that surface. A fixed projection hides exactly the property being taught.
1 Learning objectives
After completing this chapter you will be able to:
- State the Bayes classifier and the Bayes error, and explain why every method in this chapter is an attempt to approximate them.
- Compute and interpret the full family of classification metrics — sensitivity, specificity, PPV, NPV, \(F_1\), \(\kappa\), MCC, AUC, Brier score — and choose the right one for a given cost structure.
- Compare any reported accuracy against the no-information rate before concluding that a model has learned anything.
- Identify the four kinds of data leakage and design a resampling protocol that avoids them.
- Derive the Cover–Hart bound \(R_{1\text{NN}}\le 2R_{\text{Bayes}}\) and explain what it says about simple methods.
- Analyze the bias–variance tradeoff governing \(k\) in kNN, and state how the optimal \(k\) scales with \(n\) and \(d\).
- State the naive Bayes conditional-independence assumption, explain why the classifier works despite its violation, and implement it in log space to avoid underflow.
- Derive Laplace smoothing as the posterior mean of a Dirichlet–multinomial model, and tune it by held-out log-loss.
- Compare entropy, Gini, and misclassification error as splitting criteria, and explain why the last is used for pruning but not for splitting.
- Explain how bagging reduces variance, why decorrelation matters, and what out-of-bag error estimates for free.
Estimated time: 10–13 hours including exercises. Prerequisites: Chapter 2 (contingency tables, class imbalance), Chapter 3 (trees, ensembles, complexity), and Chapter 4, in particular distance concentration (§4.3.1), which is the reason kNN degrades in high dimensions.
2 PART I: THE CLASSIFICATION PROBLEM
3 Supervised and unsupervised learning
Machine learning divides along whether the algorithm is shown the answers.
Supervised learning provides input variables \(X\) and an outcome \(Y\), and seeks a mapping \(\hat f: \mathcal{X}\to\mathcal{Y}\) that generalizes to unseen data. It is “supervised” because training labels correct the learner. When \(Y\) is quantitative the task is regression (Chapter 3); when \(Y\) is categorical it is classification, the subject of this chapter.
Unsupervised learning provides only \(X\) and asks the algorithm to find structure, groupings, associations, low-dimensional coordinates, with no correct answer to compare against. Clustering and association rule mining are the canonical problems (Chapter 8, Chapter 7).
The distinction is about supervision, not about the algorithms themselves. Gaussian mixtures cluster when labels are absent and classify when they are present; the same tree machinery does regression and classification.
| Outcome type | Supervised methods | Unsupervised counterparts |
|---|---|---|
| Binary | kNN, naive Bayes, logistic regression, LDA/QDA, decision trees, SVM, boosting | k-means, Gaussian mixtures, spectral clustering |
| Multi-class | kNN, naive Bayes, multinomial logistic, trees, random forest, neural networks | Hierarchical clustering, DBSCAN, HDBSCAN |
| Ordinal | Ordinal logistic, ordinal forests, cumulative-link models | Ordinal MDS |
| Quantitative | Linear regression, regression trees, random forest, BART, SVR | PCA, factor analysis, autoencoders |
| None (structure) | – | Apriori / association rules, t-SNE, UMAP |
4 What are we actually trying to compute?
Before any algorithm, there is a target. Without it, “this classifier is 72% accurate” has no reference point.
4.1 The Bayes classifier
Suppose \((X, Y)\) are drawn from a joint distribution with \(Y\in\{1,\dots,K\}\). Under 0–1 loss, every error costs the same, the expected loss of a classifier \(g\) is its misclassification rate
\[R(g)=P\big(g(X)\ne Y\big)=E_X\Big[1-P\big(Y=g(X)\mid X\big)\Big].\]
Minimizing pointwise gives the Bayes classifier
\[\boxed{\;g^{*}(x)=\arg\max_{k\in\{1,\dots,K\}} P(Y=k\mid X=x)\;}\]
and its risk is the Bayes error
\[R^{*}=E_X\Big[1-\max_k P(Y=k\mid X=x)\Big].\]
No classifier can do better than \(R^{*}\). The Bayes error is irreducible: it is the overlap between the classes in the true joint distribution, not a deficiency of any method. If two patients have identical measured features and different outcomes, no algorithm can separate them.
Every method in this chapter is a plug-in rule: estimate \(\hat P(Y=k\mid X=x)\) somehow, then take the argmax. They differ only in how they estimate the posterior.
| Method | How it estimates \(P(Y=k\mid x)\) |
|---|---|
| kNN | Local vote fraction among the \(k\) nearest training points |
| Naive Bayes | \(\propto \hat P(k)\prod_j \hat P(x_j\mid k)\), assuming conditional independence |
| LDA / QDA | Gaussian class-conditionals with shared / separate covariance |
| Logistic regression | Direct parametric model of the log-odds |
| Decision tree | Class proportions in the leaf containing \(x\) |
| Random forest | Average of many trees’ leaf proportions |
4.2 Seeing the Bayes boundary
The Bayes error is normally unknowable, because the joint distribution is unknown. On simulated data we know it exactly, which makes simulation the only setting where “how close to optimal is this classifier?” has an answer.
# Two Gaussian classes with a KNOWN generating process
mu1 <- c(-1, -0.6); mu2 <- c(1.1, 0.7)
S1 <- matrix(c(1.0, 0.35, 0.35, 0.9), 2)
S2 <- matrix(c(1.2, -0.30, -0.30, 0.7), 2)
pi1 <- 0.45
# Posterior P(Y = 2 | x) from the true densities
post2 <- function(x1, x2) {
z <- cbind(x1, x2)
d1 <- mvtnorm::dmvnorm(z, mu1, S1)
d2 <- mvtnorm::dmvnorm(z, mu2, S2)
(1 - pi1) * d2 / (pi1 * d1 + (1 - pi1) * d2)
}
set.seed(21)
n_sim <- 1200
lab <- rbinom(n_sim, 1, 1 - pi1) # 0 -> class 1, 1 -> class 2
Xsim <- rbind(MASS::mvrnorm(sum(lab == 0), mu1, S1),
MASS::mvrnorm(sum(lab == 1), mu2, S2))
ysim <- factor(c(rep("A", sum(lab == 0)), rep("B", sum(lab == 1))))
sim <- data.frame(x1 = Xsim[, 1], x2 = Xsim[, 2], y = ysim)
# Monte-Carlo estimate of the Bayes error
set.seed(22)
mc <- 2e5
lab_mc <- rbinom(mc, 1, 1 - pi1)
Xmc <- rbind(MASS::mvrnorm(sum(lab_mc == 0), mu1, S1),
MASS::mvrnorm(sum(lab_mc == 1), mu2, S2))
p_mc <- post2(Xmc[, 1], Xmc[, 2])
bayes_error <- mean(1 - pmax(p_mc, 1 - p_mc))
c(bayes_error = round(bayes_error, 4))#> bayes_error
#> 0.1004g1 <- seq(-4.5, 4.5, length.out = 220)
g2 <- seq(-4, 4, length.out = 220)
grid_df <- expand.grid(x1 = g1, x2 = g2)
grid_df$p <- post2(grid_df$x1, grid_df$x2)
ggplot(grid_df, aes(x1, x2)) +
geom_raster(aes(fill = p)) +
geom_contour(aes(z = p), breaks = 0.5, colour = "black", linewidth = 1) +
geom_point(data = sim, aes(x1, x2, shape = y), size = 1.1, alpha = 0.55,
colour = "grey15") +
scale_fill_gradient2(low = "#3B7DD8", mid = "white", high = "#D8433B",
midpoint = 0.5, name = "P(Y = B | x)") +
scale_shape_manual(values = c(A = 1, B = 3), name = "Class") +
coord_fixed() +
labs(title = "The Bayes-optimal decision boundary",
subtitle = sprintf("Black contour is P = 0.5. Bayes error = %.3f -- no classifier can beat this",
bayes_error),
x = expression(x[1]), y = expression(x[2])) +
theme_dspa()
Because the posterior is a surface over feature space, it is worth rotating:
Zp <- matrix(grid_df$p, nrow = length(g1))
plot_ly(x = g2, y = g1, z = Zp, type = "surface", colorscale = "RdBu",
reversescale = TRUE, opacity = 0.92,
colorbar = list(title = "P(Y=B|x)")) |>
add_trace(x = sim$x2, y = sim$x1,
z = post2(sim$x1, sim$x2),
type = "scatter3d", mode = "markers", name = "Observations",
marker = list(size = 1.8, opacity = 0.45, color = "black")) |>
layout(title = "True posterior probability surface",
scene = list(xaxis = list(title = "x2"), yaxis = list(title = "x1"),
zaxis = list(title = "P(Y = B | x)", range = c(0, 1))))The quadratic boundary comes from the two classes having different covariance matrices. Had they shared one, the boundary would be a straight line — which is exactly the LDA-versus-QDA distinction in §5.16. Keep this surface in mind: every fitted classifier below produces its own version, and the question is always how well it approximates this one.
5 Evaluating a classifier
5.1 The confusion matrix
Fix a positive class. For binary predictions:
| Reference: Positive | Reference: Negative | |
|---|---|---|
| Predicted: Positive | True Positive (TP) | False Positive (FP) |
| Predicted: Negative | False Negative (FN) | True Negative (TN) |
Orientation matters, and it is easy to get backwards.
caret’s signature isconfusionMatrix(data, reference), predictions first, truth second. Passing them in the other order transposes the matrix, which exchanges sensitivity with specificity and PPV with NPV. The printed output still looks entirely plausible, which is what makes the error dangerous. Always passpositive =explicitly too, orcaretchooses the first factor level and your “sensitivity” may refer to the class you think of as negative.
5.2 The metric family
\[ \begin{aligned} \text{Accuracy} &= \frac{TP+TN}{TP+FP+FN+TN} & \text{Sensitivity (Recall, TPR)} &= \frac{TP}{TP+FN}\\[2mm] \text{Specificity (TNR)} &= \frac{TN}{TN+FP} & \text{Precision (PPV)} &= \frac{TP}{TP+FP}\\[2mm] \text{NPV} &= \frac{TN}{TN+FN} & \text{Balanced accuracy} &= \tfrac12\big(\text{Sens}+\text{Spec}\big)\\[2mm] F_1 &= \frac{2\,\text{Prec}\cdot\text{Rec}}{\text{Prec}+\text{Rec}} & \text{Prevalence} &= \frac{TP+FN}{n} \end{aligned} \]
Two more deserve their own display. Cohen’s \(\kappa\) corrects accuracy for agreement expected by chance:
\[\kappa=\frac{p_o-p_e}{1-p_e},\qquad p_o=\text{Accuracy},\quad p_e=\sum_k \hat P(\text{pred}=k)\,\hat P(\text{ref}=k).\]
Matthews correlation coefficient (MCC) is the Pearson correlation between the binary prediction and truth vectors, and is the most informative single number under class imbalance:
\[\mathrm{MCC}=\frac{TP\cdot TN - FP\cdot FN}{\sqrt{(TP+FP)(TP+FN)(TN+FP)(TN+FN)}}\ \in[-1,1].\]
binary_metrics <- function(pred, truth, positive) {
pred <- as.character(pred); truth <- as.character(truth)
TP <- sum(pred == positive & truth == positive)
FP <- sum(pred == positive & truth != positive)
FN <- sum(pred != positive & truth == positive)
TN <- sum(pred != positive & truth != positive)
n <- TP + FP + FN + TN
sens <- TP / (TP + FN); spec <- TN / (TN + FP)
ppv <- TP / (TP + FP); npv <- TN / (TN + FN)
acc <- (TP + TN) / n
p_e <- ((TP + FP) * (TP + FN) + (TN + FN) * (TN + FP)) / n^2
mcc_den <- sqrt(prod(c(TP + FP, TP + FN, TN + FP, TN + FN)))
c(n = n, TP = TP, FP = FP, FN = FN, TN = TN,
accuracy = acc,
NIR = max(mean(truth == positive), 1 - mean(truth == positive)),
sensitivity = sens, specificity = spec, PPV = ppv, NPV = npv,
balanced_accuracy = (sens + spec) / 2,
F1 = 2 * ppv * sens / (ppv + sens),
kappa = (acc - p_e) / (1 - p_e),
MCC = if (mcc_den > 0) (TP * TN - FP * FN) / mcc_den else NA_real_)
}5.3 Accuracy is not enough
Common misconception: “72% accuracy means the model learned something.” Compare it against the no-information rate, the accuracy of always predicting the majority class. If 72% of cases are positive, the NIR is 0.72 and a model scoring 0.72 has learned nothing at all.
caretprints the NIR and a one-sided binomial test of Accuracy > NIR for exactly this reason.Under imbalance, report balanced accuracy, \(\kappa\), MCC, or AUC, all of which are near zero for a majority-class predictor, and always attach a binomial confidence interval, because a 50-case test set moves 2 percentage points per case.
set.seed(31)
truth_imb <- factor(sample(c("Pos", "Neg"), 500, TRUE, prob = c(0.8, 0.2)))
always_pos <- factor(rep("Pos", 500), levels = levels(truth_imb))
round(binary_metrics(always_pos, truth_imb, positive = "Pos"), 4)#> n TP FP FN
#> 500.0000 399.0000 101.0000 0.0000
#> TN accuracy NIR sensitivity
#> 0.0000 0.7980 0.7980 1.0000
#> specificity PPV NPV balanced_accuracy
#> 0.0000 0.7980 NaN 0.5000
#> F1 kappa MCC
#> 0.8877 0.0000 NAA classifier with 80% accuracy, perfect sensitivity, and \(\kappa = 0\), \(\mathrm{MCC}\) undefined, balanced accuracy exactly 0.5. It is the majority-class rule wearing a lab coat.
#> [1] 0.52675.4 Threshold-free evaluation: ROC and precision–recall
A plug-in classifier outputs \(\hat p(x)=\hat P(Y=\text{positive}\mid x)\) and then thresholds it. Every threshold \(t\) gives a different (sensitivity, specificity) pair; the ROC curve traces them all.
\[\mathrm{AUC}=\int_0^1 \mathrm{TPR}\big(\mathrm{FPR}^{-1}(u)\big)\,du = P\big(\hat p(X^{+}) > \hat p(X^{-})\big),\]
the probability that a randomly chosen positive scores higher than a randomly chosen negative. AUC \(=0.5\) is chance; \(1\) is perfect separation. It is threshold-free and prevalence-invariant, which makes it the right summary when the operating point is not yet chosen.
The precision–recall curve is preferable when positives are rare and false positives are the expensive error, because PR does not reward the large TN count that inflates specificity.
set.seed(37)
p_true <- post2(sim$x1, sim$x2)
score_good <- p_true + rnorm(nrow(sim), sd = 0.10)
score_weak <- p_true + rnorm(nrow(sim), sd = 0.55)
roc_g <- pROC::roc(sim$y, score_good, levels = c("A", "B"), direction = "<", quiet = TRUE)
roc_w <- pROC::roc(sim$y, score_weak, levels = c("A", "B"), direction = "<", quiet = TRUE)
roc_df <- bind_rows(
data.frame(fpr = 1 - roc_g$specificities, tpr = roc_g$sensitivities,
model = sprintf("Strong score (AUC = %.3f)", as.numeric(pROC::auc(roc_g)))),
data.frame(fpr = 1 - roc_w$specificities, tpr = roc_w$sensitivities,
model = sprintf("Weak score (AUC = %.3f)", as.numeric(pROC::auc(roc_w)))))
ggplot(roc_df, aes(fpr, tpr, colour = model)) +
geom_abline(slope = 1, intercept = 0, linetype = "dashed", colour = "grey50") +
geom_line(linewidth = 1) +
scale_colour_manual(values = c("#D8433B", "#3B7DD8")) +
coord_fixed() +
labs(title = "ROC curves",
subtitle = "Dashed diagonal is chance. AUC is the probability a random positive outscores a random negative",
x = "False positive rate (1 - specificity)",
y = "True positive rate (sensitivity)", colour = NULL) +
theme_dspa()
# --- Interactive equivalent ------------------------------------------------
plot_ly() |>
add_lines(x = 1 - roc_g$specificities, y = roc_g$sensitivities,
name = sprintf("Strong (AUC=%.3f)", as.numeric(pROC::auc(roc_g)))) |>
add_lines(x = 1 - roc_w$specificities, y = roc_w$sensitivities,
name = sprintf("Weak (AUC=%.3f)", as.numeric(pROC::auc(roc_w)))) |>
add_lines(x = c(0, 1), y = c(0, 1), name = "Chance",
line = list(dash = "dash", color = "gray")) |>
layout(title = "ROC curves",
xaxis = list(title = "False positive rate"),
yaxis = list(title = "True positive rate"),
legend = list(orientation = "h"))5.5 Are the probabilities themselves any good?
AUC measures ranking. It says nothing about whether \(\hat p = 0.8\) means the event happens 80% of the time. A model can have AUC \(0.95\) and be badly miscalibrated, which matters enormously when the probability feeds a decision rule or a cost calculation.
The Brier score is the mean squared error of the probabilities, \(\mathrm{BS}=\frac1n\sum_i(\hat p_i-y_i)^2\), and a calibration curve plots observed frequency against predicted probability within bins.
calib <- function(p, y_pos, bins = 10) {
b <- cut(p, breaks = seq(0, 1, length.out = bins + 1), include.lowest = TRUE)
data.frame(p = p, y = y_pos, b = b) |>
summarise(predicted = mean(p), observed = mean(y), n = dplyr::n(), .by = b) |>
filter(n >= 5)
}
y_pos <- as.integer(sim$y == "B")
p_cal <- pmin(pmax(p_true + rnorm(nrow(sim), sd = 0.08), 0.001), 0.999)
p_overc <- pmin(pmax(plogis(3 * qlogis(p_true)), 0.001), 0.999) # over-confident
cal_df <- bind_rows(
mutate(calib(p_cal, y_pos), model = sprintf("Calibrated (Brier %.3f)", mean((p_cal - y_pos)^2))),
mutate(calib(p_overc, y_pos), model = sprintf("Over-confident (Brier %.3f)", mean((p_overc - y_pos)^2))))
ggplot(cal_df, aes(predicted, observed, colour = model)) +
geom_abline(slope = 1, intercept = 0, linetype = "dashed", colour = "grey50") +
geom_line(linewidth = 0.9) + geom_point(aes(size = n)) +
scale_colour_manual(values = c("#3B7DD8", "#D8433B")) +
scale_size_continuous(range = c(1.5, 5), guide = "none") +
coord_fixed(xlim = c(0, 1), ylim = c(0, 1)) +
labs(title = "Calibration: does a predicted 0.8 happen 80% of the time?",
subtitle = "Both models rank equally well; only one reports honest probabilities",
x = "Mean predicted probability", y = "Observed frequency", colour = NULL) +
theme_dspa()
c(AUC_calibrated = as.numeric(pROC::auc(pROC::roc(sim$y, p_cal, quiet = TRUE))),
AUC_overconfident = as.numeric(pROC::auc(pROC::roc(sim$y, p_overc, quiet = TRUE))),
Brier_calibrated = mean((p_cal - y_pos)^2),
Brier_overconfident = mean((p_overc - y_pos)^2))#> AUC_calibrated AUC_overconfident Brier_calibrated Brier_overconfident
#> 0.9563624 0.9591790 0.0762912 0.0833263Identical AUC, very different Brier scores. Ranking and calibration are separate properties; report both when the probabilities will be used, not just the labels.
6 Honest evaluation protocol
The metrics above are only as trustworthy as the data they are computed on.
6.1 Data leakage
Leakage is any flow of information from the evaluation set into the fitting process. It inflates estimated performance and is invisible in the output — the model simply looks good. Four kinds recur:
| Kind | What happens | Typical symptom |
|---|---|---|
| Preprocessing leakage | Scaling, imputation, PCA, or feature selection computed on all data before splitting | Test performance close to CV performance, both optimistic |
| Target leakage | A feature encodes the outcome (a post-diagnosis code, a treatment field) | One feature dominates; accuracy near 1 |
| Temporal leakage | Training on data collected after the test period | Model fails on genuinely prospective data |
| Group leakage | Repeated measures from the same subject split across train and test | Performance collapses on new subjects |
Common misconception: “I standardized the data first, then split, that’s the same thing.” It is not.
scale()on the full matrix computes each column’s mean and standard deviation from training and test rows together, so every training feature carries a trace of the test set. The effect is small for large \(n\) and can be substantial for small \(n\), which is exactly when people are most tempted to skip the discipline. Compute the center and scale on the training rows and apply those constants to test; inside cross-validation, recompute them within every fold.
library(recipes); library(rsample)
set.seed(41)
# A small, high-dimensional problem where preprocessing leakage bites hard:
# select the 10 features most correlated with y -- with y being PURE NOISE
n_lk <- 120; p_lk <- 2000
Xlk <- matrix(rnorm(n_lk * p_lk), n_lk, p_lk)
ylk <- factor(sample(c("A", "B"), n_lk, TRUE))
# WRONG: screen features using ALL the data, then cross-validate
r_all <- abs(apply(Xlk, 2, \(v) cor(v, as.integer(ylk))))
keep_leaky <- order(r_all, decreasing = TRUE)[1:10]
df_leaky <- data.frame(Xlk[, keep_leaky], y = ylk)
cv_leaky <- train(y ~ ., data = df_leaky, method = "knn",
trControl = trainControl(method = "cv", number = 10),
tuneGrid = data.frame(k = 5))
# RIGHT: screen INSIDE each fold, so the held-out rows never inform selection
folds <- vfold_cv(data.frame(Xlk, y = ylk), v = 10, strata = y)
acc_clean <- vapply(folds$splits, function(s) {
tr <- analysis(s); te <- assessment(s)
rr <- abs(apply(as.matrix(tr[, 1:p_lk]), 2, \(v) cor(v, as.integer(tr$y))))
kp <- order(rr, decreasing = TRUE)[1:10]
m <- knn(train = tr[, kp], test = te[, kp], cl = tr$y, k = 5)
mean(m == te$y)
}, numeric(1))
c(`selection on ALL data (leaky)` = round(max(cv_leaky$results$Accuracy), 4),
`selection inside each fold` = round(mean(acc_clean), 4),
`truth (y is pure noise)` = 0.5)#> selection on ALL data (leaky) selection inside each fold
#> 0.7500 0.4998
#> truth (y is pure noise)
#> 0.5000The leaky protocol reports well above chance on data where no signal exists. Screening 2,000 noise features against the outcome guarantees that some correlate strongly by chance (Chapter 4, §4.13.3), and doing it before the split hands those spurious features to every fold.
6.2 Splitting and resampling
Stratified splitting preserves the class proportions in each part, always use it for classification, and especially under imbalance, where a random split can leave a fold with almost no minority cases.
\(k\)-fold cross-validation partitions into \(k\) folds, trains on \(k-1\), tests on the remaining one, and averages. The tradeoff: small \(k\) (e.g. 5) means each model sees less data (pessimistic bias) but the estimates are less correlated; large \(k\) (LOOCV, \(k=n\)) is nearly unbiased but has high variance and costs \(n\) fits. \(k=5\) or \(10\) is the standard compromise. Repeated CV averages over several partitions to reduce the variance of the estimate itself.
Nested CV is required when you both tune and evaluate. The inner loop selects hyperparameters; the outer loop estimates performance. Reporting the best inner-loop score as the performance estimate is a selection bias: you have optimized over the same data you are scoring on.
set.seed(43)
dat_n <- data.frame(sim[, c("x1", "x2")], y = sim$y)
# Naive: tune and report on the SAME resamples
naive_fit <- train(y ~ ., data = dat_n, method = "knn",
trControl = trainControl(method = "cv", number = 10),
tuneGrid = data.frame(k = seq(1, 61, by = 4)))
naive_estimate <- max(naive_fit$results$Accuracy)
# Nested: inner loop tunes, outer loop scores untouched data
outer <- vfold_cv(dat_n, v = 5, strata = y)
nested <- vapply(outer$splits, function(s) {
tr <- analysis(s); te <- assessment(s)
inner <- train(y ~ ., data = tr, method = "knn",
trControl = trainControl(method = "cv", number = 5),
tuneGrid = data.frame(k = seq(1, 61, by = 4)))
mean(predict(inner, te) == te$y)
}, numeric(1))
c(naive_best_inner_score = round(naive_estimate, 4),
nested_cv_estimate = round(mean(nested), 4),
optimism = round(naive_estimate - mean(nested), 4),
bayes_ceiling = round(1 - bayes_error, 4))#> naive_best_inner_score nested_cv_estimate optimism
#> 0.8975 0.8934 0.0042
#> bayes_ceiling
#> 0.8996The naive estimate is optimistic by exactly the amount the tuning search exploited noise. The gap is small here because there is one hyperparameter and plenty of data; with dozens of hyperparameters and small \(n\), it can exceed ten percentage points.
7 PART II: INSTANCE-BASED LEARNING
8 The k-nearest-neighbour classifier
kNN makes no distributional assumptions and builds no model. It memorizes the training set and defers all computation to prediction time, which is why it is called a lazy learner. The entire algorithm:
- Store the labelled training data.
- For a new case \(x\), compute \(d(x, x_i)\) to every training point.
- Let \(A_k(x)\) be the indices of the \(k\) nearest.
- Estimate the posterior by the local vote fraction, then take the argmax.
\[\boxed{\;\hat P(Y=j\mid X=x)=\frac{1}{k}\sum_{i\in A_k(x)}\mathbb{1}\{y_i=j\}\;}\]
This is a plug-in rule in the sense of §5.2: a nonparametric estimate of the posterior, followed by the argmax.
8.1 Distance metrics and categorical features
For \(a,b\in\mathbb{R}^d\) the Minkowski family is
\[d_p(a,b)=\left(\sum_{j=1}^{d}\big|a_j-b_j\big|^p\right)^{1/p},\]
with \(p=2\) Euclidean, \(p=1\) Manhattan, \(p\to\infty\) Chebyshev. Other useful choices: Mahalanobis \(d_M(a,b)=\sqrt{(a-b)^\top\Sigma^{-1}(a-b)}\), which accounts for correlation and scale automatically; cosine distance for text; Gower for mixed numeric and categorical data; and Hamming for purely categorical features.
Nominal features must be encoded before Euclidean distance means anything. One-hot encoding creates one indicator per level:
\[\text{Sex}=\begin{cases}0 & \text{male}\\ 1 & \text{female}\end{cases} \qquad \text{Cold}=\begin{cases}0 & \text{Temp}\ge 37^\circ\mathrm{C}\\ 1 & \text{Temp} < 37^\circ\mathrm{C}\end{cases}\]
For a \(J\)-level factor, create \(J\) indicators. Note the geometric consequence: any two distinct levels are then \(\sqrt2\) apart in Euclidean distance, so all categories are equidistant, which is right for genuinely nominal data and wrong for ordinal data, where the level ordering carries information and should be encoded as a number or as ordered contrasts.
8.2 Rescaling is not optional
Euclidean distance sums squared differences across features, so a feature measured in thousands dominates one measured in units regardless of relevance.
\[\text{min–max: } X_{\text{new}}=\frac{X-\min(X)}{\max(X)-\min(X)}\in[0,1], \qquad \text{z-score: } X_{\text{new}}=\frac{X-\mu}{\sigma}.\]
Min-max bounds the range but is sensitive to outliers (one extreme value compresses everything else). Z-scoring is unbounded but robust to range, and its breakdown point is \(1/n\) (Chapter 2, §2.3.2), for heavy-tailed features, a robust version using the median and MAD is safer.
The scaling constants must come from the training set. Computing them on all the data before splitting is preprocessing leakage (§5.4.1). Inside cross-validation they must be recomputed within each fold. This is fiddly to do by hand, which is precisely why
recipesexists.
# The correct pattern: the recipe LEARNS center and scale from training data,
# then applies those same constants to any new data.
scale_recipe <- function(train_df, outcome = "y") {
recipes::recipe(as.formula(paste(outcome, "~ .")), data = train_df) |>
recipes::step_normalize(recipes::all_numeric_predictors()) |>
recipes::prep(training = train_df)
}9 Why kNN works: the Cover–Hart bound
kNN’s simplicity is disarming, and there is a striking theorem explaining why it is nevertheless competitive.
Cover–Hart theorem (1967). As \(n\to\infty\), the asymptotic risk of the 1-nearest-neighbour classifier satisfies \[R^{*}\;\le\;R_{1\mathrm{NN}}\;\le\;R^{*}\left(2-\frac{K}{K-1}R^{*}\right)\;\le\;2R^{*},\] where \(R^{*}\) is the Bayes error and \(K\) the number of classes.
Sketch (binary case). As \(n\to\infty\) the nearest neighbour of \(x\) converges to \(x\), so its label is a draw from \(\text{Bernoulli}(\eta(x))\) with \(\eta(x)=P(Y=1\mid x)\), independent of the query’s own label. The probability they disagree is \(2\eta(x)(1-\eta(x))\). Since the Bayes error at \(x\) is \(\min(\eta,1-\eta)\), and \(2\eta(1-\eta)\le 2\min(\eta,1-\eta)\), integrating over \(x\) gives the bound. \(\blacksquare\)
Read what this says. A classifier that stores the data and looks up one neighbour is, asymptotically, at worst twice as bad as the theoretically optimal rule, with no model, no assumptions, and no fitting. “In this sense it may be said that half the classification information in an infinite sample set is contained in the nearest neighbour.”
For \(k>1\) with \(k\to\infty\) and \(k/n\to0\), kNN is consistent: \(R_{k\mathrm{NN}}\to R^{*}\). Those two conditions are the whole story, \(k\) must grow (to average away noise) but more slowly than \(n\) (to keep the neighbourhood local).
set.seed(47)
ch <- function(n) {
lab <- rbinom(n, 1, 1 - pi1)
Xn <- rbind(MASS::mvrnorm(max(sum(lab == 0), 1), mu1, S1),
MASS::mvrnorm(max(sum(lab == 1), 1), mu2, S2))
yn <- factor(c(rep("A", max(sum(lab == 0), 1)), rep("B", max(sum(lab == 1), 1))))
labt <- rbinom(3000, 1, 1 - pi1)
Xt <- rbind(MASS::mvrnorm(sum(labt == 0), mu1, S1),
MASS::mvrnorm(sum(labt == 1), mu2, S2))
yt <- factor(c(rep("A", sum(labt == 0)), rep("B", sum(labt == 1))))
k_big <- max(1, round(sqrt(nrow(Xn))))
c(n = nrow(Xn),
err_1nn = mean(knn(Xn, Xt, yn, k = 1) != yt),
err_knn = mean(knn(Xn, Xt, yn, k = k_big) != yt),
k_used = k_big)
}
ch_tab <- as.data.frame(do.call(rbind, lapply(c(50, 200, 1000, 5000, 20000), ch)))
ch_tab$bayes <- bayes_error
ch_tab$upper_bound_2R <- 2 * bayes_error
round(ch_tab, 4)#> n err_1nn err_knn k_used bayes upper_bound_2R
#> 1 50 0.1647 0.1220 7 0.1004 0.2008
#> 2 200 0.1167 0.0987 14 0.1004 0.2008
#> 3 1000 0.1540 0.1053 32 0.1004 0.2008
#> 4 5000 0.1507 0.1037 71 0.1004 0.2008
#> 5 20000 0.1497 0.1090 141 0.1004 0.2008ch_tab |>
dplyr::select(n, `1-NN` = err_1nn, `k-NN (k = sqrt(n))` = err_knn) |>
pivot_longer(-n, names_to = "rule", values_to = "error") |>
ggplot(aes(n, error, colour = rule)) +
geom_hline(aes(yintercept = bayes_error, linetype = "Bayes error R*"),
colour = "black") +
geom_hline(aes(yintercept = 2 * bayes_error, linetype = "Cover-Hart bound 2R*"),
colour = "grey40") +
geom_line(linewidth = 1) + geom_point(size = 2.4) +
scale_x_log10() +
scale_colour_manual(values = c("1-NN" = "#D8433B", "k-NN (k = sqrt(n))" = "#3B7DD8")) +
scale_linetype_manual(values = c("Bayes error R*" = "solid",
"Cover-Hart bound 2R*" = "dashed")) +
labs(title = "1-NN is bounded by twice the Bayes error; k-NN converges to it",
x = "Training set size n (log scale)", y = "Test error",
colour = NULL, linetype = NULL) +
theme_dspa()
1-NN settles between \(R^{*}\) and \(2R^{*}\) exactly as the theorem promises, while \(k\)-NN with \(k\) growing converges down toward the Bayes error.
10 Choosing k: the bias–variance tradeoff
Small \(k\) gives a flexible, high-variance rule that chases individual points — \(k=1\) has zero training error and interpolates every label, including mislabelled ones. Large \(k\) averages over a wide neighbourhood, smoothing away real structure; as \(k\to n\) the classifier degenerates to the majority-class rule.
For a nonparametric regression estimate at \(x\), the local average has
\[\text{Bias}\approx C_1\left(\frac{k}{n}\right)^{2/d}, \qquad \text{Variance}\approx\frac{C_2}{k},\]
so the MSE is minimized at
\[\boxed{\;k^{\star}\asymp n^{4/(d+4)}\;}\]
Two consequences. The optimal \(k\) grows with \(n\), but slowly, like \(n^{4/(d+4)}\). And it shrinks toward a constant as \(d\) grows, because the exponent \(4/(d+4)\to0\); at the same time the achievable error rate deteriorates, which is the curse of dimensionality in its rate form.
Common misconception: “use \(k=\sqrt{n}\).” The square-root rule has no asymptotic justification. It is a serviceable starting value and nothing more: the rate-optimal \(k\) is \(n^{4/(d+4)}\), which equals \(\sqrt n\) only when \(d=4\). For \(d=20\) the exponent is \(1/6\), so with \(n=10{,}000\) the rule of thumb suggests 100 where theory suggests about 5. Cross-validate.
set.seed(53)
tr_idx <- sample(nrow(sim), 800)
tr <- sim[tr_idx, ]; te <- sim[-tr_idx, ]
ks <- c(1, 3, 5, 9, 15, 25, 41, 65, 101, 161, 251, 401, 601, 799)
ks <- ks[ks < nrow(tr) - 1] # simply drop 799
### OLD
# bv <- do.call(rbind, lapply(ks, function(k) {
# data.frame(
# k = k,
# train_error = mean(knn(tr[, 1:2], tr[, 1:2], tr$y, k = k) != tr$y),
# test_error = mean(knn(tr[, 1:2], te[, 1:2], tr$y, k = k) != te$y))
# }))
library(FNN)
bv <- do.call(rbind, lapply(ks, function(k) {
train_pred <- knn(train = tr[, 1:2], test = tr[, 1:2], cl = tr$y, k = k)
test_pred <- knn(train = tr[, 1:2], test = te[, 1:2], cl = tr$y, k = k)
data.frame(
k = k,
train_error = mean(train_pred != tr$y),
test_error = mean(test_pred != te$y)
)
}))
bv$bayes <- bayes_error
round(head(bv, 6), 4)#> k train_error test_error bayes
#> 1 1 0.0000 0.1325 0.1004
#> 2 3 0.0750 0.1025 0.1004
#> 3 5 0.0925 0.1000 0.1004
#> 4 9 0.1050 0.0850 0.1004
#> 5 15 0.1175 0.0775 0.1004
#> 6 25 0.1088 0.0850 0.1004set.seed(59)
# One model per (k, fold) pair -- the sweep genuinely varies k
cvf <- rsample::vfold_cv(tr, v = 10, strata = y)
cv_err <- vapply(ks, function(k) {
mean(vapply(cvf$splits, function(s) {
a <- rsample::analysis(s); b <- rsample::assessment(s)
ctr <- colMeans(a[, 1:2]); scl <- apply(a[, 1:2], 2, sd) # TRAINING stats only
A <- scale(a[, 1:2], ctr, scl); B <- scale(b[, 1:2], ctr, scl)
mean(knn(A, B, a$y, k = k) != b$y)
}, numeric(1)))
}, numeric(1))
bv$cv_error <- cv_err
bv |>
dplyr::select(k, Train = train_error, `10-fold CV` = cv_error, Test = test_error) |>
pivot_longer(-k, names_to = "set", values_to = "error") |>
ggplot(aes(k, error, colour = set)) +
geom_hline(yintercept = bayes_error, linetype = "dashed", colour = "black") +
ggplot2::annotate("text", x = 2, y = bayes_error - 0.012, label = "Bayes error",
hjust = 0, size = 3) +
geom_line(linewidth = 0.9) + geom_point(size = 1.8) +
geom_vline(xintercept = ks[which.min(cv_err)], colour = "grey45",
linetype = "dotted") +
scale_x_log10(breaks = ks) +
scale_colour_manual(values = c(Train = "#7FB069", `10-fold CV` = "#D8433B",
Test = "#3B7DD8")) +
labs(title = "kNN error against k: the classic U-shape",
subtitle = sprintf("CV-optimal k = %d (dotted). Train error is 0 at k = 1 and tells you nothing",
ks[which.min(cv_err)]),
x = "Number of neighbours k (log scale)", y = "Classification error",
colour = NULL) +
theme_dspa()
# --- Interactive equivalent ------------------------------------------------
dataL <- bv |> dplyr::select(K = k, Train = train_error, CV = cv_error, Test = test_error) |>
pivot_longer(-K, names_to = "variable", values_to = "value")
plot_ly(dataL, x = ~K, y = ~value, color = ~variable,
type = "scatter", mode = "markers+lines") |>
add_segments(x = ks[which.min(cv_err)], xend = ks[which.min(cv_err)],
y = 0, yend = 0.35, line = list(color = "gray", dash = "dot"),
showlegend = FALSE) |>
layout(title = "k-NN training, CV, and testing error rates against k",
xaxis = list(title = "Number of nearest neighbours (k)", type = "log"),
yaxis = list(title = "Classification error"),
legend = list(title = list(text = "<b>Sample</b>")))Note that training error is useless for choosing \(k\): it is exactly zero at \(k=1\) and increases monotonically, so minimizing it always selects the most overfit model.
10.1 The tradeoff, made visible
What \(k\) actually controls is the roughness of the posterior surface. This is a three-dimensional object, and it is worth rotating:
gx <- seq(-4, 4, length.out = 80); gy <- seq(-3.5, 3.5, length.out = 80)
gd <- expand.grid(x1 = gx, x2 = gy)
knn_surface <- function(k) {
p <- knn(tr[, 1:2], gd, tr$y, k = k, prob = TRUE)
# class::knn returns the vote share of the WINNING class -- convert to P(Y = B)
pr <- attr(p, "prob")
matrix(ifelse(p == "B", pr, 1 - pr), nrow = length(gx))
}
p3 <- plot_ly()
for (i in seq_along(c(1, 15, 100))) {
k <- c(1, 15, 100)[i]
p3 <- add_surface(p3, x = gy, y = gx, z = knn_surface(k) + (i - 1) * 1.25,
showscale = FALSE, opacity = 0.95, colorscale = "RdBu",
reversescale = TRUE, name = paste("k =", k))
}
p3 |> layout(title = "kNN posterior surfaces at k = 1 (bottom), 15 (middle), 100 (top)",
scene = list(xaxis = list(title = "x2"), yaxis = list(title = "x1"),
zaxis = list(title = "P(Y = B | x) + offset")))At \(k=1\) the surface is a jagged step function, every training point creates its own cell, and the estimate is pure variance. At \(k=100\) it is a smooth ramp that has lost the curvature of the true boundary, pure bias. At \(k=15\) it approximates the true posterior surface from §5.2.2 closely.
11 kNN in high dimensions
Chapter 4 established that pairwise distances concentrate as \(d\) grows: the ratio of their spread to their mean tends to zero. kNN is the primary casualty, because “nearest” stops being meaningfully different from “farthest”.
curse <- function(d, n = 600, n_te = 400, seed = 61) {
set.seed(seed + d)
# Signal lives in the FIRST TWO dimensions only; the rest are pure noise
mk <- function(m) {
lab <- rbinom(m, 1, 0.5)
Z <- matrix(rnorm(m * d), m, d)
Z[lab == 1, 1:2] <- Z[lab == 1, 1:2] + 1.6
list(X = Z, y = factor(ifelse(lab == 1, "B", "A")))
}
a <- mk(n); b <- mk(n_te)
dd <- as.numeric(dist(a$X))
c(d = d,
knn_error = mean(knn(a$X, b$X, a$y, k = 15) != b$y),
knn_error_signal_only = mean(knn(a$X[, 1:2], b$X[, 1:2], a$y, k = 15) != b$y),
relative_spread = sd(dd) / mean(dd))
}
as.data.frame(do.call(rbind, lapply(c(2, 5, 10, 25, 50, 100, 250, 500), curse))) |>
round(4)#> d knn_error knn_error_signal_only relative_spread
#> 1 2 0.1475 0.1475 0.5215
#> 2 5 0.1575 0.1500 0.3301
#> 3 10 0.1875 0.1650 0.2328
#> 4 25 0.1675 0.1350 0.1429
#> 5 50 0.1975 0.1250 0.1027
#> 6 100 0.2425 0.1750 0.0712
#> 7 250 0.2750 0.1325 0.0449
#> 8 500 0.3250 0.1425 0.0311cu <- as.data.frame(do.call(rbind, lapply(c(2, 5, 10, 25, 50, 100, 250, 500), curse)))
cu |>
dplyr::select(d, `All d features` = knn_error,
`Only the 2 informative features` = knn_error_signal_only) |>
pivot_longer(-d, names_to = "features", values_to = "error") |>
ggplot(aes(d, error, colour = features)) +
geom_line(linewidth = 1) + geom_point(size = 2.2) +
scale_x_log10() +
scale_colour_manual(values = c("#D8433B", "#3B7DD8")) +
labs(title = "Adding noise features destroys kNN",
subtitle = "The signal never changes -- it always lives in the first two dimensions",
x = "Ambient dimension d (log scale)", y = "Test error", colour = NULL) +
theme_dspa()
The signal is constant across the whole sweep; only irrelevant dimensions are added. kNN degrades from near-perfect to near-chance, because every noise feature contributes equally to the distance and swamps the two that matter.
The three remedies, in order of preference: reduce dimension first (PCA, UMAP, Chapter 4); select features (Chapter 11); or learn a metric that downweights irrelevant directions (Mahalanobis, large-margin nearest neighbour).
12 Computational complexity
kNN has no training cost, it stores the data. All the work is at prediction time.
| Structure | Build | Query (one point) | Practical range |
|---|---|---|---|
| Brute force | \(O(1)\) | \(O(nd)\) | Always correct; fine for \(n\lesssim10^4\) |
| k-d tree | \(O(dn\log n)\) | \(O(d\log n)\) average | Degrades to \(O(nd)\) above \(d\approx20\) |
| Ball tree | \(O(dn\log n)\) | \(O(d\log n)\) average | Better than k-d in moderate \(d\) |
| Cover tree | \(O(n\log n)\) | \(O(\log n)\) | Depends on intrinsic, not ambient, dimension |
| LSH (approximate) | \(O(nd)\) | \(O(d)\) sublinear | Returns approximate neighbours |
The k-d tree collapse at \(d\approx20\) is the practical fact to remember: in high dimensions the tree must examine nearly every leaf, so its guarantees evaporate exactly when you most wanted them. Cover trees adapt to intrinsic dimension (Chapter 4, §4.3.2), which is why they survive where k-d trees do not.
Memory is \(O(nd)\), the entire training set must be retained, which rules kNN out of memory-constrained deployment.
set.seed(67)
bench <- function(n, d) {
Xtr <- matrix(rnorm(n * d), n, d); ytr <- factor(sample(c("A", "B"), n, TRUE))
Xte <- matrix(rnorm(200 * d), 200, d)
c(n = n, d = d,
seconds = round(system.time(knn(Xtr, Xte, ytr, k = 15))[["elapsed"]], 3))
}
as.data.frame(do.call(rbind, list(bench(1000, 10), bench(10000, 10),
bench(50000, 10), bench(10000, 100))))#> n d seconds
#> 1 1000 10 0.01
#> 2 10000 10 0.04
#> 3 50000 10 0.08
#> 4 10000 100 0.1813 Case study: youth development
The Boys Town Study of Youth Development records 200 adolescents on academic, family, and behavioural measures.
Variables. ID; sex (1 =
male, 2 = female); gpa (0 = “A” average to 5 = “F”);
Alcoholuse (0 = daily to 11 = never); alcatt,
attitudes toward household drinking (0 = totally approve to 6 = totally
disapprove); dadjob, momjob (1 = yes, 2 = no);
dadclose, momclose, parental closeness (0 =
usually to 7 = never); larceny (0–4) and
vandalism (0–7), counts of delinquent acts.
13.1 Preparing the data
boystown <- dspa_read(
"https://umich.instructure.com/files/399119/download?download_frd=1",
"CaseStudy02_Boystown_Data.csv", sep = " ")
# Recode the 1/2 indicators to clean 0/1 dummies
boystown$sex <- boystown$sex - 1 # 0 = male, 1 = female
boystown$dadjob <- 2 - boystown$dadjob # 1 = has job, 0 = no
boystown$momjob <- 2 - boystown$momjob
str(boystown)#> 'data.frame': 200 obs. of 11 variables:
#> $ id : int 1 2 3 4 5 6 7 8 9 10 ...
#> $ sex : num 0 0 0 0 1 1 0 0 1 1 ...
#> $ gpa : int 5 0 3 2 3 3 1 5 1 3 ...
#> $ Alcoholuse: int 2 4 2 2 6 3 2 6 5 2 ...
#> $ alcatt : int 3 2 3 1 2 0 0 3 0 1 ...
#> $ dadjob : num 1 1 1 1 1 1 1 1 1 1 ...
#> $ momjob : num 0 0 0 0 1 0 0 0 1 1 ...
#> $ dadclose : int 1 3 2 1 2 1 3 6 3 1 ...
#> $ momclose : int 1 4 2 2 1 2 1 2 3 2 ...
#> $ larceny : int 1 0 0 3 1 0 0 0 1 1 ...
#> $ vandalism : int 3 0 2 2 2 0 5 1 4 0 ...The outcome is recidivism: two or more combined larceny and vandalism infractions.
bt_y <- factor(ifelse(boystown$larceny + boystown$vandalism > 1,
"Recidivism", "Control"),
levels = c("Control", "Recidivism"))
# Predictors EXCLUDE the ID and the two variables that define the outcome --
# leaving them in would be target leakage (Section 5.4.1)
bt_X <- boystown |> dplyr::select(-any_of(c("id", "ID", "larceny", "vandalism")))
c(n = nrow(bt_X), predictors = ncol(bt_X))#> n predictors
#> 200 8#> bt_y
#> Control Recidivism
#> 56 144#> bt_y
#> Control Recidivism
#> 0.28 0.72The classes are imbalanced at roughly 3:1. Record the no-information rate now, before fitting anything:
#> no_information_rate
#> 0.72Any accuracy near 0.75 means the classifier has matched the majority-class rule, not beaten it.
13.2 Splitting, with the scaling learned from training only
set.seed(1234)
split_bt <- rsample::initial_split(data.frame(bt_X, y = bt_y),
prop = 0.75, strata = y)
bt_train_raw <- rsample::training(split_bt)
bt_test_raw <- rsample::testing(split_bt)
# center and scale using TRAINING statistics, then apply to test
rec_bt <- scale_recipe(bt_train_raw, outcome = "y")
bt_train <- recipes::bake(rec_bt, bt_train_raw)
bt_test <- recipes::bake(rec_bt, bt_test_raw)
bt_train_X <- as.matrix(dplyr::select(bt_train, -y)); bt_train_y <- bt_train$y
bt_test_X <- as.matrix(dplyr::select(bt_test, -y)); bt_test_y <- bt_test$y
c(train = nrow(bt_train_X), test = nrow(bt_test_X))#> train test
#> 150 50#> Control Recidivism
#> train 0.28 0.72
#> test 0.28 0.72# The test columns are NOT mean 0 / sd 1 -- and that is exactly right
round(head(colMeans(bt_test_X), 4), 4)#> sex gpa Alcoholuse alcatt
#> 0.0552 -0.0497 -0.2891 -0.1674A stratified split keeps the class proportions matched. The test columns do not have mean zero, because they were transformed with the training constants — that asymmetry is the signature of a correctly built pipeline.
13.3 Extracting the posterior correctly
Common misconception: “
attr(knn_result, 'prob')is the probability of the positive class.” It is the vote share of the winning class, whichever that was. For a binary problem it is always at least 0.5, so thresholding it at 0.6 does not select high-probability positives, it selects cases where the majority for either class was decisive, relabelling confident negatives as positives. Convert it explicitly.
knn_prob <- function(train_X, test_X, train_y, k, positive) {
p <- knn(train_X, test_X, train_y, k = k, prob = TRUE)
win <- attr(p, "prob") # vote share of the WINNING class
list(class = p,
prob = ifelse(p == positive, win, 1 - win)) # P(Y = positive | x)
}demo <- knn_prob(bt_train_X, bt_test_X, bt_train_y, k = 11, positive = "Recidivism")
head(data.frame(predicted = demo$class,
raw_prob_attribute = round(attr(knn(bt_train_X, bt_test_X,
bt_train_y, k = 11,
prob = TRUE), "prob"), 3),
P_recidivism = round(demo$prob, 3)), 8)#> predicted raw_prob_attribute P_recidivism
#> 1 Recidivism 0.909 0.909
#> 2 Recidivism 0.818 0.818
#> 3 Recidivism 0.818 0.818
#> 4 Recidivism 0.636 0.636
#> 5 Recidivism 0.818 0.818
#> 6 Recidivism 0.818 0.818
#> 7 Recidivism 0.909 0.909
#> 8 Recidivism 0.727 0.727Look at the middle column: it never drops below 0.5, and it takes the same value for a confident “Control” as for a confident “Recidivism”. Only the third column is a posterior probability.
13.4 Choosing k honestly
set.seed(1234)
folds_bt <- rsample::vfold_cv(bt_train_raw, v = 10, repeats = 5, strata = y)
k_grid <- seq(1, 45, by = 2)
cv_bt <- vapply(k_grid, function(k) {
mean(vapply(folds_bt$splits, function(s) {
a <- rsample::analysis(s); b <- rsample::assessment(s)
r <- scale_recipe(a, "y") # rescale INSIDE the fold
A <- recipes::bake(r, a); B <- recipes::bake(r, b)
## OLD: errors since comparing two factors (knn predictions and the true B$y)
## with !=, but their levels do not match exactly. In a cross‑validation fold,
## it is possible that either the training set A or the assessment set B
## contains only one class (even with stratification, if the dataset is small
## or the class balance is extreme).
# mean(knn(dplyr::select(A, -y), dplyr::select(B, -y), A$y, k = k) != B$y)
## WORKING
# mean(as.character(knn(dplyr::select(A, -y), dplyr::select(B, -y), A$y, k = k))
# != as.character(B$y))
full_levels <- levels(bt_train_raw$y)
pred <- knn(dplyr::select(A, -y), dplyr::select(B, -y), A$y, k = k)
pred <- factor(pred, levels = full_levels)
B_y <- factor(B$y, levels = full_levels)
mean(pred != B_y)
}, numeric(1)))
}, numeric(1))
k_best <- k_grid[which.min(cv_bt)]
c(cv_optimal_k = k_best, cv_error = round(min(cv_bt), 4),
NIR_error = round(1 - NIR_bt, 4))#> cv_optimal_k cv_error NIR_error
#> 31.0000 0.2796 0.2800ggplot(data.frame(k = k_grid, error = cv_bt), aes(k, error)) +
geom_hline(yintercept = 1 - NIR_bt, linetype = "dashed", colour = "firebrick") +
ggplot2::annotate("text", x = 40, y = 1 - NIR_bt + 0.008,
label = "Majority-class error", size = 3, colour = "firebrick") +
geom_line(linewidth = 1, colour = "steelblue") +
geom_point(size = 2) +
geom_point(data = data.frame(k = k_best, error = min(cv_bt)),
colour = "firebrick", size = 4) +
scale_x_continuous(breaks = k_grid) +
labs(title = "Repeated stratified 10-fold CV error for kNN",
subtitle = "Scaling recomputed inside every fold; 5 repeats to stabilize the estimate",
x = "k", y = "CV classification error") +
theme_dspa()
13.5 Evaluating on held-out data
fit_bt <- knn_prob(bt_train_X, bt_test_X, bt_train_y, k = k_best,
positive = "Recidivism")
cm_bt <- confusionMatrix(data = fit_bt$class, reference = bt_test_y,
positive = "Recidivism")
cm_bt#> Confusion Matrix and Statistics
#>
#> Reference
#> Prediction Control Recidivism
#> Control 0 0
#> Recidivism 14 36
#>
#> Accuracy : 0.72
#> 95% CI : (0.575, 0.838)
#> No Information Rate : 0.72
#> P-Value [Acc > NIR] : 0.571392
#>
#> Kappa : 0
#>
#> Mcnemar's Test P-Value : 0.000512
#>
#> Sensitivity : 1.00
#> Specificity : 0.00
#> Pos Pred Value : 0.72
#> Neg Pred Value : NaN
#> Prevalence : 0.72
#> Detection Rate : 0.72
#> Detection Prevalence : 1.00
#> Balanced Accuracy : 0.50
#>
#> 'Positive' Class : Recidivism
#> #> n TP FP FN
#> 50.0000 36.0000 14.0000 0.0000
#> TN accuracy NIR sensitivity
#> 0.0000 0.7200 0.7200 1.0000
#> specificity PPV NPV balanced_accuracy
#> 0.0000 0.7200 NaN 0.5000
#> F1 kappa MCC
#> 0.8372 0.0000 NA# Is the accuracy distinguishable from the majority-class rule?
n_te <- length(bt_test_y)
## OLD: same bug as before; since fit_bt$class and bt_test_y are factors with
## different level sets (one may have only a subset of levels, or the levels are
## in a different order/name set). Comparing factors with == requires identical
## level sets
# binom.test(sum(fit_bt$class == bt_test_y), n_te, p = NIR_bt,
# alternative = "greater")
pred <- factor(fit_bt$class, levels = levels(bt_test_y))
sum(pred == bt_test_y)#> [1] 36n_correct <- sum(as.character(fit_bt$class) == as.character(bt_test_y))
binom.test(n_correct, n_te, p = NIR_bt, alternative = "greater")#>
#> Exact binomial test
#>
#> data: n_correct and n_te
#> number of successes = 36, number of trials = 50, p-value = 0.571
#> alternative hypothesis: true probability of success is greater than 0.72
#> 95 percent confidence interval:
#> 0.597383 1.000000
#> sample estimates:
#> probability of success
#> 0.72Read this honestly. The accuracy sits close to the no-information rate and the binomial test does not reject. With 50 test cases, the 95% confidence interval on accuracy spans roughly 15 percentage points, the sample is simply too small to distinguish this classifier from the majority-class rule. That is a finding about the study design, not a defect in the method, and it is exactly the kind of statement a reported point accuracy would conceal.
roc_bt <- pROC::roc(bt_test_y, fit_bt$prob, levels = c("Control", "Recidivism"),
direction = "<", quiet = TRUE)
auc_bt <- as.numeric(pROC::auc(roc_bt))
ci_bt <- as.numeric(pROC::ci.auc(roc_bt))
c(AUC = round(auc_bt, 4), CI_lower = round(ci_bt[1], 4),
CI_upper = round(ci_bt[3], 4))#> AUC CI_lower CI_upper
#> 0.5179 0.3455 0.6902ggplot(data.frame(fpr = 1 - roc_bt$specificities, tpr = roc_bt$sensitivities),
aes(fpr, tpr)) +
geom_abline(slope = 1, intercept = 0, linetype = "dashed", colour = "grey50") +
geom_line(linewidth = 1, colour = "steelblue") +
coord_fixed() +
labs(title = "kNN ROC curve on held-out Boys Town data",
subtitle = sprintf("AUC = %.3f (95%% CI %.3f - %.3f); the CI straddles chance",
auc_bt, ci_bt[1], ci_bt[3]),
x = "False positive rate", y = "True positive rate") +
theme_dspa()
13.6 Choosing an operating point
The default 0.5 threshold minimizes 0–1 loss, which assumes false positives and false negatives cost the same. They rarely do. Sweeping the threshold traces the sensitivity/specificity tradeoff explicitly.
thr_grid <- seq(0.05, 0.95, by = 0.05)
thr_res <- do.call(rbind, lapply(thr_grid, function(t) {
cl <- factor(ifelse(fit_bt$prob >= t, "Recidivism", "Control"),
levels = levels(bt_test_y))
m <- binary_metrics(cl, bt_test_y, positive = "Recidivism")
data.frame(threshold = t, sensitivity = m[["sensitivity"]],
specificity = m[["specificity"]],
balanced_accuracy = m[["balanced_accuracy"]],
accuracy = m[["accuracy"]])
}))
thr_res |>
pivot_longer(-threshold, names_to = "metric", values_to = "value") |>
ggplot(aes(threshold, value, colour = metric)) +
geom_line(linewidth = 0.9) + geom_point(size = 1.5) +
geom_vline(xintercept = 0.5, linetype = "dotted", colour = "grey40") +
scale_colour_brewer(palette = "Set1") +
labs(title = "Metrics across the decision threshold",
subtitle = "The default 0.5 (dotted) is a choice about relative costs, not a neutral default",
x = "Threshold on P(Recidivism | x)", y = NULL, colour = NULL) +
theme_dspa()
# --- Interactive equivalent ------------------------------------------------
plot_ly(thr_res, x = ~threshold, y = ~sensitivity, type = "scatter",
mode = "lines+markers", name = "Sensitivity") |>
add_trace(y = ~specificity, name = "Specificity") |>
add_trace(y = ~balanced_accuracy, name = "Balanced accuracy") |>
layout(title = "Metric tradeoffs across the decision threshold",
xaxis = list(title = "Threshold on P(Recidivism | x)"),
yaxis = list(title = "Metric value"),
legend = list(orientation = "h"))If missing a recidivism case is four times as costly as a false alarm, the threshold minimizing expected cost is well below 0.5. Choose it from the training or validation data, never by scanning the test set for the value that looks best, that is threshold selection on the test set, a subtle but real form of leakage.
13.7 A second application: galaxy spins
SOCR Case Study 22 records morphological measurements on ~51,000 galaxies with a left- or right-handed spin label.
galaxy <- dspa_read(
"https://umich.instructure.com/files/6105118/download?download_frd=1",
"galaxy_spins.csv", sep = ",")
galaxy <- galaxy[, -1] # drop the galaxy ID
dim(galaxy); str(galaxy)#> [1] 51122 11#> 'data.frame': 51122 obs. of 11 variables:
#> $ RA : num 236 237 238 238 238 ...
#> $ DEC : num -0.493 -0.482 -0.506 -0.544 -0.527 ...
#> $ HAND : chr "R" "L" "L" "L" ...
#> $ UZS : num 17.4 19.2 19.4 19.5 17.8 ...
#> $ GZS : num 16.2 18 17.5 18.3 16.6 ...
#> $ RZS : num 15.6 17.4 16.6 17.9 16.2 ...
#> $ IZS : num 15.2 17 16.2 17.6 15.8 ...
#> $ ZZS : num 14.9 16.7 15.8 17.3 15.6 ...
#> $ ELLIPS: num 0.825 0.961 0.411 0.609 0.505 ...
#> $ PHIS : num 154 100.2 29.2 109.4 39.6 ...
#> $ RSS : num 0.0547 0.0977 0.0784 0.076 0.0796 ...set.seed(1234) # seeded: the split is reproducible
label_col <- 3
gal_y <- factor(galaxy[[label_col]])
gal_X <- galaxy[, -label_col, drop = FALSE]
split_gal <- rsample::initial_split(data.frame(gal_X, y = gal_y),
prop = 0.96, strata = y)
gal_tr <- rsample::training(split_gal); gal_te <- rsample::testing(split_gal)
rec_gal <- scale_recipe(gal_tr, "y") # scaling from TRAINING only
Gtr <- recipes::bake(rec_gal, gal_tr); Gte <- recipes::bake(rec_gal, gal_te)
c(train = nrow(Gtr), test = nrow(Gte),
NIR = round(max(prop.table(table(gal_y))), 4))#> train test NIR
#> 49077.0000 2045.0000 0.5042k_gal <- c(1, 5, 11, 21, 51, 101, 201)
set.seed(71)
val <- rsample::validation_split(Gtr, prop = 0.85, strata = y)
v_tr <- rsample::analysis(val$splits[[1]]); v_va <- rsample::assessment(val$splits[[1]])
gal_err <- vapply(k_gal, function(k)
mean(knn(dplyr::select(v_tr, -y), dplyr::select(v_va, -y), v_tr$y, k = k) != v_va$y),
numeric(1))
data.frame(k = k_gal, validation_error = round(gal_err, 4))#> k validation_error
#> 1 1 0.0906
#> 2 5 0.3459
#> 3 11 0.4308
#> 4 21 0.4471
#> 5 51 0.4687
#> 6 101 0.4745
#> 7 201 0.4796k_gal_best <- k_gal[which.min(gal_err)]
gal_fit <- knn(dplyr::select(Gtr, -y), dplyr::select(Gte, -y), Gtr$y, k = k_gal_best)
confusionMatrix(data = gal_fit, reference = Gte$y)$overall[1:6] |> round(4)#> Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull
#> 0.9883 0.9765 0.9826 0.9925 0.5042
#> AccuracyPValue
#> 0.0000ggplot(data.frame(k = k_gal, error = gal_err), aes(k, error)) +
geom_line(linewidth = 1, colour = "steelblue") + geom_point(size = 2.2) +
geom_point(data = data.frame(k = k_gal_best, error = min(gal_err)),
colour = "firebrick", size = 4) +
scale_x_log10(breaks = k_gal) +
labs(title = "Galaxy spin: validation error against k",
subtitle = "Features scaled with training statistics; split seeded for reproducibility",
x = "k (log scale)", y = "Validation error") +
theme_dspa()
# --- Interactive equivalent, plus the e1071 tuner --------------------------
library(e1071)
knn_tune <- tune.knn(x = as.matrix(dplyr::select(Gtr, -y)), y = Gtr$y, k = 1:20,
tunecontrol = tune.control(sampling = "fix", fix = 10))
summary(knn_tune)
plot_ly(x = knn_tune$performance$k, y = knn_tune$performance$error,
type = "scatter", mode = "markers+lines") |>
layout(title = "Galaxy-spin kNN prediction: error rate against k",
xaxis = list(title = "Number of nearest neighbours (k)"),
yaxis = list(title = "Classification error"))Note that fix belongs inside
tune.control(); passed at the top level it is silently
absorbed and the intended fixed-split fraction never takes effect.
14 PART III: PROBABILISTIC LEARNING
15 Naive Bayes
kNN estimates the posterior by local voting. Naive Bayes estimates it by modelling how the features are generated within each class, a generative rather than discriminative approach, and then inverting with Bayes’ rule.
15.1 From Bayes’ theorem to a classifier
For events \(A\) and \(B\),
\[P(A\mid B)=\frac{P(B\mid A)\,P(A)}{P(B)} \qquad\text{i.e.}\qquad \text{posterior}=\frac{\text{likelihood}\times\text{prior}}{\text{marginal likelihood}} .\]
When \(\{B_i\}\) partition the sample space, the law of total probability expands the denominator:
\[P(A\mid B)=\frac{P(B\mid A)P(A)}{\sum_{i} P(B\mid B_i)P(B_i)} .\]
For classification with features \(F_1,\dots,F_n\) and classes \(C_1,\dots,C_K\):
\[P(C_k\mid F_1,\dots,F_n)=\frac{P(F_1,\dots,F_n\mid C_k)\,P(C_k)}{P(F_1,\dots,F_n)} .\]
The denominator does not depend on \(k\), so only the numerator matters for the argmax. Applying the chain rule to that numerator:
\[ \begin{aligned} P(F_1,\dots,F_n,C_k) &= P(F_1\mid F_2,\dots,F_n,C_k)\,P(F_2,\dots,F_n,C_k)\\ &= P(F_1\mid F_2,\dots,F_n,C_k)\,P(F_2\mid F_3,\dots,F_n,C_k)\,P(F_3,\dots,F_n,C_k)\\ &= \cdots\\ &= P(C_k)\prod_{i=1}^{n}P\big(F_i\mid F_{i+1},\dots,F_n,C_k\big). \end{aligned} \]
This is exact, and useless, because estimating \(P(F_i\mid F_{i+1},\dots,F_n,C_k)\) requires a table with \(O(2^n)\) cells. With \(n=118\) binary features there are more cells than atoms in the observable universe.
15.2 The naive assumption
Assume each feature is conditionally independent of every other, given the class:
\[\boxed{\;P\big(F_i\mid F_{i+1},\dots,F_n,\ C_k\big)=P(F_i\mid C_k)\;}\]
The joint model collapses to a product of \(n\) one-dimensional tables:
\[P(F_1,\dots,F_n,C_k)=P(C_k)\prod_{i=1}^{n}P(F_i\mid C_k),\]
and the classifier is the maximum a posteriori rule
\[\hat C=\arg\max_{k}\;\frac{P(C_k)\prod_{i=1}^{n}P(F_i\mid C_k)}{P(F_1,\dots,F_n)} =\arg\max_{k}\;P(C_k)\prod_{i=1}^{n}P(F_i\mid C_k),\]
since the denominator is a constant across \(k\).
Parameter count drops from \(O(K\cdot2^n)\) to \(O(K\cdot n)\), the reason naive Bayes trains in a single pass and works with \(n\gg N\), the regime where text classification lives.
15.3 Why it works when the assumption is false
Conditional independence is essentially never true. In clinical text, “metastatic” and “palliative” co-occur far more than independence predicts. Yet naive Bayes remains competitive. Why?
Common misconception: “naive Bayes needs the features to be independent.” It needs that assumption for its probability estimates to be accurate. It does not need it for its classifications to be correct, and classification is usually the goal.
The argmax is invariant to any monotone distortion that preserves the ordering of the class scores. Dependent features cause the same evidence to be counted repeatedly, which drives the posterior toward 0 or 1, naive Bayes is famously over-confident. But if the inflation affects the correct class most, the ranking survives and the label is right. Domingos and Pazzani (1997) showed that naive Bayes is optimal under conditions far weaker than independence, including for many functions with strong feature dependence.
The practical implication: trust naive Bayes labels more than naive Bayes probabilities. If you need calibrated probabilities, post-process with Platt scaling or isotonic regression (§5.3.5).
set.seed(79)
# Duplicate one informative feature many times: independence is maximally violated
make_dup <- function(n_dup) {
n <- 800
y <- factor(sample(c("A", "B"), n, TRUE))
z <- rnorm(n) + ifelse(y == "B", 1.4, 0)
X <- matrix(rep(z, n_dup), n, n_dup) + matrix(rnorm(n * n_dup, sd = 0.01), n, n_dup)
data.frame(X, y = y)
}
res_dup <- do.call(rbind, lapply(c(1, 3, 10, 30), function(nd) {
d <- make_dup(nd)
tr <- d[1:500, ]; te <- d[501:800, ]
m <- e1071::naiveBayes(y ~ ., data = tr)
p <- predict(m, te, type = "raw")[, "B"]
cl <- predict(m, te)
data.frame(duplicates = nd,
accuracy = mean(cl == te$y),
mean_max_prob = mean(pmax(p, 1 - p)),
brier = mean((p - as.integer(te$y == "B"))^2))
}))
round(res_dup, 4)#> duplicates accuracy mean_max_prob brier
#> 1 1 0.7667 0.7788 0.1501
#> 2 3 0.7700 0.8913 0.1810
#> 3 10 0.7600 0.9693 0.2157
#> 4 30 0.7500 0.9923 0.2462Accuracy barely moves as the same feature is duplicated thirty times, while the mean confidence climbs toward 1 and the Brier score deteriorates. The labels survive; the probabilities do not.
15.4 Smoothing: Laplace as a Dirichlet prior
If a word never appears in the training documents of class \(C_L\), then \(\hat P(F_i\mid C_L)=0\), and the entire product \(\prod_i \hat P(F_i\mid C_L)\) collapses to zero, one unseen feature vetoes the class regardless of all other evidence. A zero from a finite sample should not be treated as a zero in the population.
Lidstone smoothing adds a pseudo-count \(\alpha\) to every cell:
\[\hat P(F_i=v\mid C_k)=\frac{N_{ikv}+\alpha}{N_{k}+\alpha V_i},\]
where \(N_{ikv}\) counts occurrences, \(N_k\) the class total, and \(V_i\) the number of levels. \(\alpha=1\) is Laplace smoothing; \(\alpha=1/2\) is the Jeffreys prior.
This is not an ad-hoc fix. Placing a \(\text{Dirichlet}(\alpha,\dots,\alpha)\) prior on the multinomial cell probabilities gives a Dirichlet posterior whose mean is exactly the smoothed estimate, so \(\alpha\) is the strength of a prior belief that all levels are possible, measured in pseudo-observations.
\[\theta\sim\text{Dir}(\alpha\mathbf{1}),\quad N\mid\theta\sim\text{Mult}(N_k,\theta) \;\Longrightarrow\; \theta\mid N\sim\text{Dir}(\alpha\mathbf{1}+N),\quad E[\theta_v\mid N]=\frac{N_v+\alpha}{N_k+\alpha V}.\]
\(\alpha\) is a hyperparameter and must be tuned. Too small and zeros persist; too large and every conditional probability is dragged toward \(1/V_i\), flattening the likelihood until the classifier reduces to the prior, it predicts the majority class for everything, while accuracy rises to the no-information rate and looks like an improvement. Tune \(\alpha\) by held-out log-loss, which is sensitive to the probabilities, rather than by accuracy, which is not.
set.seed(83)
V <- 5; N_k <- 20
counts <- c(12, 5, 3, 0, 0)
alphas <- c(0, 0.1, 0.5, 1, 5, 20, 100)
smooth_tab <- do.call(rbind, lapply(alphas, function(a)
data.frame(alpha = a, level = factor(1:V),
prob = (counts + a) / (N_k + a * V))))
ggplot(smooth_tab, aes(level, prob, fill = factor(alpha))) +
geom_col(position = position_dodge(width = 0.85), width = 0.8) +
geom_hline(yintercept = 1 / V, linetype = "dashed", colour = "grey40") +
ggplot2::annotate("text", x = 4.6, y = 1/V + 0.02, label = "uniform 1/V", size = 3) +
scale_fill_viridis_d(name = expression(alpha)) +
labs(title = "Smoothing pulls every estimate toward the uniform distribution",
subtitle = "Counts (12, 5, 3, 0, 0) out of 20; large alpha erases the signal entirely",
x = "Feature level", y = expression(hat(P)(F == v ~ "|" ~ C))) +
theme_dspa()
15.5 Computing in log space
\[\hat C=\arg\max_k\;P(C_k)\prod_{i=1}^{n}P(F_i\mid C_k)\]
is correct mathematically and catastrophic numerically. With 118 features each contributing a probability around \(0.1\), the product is \(10^{-118}\); with 500 features it is \(10^{-500}\), which is exactly zero in double precision (underflow below \(\approx1.8\times10^{-308}\)). Every class scores zero and the argmax is arbitrary.
Every production implementation works with the monotone transform
\[\boxed{\;\hat C=\arg\max_k\left[\log P(C_k)+\sum_{i=1}^{n}\log P(F_i\mid C_k)\right]\;}\]
To recover normalized probabilities, use the log-sum-exp trick with a maximum subtracted for stability:
\[\log\sum_k e^{s_k}=m+\log\sum_k e^{s_k-m},\qquad m=\max_k s_k .\]
n_feat <- 500
p_per_feature <- 0.1
c(naive_product = prod(rep(p_per_feature, n_feat)), # underflows to 0
log_space_sum = sum(log(rep(p_per_feature, n_feat))), # finite and exact
recovered = exp(sum(log(rep(p_per_feature, n_feat))))) # underflows on exp#> naive_product log_space_sum recovered
#> 0.00 -1151.29 0.00log_sum_exp <- function(s) { m <- max(s); m + log(sum(exp(s - m))) }
scores <- c(-1150.2, -1148.7, -1153.9) # three class log-scores
round(exp(scores - log_sum_exp(scores)), 6) # normalized, no underflow#> [1] 0.181606 0.813904 0.004490The direct product is zero and unusable; the log-space computation is exact, and log-sum-exp recovers well-defined posteriors from log-scores near \(-1150\).
15.6 Naive Bayes variants
The naive assumption fixes the factorization; the choice of \(P(F_i\mid C_k)\) fixes the model.
| Variant | \(P(F_i\mid C_k)\) | Feature type | Typical use |
|---|---|---|---|
| Bernoulli | \(p_{ik}^{F_i}(1-p_{ik})^{1-F_i}\) | Binary presence | Short text; explicitly models absence |
| Multinomial | \(\propto\prod_i\theta_{ik}^{F_i}\) | Counts | Longer documents; uses term frequency |
| Gaussian | \(\mathcal N(\mu_{ik},\sigma_{ik}^2)\) | Continuous | Numeric features |
| Categorical | Frequency table | Nominal | Mixed survey data |
e1071::naiveBayes() dispatches on column type: factors
get frequency tables (with laplace applied), numerics get
Gaussians. laplace has no effect on numeric
columns, smoothing a Gaussian mean is not defined, so passing a
large laplace alongside continuous predictors silently does
nothing.
15.7 Relationships to LDA, QDA, and logistic regression
All four are plug-in rules; they differ in what they assume about \(P(X\mid Y)\) or model directly.
| Method | Class-conditional model | Boundary | Parameters |
|---|---|---|---|
| Gaussian NB | \(\mathcal N(\mu_k,\operatorname{diag}(\sigma_k^2))\) | Quadratic | \(2Kd\) |
| LDA | \(\mathcal N(\mu_k,\Sigma)\), shared \(\Sigma\) | Linear | \(Kd+d(d+1)/2\) |
| QDA | \(\mathcal N(\mu_k,\Sigma_k)\), separate | Quadratic | \(Kd+Kd(d+1)/2\) |
| Logistic | none — models \(\log\frac{P(Y=1\mid x)}{P(Y=0\mid x)}\) directly | Linear | \(d+1\) |
Gaussian NB is QDA with diagonal covariances: it allows different variances per class (hence a quadratic boundary) but forbids within-class correlation. LDA allows correlation but forces it to be shared, which makes the quadratic terms cancel and the boundary linear.
The bias–variance ordering follows the parameter counts: logistic and LDA are most constrained and most stable in small samples; QDA is most flexible and needs \(n_k \gg d^2/2\) per class to estimate its covariances.
sim_tr <- sim[tr_idx, ]
gd2 <- expand.grid(x1 = gx, x2 = gy)
nb_m <- e1071::naiveBayes(y ~ x1 + x2, data = sim_tr)
lda_m <- MASS::lda(y ~ x1 + x2, data = sim_tr)
qda_m <- MASS::qda(y ~ x1 + x2, data = sim_tr)
surf <- list(
`Naive Bayes` = matrix(predict(nb_m, gd2, type = "raw")[, "B"], length(gx)),
LDA = matrix(predict(lda_m, gd2)$posterior[, "B"], length(gx)),
QDA = matrix(predict(qda_m, gd2)$posterior[, "B"], length(gx)),
`True Bayes` = matrix(post2(gd2$x1, gd2$x2), length(gx)))
p4 <- plot_ly()
for (i in seq_along(surf)) {
p4 <- add_surface(p4, x = gy, y = gx, z = surf[[i]] + (i - 1) * 1.25,
showscale = FALSE, opacity = 0.95, colorscale = "RdBu",
reversescale = TRUE, name = names(surf)[i])
}
p4 |> layout(title = "Posterior surfaces: naive Bayes, LDA, QDA, and the true Bayes rule (bottom to top)",
scene = list(xaxis = list(title = "x2"), yaxis = list(title = "x1"),
zaxis = list(title = "P(Y = B | x) + offset")))sim_te <- sim[-tr_idx, ]
cmp_models <- data.frame(
model = c("Naive Bayes", "LDA", "QDA", "Logistic", "Bayes (oracle)"),
test_error = c(
mean(predict(nb_m, sim_te) != sim_te$y),
mean(predict(lda_m, sim_te)$class != sim_te$y),
mean(predict(qda_m, sim_te)$class != sim_te$y),
mean(factor(ifelse(predict(glm(y ~ x1 + x2, sim_tr, family = binomial()),
sim_te, type = "response") > 0.5, "B", "A"),
levels = levels(sim_te$y)) != sim_te$y),
bayes_error))
cmp_models |> mutate(test_error = round(test_error, 4))#> model test_error
#> 1 Naive Bayes 0.1025
#> 2 LDA 0.0925
#> 3 QDA 0.0750
#> 4 Logistic 0.0800
#> 5 Bayes (oracle) 0.1004QDA is closest to the oracle, because the data were generated with unequal covariances, precisely the structure QDA models and LDA cannot. Naive Bayes loses a little by forcing the within-class correlation to zero. When the model matches the generating process, the plug-in rule approaches the Bayes error.
16 Case study: head and neck cancer medication
16.1 The data
Inpatient head-and-neck cancer medication records: 662 encounters
with a free text MEDICATION_SUMMARY and a SEER stage
code.
Variables. PID, ENC_ID;
seer_stage (0 = in situ, 1 = localized, 2 = regional by
direct extension, 3 = regional to lymph nodes, 4 = regional both, 5 =
regional NOS, 7 = distant metastases, 8 = not applicable, 9 = unstaged),
see SEER SSM;
MEDICATION_DESC, MEDICATION_SUMMARY,
DOSE, UNIT, FREQUENCY,
TOTAL_DOSE_COUNT.
hn_med <- dspa_read(
"https://umich.instructure.com/files/1614350/download?download_frd=1",
"HeadNeck_Cancer_Medication.csv", stringsAsFactors = FALSE)
str(hn_med)#> 'data.frame': 662 obs. of 9 variables:
#> $ PID : int 10000 10008 10029 10063 10071 10103 1012 10135 10136 10143 ...
#> $ ENC_ID : int 46836 46886 47034 47240 47276 47511 3138 47739 47744 47769 ...
#> $ seer_stage : int 1 1 4 1 9 1 1 1 9 1 ...
#> $ MEDICATION_DESC : chr "ranitidine" "heparin injection" "ampicillin/sulbactam IVPB UH" "fentaNYL injection UH" ...
#> $ MEDICATION_SUMMARY: chr "(Zantac) 150 mg tablet oral two times a day" "5,000 unit subcutaneous three times a day" "(Unasyn) 15 g IV every 6 hours" "25 - 50 microgram IV every 5 minutes PRN severe pain\nMaximum dose 200 mcg Per PACU protocol" ...
#> $ DOSE : chr "150" "5000" "1.5" "50" ...
#> $ UNIT : chr "mg" "unit" "g" "microgram" ...
#> $ FREQUENCY : chr "two times a day" "three times a day" "every 6 hours" "every 5 minutes" ...
#> $ TOTAL_DOSE_COUNT : int 5 3 11 2 1 2 2 6 15 1 ...#>
#> 0 1 2 3 4 5 7 8 9
#> 21 265 53 90 46 18 87 14 6816.2 Text preprocessing
#> [1] "(Zantac) 150 mg tablet oral two times a day"# content_transformer() wraps a plain function so the corpus keeps its
# PlainTextDocument structure -- without it, downstream tm functions break.
hn_clean <- hn_corpus |>
tm_map(content_transformer(tolower)) |>
tm_map(removePunctuation) |>
tm_map(removeNumbers) |>
tm_map(removeWords, stopwords("english")) |> # drop "the", "and", "for", ...
tm_map(stemDocument) |> # collapse inflectional variants
tm_map(stripWhitespace)
hn_clean[[1]]$content#> [1] "zantac mg tablet oral two time day"Stopword removal and stemming are not cosmetic. Without them the most frequent terms, and therefore the highest-variance features, are function words that carry no clinical signal, and morphological variants of the same drug (“tablet”, “tablets”) are treated as unrelated features.
16.3 Splitting before anything else
The dictionary must be built from training documents only. Selecting vocabulary by document frequency across the whole corpus lets test documents vote on which features exist, preprocessing leakage (§5.4.1).
set.seed(12345)
n_hn <- nrow(hn_med)
train_idx <- sample(n_hn, floor(0.8 * n_hn))
hn_train_meta <- hn_med[train_idx, ]; hn_test_meta <- hn_med[-train_idx, ]
corpus_train <- hn_clean[train_idx]; corpus_test <- hn_clean[-train_idx]
# Dichotomize stage: 0-4 early / no stage; 5-9 later stage
mk_stage <- function(s) factor(s %in% 5:9, levels = c(FALSE, TRUE),
labels = c("early_stage", "later_stage"))
hn_train_y <- mk_stage(hn_train_meta$seer_stage)
hn_test_y <- mk_stage(hn_test_meta$seer_stage)
round(rbind(train = prop.table(table(hn_train_y)),
test = prop.table(table(hn_test_y))), 3)#> early_stage later_stage
#> train 0.730 0.270
#> test 0.669 0.331#> no_information_rate
#> 0.6692dtm_train_full <- DocumentTermMatrix(corpus_train)
# Vocabulary from TRAINING documents only
hn_dict <- findFreqTerms(dtm_train_full, lowfreq = 5)
c(vocabulary_size = length(hn_dict))#> vocabulary_size
#> 93hn_train_dtm <- DocumentTermMatrix(corpus_train, list(dictionary = hn_dict))
hn_test_dtm <- DocumentTermMatrix(corpus_test, list(dictionary = hn_dict))
dim(hn_train_dtm); dim(hn_test_dtm)#> [1] 529 93#> [1] 133 9316.4 Word clouds
library(wordcloud)
set.seed(101)
wordcloud(corpus_train, min.freq = 30, random.order = FALSE,
colors = RColorBrewer::brewer.pal(5, "Dark2"))
# Both clouds built from the SAME cleaned corpus, so they are comparable
## OLD - buggy
# early_txt <- sapply(corpus_train[hn_train_y == "early_stage"], \(d) d$content)
# later_txt <- sapply(corpus_train[hn_train_y == "later_stage"], \(d) d$content)
early_txt <- corpus_train$content[hn_train_y == "early_stage"]
later_txt <- corpus_train$content[hn_train_y == "later_stage"]
op <- par(mfrow = c(1, 2), mar = c(0, 0, 2, 0))
set.seed(102); wordcloud(early_txt, max.words = 25,
colors = RColorBrewer::brewer.pal(3, "Dark2"))
title("Early stage")
set.seed(103); wordcloud(later_txt, max.words = 25,
colors = RColorBrewer::brewer.pal(3, "Dark2"))
title("Later stage")
16.5 Fitting and tuning
# Binary presence indicators -- the Bernoulli NB representation
to_indicator <- function(dtm) {
m <- as.matrix(dtm)
storage.mode(m) <- "integer"
as.data.frame(lapply(as.data.frame(m > 0),
\(v) factor(v, levels = c(FALSE, TRUE), labels = c("No", "Yes"))))
}
hn_train_X <- to_indicator(hn_train_dtm)
hn_test_X <- to_indicator(hn_test_dtm)
dim(hn_train_X)#> [1] 529 93Hold the feature representation fixed while tuning \(\alpha\). Changing both the smoothing and the encoding at once makes any accuracy change uninterpretable, the classic confounded comparison.
set.seed(105)
alpha_grid <- c(0, 0.1, 0.5, 1, 2, 5, 10, 20, 50)
inner <- rsample::vfold_cv(data.frame(hn_train_X, y = hn_train_y), v = 5, strata = y)
logloss <- function(p, y_pos, eps = 1e-15) {
p <- pmin(pmax(p, eps), 1 - eps)
-mean(y_pos * log(p) + (1 - y_pos) * log(1 - p))
}
tune_alpha <- do.call(rbind, lapply(alpha_grid, function(a) {
sc <- vapply(inner$splits, function(s) {
A <- rsample::analysis(s); B <- rsample::assessment(s)
m <- e1071::naiveBayes(dplyr::select(A, -y), A$y, laplace = a)
p <- predict(m, dplyr::select(B, -y), type = "raw")[, "later_stage"]
cl <- predict(m, dplyr::select(B, -y))
c(logloss(p, as.integer(B$y == "later_stage")), mean(cl == B$y))
}, numeric(2))
data.frame(alpha = a, cv_logloss = mean(sc[1, ]), cv_accuracy = mean(sc[2, ]))
}))
tune_alpha |> mutate(across(where(is.numeric), \(z) round(z, 4)))#> alpha cv_logloss cv_accuracy
#> 1 0.0 1.0459 0.6350
#> 2 0.1 1.0419 0.6294
#> 3 0.5 1.0092 0.6293
#> 4 1.0 0.9853 0.6274
#> 5 2.0 0.9530 0.6481
#> 6 5.0 0.9302 0.6633
#> 7 10.0 1.0759 0.7108
#> 8 20.0 1.5717 0.7277
#> 9 50.0 2.6631 0.7297a_best <- tune_alpha$alpha[which.min(tune_alpha$cv_logloss)]
c(alpha_by_logloss = a_best,
alpha_by_accuracy = tune_alpha$alpha[which.max(tune_alpha$cv_accuracy)])#> alpha_by_logloss alpha_by_accuracy
#> 5 50tune_alpha |>
mutate(cv_error = 1 - cv_accuracy) |>
dplyr::select(alpha, `CV log-loss` = cv_logloss, `CV error` = cv_error) |>
pivot_longer(-alpha, names_to = "metric", values_to = "value") |>
ggplot(aes(alpha + 0.05, value, colour = metric)) +
geom_line(linewidth = 1) + geom_point(size = 2.2) +
scale_x_log10() +
facet_wrap(~ metric, scales = "free_y") +
scale_colour_manual(values = c("#D8433B", "#3B7DD8"), guide = "none") +
labs(title = "Tuning the smoothing parameter",
subtitle = "Log-loss responds to over-smoothing; accuracy is nearly flat because it ignores the probabilities",
x = expression(alpha + 0.05 ~ "(log scale)"), y = NULL) +
theme_dspa(10)
Log-loss has a clear minimum; accuracy is nearly flat and then rises at large \(\alpha\), because over-smoothing pushes every prediction to the majority class, whose accuracy is the no-information rate. Tuning on accuracy would select the degenerate model.
hn_nb <- e1071::naiveBayes(hn_train_X, hn_train_y, laplace = a_best)
hn_pred <- predict(hn_nb, hn_test_X)
hn_prob <- predict(hn_nb, hn_test_X, type = "raw")[, "later_stage"]
cm_hn <- confusionMatrix(data = hn_pred, reference = hn_test_y,
positive = "later_stage")
cm_hn#> Confusion Matrix and Statistics
#>
#> Reference
#> Prediction early_stage later_stage
#> early_stage 67 34
#> later_stage 22 10
#>
#> Accuracy : 0.579
#> 95% CI : (0.49, 0.664)
#> No Information Rate : 0.669
#> P-Value [Acc > NIR] : 0.988
#>
#> Kappa : -0.021
#>
#> Mcnemar's Test P-Value : 0.142
#>
#> Sensitivity : 0.2273
#> Specificity : 0.7528
#> Pos Pred Value : 0.3125
#> Neg Pred Value : 0.6634
#> Prevalence : 0.3308
#> Detection Rate : 0.0752
#> Detection Prevalence : 0.2406
#> Balanced Accuracy : 0.4900
#>
#> 'Positive' Class : later_stage
#> #> n TP FP FN
#> 133.0000 10.0000 22.0000 34.0000
#> TN accuracy NIR sensitivity
#> 67.0000 0.5789 0.6692 0.2273
#> specificity PPV NPV balanced_accuracy
#> 0.7528 0.3125 0.6634 0.4900
#> F1 kappa MCC
#> 0.2632 -0.0214 -0.0219cm_df <- as.data.frame(cm_hn$table)
names(cm_df) <- c("Predicted", "Reference", "Count")
ggplot(cm_df, aes(Reference, Predicted, fill = Count)) +
geom_tile(colour = "white", linewidth = 1) +
geom_text(aes(label = Count), size = 6, fontface = "bold") +
scale_fill_gradient(low = "#EAF2FB", high = "#3B7DD8", guide = "none") +
labs(title = "Naive Bayes confusion matrix (held-out data)",
subtitle = sprintf("Accuracy %.3f vs. no-information rate %.3f",
cm_hn$overall[["Accuracy"]], NIR_hn)) +
theme_dspa()
# --- Interactive equivalent ------------------------------------------------
m <- binary_metrics(hn_pred, hn_test_y, positive = "later_stage")
plot_ly(x = c("TN", "FN", "FP", "TP"),
y = c(m[["TN"]], m[["FN"]], m[["FP"]], m[["TP"]]),
type = "bar", color = c("TN", "FN", "FP", "TP")) |>
layout(title = "Confusion matrix counts",
legend = list(title = list(text = "<b>Cell</b>")),
xaxis = list(title = "Cell"), yaxis = list(title = "Count"))16.6 Comparing generative models on the same features
# Multinomial NB from counts (rather than binary presence)
counts_train <- as.matrix(hn_train_dtm); counts_test <- as.matrix(hn_test_dtm)
multinomial_nb <- function(Xtr, ytr, Xte, alpha = 1) {
cls <- levels(ytr)
logprior <- log(table(ytr) / length(ytr))
loglik <- sapply(cls, function(k) {
tot <- colSums(Xtr[ytr == k, , drop = FALSE])
log((tot + alpha) / (sum(tot) + alpha * ncol(Xtr)))
})
scores <- Xte %*% loglik + matrix(logprior, nrow(Xte), length(cls), byrow = TRUE)
factor(cls[max.col(scores)], levels = cls)
}
mn_pred <- multinomial_nb(counts_train, hn_train_y, counts_test, alpha = a_best)
# LDA on the leading principal components: reduces dimension AND decorrelates,
# which is what LDA's shared-covariance assumption needs.
pc <- prcomp(counts_train, center = TRUE, scale. = FALSE)
n_pc <- 25
Ztr <- data.frame(pc$x[, 1:n_pc], y = hn_train_y)
Zte <- data.frame(scale(counts_test, pc$center, FALSE) %*% pc$rotation[, 1:n_pc])
lda_hn <- MASS::lda(y ~ ., data = Ztr)
lda_pred <- predict(lda_hn, Zte)$class
data.frame(
model = c("Majority class", "Bernoulli NB", "Multinomial NB", paste0("LDA on ", n_pc, " PCs")),
accuracy = round(c(NIR_hn,
mean(hn_pred == hn_test_y),
mean(mn_pred == hn_test_y),
mean(lda_pred == hn_test_y)), 4),
balanced_accuracy = round(c(0.5,
binary_metrics(hn_pred, hn_test_y, "later_stage")[["balanced_accuracy"]],
binary_metrics(mn_pred, hn_test_y, "later_stage")[["balanced_accuracy"]],
binary_metrics(lda_pred, hn_test_y, "later_stage")[["balanced_accuracy"]]), 4))#> model accuracy balanced_accuracy
#> 1 Majority class 0.6692 0.500
#> 2 Bernoulli NB 0.5789 0.490
#> 3 Multinomial NB 0.6541 0.506
#> 4 LDA on 25 PCs 0.6692 0.500Note the change from a naive application of LDA. Fitting LDA directly on hundreds of 0/1 indicators violates its multivariate-normality and shared-covariance assumptions comprehensively, and with \(p>n\) the pooled covariance is singular. Projecting onto principal components first (Chapter 4) reduces dimension, decorrelates the features, and produces approximately continuous scores, which is the setting LDA was designed for.
Read the comparison against the baseline. Balanced accuracy is the metric to watch: it is 0.5 for the majority-class rule by construction, so any value near 0.5 means the model is not separating the classes regardless of what its raw accuracy says.
17 PART IV: TREES, RULES, AND ENSEMBLES
kNN and naive Bayes give no explicit rules. When the classification logic must be inspectable, a credit-scoring rubric, a clinical triage protocol, a regulatory audit, a decision tree states its reasoning as a sequence of conditions any reader can follow.
18 Recursive partitioning
A decision tree partitions feature space by axis-aligned splits, applied recursively, and predicts the majority class in each resulting region:
\[\hat f(x)=\sum_{m=1}^{M}\hat c_m\,\mathbb{1}\{x\in R_m\}, \qquad \hat c_m=\arg\max_k \hat p_{mk},\]
with \(\hat p_{mk}\) the proportion of class \(k\) among training points in region \(R_m\). Prediction is a walk from the root to a leaf.
Splitting stops when a node is pure, when the majority class exceeds a threshold, when no attribute remains, or when a size constraint binds.
library(partykit) # the maintained successor to `party`
data(iris)
iris_ctree <- partykit::ctree(Species ~ ., data = iris)
print(iris_ctree)#>
#> Model formula:
#> Species ~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width
#>
#> Fitted party:
#> [1] root
#> | [2] Petal.Length <= 1.9: setosa (n = 50, err = 0%)
#> | [3] Petal.Length > 1.9
#> | | [4] Petal.Width <= 1.7
#> | | | [5] Petal.Length <= 4.8: versicolor (n = 46, err = 2%)
#> | | | [6] Petal.Length > 4.8: versicolor (n = 8, err = 50%)
#> | | [7] Petal.Width > 1.7: virginica (n = 46, err = 2%)
#>
#> Number of inner nodes: 3
#> Number of terminal nodes: 4
ctree builds a conditional-inference
tree: at each node it runs a permutation test of independence between
each feature and the response, splitting on the most significant one
only if it survives a multiplicity-adjusted threshold. That gives a
principled stopping rule and removes the bias toward high-cardinality
features that afflicts information-gain splitting (§5.19).
library(rpart); library(rpart.plot)
iris_rpart <- rpart(Species ~ ., data = iris, method = "class")
rpart.plot(iris_rpart, type = 4, extra = 104, box.palette = "BuGn",
main = "CART classification of iris taxa")
# A Sankey rendering of the same tree: flow width is the number of cases
tf <- iris_rpart$frame
is_leaf <- tf$var == "<leaf>"
ylev <- attr(iris_rpart, "ylevels")
node_lab <- character(nrow(tf))
node_lab[is_leaf] <- ylev[tf$yval][is_leaf]
node_lab[!is_leaf] <- labels(iris_rpart)[-1][!is_leaf[-length(is_leaf)]]
node_ids <- as.numeric(rownames(tf))
parent_of <- function(id) if (id == 1) NA_integer_ else floor(id / 2)
src <- match(vapply(node_ids, parent_of, numeric(1)), node_ids) - 1L
tgt <- seq_along(node_ids) - 1L
keep <- !is.na(src)
plot_ly(type = "sankey", orientation = "h",
node = list(label = node_lab, pad = 15, thickness = 28,
line = list(color = "black", width = 1)),
link = list(source = src[keep], target = tgt[keep],
value = tf$n[keep])) |>
layout(title = "Iris CART decision tree as a flow diagram",
font = list(size = 12))19 Impurity measures
A split is good if the children are purer than the parent. Three measures quantify impurity for a node with class proportions \(p_1,\dots,p_K\):
\[ \begin{aligned} \textbf{Entropy: }\quad & H(p)=-\sum_{k=1}^{K}p_k\log_2 p_k &&\in\big[0,\ \log_2 K\big]\\[2mm] \textbf{Gini: }\quad & G(p)=\sum_{k=1}^{K}p_k(1-p_k)=1-\sum_k p_k^2 &&\in\big[0,\ 1-\tfrac1K\big]\\[2mm] \textbf{Misclassification: }\quad & E(p)=1-\max_k p_k &&\in\big[0,\ 1-\tfrac1K\big] \end{aligned} \]
All three are zero for a pure node and maximal at \(p_k=1/K\). Note that \(\lim_{p\to0}p\log p=0\) by L’Hôpital, so the entropy of a pure node is \(-1\cdot\log_2 1 - 0 = 0\) with no undefined term. The maximum \(\log_2 K\) grows with \(K\), it is not bounded, and the base of the logarithm only rescales, since \(\log_b x=\log_2 x/\log_2 b\).
pp <- seq(0.001, 0.999, length.out = 400)
imp <- data.frame(
p = rep(pp, 3),
value = c(-pp * log2(pp) - (1 - pp) * log2(1 - pp),
2 * pp * (1 - pp),
pmin(pp, 1 - pp)),
measure = rep(c("Entropy", "Gini (x2)", "Misclassification"), each = length(pp)))
ggplot(imp, aes(p, value, colour = measure)) +
geom_line(linewidth = 1) +
geom_point(data = data.frame(p = 0.8, value = -0.8 * log2(0.8) - 0.2 * log2(0.2),
measure = "Entropy"),
size = 4, colour = "darkgreen") +
ggplot2::annotate("text", x = 0.8, y = 0.79, label = "(0.8, 0.722)", size = 3.2) +
scale_colour_manual(values = c(Entropy = "#3B7DD8", `Gini (x2)` = "#D8433B",
Misclassification = "#7FB069")) +
labs(title = "Impurity measures for a two-class node",
subtitle = "Gini is doubled to share the scale. All peak at p = 0.5 and vanish at 0 and 1",
x = "Proportion of class 1", y = "Impurity", colour = NULL) +
theme_dspa()
# --- Interactive equivalent ------------------------------------------------
labs_i <- c("High structure / low entropy", "Low structure / high entropy",
"High structure / low entropy")
plot_ly() |>
add_lines(x = pp, y = -pp * log2(pp) - (1 - pp) * log2(1 - pp),
name = "Entropy", line = list(width = 4)) |>
add_segments(x = 0.8, xend = 0.8, y = 0.001, yend = 0.7219281, showlegend = FALSE) |>
add_markers(x = 0.8, y = 0.7219281, name = "(0.8, 0.722)",
marker = list(size = 18, color = "green",
line = list(color = "yellow", width = 2))) |>
layout(title = "Binary class entropy",
xaxis = list(title = "Proportion of class 1 observations"),
yaxis = list(title = "Entropy", range = c(-0.1, 1.1)),
annotations = list(text = labs_i, x = c(0.03, 0.5, 0.97),
y = c(0.5, 0.7, 0.5), textangle = 90,
showarrow = FALSE))Common misconception: “lower entropy means more information.” These are two different quantities and the sentence conflates them. In the Shannon sense, high entropy = high information content: a fair coin carries one bit, a two-headed coin carries none. What a good split maximizes is the information gain, the reduction in entropy, which is large precisely when the children are low-entropy. Low entropy in a node means the node is predictable, i.e. carries little information and is therefore useful for classification.
19.1 Why misclassification error is a poor splitting criterion
Entropy and Gini are strictly concave; misclassification error is piecewise linear. Strict concavity guarantees that any split producing children with different class proportions strictly reduces the weighted impurity. Misclassification error can give exactly zero improvement for a split that is clearly informative, so a greedy search using it stalls.
# Parent: 400 cases, 50/50. A split that is obviously useful:
# left = 200 cases, 75/25; right = 200 cases, 25/75
parent <- c(200, 200); left <- c(150, 50); right <- c(50, 150)
imp_fns <- list(
Entropy = \(v) { p <- v / sum(v); p <- p[p > 0]; -sum(p * log2(p)) },
Gini = \(v) { p <- v / sum(v); 1 - sum(p^2) },
Misclassification = \(v) { p <- v / sum(v); 1 - max(p) })
sapply(imp_fns, function(f) {
wl <- sum(left) / sum(parent); wr <- sum(right) / sum(parent)
c(parent = f(parent), children = wl * f(left) + wr * f(right),
gain = f(parent) - (wl * f(left) + wr * f(right)))
}) |> round(4)#> Entropy Gini Misclassification
#> parent 1.0000 0.500 0.50
#> children 0.8113 0.375 0.25
#> gain 0.1887 0.125 0.25Entropy and Gini both register a substantial gain. Misclassification error registers exactly zero, the parent misclassifies 50%, and so do the weighted children, even though the split has separated the classes considerably.
The standard practice follows: split on Gini or entropy, prune on misclassification error (or on cross-validated error). Splitting needs a criterion sensitive to probability changes; pruning wants the quantity you actually care about.
20 Information gain and gain ratio
Splitting node \(S\) on feature \(F\) into segments \(S_1,\dots,S_m\) produces the weighted child impurity
\[H(S\mid F)=\sum_{i=1}^{m}w_i\,H(S_i) =\sum_{i=1}^{m}w_i\left(-\sum_{j=1}^{K}p_{j\mid i}\log_2 p_{j\mid i}\right),\]
where \(w_i=|S_i|/|S|\) and \(p_{j\mid i}\) is the proportion of class \(j\) within segment \(i\). The information gain is
\[\mathrm{Gain}(F)=H(S)-H(S\mid F)\ \ge 0 .\]
It equals \(H(S)\) when the children are perfectly pure, and zero when the split is uninformative.
Information gain has a systematic bias toward high-cardinality features. A patient ID splits every node into singletons, giving maximum gain and zero predictive value. C4.5 corrects this with the gain ratio, normalizing by the entropy of the split itself:
\[\mathrm{SplitInfo}(F)=-\sum_{i=1}^{m}w_i\log_2 w_i, \qquad \mathrm{GainRatio}(F)=\frac{\mathrm{Gain}(F)}{\mathrm{SplitInfo}(F)} .\]
A feature with many levels has large \(\mathrm{SplitInfo}\), which penalizes it.
set.seed(107)
n_g <- 600
y_g <- factor(sample(c("A", "B"), n_g, TRUE))
useful <- factor(ifelse(runif(n_g) < ifelse(y_g == "A", 0.8, 0.25), "hi", "lo"))
id_col <- factor(seq_len(n_g)) # a unique ID: zero signal
ent <- function(v) { p <- table(v) / length(v); p <- p[p > 0]; -sum(p * log2(p)) }
cond_ent <- function(x, y) {
sum(vapply(split(y, x), \(g) length(g) / length(y) * ent(g), numeric(1)))
}
split_info <- function(x) { w <- table(x) / length(x); -sum(w * log2(w)) }
rbind(
`Useful binary feature` = c(gain = ent(y_g) - cond_ent(useful, y_g),
split_info = split_info(useful),
gain_ratio = (ent(y_g) - cond_ent(useful, y_g)) / split_info(useful)),
`Patient ID (no signal)` = c(gain = ent(y_g) - cond_ent(id_col, y_g),
split_info = split_info(id_col),
gain_ratio = (ent(y_g) - cond_ent(id_col, y_g)) / split_info(id_col))
) |> round(4)#> gain split_info gain_ratio
#> Useful binary feature 0.2259 0.9937 0.2274
#> Patient ID (no signal) 0.9990 9.2288 0.1083The ID achieves the maximum possible information gain, it perfectly “predicts” every training label, while the useful feature scores far lower. The gain ratio reverses the ranking. This is not a hypothetical: administrative identifiers, timestamps, and record numbers routinely top variable-importance lists for exactly this reason (§5.23.5).
21 Pruning
A tree grown until every leaf is pure has zero training error and generalizes badly. Two controls:
Pre-pruning (early stopping), minsplit,
minbucket, maxdepth, or a significance
threshold. Cheap, but it can stop before a split that would have enabled
a valuable one below it (the horizon effect).
Post-pruning, grow a large tree, then collapse branches that do not earn their complexity. CART uses cost-complexity pruning: for \(\alpha\ge0\), minimize
\[C_\alpha(T)=\sum_{m=1}^{|T|}n_m\,E(T_m)+\alpha|T| .\]
As \(\alpha\) increases, the minimizing subtree shrinks along a nested sequence, so cross-validating \(\alpha\) is a one-dimensional search. C5.0 uses a confidence-based post-pruning rule with a similar effect.
The 1-SE rule picks the smallest tree whose cross-validated error is within one standard error of the minimum, trading a little accuracy for a model that is smaller, more stable, and easier to defend.
22 Case study: quality of life and chronic disease
22.1 The data
Case
Study 06, Case06_QoL_Symptom_ChronicIllness.csv: 41
variables on ~2,200 patients.
- Charlson Comorbidity Index: 0–10; higher indicates more comorbidity.
- Chronic Disease Score: severity summary from prescription complexity; \(-9\) codes a missing value.
qol <- dspa_read("https://umich.instructure.com/files/481332/download?download_frd=1",
"Case06_QoL_Symptom_ChronicIllness.csv")
dim(qol)#> [1] 2356 41#>
#> 1 2 3 4 5 6
#> 44 213 801 900 263 135#> missing_coded_as_neg9
#> 142#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 0.00 0.88 1.40 1.50 1.97 4.7622.2 Constructing the outcome, and its consequences
# The threshold is COMPUTED from the data, not hard-coded
cut_point <- switch(CDS_SPLIT,
median = median(qol$CHRONICDISEASESCORE),
mean = mean(qol$CHRONICDISEASESCORE))
c(split_rule = CDS_SPLIT, cut_point = round(cut_point, 4),
mean = round(mean(qol$CHRONICDISEASESCORE), 4),
median = round(median(qol$CHRONICDISEASESCORE), 4))#> split_rule cut_point mean median
#> "median" "1.395" "1.4975" "1.395"qol$cd <- factor(qol$CHRONICDISEASESCORE > cut_point,
levels = c(FALSE, TRUE),
labels = c("minor_disease", "severe_disease"))
round(prop.table(table(qol$cd)), 4)#>
#> minor_disease severe_disease
#> 0.5 0.5# What the alternative would have given
data.frame(
rule = c("median", "mean"),
cut = round(c(median(qol$CHRONICDISEASESCORE), mean(qol$CHRONICDISEASESCORE)), 3),
pct_severe = round(c(mean(qol$CHRONICDISEASESCORE > median(qol$CHRONICDISEASESCORE)),
mean(qol$CHRONICDISEASESCORE > mean(qol$CHRONICDISEASESCORE))), 4),
no_information_rate = round(c(
max(mean(qol$CHRONICDISEASESCORE > median(qol$CHRONICDISEASESCORE)),
1 - mean(qol$CHRONICDISEASESCORE > median(qol$CHRONICDISEASESCORE))),
max(mean(qol$CHRONICDISEASESCORE > mean(qol$CHRONICDISEASESCORE)),
1 - mean(qol$CHRONICDISEASESCORE > mean(qol$CHRONICDISEASESCORE)))), 4))#> rule cut pct_severe no_information_rate
#> 1 median 1.395 0.500 0.500
#> 2 mean 1.497 0.477 0.523The choice of cut-point is a modelling decision with a measurable consequence. The median gives balanced classes and a no-information rate near 0.5, so accuracy is directly interpretable. The mean gives an imbalanced split and a higher NIR, so a model can look better while learning less. Neither is “correct”, but the choice must be made deliberately, and a clinically motivated threshold would be better than either.
22.3 Excluding leakage-prone features
# CHRONICDISEASESCORE defines the outcome -> target leakage.
# INTERVIEWDATE and ID are administrative artifacts -> leakage by proxy:
# they encode collection waves and protocol changes, not clinical signal.
drop_cols <- c("ID", "CHRONICDISEASESCORE", "INTERVIEWDATE")
qol_model_df <- qol |> dplyr::select(-any_of(drop_cols))
c(dropped = paste(intersect(drop_cols, names(qol)), collapse = ", "),
predictors_remaining = ncol(qol_model_df) - 1)#> dropped
#> "ID, CHRONICDISEASESCORE, INTERVIEWDATE"
#> predictors_remaining
#> "38"Common misconception: “if a variable predicts well, keep it.” A rule keyed to
INTERVIEWDATEwill score well and generalize to nothing — it has detected a data-collection wave, not a disease mechanism. Ask of every strong predictor: would this value be available, and mean the same thing, at the moment a prediction is actually needed? Identifiers, timestamps, record numbers, and post-outcome codes almost always fail that test. This is why variable-importance rankings must be read with the data-generating process in mind, never on their own.
set.seed(1234)
split_qol <- rsample::initial_split(qol_model_df, prop = 0.8, strata = cd)
qol_train <- rsample::training(split_qol)
qol_test <- rsample::testing(split_qol)
round(rbind(train = prop.table(table(qol_train$cd)),
test = prop.table(table(qol_test$cd))), 4)#> minor_disease severe_disease
#> train 0.5 0.5
#> test 0.5 0.5NIR_qol <- max(prop.table(table(qol_test$cd)))
c(n_train = nrow(qol_train), n_test = nrow(qol_test),
no_information_rate = round(NIR_qol, 4))#> n_train n_test no_information_rate
#> 1770.0 444.0 0.522.4 C5.0
#>
#> Call:
#> C5.0.formula(formula = cd ~ ., data = qol_train)
#>
#> Classification Tree
#> Number of samples: 1770
#> Number of predictors: 38
#>
#> Tree size: 22
#>
#> Non-standard options: attempt to group attributesc50_pred <- predict(qol_c50, qol_test)
cm_c50 <- confusionMatrix(data = c50_pred, reference = qol_test$cd,
positive = "severe_disease")
cm_c50$overall[1:6] |> round(4)#> Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull
#> 0.6126 0.2252 0.5655 0.6582 0.5000
#> AccuracyPValue
#> 0.0000#> n TP FP FN
#> 444.0000 177.0000 127.0000 45.0000
#> TN accuracy NIR sensitivity
#> 95.0000 0.6126 0.5000 0.7973
#> specificity PPV NPV balanced_accuracy
#> 0.4279 0.5822 0.6786 0.6126
#> F1 kappa MCC
#> 0.6730 0.2252 0.2424Boosting grows a sequence of trees, each weighting the cases its predecessors got wrong:
set.seed(1234)
qol_boost <- C5.0(cd ~ ., data = qol_train, trials = 10)
boost_pred <- predict(qol_boost, qol_test)
data.frame(
model = c("Majority class", "C5.0 single tree", "C5.0 boosted (10 trials)"),
accuracy = round(c(NIR_qol, mean(c50_pred == qol_test$cd),
mean(boost_pred == qol_test$cd)), 4),
kappa = round(c(0,
binary_metrics(c50_pred, qol_test$cd, "severe_disease")[["kappa"]],
binary_metrics(boost_pred, qol_test$cd, "severe_disease")[["kappa"]]), 4))#> model accuracy kappa
#> 1 Majority class 0.5000 0.0000
#> 2 C5.0 single tree 0.6126 0.2252
#> 3 C5.0 boosted (10 trials) 0.6464 0.292822.5 Cost-sensitive learning
Missing a severe case is usually costlier than a false alarm. A cost matrix makes that asymmetry explicit rather than leaving it implicit in the 0.5 threshold.
# Rows = predicted, columns = actual. Penalize predicting "minor" when actually
# "severe" (a missed severe case) four times as heavily as the reverse.
error_cost <- matrix(c(0, 1, 4, 0), nrow = 2,
dimnames = list(Predicted = levels(qol_train$cd),
Actual = levels(qol_train$cd)))
error_cost#> Actual
#> Predicted minor_disease severe_disease
#> minor_disease 0 4
#> severe_disease 1 0set.seed(1234)
qol_costly <- C5.0(cd ~ ., data = qol_train, costs = error_cost)
cost_pred <- predict(qol_costly, qol_test)
comparison_cost <- rbind(
`Default (equal costs)` = binary_metrics(c50_pred, qol_test$cd, "severe_disease"),
`Cost matrix (4:1)` = binary_metrics(cost_pred, qol_test$cd, "severe_disease"))
round(comparison_cost[, c("accuracy", "sensitivity", "specificity",
"balanced_accuracy", "FN", "FP")], 4)#> accuracy sensitivity specificity balanced_accuracy FN FP
#> Default (equal costs) 0.6126 0.7973 0.4279 0.6126 45 127
#> Cost matrix (4:1) 0.5991 0.9505 0.2477 0.5991 11 167Overall accuracy falls while false negatives drop sharply, precisely the intended trade. Whether it is a good trade depends on the real cost ratio, which is a clinical question, not a statistical one. Note that a cost matrix and a shifted decision threshold achieve the same thing by different routes; the cost matrix additionally influences which splits the tree chooses.
22.6 Cost-complexity pruning with rpart
set.seed(1234)
qol_full <- rpart(cd ~ ., data = qol_train, method = "class",
control = rpart.control(cp = 0, minsplit = 2, xval = 10))
cp_tab <- as.data.frame(qol_full$cptable)
cp_min <- cp_tab$CP[which.min(cp_tab$xerror)]
se_thr <- min(cp_tab$xerror) + cp_tab$xstd[which.min(cp_tab$xerror)]
cp_1se <- cp_tab$CP[which(cp_tab$xerror <= se_thr)[1]]
c(cp_min_error = signif(cp_min, 4), cp_1se_rule = signif(cp_1se, 4),
size_min = cp_tab$nsplit[which.min(cp_tab$xerror)] + 1,
size_1se = cp_tab$nsplit[which(cp_tab$xerror <= se_thr)[1]] + 1)#> cp_min_error cp_1se_rule size_min size_1se
#> 0.01469 0.01469 2.00000 2.00000ggplot(cp_tab, aes(nsplit + 1, xerror)) +
geom_ribbon(aes(ymin = xerror - xstd, ymax = xerror + xstd), fill = "grey88") +
geom_line(linewidth = 0.9, colour = "steelblue") +
geom_hline(yintercept = se_thr, linetype = "dashed", colour = "firebrick") +
geom_point(data = cp_tab[which.min(cp_tab$xerror), ], colour = "firebrick", size = 3) +
scale_x_log10() +
labs(title = "Cost-complexity pruning path",
subtitle = "Dashed line: 1-SE threshold. Pick the smallest tree beneath it",
x = "Tree size (leaves, log scale)", y = "Cross-validated relative error") +
theme_dspa()
# --- Interactive equivalent ------------------------------------------------
plot_ly(cp_tab, x = ~factor(signif(CP, 3)), y = ~xerror, type = "scatter",
mode = "lines+markers", name = "CV error",
error_y = ~list(array = xstd, color = "gray")) |>
layout(title = "Complexity parameter vs. cross-validated error rate",
xaxis = list(title = "cp"), yaxis = list(title = "xerror"))qol_pruned <- prune(qol_full, cp = cp_1se)
rpart.plot(qol_pruned, type = 4, extra = 104, box.palette = "BuGn",
main = "Pruned tree (1-SE rule)")
pruned_pred <- predict(qol_pruned, qol_test, type = "class")
data.frame(
tree = c("Unpruned (cp = 0)", "Pruned at min CV error", "Pruned by 1-SE rule"),
leaves = c(sum(qol_full$frame$var == "<leaf>"),
sum(prune(qol_full, cp = cp_min)$frame$var == "<leaf>"),
sum(qol_pruned$frame$var == "<leaf>")),
test_accuracy = round(c(
mean(predict(qol_full, qol_test, type = "class") == qol_test$cd),
mean(predict(prune(qol_full, cp = cp_min), qol_test, type = "class") == qol_test$cd),
mean(pruned_pred == qol_test$cd)), 4))#> tree leaves test_accuracy
#> 1 Unpruned (cp = 0) 394 0.5766
#> 2 Pruned at min CV error 2 0.6149
#> 3 Pruned by 1-SE rule 2 0.6149The unpruned tree has hundreds of leaves and no better test accuracy. Pruning buys interpretability at no predictive cost, which is the usual outcome and the reason to do it.
22.7 The tree as a surface
tree_sim <- rpart(y ~ x1 + x2, data = sim_tr, method = "class",
control = rpart.control(cp = 0.005))
z_tree <- matrix(predict(tree_sim, gd2, type = "prob")[, "B"], length(gx))
plot_ly() |>
add_surface(x = gy, y = gx, z = z_tree, showscale = FALSE, opacity = 0.95,
colorscale = "RdBu", reversescale = TRUE) |>
add_surface(x = gy, y = gx, z = matrix(post2(gd2$x1, gd2$x2), length(gx)) + 1.3,
showscale = FALSE, opacity = 0.9, colorscale = "RdBu",
reversescale = TRUE) |>
layout(title = "Decision tree posterior (bottom) vs. the true Bayes posterior (top)",
scene = list(xaxis = list(title = "x2"), yaxis = list(title = "x1"),
zaxis = list(title = "P(Y = B | x) + offset")))The tree surface is piecewise constant with axis-aligned cliffs, its defining feature. That structure makes trees invariant to monotone transformations of any feature and able to represent interactions automatically. It also makes them poor at representing a smooth diagonal boundary, which they can only approximate with a staircase.
23 Classification rules
Separate-and-conquer repeatedly finds a rule covering some cases, removes them, and repeats on the remainder. Unlike a tree, where every node inherits all ancestors’ conditions, each rule is independent and readable on its own.
23.1 OneR
ZeroR predicts the mode for everything, the
no-information baseline (§5.3.3).
OneR improves on it by using a single feature: for each
candidate, bin its values, assign each bin its majority class, compute
the error, and keep the feature with the lowest.
Trivially simple, and a valuable calibration point: if an elaborate model cannot beat a single rule on one variable, its complexity is not earning its keep.
library(OneR) # pure R -- no Java dependency
set.seed(1234)
# Fit on TRAINING data only, then evaluate on untouched test data
qol_1r <- OneR::OneR(cd ~ ., data = OneR::optbin(qol_train), verbose = TRUE)#>
#> Attribute Accuracy
#> 1 * CHARLSONSCORE 65.82%
#> 2 QOL_Q_06 54.86%
#> 3 MSA_Q_06 54.63%
#> 4 QOL_Q_03 54.24%
#> 4 QOL_Q_07 54.24%
#> 6 MSA_Q_02 54.18%
#> 7 PH2_Q_01 54.07%
#> 8 QOL_Q_05 54.01%
#> 8 MSA_Q_04 54.01%
#> 10 QOL_Q_08 53.95%
#> 11 QOL_Q_02 53.73%
#> 12 QOL_Q_01 53.67%
#> 13 QOL_Q_04 53.62%
#> 14 MSA_Q_10 53.45%
#> 15 PH2_Q_02 53.39%
#> 16 MSA_Q_09 53.33%
#> 17 MSA_Q_01 53.16%
#> 18 QOL_Q_09 52.94%
#> 19 RACE_ETHNICITY 52.82%
#> 20 MSA_Q_15 52.77%
#> 21 MSA_Q_17 52.71%
#> 22 MSA_Q_07 52.6%
#> 22 MSA_Q_08 52.6%
#> 24 AGE 52.49%
#> 25 MSA_Q_05 52.32%
#> 26 QOL_Q_10 52.15%
#> 27 MSA_Q_11 52.03%
#> 28 MSA_Q_14 51.98%
#> 29 MSA_Q_03 51.86%
#> 30 MSA_Q_12 51.81%
#> 31 MSA_Q_16 51.69%
#> 32 LANGUAGE 51.36%
#> 33 TOS_Q_01 51.19%
#> 34 MSA_Q_13 50.96%
#> 35 SEX 50.28%
#> 36 TOS_Q_02 50.23%
#> 36 TOS_Q_03 50.23%
#> 36 TOS_Q_04 50.23%
#> ---
#> Chosen attribute due to accuracy
#> and ties method (if applicable): '*'#>
#> Call:
#> OneR.formula(formula = cd ~ ., data = OneR::optbin(qol_train),
#> verbose = TRUE)
#>
#> Rules:
#> If CHARLSONSCORE = (-9.02,0.909] then cd = minor_disease
#> If CHARLSONSCORE = (0.909,10] then cd = severe_disease
#>
#> Accuracy:
#> 1165 of 1770 instances classified correctly (65.82%)oner_pred <- predict(qol_1r, qol_test)
c(train_accuracy = round(mean(predict(qol_1r, qol_train) == qol_train$cd), 4),
test_accuracy = round(mean(oner_pred == qol_test$cd, na.rm = TRUE), 4),
no_information_rate = round(NIR_qol, 4))#> train_accuracy test_accuracy no_information_rate
#> 0.6582 0.6149 0.5000The training and test sets must be disjoint. Fitting a rule learner on the full dataset and then “validating” on a subset of those same rows reports the resubstitution error twice. The tell is that internal and external accuracy agree exactly, which is not reassuring evidence of generalization but proof that no generalization was tested.
23.2 RIPPER and rule induction
RIPPER (Repeated Incremental Pruning to Produce Error Reduction) runs three phases:
- Grow, add conditions greedily until the rule stops improving on the growing set;
- Prune, remove trailing conditions that raise error on a held-out pruning set;
- Optimize, revisit each rule, considering replacements and revisions, using an MDL criterion to decide when to stop.
C5.0 can also emit a rule set directly, converting each root-to-leaf path into a rule and then simplifying globally.
#>
#> Call:
#> C5.0.formula(formula = cd ~ ., data = qol_train, rules = TRUE)
#>
#> Rule-Based Model
#> Number of samples: 1770
#> Number of predictors: 38
#>
#> Number of Rules: 12
#>
#> Non-standard options: attempt to group attributesrule_pred <- predict(qol_rules, qol_test)
c(test_accuracy = round(mean(rule_pred == qol_test$cd), 4),
n_rules = qol_rules$size)#> test_accuracy n_rules
#> 0.6149 12.0000# --- The Weka RIPPER implementation, for readers with a working rJava ------
library(RWeka)
set.seed(1234)
qol_jrip <- JRip(cd ~ ., data = qol_train) # TRAINING data only
qol_jrip
summary(qol_jrip)
jrip_pred <- predict(qol_jrip, qol_test)
caret::confusionMatrix(data = jrip_pred, reference = qol_test$cd,
positive = "severe_disease")24 Ensembles
A single tree is high-variance: perturb the data slightly and the top split can change, restructuring everything below it. Ensembles exploit that instability.
24.1 Why averaging works
Average \(B\) predictors each with variance \(\sigma^2\) and pairwise correlation \(\rho\). The variance of the mean is
\[\boxed{\;\operatorname{Var}\!\left(\frac1B\sum_{b=1}^{B}\hat f_b\right) =\rho\sigma^2+\frac{1-\rho}{B}\sigma^2\;}\]
Read both terms. The second vanishes as \(B\to\infty\), more trees never hurt, they only stop helping. The first does not: it is a floor set by how correlated the trees are. Reducing \(\rho\) is therefore the only way to push the variance lower, and it is exactly what random forests do.
Bagging fits each tree on a bootstrap resample,
which decorrelates them somewhat. Random forests add a
second source: at every split, only a random subset of \(m\) features (mtry) is
considered, so a dominant predictor cannot appear at the top of every
tree. Extremely randomized trees go further, choosing
split points at random rather than optimizing them.
set.seed(113)
B_grid <- c(1, 2, 5, 10, 25, 50, 100, 250)
n_rep <- 60
ens_var <- function(B, mtry) {
preds <- replicate(n_rep, {
idx <- sample(nrow(sim_tr), replace = TRUE)
f <- ranger::ranger(y ~ ., data = sim_tr[idx, ], num.trees = B,
mtry = mtry, probability = TRUE, num.threads = 1)
predict(f, sim_te)$predictions[, "B"]
})
mean(apply(preds, 1, var)) # variance across resamples, per test point
}
ev <- data.frame(
B = B_grid,
`mtry = 2 (all features)` = vapply(B_grid, ens_var, numeric(1), mtry = 2),
`mtry = 1 (decorrelated)` = vapply(B_grid, ens_var, numeric(1), mtry = 1),
check.names = FALSE)
round(ev, 5)#> B mtry = 2 (all features) mtry = 1 (decorrelated)
#> 1 1 0.05251 0.05954
#> 2 2 0.03229 0.03290
#> 3 5 0.01865 0.01756
#> 4 10 0.01490 0.01185
#> 5 25 0.01206 0.00863
#> 6 50 0.01172 0.00802
#> 7 100 0.01090 0.00726
#> 8 250 0.01057 0.00719ev |> pivot_longer(-B, names_to = "setting", values_to = "variance") |>
ggplot(aes(B, variance, colour = setting)) +
geom_line(linewidth = 1) + geom_point(size = 2.2) +
scale_x_log10() + scale_y_log10() +
scale_colour_manual(values = c("#D8433B", "#3B7DD8")) +
labs(title = "Prediction variance falls with B, then flattens at a floor",
subtitle = "The floor is rho*sigma^2 -- lower mtry decorrelates the trees and lowers it",
x = "Number of trees B (log scale)", y = "Average prediction variance",
colour = NULL) +
theme_dspa()
Both curves fall like \(1/B\) and
then flatten. The lower-mtry forest flattens
lower, because its trees are less correlated, the \(\rho\sigma^2\) term made visible.
24.2 Random forests in practice
library(ranger)
set.seed(12)
qol_rf <- ranger(cd ~ ., data = qol_train, num.trees = 500,
importance = "permutation", probability = FALSE,
respect.unordered.factors = "order")
qol_rf#> Ranger result
#>
#> Call:
#> ranger(cd ~ ., data = qol_train, num.trees = 500, importance = "permutation", probability = FALSE, respect.unordered.factors = "order")
#>
#> Type: Classification
#> Number of trees: 500
#> Sample size: 1770
#> Number of independent variables: 38
#> Mtry: 6
#> Target node size: 1
#> Variable importance mode: permutation
#> Splitrule: gini
#> OOB prediction error: 36.44 %rf_pred <- predict(qol_rf, qol_test)$predictions
cm_rf <- confusionMatrix(data = rf_pred, reference = qol_test$cd,
positive = "severe_disease")
cm_rf$overall[1:6] |> round(4)#> Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull
#> 0.6216 0.2432 0.5747 0.6669 0.5000
#> AccuracyPValue
#> 0.0000The out-of-bag error is a free cross-validation estimate: each bootstrap sample omits about \(1-(1-1/n)^n\to e^{-1}\approx36.8\%\) of the data, so every observation is out-of-bag for roughly a third of the trees and can be predicted by them.
oob_curve <- function(mtry) {
vapply(c(10, 25, 50, 100, 200, 350, 500), function(B) {
set.seed(12)
ranger(cd ~ ., data = qol_train, num.trees = B, mtry = mtry,
num.threads = 1)$prediction.error
}, numeric(1))
}
p_qol <- ncol(qol_train) - 1
oob_df <- data.frame(
B = c(10, 25, 50, 100, 200, 350, 500),
`mtry = sqrt(p)` = oob_curve(floor(sqrt(p_qol))),
`mtry = p/3` = oob_curve(max(1, floor(p_qol / 3))),
`mtry = p (bagging)` = oob_curve(p_qol),
check.names = FALSE)
oob_df |> pivot_longer(-B, names_to = "mtry", values_to = "oob") |>
ggplot(aes(B, oob, colour = mtry)) +
geom_line(linewidth = 1) + geom_point(size = 2) +
scale_colour_brewer(palette = "Set1") +
labs(title = "Out-of-bag error against the number of trees",
subtitle = "Curves flatten: adding trees stops helping but never hurts",
x = "Number of trees", y = "OOB error", colour = NULL) +
theme_dspa()
# --- Interactive equivalent ------------------------------------------------
plot_ly(oob_df, x = ~B, y = ~`mtry = sqrt(p)`, type = "scatter", mode = "lines",
name = "mtry = sqrt(p)", line = list(width = 2)) |>
add_trace(y = ~`mtry = p/3`, name = "mtry = p/3") |>
add_trace(y = ~`mtry = p (bagging)`, name = "mtry = p (bagging)") |>
layout(title = "Out-of-bag error rates for three forests",
xaxis = list(title = "Number of trees"),
yaxis = list(title = "OOB error"))24.3 Variable importance, read carefully
Impurity importance (mean decrease in Gini) sums the improvement each feature produces across all splits. It is fast and biased toward high-cardinality and continuous features, the same bias that afflicts information gain (§5.19).
Permutation importance measures the OOB accuracy drop when a feature’s values are shuffled. It is slower, model-agnostic, and far less biased. Prefer it, and note that both are unreliable when features are strongly correlated, because importance is split arbitrarily among the correlated group.
imp_perm <- sort(ranger::importance(qol_rf), decreasing = TRUE)[1:15]
ggplot(data.frame(var = factor(names(imp_perm), levels = rev(names(imp_perm))),
imp = as.numeric(imp_perm)),
aes(imp, var)) +
geom_col(fill = "steelblue") +
labs(title = "Permutation variable importance (top 15)",
subtitle = "OOB accuracy loss when each feature is shuffled",
x = "Mean decrease in accuracy", y = NULL) +
theme_dspa(10)
# --- Interactive equivalent (both importance types) ------------------------
plot_ly(x = ~as.numeric(imp_perm),
y = ~reorder(names(imp_perm), as.numeric(imp_perm)),
type = "bar", name = "Permutation") |>
layout(title = "RF variable importance (permutation)",
xaxis = list(title = "Mean decrease in accuracy"),
yaxis = list(title = "Variable"))24.4 Partial dependence
A partial dependence plot marginalizes over all other features:
\[\mathrm{PD}_j(v)=\frac1n\sum_{i=1}^{n}\hat f\big(x_{i,-j},\ X_j=v\big).\]
It shows the average modelled effect of \(X_j\). Its blind spot: by averaging, it can conceal opposite effects in different subgroups. Individual conditional expectation curves plot one line per observation and reveal that heterogeneity.
top_vars <- names(imp_perm)[1:4]
pdp_df <- do.call(rbind, lapply(top_vars, function(v) {
vals <- sort(unique(qol_train[[v]]))
if (length(vals) > 12) vals <- quantile(qol_train[[v]], seq(0, 1, length.out = 12))
do.call(rbind, lapply(vals, function(z) {
d <- qol_train; d[[v]] <- z
p <- predict(qol_rf, d, predict.all = FALSE)$predictions
data.frame(variable = v, value = z, pd = mean(p == "severe_disease"))
}))
}))
ggplot(pdp_df, aes(value, pd)) +
geom_line(linewidth = 1, colour = "steelblue") +
geom_point(size = 1.6) +
facet_wrap(~ variable, scales = "free_x") +
labs(title = "Partial dependence of severe-disease probability",
subtitle = "Averaged over all other features; subgroup heterogeneity is hidden by construction",
x = NULL, y = "P(severe disease)") +
theme_dspa(10)
# --- Interactive equivalent ------------------------------------------------
plot_ly(pdp_df, x = ~value, y = ~pd, color = ~variable,
type = "scatter", mode = "lines+markers") |>
layout(title = "Partial dependence on the top salient features",
xaxis = list(title = "Feature value"),
yaxis = list(title = "P(severe disease)"))24.5 The forest as a surface
rf_sim <- ranger(y ~ ., data = sim_tr, num.trees = 500, probability = TRUE,
num.threads = 1)
z_rf <- matrix(predict(rf_sim, gd2)$predictions[, "B"], length(gx))
plot_ly() |>
add_surface(x = gy, y = gx, z = z_tree, showscale = FALSE, opacity = 0.95,
colorscale = "RdBu", reversescale = TRUE) |>
add_surface(x = gy, y = gx, z = z_rf + 1.3, showscale = FALSE, opacity = 0.95,
colorscale = "RdBu", reversescale = TRUE) |>
add_surface(x = gy, y = gx, z = matrix(post2(gd2$x1, gd2$x2), length(gx)) + 2.6,
showscale = FALSE, opacity = 0.9, colorscale = "RdBu",
reversescale = TRUE) |>
layout(title = "Single tree (bottom), random forest (middle), true Bayes posterior (top)",
scene = list(xaxis = list(title = "x2"), yaxis = list(title = "x1"),
zaxis = list(title = "P(Y = B | x) + offset")))Averaging 500 trees turns the staircase into a smooth ramp that tracks the true posterior closely. Ensembling buys smoothness, and smoothness is what the single tree lacked. The cost is interpretability: the forest has no rule set to read.
25 PART V: COMPARING CLASSIFIERS
26 An honest comparison
Comparing methods is where the protocol of §5.4 earns its keep. Three requirements:
- Every model sees the same resamples. Comparing models scored on different splits confounds the method with the split.
- All models are fitted on training data only. Fitting on the test set makes the comparison meaningless, and the mistake is invisible in the output, because every model looks good.
- Differences carry uncertainty. Two accuracies differing by one point on a 400-case test set are indistinguishable. Use paired tests across common resamples, where the pairing removes split-to-split variation.
set.seed(1234)
ctrl <- trainControl(method = "repeatedcv", number = 10, repeats = 3,
classProbs = TRUE, summaryFunction = twoClassSummary,
savePredictions = "final", index = NULL)
# NOTE the data argument: qol_train throughout. Fitting on qol_test would
# produce excellent-looking numbers that mean nothing.
fit_rf <- train(cd ~ ., data = qol_train, method = "ranger", metric = "ROC",
trControl = ctrl, tuneLength = 3,
num.trees = 300, num.threads = 1)
fit_knn <- train(cd ~ ., data = qol_train, method = "knn", metric = "ROC",
trControl = ctrl, tuneLength = 8,
preProcess = c("center", "scale")) # scaling INSIDE each fold
fit_nb <- train(cd ~ ., data = qol_train, method = "naive_bayes",
metric = "ROC", trControl = ctrl, tuneLength = 3)
fit_c50 <- train(cd ~ ., data = qol_train, method = "C5.0", metric = "ROC",
trControl = ctrl, tuneLength = 3)The preProcess argument on the kNN fit is doing real
work: caret recomputes the centring and scaling
within every resample, using only that resample’s
analysis set. Standardizing the whole frame beforehand would leak (§5.4.1).
res <- resamples(list(RandomForest = fit_rf, kNN = fit_knn,
NaiveBayes = fit_nb, C5.0 = fit_c50))
summary(res)$statistics$ROC |> round(4)#> Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
#> RandomForest 0.6304 0.6620 0.6983 0.6949 0.7261 0.7644 0
#> kNN 0.4756 0.5365 0.5713 0.5681 0.6007 0.6593 0
#> NaiveBayes 0.5211 0.6008 0.6297 0.6281 0.6657 0.6951 0
#> C5.0 0.6293 0.6658 0.6823 0.6924 0.7241 0.7622 0res_long <- res$values |>
dplyr::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", outlier.alpha = 0.4, width = 0.55) +
geom_jitter(width = 0.08, alpha = 0.25, size = 0.8) +
coord_flip() +
labs(title = "Cross-validated AUC across 30 common resamples",
subtitle = "Same folds for every model, so the comparison is paired",
x = NULL, y = "AUC") +
theme_dspa()
# --- Interactive equivalent, plus caret's built-in lattice displays --------
plot_ly(res_long, x = ~ROC, color = ~model, type = "box") |>
layout(title = "Cross-validated AUC by model",
xaxis = list(title = "AUC"))
scales <- list(x = list(relation = "free"), y = list(relation = "free"))
bwplot(res, scales = scales)
densityplot(res, scales = scales, pch = "|")
dotplot(res, scales = scales)
splom(res)#>
#> Call:
#> summary.diff.resamples(object = diffs)
#>
#> p-value adjustment: bonferroni
#> Upper diagonal: estimates of the difference
#> Lower diagonal: p-value for H0: difference = 0
#>
#> ROC
#> RandomForest kNN NaiveBayes C5.0
#> RandomForest 0.12673 0.06679 0.00247
#> kNN 1.28e-11 -0.05994 -0.12426
#> NaiveBayes 2.32e-06 3.58e-05 -0.06431
#> C5.0 1 1.06e-11 1.82e-05
#>
#> Sens
#> RandomForest kNN NaiveBayes C5.0
#> RandomForest 0.0319 -0.2487 -0.0019
#> kNN 0.342 -0.2806 -0.0338
#> NaiveBayes 3.60e-15 < 2e-16 0.2468
#> C5.0 1.000 0.192 1.17e-13
#>
#> Spec
#> RandomForest kNN NaiveBayes C5.0
#> RandomForest 0.15358 0.39357 -0.00614
#> kNN 5.19e-11 0.23999 -0.15972
#> NaiveBayes < 2e-16 1.85e-13 -0.39971
#> C5.0 1 7.38e-12 < 2e-16The paired test removes fold-to-fold variation, which is usually larger than the between-model difference. Without pairing, real differences vanish into noise.
# Final check on data no model has touched
holdout <- data.frame(
model = c("Majority class", "Random forest", "kNN", "Naive Bayes", "C5.0"),
accuracy = round(c(
NIR_qol,
mean(predict(fit_rf, qol_test) == qol_test$cd),
mean(predict(fit_knn, qol_test) == qol_test$cd),
mean(predict(fit_nb, qol_test) == qol_test$cd),
mean(predict(fit_c50, qol_test) == qol_test$cd)), 4))
holdout$kappa <- round(c(0, vapply(list(fit_rf, fit_knn, fit_nb, fit_c50),
\(m) binary_metrics(predict(m, qol_test), qol_test$cd,
"severe_disease")[["kappa"]], numeric(1))), 4)
holdout$AUC <- round(c(0.5, vapply(list(fit_rf, fit_knn, fit_nb, fit_c50),
\(m) as.numeric(pROC::auc(pROC::roc(qol_test$cd,
predict(m, qol_test, type = "prob")[, "severe_disease"], quiet = TRUE))),
numeric(1))), 4)
# Binomial CI on each accuracy: differences smaller than the CI width are noise
holdout$acc_CI <- vapply(list(NULL, fit_rf, fit_knn, fit_nb, fit_c50), function(m) {
k <- if (is.null(m)) round(NIR_qol * nrow(qol_test)) else
sum(predict(m, qol_test) == qol_test$cd)
ci <- binom.test(k, nrow(qol_test))$conf.int
sprintf("[%.3f, %.3f]", ci[1], ci[2])
}, character(1))
holdout#> model accuracy kappa AUC acc_CI
#> 1 Majority class 0.5000 0.0000 0.5000 [0.452, 0.548]
#> 2 Random forest 0.6306 0.2613 0.6766 [0.584, 0.676]
#> 3 kNN 0.5721 0.1441 0.5782 [0.527, 0.621]
#> 4 Naive Bayes 0.5676 0.1351 0.6382 [0.520, 0.614]
#> 5 C5.0 0.6081 0.2162 0.6779 [0.561, 0.654]Read the confidence intervals before the point estimates. If they overlap substantially, the models are not distinguishable on this test set, and ranking them is over-reading. With ~440 test cases, the CI half-width is about 4.5 percentage points.
27 Further applications
27.1 Iris species
set.seed(1234)
iris_split <- rsample::initial_split(iris, prop = 0.75, strata = Species)
iris_tr <- rsample::training(iris_split); iris_te <- rsample::testing(iris_split)
iris_nb <- e1071::naiveBayes(Species ~ ., data = iris_tr)
iris_pred <- predict(iris_nb, iris_te)
confusionMatrix(data = iris_pred, reference = iris_te$Species)$table#> Reference
#> Prediction setosa versicolor virginica
#> setosa 13 0 0
#> versicolor 0 13 2
#> virginica 0 0 11c(test_accuracy = round(mean(iris_pred == iris_te$Species), 4),
no_information_rate = round(max(prop.table(table(iris_te$Species))), 4))#> test_accuracy no_information_rate
#> 0.9487 0.3333Note the split: evaluating on the training data, a common shortcut
with iris — reports resubstitution accuracy, which is
optimistic for every method and uninformative for comparison.
27.2 Refining the cancer-stage definition
SEER codes 8 (“not applicable”) and 9 (“unstaged”) are missing-data codes, not early stages. Grouping them with early-stage cancers mixes “we know this is early” with “we do not know”. Removing them makes the outcome mean one thing.
hn2 <- hn_med[!hn_med$seer_stage %in% c(8, 9), ]
c(n_original = nrow(hn_med), n_after_removing_8_9 = nrow(hn2))#> n_original n_after_removing_8_9
#> 662 580corpus2 <- Corpus(VectorSource(hn2$MEDICATION_SUMMARY)) |>
tm_map(content_transformer(tolower)) |> tm_map(removePunctuation) |>
tm_map(removeNumbers) |> tm_map(removeWords, stopwords("english")) |>
tm_map(stemDocument) |> tm_map(stripWhitespace)
set.seed(11)
idx2 <- sample(nrow(hn2), floor(0.8 * nrow(hn2)))
y2_tr <- factor(hn2$seer_stage[idx2] %in% c(4, 5, 7), levels = c(FALSE, TRUE),
labels = c("early_stage", "later_stage"))
y2_te <- factor(hn2$seer_stage[-idx2] %in% c(4, 5, 7), levels = c(FALSE, TRUE),
labels = c("early_stage", "later_stage"))
d2_tr_full <- DocumentTermMatrix(corpus2[idx2])
dict2 <- findFreqTerms(d2_tr_full, lowfreq = 5) # TRAINING vocabulary only
X2_tr <- to_indicator(DocumentTermMatrix(corpus2[idx2], list(dictionary = dict2)))
X2_te <- to_indicator(DocumentTermMatrix(corpus2[-idx2], list(dictionary = dict2)))
nb2 <- e1071::naiveBayes(X2_tr, y2_tr, laplace = a_best)
pred2 <- predict(nb2, X2_te)
round(binary_metrics(pred2, y2_te, positive = "later_stage"), 4)#> n TP FP FN
#> 116.0000 1.0000 3.0000 29.0000
#> TN accuracy NIR sensitivity
#> 83.0000 0.7241 0.7414 0.0333
#> specificity PPV NPV balanced_accuracy
#> 0.9651 0.2500 0.7411 0.4992
#> F1 kappa MCC
#> 0.0588 -0.0022 -0.0037Compare accuracy against NIR and
balanced_accuracy against 0.5. A high accuracy paired with
a balanced accuracy near 0.5 is the signature of a classifier that has
collapsed onto the majority class, which is what over-smoothing or a
strongly imbalanced outcome produces, and which raw accuracy alone would
hide.
27.3 Multi-class: predicting MLB player position
mlb <- dspa_read("https://umich.instructure.com/files/330381/download?download_frd=1",
"01a_data.txt", reader = read.table, header = TRUE)
set.seed(123)
mlb_split <- rsample::initial_split(mlb, prop = 0.75, strata = Position)
mlb_tr <- rsample::training(mlb_split); mlb_te <- rsample::testing(mlb_split)
# Numeric predictors -> Gaussian naive Bayes. Note: `laplace` applies only to
# CATEGORICAL features and is silently ignored for numeric ones.
mlb_nb <- e1071::naiveBayes(mlb_tr[, c("Weight", "Height", "Age")],
factor(mlb_tr$Position))
mlb_pred <- predict(mlb_nb, mlb_te[, c("Weight", "Height", "Age")])
tab_mlb <- table(Predicted = mlb_pred, Actual = factor(mlb_te$Position))
c(accuracy = round(sum(diag(tab_mlb)) / sum(tab_mlb), 4),
no_information_rate = round(max(prop.table(table(mlb_te$Position))), 4),
n_classes = nlevels(factor(mlb$Position)))#> accuracy no_information_rate n_classes
#> 0.3012 0.3243 9.0000as.data.frame(tab_mlb) |>
ggplot(aes(Actual, Predicted, fill = Freq)) +
geom_tile(colour = "white") +
scale_fill_gradient(low = "white", high = "#3B7DD8", name = "Count") +
labs(title = "Multi-class naive Bayes: MLB player position",
subtitle = "Three numeric predictors cannot separate nine positions",
x = "Actual position", y = "Predicted position") +
theme_dspa(9) +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
# --- Interactive equivalent ------------------------------------------------
tab_df <- tidyr::spread(as.data.frame(tab_mlb), key = Actual, value = Freq)
plot_ly(x = colnames(tab_mlb), y = rownames(tab_mlb),
z = as.matrix(tab_df[, -1]), type = "heatmap") |>
layout(title = "Actual vs. predicted MLB positions",
xaxis = list(title = "Actual"), yaxis = list(title = "Predicted"))Height, weight, and age carry little positional information, and the confusion matrix shows predictions collapsing onto the most common positions. This is a feature problem, not a method problem, no classifier recovers information the features do not contain.
27.4 Clinical text: 40 medical specialties
dataCT <- dspa_read("https://umich.instructure.com/files/21152999/download?download_frd=1",
"MedicalSpecialty_Notes.csv", header = TRUE)
str(dataCT, max.level = 1)#> 'data.frame': 4999 obs. of 6 variables:
#> $ Index : int 0 1 2 3 4 5 6 7 8 9 ...
#> $ description : chr " A 23-year-old white female presents with complaint of allergies." " Consult for laparoscopic gastric bypass." " Consult for laparoscopic gastric bypass." " 2-D M-Mode. Doppler. " ...
#> $ medical_specialty: chr " Allergy / Immunology" " Bariatrics" " Bariatrics" " Cardiovascular / Pulmonary" ...
#> $ sample_name : chr " Allergic Rhinitis " " Laparoscopic Gastric Bypass Consult - 2 " " Laparoscopic Gastric Bypass Consult - 1 " " 2-D Echocardiogram - 1 " ...
#> $ transcription : chr "SUBJECTIVE:, This 23-year-old white female presents with complaint of allergies. She used to have allergies w"| __truncated__ "PAST MEDICAL HISTORY:, He has difficulty climbing stairs, difficulty with airline seats, tying shoes, used to p"| __truncated__ "HISTORY OF PRESENT ILLNESS: , I have seen ABC today. He is a very pleasant gentleman who is 42 years old, 344 "| __truncated__ "2-D M-MODE: , ,1. Left atrial enlargement with left atrial diameter of 4.7 cm.,2. Normal size right and left "| __truncated__ ...
#> $ keywords : chr "allergy / immunology, allergic rhinitis, allergies, asthma, nasal sprays, rhinitis, nasal, erythematous, allegr"| __truncated__ "bariatrics, laparoscopic gastric bypass, weight loss programs, gastric bypass, atkin's diet, weight watcher's, "| __truncated__ "bariatrics, laparoscopic gastric bypass, heart attacks, body weight, pulmonary embolism, potential complication"| __truncated__ "cardiovascular / pulmonary, 2-d m-mode, doppler, aortic valve, atrial enlargement, diastolic function, ejection"| __truncated__ ...#> medical_specialty n
#> 1 Surgery 1103
#> 2 Consult - History and Phy. 516
#> 3 Cardiovascular / Pulmonary 372
#> 4 Orthopedic 355
#> 5 Radiology 273
#> 6 General Medicine 259
#> 7 Gastroenterology 230
#> 8 Neurology 223
#> 9 SOAP / Chart / Progress Notes 166
#> 10 Obstetrics / Gynecology 160c(n_documents = nrow(dataCT), n_specialties = nrow(spec_counts),
largest_class_share = round(max(spec_counts$n) / nrow(dataCT), 4))#> n_documents n_specialties largest_class_share
#> 4999.0000 40.0000 0.2206spec_counts |>
slice_head(n = 20) |>
ggplot(aes(n, reorder(medical_specialty, n))) +
geom_col(fill = "steelblue") +
labs(title = "Document counts by medical specialty (top 20)",
subtitle = "Severely imbalanced: any accuracy must be read against the largest-class share",
x = "Number of transcriptions", y = NULL) +
theme_dspa(9)
# --- Interactive equivalent ------------------------------------------------
plot_ly(dataCT, x = ~medical_specialty) |> add_histogram() |>
layout(title = "Distribution of medical specialties",
xaxis = list(title = ""), yaxis = list(title = "Count"))ct_corpus <- Corpus(VectorSource(dataCT$transcription)) |>
tm_map(content_transformer(tolower)) |> tm_map(removePunctuation) |>
tm_map(removeNumbers) |> tm_map(removeWords, stopwords("english")) |>
tm_map(stripWhitespace)
set.seed(1234)
ct_idx <- sample(nrow(dataCT), floor(0.8 * nrow(dataCT)))
ct_y_tr <- factor(dataCT$medical_specialty[ct_idx])
ct_y_te <- factor(dataCT$medical_specialty[-ct_idx], levels = levels(ct_y_tr))
# Vocabulary from TRAINING documents; capped by document frequency to keep the
# matrix manageable at ~5,000 documents
dtm_ct_tr <- DocumentTermMatrix(ct_corpus[ct_idx],
control = list(bounds = list(global = c(20, Inf))))
ct_dict <- Terms(dtm_ct_tr)
c(vocabulary_size = length(ct_dict))#> vocabulary_size
#> 4708ct_tr <- to_indicator(DocumentTermMatrix(ct_corpus[ct_idx],
list(dictionary = ct_dict)))
ct_te <- to_indicator(DocumentTermMatrix(ct_corpus[-ct_idx],
list(dictionary = ct_dict)))
ct_nb <- e1071::naiveBayes(ct_tr, ct_y_tr, laplace = 1)
ct_pred <- predict(ct_nb, ct_te)
c(accuracy = round(mean(ct_pred == ct_y_te, na.rm = TRUE), 4),
no_information_rate = round(max(prop.table(table(ct_y_te))), 4),
n_classes = nlevels(ct_y_tr))#> accuracy no_information_rate n_classes
#> 0.325 0.221 40.000top_spec <- head(spec_counts$medical_specialty, 12)
cm_ct <- table(Predicted = ct_pred, Actual = ct_y_te)
cm_sub <- cm_ct[rownames(cm_ct) %in% top_spec, colnames(cm_ct) %in% top_spec]
as.data.frame(cm_sub) |>
ggplot(aes(Actual, Predicted, fill = log1p(Freq))) +
geom_tile(colour = "white") +
scale_fill_viridis_c(option = "magma", direction = -1,
name = "log(1 + count)") +
labs(title = "Naive Bayes classification of clinical transcriptions",
subtitle = "Twelve most common specialties; a strong diagonal would indicate success",
x = "True specialty", y = "Predicted specialty") +
theme_dspa(8) +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
# --- Interactive equivalent ------------------------------------------------
tab_df <- tidyr::spread(as.data.frame(cm_ct), key = Actual, value = Freq)
plot_ly(x = colnames(cm_ct), y = rownames(cm_ct),
z = as.matrix(tab_df[, -1]), type = "heatmap") |>
layout(title = "Naive Bayes classification of 40 medical specialties",
xaxis = list(title = "True specialty"),
yaxis = list(title = "Predicted specialty"))Forty classes with a heavily skewed distribution is a hard problem for a bag-of-words model, and the diagonal is faint. Two extensions help: TF-IDF weighting rather than binary presence, and contextual embeddings, which Chapter 14 develops.
27.5 Three-class chronic disease
# Tertiles COMPUTED from the data, via cut() -- no hard-coded row counts or breaks
qol$cdthree <- cut(qol$CHRONICDISEASESCORE,
breaks = c(-Inf, quantile(qol$CHRONICDISEASESCORE, c(1/3, 2/3)), Inf),
labels = c("minor_disease", "mild_disease", "severe_disease"),
include.lowest = TRUE)
table(qol$cdthree)#>
#> minor_disease mild_disease severe_disease
#> 803 680 731#> 33.33333% 66.66667%
#> 1.06 1.80qol3 <- qol |> dplyr::select(-any_of(c(drop_cols, "cd")))
set.seed(1234)
sp3 <- rsample::initial_split(qol3, prop = 0.8, strata = cdthree)
tr3 <- rsample::training(sp3); te3 <- rsample::testing(sp3)
NIR3 <- max(prop.table(table(te3$cdthree)))
set.seed(1234)
m_c50 <- C5.0(cdthree ~ ., data = tr3, trials = 10)
m_1r <- OneR::OneR(cdthree ~ ., data = OneR::optbin(tr3))
m_rul <- C5.0(cdthree ~ ., data = tr3, rules = TRUE)
m_rf <- ranger(cdthree ~ ., data = tr3, num.trees = 500, num.threads = 1)
## OLD: need to convert both the predicted and true values to character inside
## the accuracy function. This bypasses the factor level mismatch entirely.
# acc <- function(p) mean(p == te3$cdthree, na.rm = TRUE)
acc <- function(p) {
mean(as.character(p) == as.character(te3$cdthree), na.rm = TRUE)
}
data.frame(
model = c("Majority class", "OneR", "C5.0 rules", "C5.0 boosted", "Random forest"),
test_accuracy = round(c(NIR3,
acc(predict(m_1r, te3)),
acc(predict(m_rul, te3)),
acc(predict(m_c50, te3)),
acc(predict(m_rf, te3)$predictions)), 4))#> model test_accuracy
#> 1 Majority class 0.3626
#> 2 OneR 0.4572
#> 3 C5.0 rules 0.4234
#> 4 C5.0 boosted 0.4212
#> 5 Random forest 0.4505confusionMatrix(data = predict(m_c50, te3), reference = te3$cdthree)$byClass[, c(1, 2, 11)] |>
round(4)#> Sensitivity Specificity Balanced Accuracy
#> Class: minor_disease 0.4845 0.7527 0.6186
#> Class: mild_disease 0.2500 0.6981 0.4740
#> Class: severe_disease 0.5102 0.6835 0.5969Every model is fitted on tr3 and scored on untouched
te3. Accuracy is lower than the two-class problem for the
obvious reason, three classes make the task harder and drop the
no-information rate to about \(1/3\),
but the models now beat their baseline by a clear margin, which the
two-class version did not. A lower accuracy against a lower
baseline can represent more actual learning, which is why the
comparison must always be to the baseline rather than across
problems.
Per-class sensitivity is the diagnostic to read: a model can score
well overall while being blind to the middle class, and the
byClass table exposes that where a single accuracy figure
would not.
28 Computational complexity summary
\(n\) = training cases, \(d\) = features, \(K\) = classes, \(B\) = trees, \(T\) = boosting iterations, \(L\) = leaves.
| Method | Training | Prediction (one case) | Memory | Notes |
|---|---|---|---|---|
| Majority class | \(O(n)\) | \(O(1)\) | \(O(1)\) | The baseline every model must beat |
| kNN (brute force) | \(O(1)\) | \(O(nd)\) | \(\mathbf{O(nd)}\) | No training; stores everything |
| kNN (k-d tree) | \(O(dn\log n)\) | \(O(d\log n)\) avg | \(O(nd)\) | Degrades to \(O(nd)\) above \(d\approx20\) |
| kNN (ball / cover tree) | \(O(dn\log n)\) | \(O(\log n)\) | \(O(nd)\) | Cover tree adapts to intrinsic dimension |
| kNN (LSH, approximate) | \(O(nd)\) | sublinear in \(n\) | \(O(nd)\) | Trades exactness for speed |
| Naive Bayes (categorical) | \(O(nd)\) | \(O(Kd)\) | \(O(Kd)\) | One pass; the fastest classifier here |
| Naive Bayes (Gaussian) | \(O(nd)\) | \(O(Kd)\) | \(O(Kd)\) | Two sufficient statistics per cell |
| LDA | \(O(nd^2+d^3)\) | \(O(Kd)\) | \(O(d^2)\) | Pooled covariance inversion |
| QDA | \(O(nd^2+Kd^3)\) | \(O(Kd^2)\) | \(O(Kd^2)\) | Needs \(n_k \gg d^2/2\) per class |
| Logistic regression | \(O(Tnd^2)\) IRLS | \(O(d)\) | \(O(d^2)\) | \(O(Tnd)\) with gradient methods |
| Decision tree (CART) | \(O(dn\log n)\) | \(O(\text{depth})\) | \(O(L)\) | Sort once per feature, then running sums |
| Cost-complexity pruning + CV | \(O(k\cdot dn\log n)\) | — | \(O(L)\) | \(k\)-fold over a nested subtree sequence |
| C5.0 boosting | \(O(T\,dn\log n)\) | \(O(T\cdot\text{depth})\) | \(O(TL)\) | Sequential; cannot parallelize over \(T\) |
| Random forest | \(O(B\,m\,n\log n)\), \(m=\)mtry |
\(O(B\cdot\text{depth})\) | \(O(BL)\) | Embarrassingly parallel over \(B\) |
| Extra-trees | \(O(B\,m\,n)\) | \(O(B\cdot\text{depth})\) | \(O(BL)\) | Random cut points remove the sort |
| OneR | \(O(nd)\) | \(O(1)\) | \(O(\text{bins})\) | The interpretability baseline |
| RIPPER | \(O(dn\log^2 n)\) | \(O(\text{rules})\) | \(O(\text{rules})\) | Grow, prune, optimize |
Three practical rules follow.
kNN’s cost is at prediction time, and it never goes away. Every query touches the whole training set unless an index helps, and indices stop helping above \(d\approx20\). If predictions must be fast or memory is bounded, kNN is the wrong choice regardless of its accuracy.
Naive Bayes is the cheapest thing that works. One pass to train, \(O(Kd)\) to predict, \(O(Kd)\) to store. It is the right first model on any text problem and the right sanity baseline everywhere else.
Random forests parallelize; boosting does not. Trees in a forest are independent, so \(B\) divides across cores. Boosting fits each tree to the previous residuals, so \(T\) is inherently sequential, which is why forests win on wall-clock time even when boosting wins on accuracy.
29 Common pitfalls
| # | Pitfall | Consequence | Fix |
|---|---|---|---|
| 1 | Reporting accuracy without the no-information rate | A majority-class rule looks like a model | Report NIR, \(\kappa\), balanced accuracy, MCC |
| 2 | Scaling or imputing before the split | Preprocessing leakage; optimistic everywhere | Learn constants on training; recompute per fold |
| 3 | Selecting features on all the data, then cross-validating | Above-chance results on pure noise | Screen inside every fold |
| 4 | Training on the test set | Meaningless comparison, invisible in output | Check the data = argument of every fit |
| 5 | Fitting a rule learner on all rows, “validating” on a subset | Resubstitution error reported twice | Disjoint train and test |
| 6 | confusionMatrix(truth, pred) |
Transposes the matrix; swaps sensitivity and specificity | confusionMatrix(data = pred, reference = truth, positive = ...) |
| 7 | Treating attr(knn(...), "prob") as \(P(\text{positive})\) |
It is the winning class’s vote share | Convert: ifelse(pred == pos, p, 1 - p) |
| 8 | Confusing sensitivity with specificity | Inverted clinical conclusions | \(\text{Sens}=TP/(TP+FN)\), \(\text{Spec}=TN/(TN+FP)\) |
| 9 | Computing rates from prop.row |
Denominators are identically 1 | Compute from raw counts |
| 10 | Choosing \(k\) by training error | Always selects \(k=1\) | Cross-validate |
| 11 | Using \(k=\sqrt n\) as if it were theory | Wrong by orders of magnitude in high \(d\) | \(k^\star \asymp n^{4/(d+4)}\); cross-validate |
| 12 | kNN on unscaled features | Large-scale features dominate the distance | Standardize with training statistics |
| 13 | kNN with many irrelevant features | Distance concentration; accuracy collapses to chance | Reduce dimension or select features first |
| 14 | Describing kNN as stochastic | Masks an unset seed or an unstable split | Only tie-breaking is random |
| 15 | Trusting naive Bayes probabilities | Systematically over-confident under dependence | Trust the labels; recalibrate the probabilities |
| 16 | Computing \(\prod_i P(F_i \mid C)\) directly | Underflows to zero past ~300 features | Sum logs; log-sum-exp to normalize |
| 17 | Large Laplace \(\alpha\) without tuning | Collapses to the prior; accuracy rises to the NIR | Tune \(\alpha\) by held-out log-loss |
| 18 | Applying laplace to numeric predictors |
Silently ignored | Smoothing applies to categorical features only |
| 19 | Changing two things and attributing to one | Confounded comparison | Vary one factor at a time |
| 20 | Splitting on misclassification error | Zero gain for informative splits | Split on Gini/entropy, prune on error |
| 21 | Reading raw information gain as importance | High-cardinality features win by construction | Gain ratio, or permutation importance |
| 22 | Keeping IDs, dates, or record numbers as features | Leakage by proxy; will not generalize | Ask: available and meaningful at prediction time? |
| 23 | Reading impurity importance with correlated features | Importance split arbitrarily within groups | Permutation importance; group-aware methods |
| 24 | Ranking models by point accuracy | Differences smaller than the CI are noise | Paired tests on common resamples; report CIs |
30 Practice problems
30.1 Problem 1: Verify the Cover–Hart bound
Construct a binary problem whose Bayes error you can compute exactly, and confirm empirically that 1-NN error lies in \([R^{*}, 2R^{*}]\) as \(n\) grows.
Solution
# One-dimensional problem: two Gaussians with equal variance and equal priors.
# The Bayes rule cuts at the midpoint, so R* has a closed form.
m1 <- 0; m2 <- 1.5; s <- 1
R_star <- pnorm((m1 + m2) / 2, m2, s) # P(wrong | class 2), by symmetry
sim1d <- function(n, n_te = 4000, seed = 301) {
set.seed(seed + n)
y <- rbinom(n, 1, 0.5); x <- rnorm(n, ifelse(y == 1, m2, m1), s)
yt <- rbinom(n_te, 1, 0.5); xt <- rnorm(n_te, ifelse(yt == 1, m2, m1), s)
c(n = n,
err_1nn = mean(knn(matrix(x), matrix(xt), factor(y), k = 1) != factor(yt)))
}
tab1 <- as.data.frame(do.call(rbind, lapply(c(25, 100, 500, 2500, 12500), sim1d)))
tab1$R_star <- R_star; tab1$upper <- 2 * R_star * (1 - R_star)
tab1$within_bound <- tab1$err_1nn >= R_star & tab1$err_1nn <= 2 * R_star
round(tab1, 4)#> n err_1nn R_star upper within_bound
#> 1 25 0.2688 0.2266 0.3505 1
#> 2 100 0.3142 0.2266 0.3505 1
#> 3 500 0.2978 0.2266 0.3505 1
#> 4 2500 0.2948 0.2266 0.3505 1
#> 5 12500 0.3172 0.2266 0.3505 1ggplot(tab1, aes(n, err_1nn)) +
geom_ribbon(aes(ymin = R_star, ymax = 2 * R_star), fill = "grey88") +
geom_line(linewidth = 1, colour = "steelblue") + geom_point(size = 2.4) +
geom_hline(yintercept = R_star, linetype = "dashed") +
scale_x_log10() +
labs(title = "1-NN error converges into the Cover-Hart band",
subtitle = "Shaded region is [R*, 2R*]; dashed line is the Bayes error",
x = "n (log scale)", y = "1-NN test error") +
theme_dspa()
30.2 Problem 2: Quantify preprocessing leakage
Measure how much optimism preprocessing leakage introduces, as a function of sample size.
Solution
leak_gap <- function(n, d = 400, reps = 12, seed = 311) {
set.seed(seed + n)
mean(replicate(reps, {
X <- matrix(rnorm(n * d), n, d)
y <- factor(sample(c("A", "B"), n, TRUE)) # NO signal at all
# Leaky: standardize the whole matrix, then split
Xs <- scale(X)
tr <- sample(n, floor(0.7 * n))
leaky <- mean(knn(Xs[tr, ], Xs[-tr, ], y[tr], k = 5) != y[-tr])
# Clean: standardize using TRAINING statistics only
ctr <- colMeans(X[tr, ]); scl <- apply(X[tr, ], 2, sd)
clean <- mean(knn(scale(X[tr, ], ctr, scl), scale(X[-tr, ], ctr, scl),
y[tr], k = 5) != y[-tr])
clean - leaky # positive = leaky looks better
}))
}
gaps <- data.frame(n = c(40, 80, 160, 320, 640),
optimism = vapply(c(40, 80, 160, 320, 640), leak_gap, numeric(1)))
round(gaps, 4)#> n optimism
#> 1 40 -0.0278
#> 2 80 -0.0035
#> 3 160 0.0087
#> 4 320 0.0113
#> 5 640 0.0026ggplot(gaps, aes(n, optimism)) +
geom_hline(yintercept = 0, colour = "grey50") +
geom_line(linewidth = 1, colour = "firebrick") + geom_point(size = 2.4) +
scale_x_log10() +
labs(title = "Optimism from standardizing before splitting",
subtitle = "Error difference (clean - leaky) on data with no signal; positive means leakage helped",
x = "n (log scale)", y = "Error difference") +
theme_dspa()
30.3 Problem 3: The
knn(prob = TRUE) trap
Show concretely that thresholding the raw prob attribute
produces different — and wrong, labels compared to thresholding the
positive-class posterior.
Solution
set.seed(317)
n3 <- 400
X3 <- matrix(rnorm(n3 * 2), n3, 2)
y3 <- factor(ifelse(X3[, 1] + X3[, 2] + rnorm(n3, sd = 0.7) > 0.6, "Pos", "Neg"))
tr3i <- 1:300; te3i <- 301:n3
p <- knn(X3[tr3i, ], X3[te3i, ], y3[tr3i], k = 11, prob = TRUE)
raw <- attr(p, "prob")
correct <- ifelse(p == "Pos", raw, 1 - raw)
thresh_raw <- factor(ifelse(raw > 0.6, "Pos", "Neg"), levels = levels(y3))
thresh_ok <- factor(ifelse(correct > 0.6, "Pos", "Neg"), levels = levels(y3))
c(range_of_raw_attribute = paste(round(range(raw), 3), collapse = " - "),
labels_that_differ = sum(thresh_raw != thresh_ok),
accuracy_raw_threshold = round(mean(thresh_raw == y3[te3i]), 4),
accuracy_correct_threshold = round(mean(thresh_ok == y3[te3i]), 4))#> range_of_raw_attribute labels_that_differ
#> "0.545 - 1" "63"
#> accuracy_raw_threshold accuracy_correct_threshold
#> "0.38" "0.89"# The decisive check: cases the classifier confidently called Neg
confident_neg <- p == "Neg" & raw > 0.6
c(confidently_negative_cases = sum(confident_neg),
relabelled_Pos_by_raw_rule = sum(thresh_raw[confident_neg] == "Pos"))#> confidently_negative_cases relabelled_Pos_by_raw_rule
#> 63 6330.4 Problem 4: Locate the bias–variance minimum
Show that the CV-optimal \(k\) grows with \(n\), and compare its growth against \(\sqrt n\) and \(n^{4/(d+4)}\).
Solution
## OLD: error in comparing factors with different level sets inside the
## cross‑validation loop. In some folds, the training or assessment set may
## contain only one class (especially for small n or imbalanced classes),
## causing knn() to return a factor with only one level, while b$y may have both
## original levels. Converting both to character before comparing
# best_k <- function(n, d = 2, seed = 331) {
# set.seed(seed + n)
# lab <- rbinom(n, 1, 1 - pi1)
# X <- rbind(MASS::mvrnorm(max(sum(lab == 0), 2), mu1, S1),
# MASS::mvrnorm(max(sum(lab == 1), 2), mu2, S2))
# y <- factor(c(rep("A", max(sum(lab == 0), 2)), rep("B", max(sum(lab == 1), 2))))
# ks <- unique(pmax(1, round(seq(1, min(201, nrow(X) - 1), length.out = 22))))
# cvf <- rsample::vfold_cv(data.frame(X, y = y), v = 5, strata = y)
# e <- vapply(ks, function(k) mean(vapply(cvf$splits, function(s) {
# a <- rsample::analysis(s); b <- rsample::assessment(s)
# mean(knn(a[, 1:2], b[, 1:2], a$y, k = k) != b$y)
# }, numeric(1))), numeric(1))
# c(n = nrow(X), k_cv = ks[which.min(e)])
# }
# kk <- as.data.frame(do.call(rbind, lapply(c(50, 100, 250, 500, 1000, 2500), best_k)))
best_k <- function(n, d = 2, seed = 331) {
set.seed(seed + n)
lab <- rbinom(n, 1, 1 - pi1) # pi1 assumed defined in environment
X <- rbind(MASS::mvrnorm(max(sum(lab == 0), 2), mu1, S1),
MASS::mvrnorm(max(sum(lab == 1), 2), mu2, S2))
y <- factor(c(rep("A", max(sum(lab == 0), 2)), rep("B", max(sum(lab == 1), 2))))
# Create cross-validation folds
cvf <- rsample::vfold_cv(data.frame(X, y = y), v = 5, strata = y)
# Determine the smallest training set size across all folds
min_train_size <- min(vapply(cvf$splits, function(s) nrow(rsample::analysis(s)), integer(1)))
# Candidate k values that are safe (k < min_train_size) and also cap at 201
ks <- unique(pmax(1, round(seq(1, min(201, min_train_size - 1), length.out = 22))))
# Compute cross-validated error for each k, using character comparison
e <- vapply(ks, function(k) {
mean(vapply(cvf$splits, function(s) {
a <- rsample::analysis(s)
b <- rsample::assessment(s)
pred <- class::knn(a[, 1:2, drop = FALSE], b[, 1:2, drop = FALSE], a$y, k = k, use.all = FALSE)
mean(as.character(pred) != as.character(b$y))
}, numeric(1)))
}, numeric(1))
c(n = nrow(X), k_cv = ks[which.min(e)])
}
kk <- as.data.frame(do.call(rbind, lapply(c(50, 100, 250, 500, 1000, 2500), best_k)))
kk$sqrt_n <- round(sqrt(kk$n))
kk$rate_d2 <- round(kk$n^(4 / (2 + 4))) # d = 2
kk#> n k_cv sqrt_n rate_d2
#> 1 50 3 7 14
#> 2 100 19 10 22
#> 3 250 20 16 40
#> 4 500 11 22 63
#> 5 1000 49 32 100
#> 6 2500 191 50 184kk |> pivot_longer(-n, names_to = "rule", values_to = "k") |>
ggplot(aes(n, k, colour = rule)) +
geom_line(linewidth = 1) + geom_point(size = 2.2) +
scale_x_log10() + scale_y_log10() +
scale_colour_brewer(palette = "Set1") +
labs(title = "How the optimal k grows with n",
subtitle = "d = 2, so the theoretical rate is n^(4/6) = n^(2/3)",
x = "n (log scale)", y = "k (log scale)", colour = NULL) +
theme_dspa()
30.5 Problem 5: Naive Bayes under dependence
Vary the correlation between two informative features and track accuracy, confidence, and Brier score.
Solution
nb_dep <- function(rho, n = 1200, seed = 337) {
set.seed(seed)
Sg <- matrix(c(1, rho, rho, 1), 2)
y <- factor(sample(c("A", "B"), n, TRUE))
Z <- MASS::mvrnorm(n, c(0, 0), Sg)
Z[y == "B", ] <- Z[y == "B", ] + 1.1
d <- data.frame(Z, y = y); tr <- 1:800; te <- 801:n
m <- e1071::naiveBayes(y ~ ., data = d[tr, ])
p <- predict(m, d[te, ], type = "raw")[, "B"]
c(rho = rho, accuracy = mean(predict(m, d[te, ]) == d$y[te]),
mean_confidence = mean(pmax(p, 1 - p)),
brier = mean((p - as.integer(d$y[te] == "B"))^2))
}
dep <- as.data.frame(do.call(rbind, lapply(c(0, 0.3, 0.6, 0.85, 0.95, 0.99), nb_dep)))
round(dep, 4)#> rho accuracy mean_confidence brier
#> 1 0.00 0.7575 0.7743 0.1530
#> 2 0.30 0.7800 0.7885 0.1593
#> 3 0.60 0.7600 0.7981 0.1737
#> 4 0.85 0.7475 0.8051 0.1844
#> 5 0.95 0.7475 0.8077 0.1885
#> 6 0.99 0.7450 0.8087 0.1900dep |> dplyr::select(rho, Accuracy = accuracy, `Mean confidence` = mean_confidence,
`Brier score` = brier) |>
pivot_longer(-rho, names_to = "metric", values_to = "value") |>
ggplot(aes(rho, value, colour = metric)) +
geom_line(linewidth = 1) + geom_point(size = 2.2) +
facet_wrap(~ metric, scales = "free_y") +
scale_colour_brewer(palette = "Set1", guide = "none") +
labs(title = "Naive Bayes as the independence assumption is violated",
x = expression(rho), y = NULL) +
theme_dspa(10)
30.6 Problem 6: Impurity measures compared
Fit trees using Gini and entropy, compare their structures and accuracies, and construct a split where misclassification error gives zero gain.
Solution
set.seed(1234)
t_gini <- rpart(cd ~ ., data = qol_train, method = "class",
parms = list(split = "gini"))
t_info <- rpart(cd ~ ., data = qol_train, method = "class",
parms = list(split = "information"))
data.frame(
criterion = c("Gini", "Information (entropy)"),
leaves = c(sum(t_gini$frame$var == "<leaf>"), sum(t_info$frame$var == "<leaf>")),
first_split = c(as.character(t_gini$frame$var[1]), as.character(t_info$frame$var[1])),
test_accuracy = round(c(
mean(predict(t_gini, qol_test, type = "class") == qol_test$cd),
mean(predict(t_info, qol_test, type = "class") == qol_test$cd)), 4))#> criterion leaves first_split test_accuracy
#> 1 Gini 3 CHARLSONSCORE 0.6194
#> 2 Information (entropy) 3 CHARLSONSCORE 0.6194# A split where misclassification error registers ZERO gain
parent <- c(300, 300); L <- c(240, 160); R <- c(60, 140)
w <- c(sum(L), sum(R)) / sum(parent)
sapply(imp_fns, function(f)
c(gain = f(parent) - (w[1] * f(L) + w[2] * f(R)))) |> round(5)#> Entropy.gain Gini.gain Misclassification.gain
#> 0.05894 0.04000 0.1333330.7 Problem 7: Verify the ensemble variance identity
Confirm \(\operatorname{Var}(\bar
f)=\rho\sigma^2+\frac{1-\rho}{B}\sigma^2\) by simulation, and
show that lowering mtry lowers \(\rho\).
Solution
set.seed(347)
tree_preds <- function(mtry, B = 40, reps = 50) {
P <- replicate(reps, {
idx <- sample(nrow(sim_tr), replace = TRUE)
vapply(seq_len(B), function(b) {
bb <- sample(idx, replace = TRUE)
f <- ranger::ranger(y ~ ., data = sim_tr[bb, ], num.trees = 1,
mtry = mtry, probability = TRUE, num.threads = 1)
mean(predict(f, sim_te[1:50, ])$predictions[, "B"])
}, numeric(1))
})
sig2 <- mean(apply(P, 2, var)) # variance of a single tree
rho <- mean(cor(t(P))[upper.tri(diag(reps))]) # correlation between trees
c(mtry = mtry, sigma2 = sig2, rho = rho,
predicted_var = rho * sig2 + (1 - rho) * sig2 / B,
observed_var = var(colMeans(P)))
}
rbind(tree_preds(2), tree_preds(1)) |> round(6)#> mtry sigma2 rho predicted_var observed_var
#> [1,] 2 0.001141 NA NA 0.000188
#> [2,] 1 0.001281 NA NA 0.000215mtry from 2 to 1 lowers \(\rho\), and the ensemble variance falls
with it, the decorrelation mechanism made quantitative.
30.8 Problem 8: Detect leakage by proxy
Build a dataset with a leaky administrative feature, show that it dominates importance rankings, and demonstrate that it does not generalize.
Solution
set.seed(353)
n8 <- 1500
# Cohort 1 (first 750 rows) collected under a protocol biased toward class B
cohort <- rep(1:2, each = n8 / 2)
y8 <- factor(ifelse(runif(n8) < ifelse(cohort == 1, 0.75, 0.25), "B", "A"))
X8 <- matrix(rnorm(n8 * 4), n8, 4)
X8[, 1] <- X8[, 1] + ifelse(y8 == "B", 0.55, 0) # a weak REAL signal
intake_date <- as.numeric(cohort) * 1000 + seq_len(n8) # an administrative field
d8 <- data.frame(X8, intake_date = intake_date, y = y8)
# Random split: cohorts are mixed, so intake_date "works"
set.seed(354)
mix <- sample(n8, 1000)
f_mix <- ranger(y ~ ., data = d8[mix, ], num.trees = 300,
importance = "permutation", num.threads = 1)
acc_mix <- mean(predict(f_mix, d8[-mix, ])$predictions == d8$y[-mix])
# Temporal split: train on cohort 1, test on cohort 2 -- the realistic setting
f_tmp <- ranger(y ~ ., data = d8[cohort == 1, ], num.trees = 300,
num.threads = 1)
acc_tmp <- mean(predict(f_tmp, d8[cohort == 2, ])$predictions == d8$y[cohort == 2])
# Same temporal split, with the administrative field removed
f_clean <- ranger(y ~ ., data = d8[cohort == 1, setdiff(names(d8), "intake_date")],
num.trees = 300, num.threads = 1)
acc_clean <- mean(predict(f_clean, d8[cohort == 2, ])$predictions == d8$y[cohort == 2])
round(sort(ranger::importance(f_mix), decreasing = TRUE), 4)#> intake_date X1 X4 X3 X2
#> 0.1333 0.0300 0.0012 0.0008 -0.0002c(random_split_with_intake_date = round(acc_mix, 4),
temporal_split_with_intake_date = round(acc_tmp, 4),
temporal_split_without_it = round(acc_clean, 4),
majority_class_in_cohort2 = round(max(prop.table(table(d8$y[cohort == 2]))), 4))#> random_split_with_intake_date temporal_split_with_intake_date
#> 0.7320 0.3573
#> temporal_split_without_it majority_class_in_cohort2
#> 0.3240 0.7853intake_date tops the importance ranking and gives
excellent accuracy under a random split. Under the
temporal split, training on one cohort, testing on the
next, which is what deployment looks like, it collapses to
below the majority-class rule, because the model learned a
mapping from date to class that inverts in the new cohort. Removing it
restores sensible performance.
31 Checkpoint
- A colleague reports 91% accuracy on a dataset where 90% of cases are negative. What do you ask for?
attr(knn(...), "prob")returns 0.82 for a case the classifier labelled “Control”. What is \(\hat P(\text{Recidivism}\mid x)\)?- Why does training error always select \(k=1\), and what should be used instead?
- Naive Bayes gives 88% accuracy but its predicted probabilities are almost all above 0.99. What is happening, and does it matter?
- Misclassification error registers zero gain for a split that entropy scores highly. Which measure is wrong, and what does each one get used for?
- Adding trees to a random forest stops improving accuracy at \(B=200\). What would you change to improve it further?
Answers
- The no-information rate, which is 0.90, so 91% is one point above always predicting the majority class. Ask for \(\kappa\), balanced accuracy, MCC, or AUC, all of which are near their null value for a majority-class rule, and for a binomial confidence interval on the accuracy plus the one-sided test of Accuracy > NIR. On a 500-case test set the CI half-width is about 2.5 points, so 91% versus 90% is not distinguishable.
- \(1 - 0.82 = \mathbf{0.18}\). The attribute is the vote share of the winning class, and the winner here was “Control”, so 0.82 is \(\hat P(\text{Control}\mid x)\). Thresholding the raw attribute at 0.6 would relabel this confidently-Control case as “Recidivism”.
- At \(k=1\) every training point is its own nearest neighbour, so the training error is exactly zero, kNN interpolates every label including mislabelled ones. Training error decreases monotonically as \(k\) decreases, so minimizing it always picks the most overfit model. Use cross-validated error, with any scaling recomputed inside each fold.
- The features are dependent, so correlated evidence is counted repeatedly and the posterior is pushed toward 0 or 1. It matters if you use the probabilities, for expected-cost decisions, risk stratification, or downstream combination, because they are badly calibrated (check the Brier score and a calibration curve). It does not matter much if you only use the labels, since the argmax is invariant to monotone distortion. Recalibrate with Platt scaling or isotonic regression if you need honest probabilities.
- Neither is wrong; they measure different things. Misclassification error is piecewise linear, so a split that changes the class proportions without changing which class is the majority in either child produces zero improvement. Entropy and Gini are strictly concave, so any change in proportions registers. Standard practice: split on Gini or entropy, prune on misclassification error or cross-validated error, splitting needs sensitivity to probabilities, pruning wants the quantity you actually care about.
- Nothing about \(B\), the variance
identity \(\rho\sigma^2 +
\frac{1-\rho}{B}\sigma^2\) says the second term is already
negligible and more trees will not help. The floor is \(\rho\sigma^2\), so improvement must come
from lowering \(\rho\)
(reduce
mtry, or switch to extremely randomized trees) or from lowering \(\sigma^2\) (deeper or better-tuned individual trees, better features). If the floor is set by irreducible noise, no ensemble change helps and the problem is the features.
32 Summary
The problem
- The Bayes classifier \(\arg\max_k P(Y=k\mid x)\) minimizes 0–1 loss, and the Bayes error is an irreducible floor. Every method here is a plug-in estimate of that posterior.
- Accuracy alone is uninterpretable. Compare against the no-information rate; report \(\kappa\), balanced accuracy, MCC, AUC, and a confidence interval.
- Ranking and calibration are different properties. AUC measures the first, the Brier score and calibration curves the second.
- Leakage, preprocessing, target, temporal, group, inflates every estimate and is invisible in the output. Learn all preprocessing constants from training data, inside every fold.
Instance-based learning
- kNN is a local vote. Cover–Hart: asymptotic 1-NN error is at most \(2R^{*}\), with no model and no fitting.
- \(k\) trades bias against variance; the rate-optimal choice is \(k^\star \asymp n^{4/(d+4)}\), not \(\sqrt n\). Cross-validate.
- Distance concentration makes kNN fail in high dimensions. Reduce dimension first.
- No training cost; \(O(nd)\) per query and \(O(nd)\) memory, and spatial indices stop helping above \(d\approx20\).
Probabilistic learning
- Naive Bayes assumes conditional independence, which is essentially always false, and classifies well anyway, because the argmax survives distortions that ruin the probabilities.
- Trust naive Bayes labels, not its probabilities; recalibrate if you need them.
- Laplace smoothing is the posterior mean under a Dirichlet prior. Tune \(\alpha\) by held-out log-loss; tuning on accuracy selects the degenerate model.
- Compute in log space. The direct product underflows past a few hundred features.
- Gaussian NB is diagonal-covariance QDA; LDA shares one covariance and gets a linear boundary.
Trees, rules, ensembles
- Split on Gini or entropy (strictly concave, so informative splits always register); prune on misclassification or CV error.
- Information gain favours high-cardinality features; gain ratio corrects it. The same bias afflicts impurity-based variable importance.
- Pruning buys interpretability at essentially no predictive cost.
- Averaging gives \(\operatorname{Var}=\rho\sigma^2+\frac{1-\rho}{B}\sigma^2\): more trees remove the second term, decorrelation lowers the first.
- OOB error is a free cross-validation estimate. Prefer permutation importance, and read it against the data-generating process.
Where these threads continue
| Thread | Continues in |
|---|---|
| Kernels and margins; boosting in depth | Black-box methods |
| Association rules and text mining | NLP and rule learning |
| Clustering without labels | Unsupervised clustering |
| Resampling, tuning, and calibration in depth | Model assessment |
| Regularized selection; FDR control | Feature selection |
| Classifiers over time; concept drift | Longitudinal analysis |
| Representation learning for text and images | Deep learning |
33 Chapter roadmap
- Chapter 1:
Foundations. R toolchain, reproducibility conventions,
dspa_read(), simulation. - Chapter 2: Data quality and exploratory visual analytics. Missingness, contingency tables, class imbalance, SMOTE.
- Chapter 3: Linear algebra, matrix computing, and regression. Projection, SVD, conditioning, trees, BART.
- Chapter 4: Dimensionality reduction. PCA, ICA, factor analysis, distance concentration, t-SNE, UMAP.
- Black-box methods. Neural networks, SVM and the kernel trick, gradient boosting, XGBoost.
- Text mining and association rules. NLP pipelines, Apriori, qualitative learning.
- Unsupervised clustering. k-means, hierarchical, spectral, Gaussian mixtures.
- Model assessment, validation, improvement. Nested resampling, calibration, hyperparameter optimization, cost-sensitive learning.
- Variable importance and feature selection. Ridge, LASSO, elastic net, stability selection, FDR control.
- Longitudinal and time-series analysis. Mixed models, GEE, concept drift.
- Deep learning. Convolutional and recurrent architectures, embeddings, transfer learning.
34 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] grid stats graphics grDevices utils datasets methods
#> [8] base
#>
#> other attached packages:
#> [1] ranger_0.16.0 OneR_2.2 C50_0.1.8 rpart.plot_3.1.2
#> [5] rpart_4.1.23 partykit_1.2-20 mvtnorm_1.2-4 libcoin_1.0-10
#> [9] wordcloud_2.6 RColorBrewer_1.1-3 tm_0.7-13 NLP_0.2-1
#> [13] FNN_1.1.4 rsample_1.2.1 recipes_1.4.0 pROC_1.18.5
#> [17] caret_6.0-94 lattice_0.22-6 class_7.3-22 plotly_4.12.0
#> [21] patchwork_1.3.0 tidyr_1.3.1 dplyr_1.1.4 ggplot2_4.0.1
#>
#> loaded via a namespace (and not attached):
#> [1] rlang_1.1.5 magrittr_2.0.3 furrr_0.3.1
#> [4] otel_0.2.0 e1071_1.7-14 compiler_4.3.3
#> [7] vctrs_0.6.5 reshape2_1.4.4 stringr_1.5.1
#> [10] pkgconfig_2.0.3 fastmap_1.2.0 inum_1.0-5
#> [13] labeling_0.4.3 rmarkdown_2.31 prodlim_2024.06.25
#> [16] purrr_1.0.2 xfun_0.52 cachem_1.1.0
#> [19] jsonlite_1.8.9 SnowballC_0.7.1 parallel_4.3.3
#> [22] R6_2.6.1 bslib_0.9.0 stringi_1.8.4
#> [25] parallelly_1.37.1 lubridate_1.9.3 jquerylib_0.1.4
#> [28] Rcpp_1.0.14 iterators_1.0.14 knitr_1.51
#> [31] future.apply_1.11.2 Matrix_1.6-5 splines_4.3.3
#> [34] nnet_7.3-19 timechange_0.3.0 tidyselect_1.2.1
#> [37] rstudioapi_0.18.0 yaml_2.3.10 timeDate_4032.109
#> [40] codetools_0.2-20 listenv_0.9.1 tibble_3.2.1
#> [43] plyr_1.8.9 withr_3.0.2 S7_0.2.1
#> [46] evaluate_1.0.3 future_1.33.2 survival_3.7-0
#> [49] proxy_0.4-27 isoband_0.2.7 xml2_1.3.6
#> [52] pillar_1.10.1 foreach_1.5.2 stats4_4.3.3
#> [55] generics_0.1.3 scales_1.4.0 globals_0.16.3
#> [58] glue_1.8.0 slam_0.1-50 lazyeval_0.2.2
#> [61] tools_4.3.3 data.table_1.16.4 ModelMetrics_1.2.2.2
#> [64] gower_1.0.1 crosstalk_1.2.1 ipred_0.9-14
#> [67] nlme_3.1-165 Cubist_0.4.4 Formula_1.2-5
#> [70] cli_3.6.3 viridisLite_0.4.2 lava_1.8.0
#> [73] gtable_0.3.6 sass_0.4.9 digest_0.6.37
#> [76] naivebayes_1.0.0 htmlwidgets_1.6.4 farver_2.1.2
#> [79] htmltools_0.5.8.1 lifecycle_1.0.5 hardhat_1.4.3
#> [82] httr_1.4.7 MASS_7.3-60.0.1