Chapter 8, Part 2: Sections 8.9–8.18

Biomedical Physics with Applications to Disease (BPAD)

Part 1 asked what a model can legitimately learn from a real cohort. Part 2 asks what happens when the model is given the image itself, and then asks the harder question: whether anything it learned survives contact with a different scanner, a different hospital, and a different patient.

How to Use Part 2

This is part 2 of BPAD Chapter 8, part 1 of Chapter 8 is available here.

Part 1 (Sections 8.1–8.8) developed the framing, the mathematics, study design, the real KiTS19 case study, image-to-feature extraction, preprocessing, supervised learning, and time-to-event modelling. Part 2 continues directly:

8.9 Unsupervised learning, segmentation, and spatial structure
8.10 Deep learning and image-native AI
8.11 Evaluation, calibration, uncertainty, and clinical utility
8.12 Translation: shift, fairness, interpretability, privacy, deployment
8.13 Reproducible computational practice and the capstone project
8.14–8.18 Chapter summary, key equations, glossary, practice problems, and solutions

This file is self-contained and rebuilds the same real data layer, so it can be knitted independently of Part 1. Section numbering, notation, the colour palette, and the callout conventions are identical, and cross-references to Sections 8.1–8.8 point into Part 1.

Part 2 clinical puzzle. The malignancy classifier of Part 1 was honest and unimpressive. Three responses are now on the table, and only one of them is scientifically defensible. Response A: cluster the patients and look for a phenotype the supervised model missed. Response B: abandon hand-crafted features and train a convolutional network directly on the images. Response C: accept the result, quantify the uncertainty properly, and ask whether the limiting factor is the sample size, the labels, or the physics. Sections 8.9 and 8.10 test the first two responses on the real cohort, and Sections 8.11 through 8.13 pursue the third.

Expectations, restated. Part 1 ended with a malignancy model of modest performance, and Section 8.11.1 will show that its bootstrap confidence interval is about 0.30 wide, spanning everything from barely-better-than-chance to moderately useful. Nothing in Part 2 narrows it, and the reader should not expect it to. Unsupervised clustering will recover the surgical decision rather than the disease; a neural network will improve its training discrimination monotonically with capacity while its out-of-fold discrimination stays flat; and the subgroup analysis will produce confidence intervals wider than the entire range of the metric. Every one of these is the correct result, and learning to recognize them is the point.

Learning Objectives for Part 2

#	Objective	Section
28	State what unsupervised structure can and cannot establish, and identify the circularity in post-hoc phenotype naming.	8.9
29	Apply k-means, assess cluster quality with silhouette width and the gap statistic, and test partition stability by resampling.	8.9
30	Distinguish overlap metrics (Dice, Jaccard) from boundary metrics (Hausdorff, ASSD) and show that they fail in different directions.	8.9
31	Construct a graph Laplacian from an image and use the Fiedler vector to obtain a spectral bipartition.	8.9
32	Derive and implement backpropagation for a one-hidden-layer network, and verify the gradients numerically.	8.10
33	Explain convolution as a linear, weight-shared operator and compute the parameter saving over a dense layer.	8.10
34	Compute receptive field and parameter counts for a CNN, and relate the receptive field to the physical size of the target lesion.	8.10
35	Describe encoder-decoder and U-Net architectures and the role of skip connections.	8.10
36	Design augmentation that reflects physically plausible acquisition variability.	8.10
37	Explain transfer learning, self-supervision, attention, and multiple-instance learning, and state when each is appropriate.	8.10
38	Report discrimination with confidence intervals, and diagnose miscalibration using calibration-in-the-large and the calibration slope.	8.11
39	Construct split-conformal prediction intervals and verify their coverage empirically.	8.11
40	Evaluate subgroup performance with uncertainty, and recognize when subgroup claims are unsupportable.	8.11
41	Classify dataset shift as covariate, prior-probability, or concept shift, and choose an appropriate response.	8.12
42	Distinguish fairness as measurement, as performance, and as consequence, and test each.	8.12
43	Detect shortcut learning using metadata-only and label-permutation controls.	8.12
44	Describe the privacy-utility trade-off quantitatively using the Laplace mechanism.	8.12
45	Specify deployment monitoring, model updating, and the contents of a model card.	8.12
46	Build a reproducible project with configuration, seeds, provenance, assertions, and unit tests.	8.13
47	Interpret a learning curve to decide whether more data, better labels, or better features is the binding constraint.	8.13

Rebuilding the Real Data Layer

num_with_flag <- function(s) {
  isnum <- grepl("^-?[0-9]+(\\.[0-9]+)?$", s) & !is.na(s)
  list(value = ifelse(isnum, suppressWarnings(as.numeric(s)), NA_real_),
       flag  = ifelse(is.na(s), "missing", ifelse(isnum, "numeric", s)))
}
kits_url  <- "https://raw.githubusercontent.com/neheller/kits19/master/data/kits.json"
kits_path <- if (has_pkg("jsonlite")) bpad_fetch(kits_url, "kits.json") else NA_character_

build_real_cohort <- function(path) {
  j <- jsonlite::fromJSON(path, simplifyDataFrame = FALSE)
  pull <- function(...) {
    pth <- c(...)
    vapply(j, function(z) { v <- z
      for (p in pth) { v <- v[[p]]; if (is.null(v)) return(NA_character_) }
      if (length(v) != 1) return(NA_character_); as.character(v) }, character(1))
  }
  ep <- num_with_flag(pull("last_preop_egfr", "value"))
  eo <- num_with_flag(pull("last_postop_egfr", "value"))
  data.frame(
    case_id = pull("case_id"),
    age = suppressWarnings(as.numeric(pull("age_at_nephrectomy"))),
    gender = pull("gender"),
    bmi = suppressWarnings(as.numeric(pull("body_mass_index"))),
    radiographic_size_cm = suppressWarnings(as.numeric(pull("radiographic_size"))),
    pathologic_size_cm   = suppressWarnings(as.numeric(pull("pathologic_size"))),
    malignant = pull("malignant") == "TRUE",
    procedure = pull("surgical_procedure"),   # extent: partial vs radical
    access    = pull("surgery_type"),         # open / laparoscopic / robotic
    egfr_pre  = ep$value, egfr_pre_flag = ep$flag,
    egfr_post = eo$value, egfr_post_flag = eo$flag,
    vital_status = pull("vital_status"),
    vital_days   = suppressWarnings(as.numeric(pull("vital_days_after_surgery"))),
    spacing_x_mm = suppressWarnings(as.numeric(pull("voxel_spacing", "x_spacing"))),
    spacing_z_mm = suppressWarnings(as.numeric(pull("voxel_spacing", "z_spacing"))),
    data_source = "KiTS19 (real, publicly released)", stringsAsFactors = FALSE)
}

build_surrogate_cohort <- function(n = 210, seed = 19) {
  set.seed(seed)
  procedure <- sample(c("partial_nephrectomy","radical_nephrectomy"), n, TRUE, c(140,70)/n)
  size <- pmax(0.6, rgamma(n, 3.1, scale = 1.45) + 2.4*(procedure == "radical_nephrectomy"))
  age  <- pmin(90, pmax(6, round(rnorm(n, 60, 13))))
  egfr <- pmin(120, pmax(15, round(rnorm(n, 76, 22) - 0.25*(age - 60))))
  drop <- pmax(-30, round(0.35*egfr - 20 + 21*(procedure=="radical_nephrectomy") + rnorm(n,0,12)))
  ep <- num_with_flag(ifelse(egfr >= 90, ">=90", as.character(egfr)))
  eo <- num_with_flag(ifelse(pmax(5, egfr-drop) >= 90, ">=90", as.character(pmax(5, egfr-drop))))
  tt <- rexp(n, 0.02*exp(0.02*(age-60) + 0.16*size)); cc <- runif(n, 0.05, 9)
  data.frame(case_id = sprintf("case_%05d", seq_len(n)-1L), age = age,
    gender = sample(c("male","female"), n, TRUE, c(.62,.38)),
    bmi = round(pmax(16, rnorm(n, 29.5, 6.4)), 1),
    radiographic_size_cm = round(size, 1),
    pathologic_size_cm = round(pmax(0.3, size + rnorm(n, -0.05, 1.17)), 1),
    malignant = runif(n) < inv_logit(-0.4 + 0.55*size),
    procedure = procedure,
    access = sample(c("robotic","open","laparoscopic"), n, TRUE, c(122,60,28)/n),
    egfr_pre = ep$value, egfr_pre_flag = ep$flag,
    egfr_post = eo$value, egfr_post_flag = eo$flag,
    vital_status = ifelse(tt <= cc, "dead", "censored"),
    vital_days = round(pmin(tt, cc)*365.25),
    spacing_x_mm = round(runif(n, .6, 1), 4),
    spacing_z_mm = sample(c(.5,1,1.25,2,3,5), n, TRUE, c(20,29,60,55,30,16)/n),
    data_source = "SURROGATE (network unavailable; schema and pathologies preserved)",
    stringsAsFactors = FALSE)
}

kits <- if (!is.na(kits_path)) build_real_cohort(kits_path) else build_surrogate_cohort()
kits$event     <- as.integer(kits$vital_status == "dead")
kits$time_yr   <- kits$vital_days/365.25
kits$radical   <- as.integer(kits$procedure == "radical_nephrectomy")
kits$open      <- as.integer(kits$access == "open")
kits$male      <- as.integer(kits$gender == "male")
kits$log_size  <- log(pmax(kits$radiographic_size_cm, 0.1))
kits$egfr_drop <- kits$egfr_pre - kits$egfr_post
kits$acq_group <- ifelse(kits$spacing_z_mm <= 1, "thin_slice", "thick_slice")

cat("data source:", unique(kits$data_source), "\n")

## data source: KiTS19 (real, publicly released)

cat(sprintf("n = %d, malignant prevalence = %.3f (%d benign), deaths = %d\n",
            nrow(kits), mean(kits$malignant), sum(!kits$malignant), sum(kits$event)))

## n = 210, malignant prevalence = 0.914 (18 benign), deaths = 21

## The Part 1 malignancy model, refitted here so Part 2 is self-contained.
mvars <- c("log_size", "age", "bmi", "radical", "male")
D  <- kits[complete.cases(kits[, c(mvars, "malignant")]), ]
y  <- as.numeric(D$malignant)
folds <- make_folds(nrow(D), 5, strata = D$malignant, seed = 11)
oof <- rep(NA_real_, nrow(D))
for (k in seq_len(5)) {
  tr <- folds != k
  f  <- suppressWarnings(glm(malignant ~ log_size + age + bmi + radical + male,
                             D[tr, ], family = binomial()))
  oof[!tr] <- predict(f, D[!tr, ], type = "response")
}
cat(sprintf("reference model: n = %d, benign = %d, out-of-fold AUC = %.3f, Brier = %.4f\n",
            nrow(D), sum(!y), auc_rank(oof, y), brier(oof, y)))

## reference model: n = 210, benign = 18, out-of-fold AUC = 0.698, Brier = 0.0744

8.9 Unsupervised Learning, Segmentation, and Spatial Structure

8.9.1 What unsupervised structure can and cannot establish

Unsupervised methods look for structure in $X$ with no outcome label. They can compress data, reveal acquisition artefacts, propose candidate phenotypes, and initialize a segmentation. They cannot, on their own, establish that a cluster is a disease subtype or that a component is biologically meaningful.

For $k$-means the objective is

\[\begin{equation} \min_{C_1,\ldots,C_K}\ \sum_{k=1}^{K}\sum_{i\in C_k}\lVert x_i - \mu_k\rVert_2^2, \qquad \mu_k = \frac{1}{|C_k|}\sum_{i\in C_k} x_i , \tag{1} \end{equation}\]

a non-convex problem solved by alternating assignment and update steps (Lloyd’s algorithm), which converges to a local optimum that depends on initialization. Equation (1) also encodes three assumptions that are rarely stated: clusters are convex and roughly spherical, they are similarly sized, and distance is Euclidean, which makes the result entirely dependent on feature scaling. If tumour volume is in mm$^3$ and age in years, volume alone determines the partition.

The circular trap of unsupervised phenotyping. If you cluster patients without the outcome and then name a cluster “aggressive phenotype” because it happens to contain more events, you have used the outcome after all, in the naming step. To be scientifically defensible a cluster must be stable (reproducible under resampling or in a held-out sample), upstream (defined from variables that are causally or clinically prior to the outcome), and externally validated (its association with the outcome demonstrated in an independent cohort). Section 8.9.3 tests the first of these on the real cohort and Section 8.9.2 shows why the third matters.

8.9.2 Patient phenotyping in standardized feature space

cvars <- c("age", "bmi", "log_size", "radical", "egfr_pre")
Zc <- kits[complete.cases(kits[, cvars]), ]
X  <- scale(as.matrix(Zc[, cvars]))
cat(sprintf("clustering matrix: %d patients x %d standardized features\n", nrow(X), ncol(X)))

## clustering matrix: 120 patients x 5 standardized features

set.seed(8)
km3 <- kmeans(X, centers = 3, nstart = 25)
pc  <- prcomp(X)

p1 <- ggplot(data.frame(PC1 = pc$x[,1], PC2 = pc$x[,2],
                        cl = factor(km3$cluster)), aes(PC1, PC2, colour = cl)) +
  geom_point(alpha = 0.75, size = 1.9) +
  scale_colour_manual(values = bpad_pal[1:3]) +
  labs(title = "clusters in PC space")

prof <- do.call(rbind, lapply(cvars, function(v)
  data.frame(feature = v, cluster = factor(1:3),
             value = as.numeric(tapply(Zc[[v]], km3$cluster, median)))))
prof$scaled <- ave(prof$value, prof$feature,
                   FUN = function(z) (z - min(z))/(max(z) - min(z) + 1e-9))
p2 <- ggplot(prof, aes(cluster, feature, fill = scaled)) +
  geom_tile(colour = "white") +
  geom_text(aes(label = round(value, 1)), size = 2.8) +
  scale_fill_gradient(low = "white", high = bpad_pal[1]) +
  theme(legend.position = "none") +
  labs(title = "cluster medians (original scale)", x = "cluster", y = NULL)

comp <- data.frame(cluster = factor(km3$cluster),
                   status = ifelse(Zc$malignant, "malignant", "benign"))
p3 <- ggplot(comp, aes(cluster, fill = status)) +
  geom_bar(position = "fill", alpha = 0.85) +
  scale_fill_manual(values = c(malignant = bpad_pal[2], benign = bpad_pal[1])) +
  labs(title = "composition by malignancy", x = "cluster", y = "proportion")
bpad_grid(p1, p2, p3, ncol = 3)

Unsupervised phenotyping of the real cohort in standardized feature space. Left: patients in the first two principal components, coloured by k-means cluster. Centre: cluster profiles on the original scale. Right: composition of each cluster by malignancy. The clusters are well separated in feature space and carry essentially no information about the outcome: the malignant proportion is nearly identical in all three. What they do separate, almost perfectly, is the surgical decision.

Figure 1: Unsupervised phenotyping of the real cohort in standardized feature space. Left: patients in the first two principal components, coloured by k-means cluster. Centre: cluster profiles on the original scale. Right: composition of each cluster by malignancy. The clusters are well separated in feature space and carry essentially no information about the outcome: the malignant proportion is nearly identical in all three. What they do separate, almost perfectly, is the surgical decision.

print(knitr::kable(data.frame(
  cluster              = 1:3,
  n                    = as.vector(table(km3$cluster)),
  malignant_proportion = round(as.vector(tapply(Zc$malignant, km3$cluster, mean)), 3),
  median_size_cm       = round(as.vector(tapply(Zc$radiographic_size_cm, km3$cluster, median)), 2),
  radical_fraction     = round(as.vector(tapply(Zc$radical, km3$cluster, mean)), 3),
  median_age           = round(as.vector(tapply(Zc$age, km3$cluster, median)), 0)),
  row.names = FALSE,
  caption = "What the clusters actually separate. Malignancy is nearly constant across clusters; operative extent and tumour size are almost perfectly separated."))

## 
## 
## Table: (\#tab:patient-clustering)What the clusters actually separate. Malignancy is nearly constant across clusters; operative extent and tumour size are almost perfectly separated.
## 
## | cluster|  n| malignant_proportion| median_size_cm| radical_fraction| median_age|
## |-------:|--:|--------------------:|--------------:|----------------:|----------:|
## |       1| 45|                0.889|            6.3|            1.000|         66|
## |       2| 30|                0.833|            3.0|            0.067|         68|
## |       3| 45|                0.889|            2.7|            0.022|         56|

The clusters found the surgeon, not the disease. Malignant proportion is essentially identical across the three clusters, so a supervised model would gain nothing from the cluster label. What the partition recovers almost perfectly is operative extent, and with it tumour size. This is the same lesson as the PCA figure in Part 1, Section 8.6.4, in a different algorithm: the dominant directions of variance in a clinical cohort are usually the directions of clinical decision-making, not of biology. A phenotype paper that stopped one figure earlier would have reported three “novel renal-mass phenotypes” that are, on inspection, partial nephrectomy, small radical, and large radical.

8.9.3 Internal validation: silhouette, gap statistic, and stability

Three questions must be answered before a partition is reported: is it cohesive, is the number of clusters supported, and is it reproducible?

The silhouette width for observation $i$ compares its mean distance to its own cluster with its mean distance to the nearest other cluster:

\[\begin{equation} s(i) = \frac{b(i) - a(i)}{\max\{a(i),\,b(i)\}} \in [-1, 1], \qquad a(i) = \frac{1}{|C_{k(i)}|-1}\sum_{j\in C_{k(i)},\, j\ne i} d_{ij}, \quad b(i) = \min_{k\ne k(i)}\frac{1}{|C_k|}\sum_{j\in C_k} d_{ij}. \tag{2} \end{equation}\]

The gap statistic compares the observed within-cluster dispersion $W_k$ with its expectation under a null reference distribution of uniform points on the data’s bounding box:

\[\begin{equation} \mathrm{Gap}(k) = \mathbb{E}^{\ast}\{\log W_k\} - \log W_k , \tag{3} \end{equation}\]

and the conventional choice is the smallest $k$ satisfying $\mathrm{Gap}(k) \ge \mathrm{Gap}(k+1) - s_{k+1}$, the one-standard-error rule.

Stability asks a different question again: does the same partition reappear when the patients are resampled? We measure this with the adjusted Rand index between the reference partition and a partition fitted to a bootstrap resample.

sil_width <- function(X, cl) {
  Dm <- as.matrix(dist(X)); n <- nrow(X); s <- numeric(n)
  for (i in seq_len(n)) {
    own <- cl == cl[i]; own[i] <- FALSE
    a <- if (sum(own)) mean(Dm[i, own]) else 0
    others <- setdiff(unique(cl), cl[i])
    b <- min(vapply(others, function(k) mean(Dm[i, cl == k]), numeric(1)))
    s[i] <- if (max(a, b) > 0) (b - a)/max(a, b) else 0
  }
  s
}
adj_rand <- function(a, b) {
  tb <- table(a, b); n <- sum(tb)
  si <- sum(choose(tb, 2)); sa <- sum(choose(rowSums(tb), 2)); sb <- sum(choose(colSums(tb), 2))
  ex <- sa*sb/choose(n, 2); mx <- (sa + sb)/2
  (si - ex)/(mx - ex)
}

set.seed(8)
ks <- 2:6
val <- t(vapply(ks, function(k) {
  km <- kmeans(X, k, nstart = 25)
  c(wss = km$tot.withinss, sil = mean(sil_width(X, km$cluster)))
}, numeric(2)))

gap_statistic <- function(X, kmax = 6, B = 30) {
  set.seed(3); rng <- apply(X, 2, range)
  logW <- vapply(seq_len(kmax), function(k) log(kmeans(X, k, nstart = 10)$tot.withinss),
                 numeric(1))
  Wb <- vapply(seq_len(B), function(b) {
    Xb <- vapply(seq_len(ncol(X)), function(j) runif(nrow(X), rng[1,j], rng[2,j]),
                 numeric(nrow(X)))
    vapply(seq_len(kmax), function(k) log(kmeans(Xb, k, nstart = 10)$tot.withinss), numeric(1))
  }, numeric(kmax))
  data.frame(k = seq_len(kmax), gap = rowMeans(Wb) - logW,
             se = apply(Wb, 1, sd)*sqrt(1 + 1/B))
}
g <- gap_statistic(X)
k_gap <- which(g$gap[-nrow(g)] >= g$gap[-1] - g$se[-1])[1]

stability <- vapply(ks, function(k) {
  base <- kmeans(X, k, nstart = 25)$cluster
  r <- replicate(60, {
    idx <- sample(nrow(X), replace = TRUE)
    adj_rand(base[idx], kmeans(X[idx, ], k, nstart = 10)$cluster)
  })
  c(mean = mean(r), sd = sd(r))
}, numeric(2))

print(knitr::kable(data.frame(
  k = ks, within_SS = round(val[, "wss"], 1),
  mean_silhouette = round(val[, "sil"], 3),
  gap = round(g$gap[ks], 3),
  stability_adjRand = round(stability["mean", ], 3),
  stability_sd = round(stability["sd", ], 3)),
  row.names = FALSE,
  caption = "Three internal criteria and a stability criterion. Silhouette is maximized at the smallest k considered, the gap one-standard-error rule selects the same value, and bootstrap stability falls sharply thereafter."))

## 
## 
## Table: (\#tab:cluster-validation)Three internal criteria and a stability criterion. Silhouette is maximized at the smallest k considered, the gap one-standard-error rule selects the same value, and bootstrap stability falls sharply thereafter.
## 
## |  k| within_SS| mean_silhouette|   gap| stability_adjRand| stability_sd|
## |--:|---------:|---------------:|-----:|-----------------:|------------:|
## |  2|     415.6|           0.307| 0.680|             0.946|        0.152|
## |  3|     340.5|           0.252| 0.682|             0.712|        0.194|
## |  4|     285.0|           0.255| 0.693|             0.741|        0.147|
## |  5|     252.2|           0.232| 0.683|             0.713|        0.146|
## |  6|     225.9|           0.236| 0.680|             0.697|        0.116|

cat(sprintf("gap one-standard-error choice: k = %d\n", g$k[k_gap]))

## gap one-standard-error choice: k = 2

q1 <- ggplot(data.frame(k = ks, wss = val[, "wss"], sil = val[, "sil"]), aes(k)) +
  geom_line(aes(y = wss/max(wss), colour = "within-cluster SS (scaled)"), linewidth = 0.9) +
  geom_line(aes(y = sil/max(sil), colour = "mean silhouette (scaled)"), linewidth = 0.9) +
  geom_point(aes(y = wss/max(wss)), size = 1.7) + geom_point(aes(y = sil/max(sil)), size = 1.7) +
  scale_colour_manual(values = bpad_pal[1:2]) +
  labs(title = "elbow and silhouette", x = "k", y = "scaled criterion")
q2 <- ggplot(g, aes(k, gap)) +
  geom_errorbar(aes(ymin = gap - se, ymax = gap + se), width = 0.14, colour = bpad_pal[1]) +
  geom_line(colour = bpad_pal[1], linewidth = 0.9) + geom_point(size = 2) +
  labs(title = "gap statistic", x = "k", y = "Gap(k)")
q3 <- ggplot(data.frame(k = ks, m = stability["mean", ], s = stability["sd", ]), aes(k, m)) +
  geom_errorbar(aes(ymin = m - s, ymax = m + s), width = 0.14, colour = bpad_pal[3]) +
  geom_line(colour = bpad_pal[3], linewidth = 0.9) + geom_point(size = 2) +
  coord_cartesian(ylim = c(0, 1.05)) +
  labs(title = "bootstrap stability", x = "k", y = "adjusted Rand index")
bpad_grid(q1, q2, q3, ncol = 3)

$Internal validation of the patient partition. Left: within-cluster sum of squares (the elbow) and mean silhouette width against k. Centre: the gap statistic with one-standard-error bars, Eq. \@ref(eq:gap). Right: partition stability under bootstrap resampling, measured by the adjusted Rand index. All three criteria agree that only k = 2 is well supported, and stability collapses sharply beyond it.$

Figure 2: Internal validation of the patient partition. Left: within-cluster sum of squares (the elbow) and mean silhouette width against k. Centre: the gap statistic with one-standard-error bars, Eq. (3). Right: partition stability under bootstrap resampling, measured by the adjusted Rand index. All three criteria agree that only k = 2 is well supported, and stability collapses sharply beyond it.

8.9.4 Unsupervised intensity segmentation

The same objective, Eq. (1), applied to voxel intensities rather than patients, gives the simplest possible segmentation. It has no notion of spatial contiguity, which is both its weakness and a useful demonstration: any spatial coherence in the result comes from the anatomy, not the algorithm.

set.seed(21)
n <- 110
gx <- outer(1:n, 1:n, function(i, j) i); gy <- outer(1:n, 1:n, function(i, j) j)
img <- matrix(20, n, n)                                    # background soft tissue
img[(gx - 45)^2/1.4 + (gy - 55)^2 < 620] <- 110            # enhancing parenchyma
img[(gx - 68)^2 + (gy - 62)^2 < 150]     <- 65             # lesion
img <- img + matrix(rnorm(n*n, 0, 7), n)                   # detector noise
bias <- 1 + 1.20*(gx/n)                                    # smooth multiplicative gain
img_b <- img*bias

seg_km <- function(v, k, seed = 1) { set.seed(seed)
  km <- kmeans(v, centers = k, nstart = 20)
  ord <- order(km$centers); matrix(match(km$cluster, ord), sqrt(length(v)))
}
S  <- seg_km(as.numeric(img),   3)
Sb <- seg_km(as.numeric(img_b), 3)
km_plain <- kmeans(as.numeric(img), 3, nstart = 20)
cuts <- sort(km_plain$centers)
bnds <- (cuts[-1] + cuts[-length(cuts)])/2

op <- par(mfrow = c(2, 3), mar = c(4, 3.4, 2.6, 1))
image(t(img[n:1, ]), col = grey.colors(64), axes = FALSE, asp = 1,
      main = "synthetic CT slice"); box(col = "grey70")
hist(as.numeric(img), breaks = 60, col = "grey85", border = "grey60",
     main = "intensity histogram", xlab = "intensity")
abline(v = bnds, col = bpad_pal[2], lwd = 2, lty = 2)
image(t(S[n:1, ]), col = bpad_pal[c(7,1,2)], axes = FALSE, asp = 1,
      main = "3-class k-means"); box(col = "grey70")
image(t(img_b[n:1, ]), col = grey.colors(64), axes = FALSE, asp = 1,
      main = "with a smooth bias field"); box(col = "grey70")
hist(as.numeric(img_b), breaks = 60, col = "grey85", border = "grey60",
     main = "histogram after bias", xlab = "intensity")
image(t(Sb[n:1, ]), col = bpad_pal[c(7,1,2)], axes = FALSE, asp = 1,
      main = "segmentation now fails"); box(col = "grey70")

Unsupervised intensity segmentation of a synthetic image built to mimic contrast-enhanced abdominal CT, with Poisson-like noise and a smooth intensity gradient. Top: the image, its intensity histogram with the k-means boundaries, and the three-class segmentation. Bottom: the same segmentation after a bias field is added, showing how a smoothly varying gain, which the eye ignores, destroys a purely intensity-based partition. Spatial regularization or bias-field correction, both covered in Chapter 7, are required before intensity clustering is usable.

Figure 3: Unsupervised intensity segmentation of a synthetic image built to mimic contrast-enhanced abdominal CT, with Poisson-like noise and a smooth intensity gradient. Top: the image, its intensity histogram with the k-means boundaries, and the three-class segmentation. Bottom: the same segmentation after a bias field is added, showing how a smoothly varying gain, which the eye ignores, destroys a purely intensity-based partition. Spatial regularization or bias-field correction, both covered in Chapter 7, are required before intensity clustering is usable.

par(op)

agree <- mean(S == Sb)
cat(sprintf("voxel agreement between the two segmentations: %.3f\n", agree))

## voxel agreement between the two segmentations: 0.665

cat(sprintf("bias field spans a factor of %.2f across the image (a %.0f%% gain change)\n",
            max(bias)/min(bias), 100*(max(bias)/min(bias) - 1)))

## bias field spans a factor of 2.18 across the image (a 118% gain change)

cat(sprintf("%.1f%% of voxels are reassigned to a different class by the gain alone.\n",
            100*(1 - agree)))

## 33.5% of voxels are reassigned to a different class by the gain alone.

cat("The eye compensates for a smooth gain automatically; an intensity histogram does not.\n")

## The eye compensates for a smooth gain automatically; an intensity histogram does not.

8.9.5 Overlap and boundary metrics: two failures in opposite directions

For a predicted mask $P$ and reference mask $G$,

\[\begin{equation} \mathrm{Dice}(P,G) = \frac{2\lvert P\cap G\rvert}{\lvert P\rvert + \lvert G\rvert}, \qquad \mathrm{Jaccard}(P,G) = \frac{\lvert P\cap G\rvert}{\lvert P\cup G\rvert}, \qquad \mathrm{Dice} = \frac{2J}{1+J}. \tag{4} \end{equation}\]

These are overlap metrics: they count voxels. Boundary agreement is measured instead by surface distances. With $\partial P$ and $\partial G$ the boundary sets and $d(a, S) = \min_{b\in S}\lVert a-b\rVert$,

\[\begin{equation} \mathrm{HD}(P,G) = \max\Big\{\max_{a\in\partial P} d(a,\partial G),\ \max_{b\in\partial G} d(b,\partial P)\Big\}, \tag{5} \end{equation}\]

\[\begin{equation} \mathrm{ASSD}(P,G) = \frac{\sum_{a\in\partial P} d(a,\partial G) + \sum_{b\in\partial G} d(b,\partial P)} {\lvert\partial P\rvert + \lvert\partial G\rvert}. \tag{6} \end{equation}\]

The Hausdorff distance is a maximum and is therefore dominated by the single worst point, which is why the 95th-percentile variant $\mathrm{HD}_{95}$ is normally reported alongside it.

The demonstration below takes a real expert kidney mask from one axial plane of a KiTS19 reference segmentation and perturbs it in two structurally different ways: a dilation, which inflates volume while keeping boundaries close, and a translation, which preserves volume exactly while displacing every boundary point.

## ---------------------------------------------------------------------------
## Obtain a real axial label plane from a KiTS19 expert segmentation by
## streaming, so no large volume is held in memory. Falls back to a synthetic
## kidney-shaped mask if the volume cannot be fetched.
## ---------------------------------------------------------------------------
read_nifti_header <- function(path) {
  con <- if (grepl("\\.gz$", path)) gzfile(path, "rb") else file(path, "rb")
  on.exit(close(con)); hdr <- readBin(con, "raw", n = 348)
  endian <- "little"; sz <- readBin(hdr[1:4], "integer", size = 4, endian = "little")
  if (!identical(sz, 348L)) { sz <- readBin(hdr[1:4], "integer", size = 4, endian = "big")
                              endian <- "big" }
  if (!identical(sz, 348L)) stop("not a NIfTI-1 file")
  gi <- function(o, n, s) readBin(hdr[(o+1):(o+n*s)], "integer", n = n, size = s, endian = endian)
  gd <- function(o, n)    readBin(hdr[(o+1):(o+n*4)], "double",  n = n, size = 4, endian = endian)
  d <- gi(40, 8, 2); pd <- gd(76, 8)
  list(dim = d[2:(1+d[1])], datatype = gi(70,1,2), pixdim = pd[2:(1+d[1])],
       vox_offset = gd(108,1), endian = endian)
}
nifti_type <- function(dt) switch(as.character(dt),
  "2" = list(what="integer", size=1, signed=FALSE), "4" = list(what="integer", size=2, signed=TRUE),
  "8" = list(what="integer", size=4, signed=TRUE), "16"= list(what="double", size=4, signed=TRUE),
  "64"= list(what="double", size=8, signed=TRUE), "512"=list(what="integer",size=2,signed=FALSE),
  stop("unsupported datatype ", dt))

axial_label_plane <- function(path) {
  h <- read_nifti_header(path); tp <- nifti_type(h$datatype)
  d <- h$dim[1:3]; nx <- d[1]; ny <- d[2]; nz <- d[3]
  open_it <- function() { cn <- if (grepl("\\.gz$", path)) gzfile(path,"rb") else file(path,"rb")
                          invisible(readBin(cn, "raw", n = h$vox_offset)); cn }
  cn <- open_it(); prof <- numeric(nx)
  for (k in seq_len(nz)) {
    m <- matrix(readBin(cn, tp$what, n = nx*ny, size = tp$size, signed = tp$signed), nx, ny)
    prof <- prof + rowSums(m == 2)                     # tumour count per axial index
  }
  close(cn); i0 <- which.max(prof)
  cn <- open_it(); A <- matrix(0L, ny, nz)
  for (k in seq_len(nz)) {
    m <- matrix(readBin(cn, tp$what, n = nx*ny, size = tp$size, signed = tp$signed), nx, ny)
    A[, k] <- m[i0, ]
  }
  close(cn)
  list(plane = A, spacing = h$pixdim[2:3], index = i0)
}

seg_url  <- "https://raw.githubusercontent.com/neheller/kits19/master/data/case_00000/segmentation.nii.gz"
plane_rds <- file.path(bpad_cache_dir, "axial_plane.rds")
if (file.exists(plane_rds)) {
  ax <- readRDS(plane_rds)
} else {
  sp_path <- bpad_fetch(seg_url, "seg_00000.nii.gz")
  ax <- if (!is.na(sp_path)) {
    a <- axial_label_plane(sp_path)
    r1 <- range(which(rowSums(a$plane > 0) > 0)); r2 <- range(which(colSums(a$plane > 0) > 0))
    pad <- 12
    a$plane <- a$plane[max(1, r1[1]-pad):min(nrow(a$plane), r1[2]+pad),
                       max(1, r2[1]-pad):min(ncol(a$plane), r2[2]+pad)]
    a$source <- "real KiTS19 expert segmentation, case_00000"
    a
  } else {
    g <- expand.grid(i = 1:90, j = 1:80)
    pl <- matrix(0L, 90, 80)
    pl[((g$i-45)^2/2.2 + (g$j-40)^2/1.1) < 520] <- 1L
    pl[((g$i-58)^2 + (g$j-50)^2) < 130] <- 2L
    list(plane = pl, spacing = c(0.92, 0.92), index = NA,
         source = "synthetic kidney mask (real volume unavailable)")
  }
  saveRDS(ax, plane_rds)
}
cat("mask source:", ax$source, "\n")

## mask source: real KiTS19 expert segmentation, case_00000

ref <- ax$plane == 1; sp <- ax$spacing
cat(sprintf("plane %d x %d, in-plane spacing %.3f x %.3f mm, kidney pixels %d, tumour pixels %d\n",
            nrow(ax$plane), ncol(ax$plane), sp[1], sp[2], sum(ref), sum(ax$plane == 2)))

## plane 87 x 77, in-plane spacing 0.920 x 0.920 mm, kidney pixels 1905, tumour pixels 537

surface_distances <- function(A, B, spacing) {
  ia <- which(boundary2d(A), arr.ind = TRUE); ib <- which(boundary2d(B), arr.ind = TRUE)
  if (!nrow(ia) || !nrow(ib)) return(c(hd = NA, hd95 = NA, assd = NA))
  ia <- sweep(ia, 2, spacing, "*"); ib <- sweep(ib, 2, spacing, "*")
  dab <- apply(ia, 1, function(p) sqrt(min(colSums((t(ib) - p)^2))))
  dba <- apply(ib, 1, function(p) sqrt(min(colSums((t(ia) - p)^2))))
  c(hd   = max(max(dab), max(dba)),
    hd95 = max(stats::quantile(dab, 0.95), stats::quantile(dba, 0.95)),
    assd = (sum(dab) + sum(dba))/(length(dab) + length(dba)))
}
shift_mask <- function(B, k) { s <- B; s[] <- FALSE
  s[(k+1):nrow(B), ] <- B[1:(nrow(B)-k), ]; s }

variants <- list("dilate 1 px" = dilate2d(ref),
                 "erode 1 px"  = erode2d(ref),
                 "dilate 2 px" = dilate2d(dilate2d(ref)),
                 "shift 2 px"  = shift_mask(ref, 2),
                 "shift 4 px"  = shift_mask(ref, 4))
seg_tab <- do.call(rbind, lapply(names(variants), function(nm) {
  v <- variants[[nm]]; s <- surface_distances(ref, v, sp)
  data.frame(perturbation = nm,
             Dice = round(dice_coef(ref, v), 3), Jaccard = round(jaccard(ref, v), 3),
             HD_mm = round(s["hd"], 2), HD95_mm = round(s["hd95"], 2),
             ASSD_mm = round(s["assd"], 3),
             area_change_pct = round(100*(sum(v) - sum(ref))/sum(ref), 1))
}))
print(knitr::kable(seg_tab, row.names = FALSE,
  caption = paste("Overlap and boundary metrics on a real expert kidney mask.",
                  "The translations change no area at all yet displace every boundary",
                  "point; the dilations inflate area while keeping boundaries close.")))

## 
## 
## Table: (\#tab:segmentation-metrics)Overlap and boundary metrics on a real expert kidney mask. The translations change no area at all yet displace every boundary point; the dilations inflate area while keeping boundaries close.
## 
## |perturbation |  Dice| Jaccard| HD_mm| HD95_mm| ASSD_mm| area_change_pct|
## |:------------|-----:|-------:|-----:|-------:|-------:|---------------:|
## |dilate 1 px  | 0.953|   0.909|  0.92|    0.92|   0.920|            10.0|
## |erode 1 px   | 0.949|   0.903|  3.32|    0.92|   0.946|            -9.7|
## |dilate 2 px  | 0.908|   0.832|  2.06|    1.84|   1.613|            20.2|
## |shift 2 px   | 0.936|   0.881|  1.84|    1.84|   1.029|             0.0|
## |shift 4 px   | 0.876|   0.780|  3.68|    3.68|   1.983|             0.0|

## visual: boundary overlays plus the metric plane
bnd_df <- function(M, lab) { idx <- which(boundary2d(M), arr.ind = TRUE)
  data.frame(x = idx[,1], y = idx[,2], what = lab) }
ov1 <- rbind(bnd_df(ref, "reference"), bnd_df(variants[["dilate 2 px"]], "dilate 2 px"))
ov2 <- rbind(bnd_df(ref, "reference"), bnd_df(variants[["shift 4 px"]],  "shift 4 px"))
v1 <- ggplot(ov1, aes(x, y, colour = what)) + geom_point(size = 0.55) + coord_equal() +
  scale_colour_manual(values = c(reference = "grey55", `dilate 2 px` = bpad_pal[2])) +
  labs(title = "dilation: area grows, boundary stays near", x = NULL, y = NULL) +
  theme(axis.text = element_blank())
v2 <- ggplot(ov2, aes(x, y, colour = what)) + geom_point(size = 0.55) + coord_equal() +
  scale_colour_manual(values = c(reference = "grey55", `shift 4 px` = bpad_pal[1])) +
  labs(title = "translation: area identical, boundary moves", x = NULL, y = NULL) +
  theme(axis.text = element_blank())
v3 <- ggplot(seg_tab, aes(ASSD_mm, Dice)) +
  geom_point(aes(size = abs(area_change_pct), colour = perturbation)) +
  geom_text(aes(label = perturbation), size = 2.5, vjust = -1.1) +
  scale_colour_manual(values = bpad_pal[1:5], guide = "none") +
  scale_size_continuous(range = c(2, 7), name = "|area change| %") +
  coord_cartesian(ylim = c(0.83, 1.0)) +
  labs(title = "Dice versus boundary distance",
       x = "mean symmetric surface distance (mm)", y = "Dice")
bpad_grid(v1, v2, v3, ncol = 3)

Overlap and boundary metrics disagree by construction. A real expert kidney mask from one axial plane is perturbed by dilation, erosion, and rigid translation. Left: the reference boundary in grey with the perturbed boundary overlaid, for a dilation and for a four-pixel shift. Right: Dice against mean symmetric surface distance for all perturbations, with point size proportional to the absolute area change. The dilation and the shift have similar Dice and opposite failure modes: one changes volume without moving the boundary far, the other moves the whole boundary without changing volume at all.

Figure 4: Overlap and boundary metrics disagree by construction. A real expert kidney mask from one axial plane is perturbed by dilation, erosion, and rigid translation. Left: the reference boundary in grey with the perturbed boundary overlaid, for a dilation and for a four-pixel shift. Right: Dice against mean symmetric surface distance for all perturbations, with point size proportional to the absolute area change. The dilation and the shift have similar Dice and opposite failure modes: one changes volume without moving the boundary far, the other moves the whole boundary without changing volume at all.

Why a single segmentation metric is never enough. In radiation-therapy planning, a two-pixel boundary displacement can preserve Dice above 0.93 while shifting the dose gradient outside the planning target volume; a uniform dilation can produce the same Dice while overdosing an adjacent organ at risk. The two errors have the same overlap score and opposite clinical consequences. Contemporary imaging-biomarker guidance therefore requires complementary metrics: an overlap measure (Dice or Jaccard), a boundary measure ($\mathrm{HD}_{95}$ or ASSD), and a clinically anchored measure (volumetric agreement, lesion detection rate). Part 1, Section 8.5.6 supplies the matching quantity for a volumetric biomarker: the one-voxel perturbation band.

8.9.6 Hierarchical clustering and mixture models

Hierarchical clustering builds a dendrogram rather than a flat partition, and it does not require $k$ in advance. The linkage rule determines the geometry it can recover: single linkage (nearest neighbour) can follow elongated, chained structures but is sensitive to noise bridges; complete linkage (farthest neighbour) produces compact, roughly spherical clusters; Ward’s method minimizes the increase in within-cluster variance and behaves most like $k$-means.

Gaussian mixture models replace hard assignment with a probabilistic one,

\[\begin{equation} p(x) = \sum_{k=1}^{K} \pi_k\, \mathcal{N}(x \mid \mu_k, \Sigma_k), \qquad \sum_k \pi_k = 1, \tag{7} \end{equation}\]

fitted by expectation-maximization. Because each component carries its own covariance $\Sigma_k$, a GMM can represent elongated and differently oriented clusters that $k$-means cannot, and it returns a posterior probability of membership rather than a hard label, which is far more honest for a clinical phenotype.

Dm <- dist(X)
links <- c("single", "complete", "average", "ward.D2")
hc_res <- lapply(links, function(L) cutree(hclust(Dm, method = L), k = 3))
names(hc_res) <- links

size_tab <- do.call(rbind, lapply(links, function(L) {
  s <- sort(as.vector(table(hc_res[[L]])), decreasing = TRUE)
  data.frame(linkage = L, largest = s[1], middle = s[2], smallest = s[length(s)],
             agreement_with_kmeans = round(adj_rand(hc_res[[L]], km3$cluster), 3))
}))
print(knitr::kable(size_tab, row.names = FALSE,
  caption = "Cluster sizes at k = 3 under four linkage rules, and agreement with the k-means partition."))

## 
## 
## Table: (\#tab:hierarchical-clustering)Cluster sizes at k = 3 under four linkage rules, and agreement with the k-means partition.
## 
## |linkage  | largest| middle| smallest| agreement_with_kmeans|
## |:--------|-------:|------:|--------:|---------------------:|
## |single   |     118|      1|        1|                 0.003|
## |complete |      66|     43|       11|                 0.147|
## |average  |     118|      1|        1|                 0.010|
## |ward.D2  |      48|     41|       31|                 0.821|

sz <- do.call(rbind, lapply(links, function(L)
  data.frame(linkage = L, cluster = factor(1:3),
             n = as.vector(table(factor(hc_res[[L]], levels = 1:3))))))
h1 <- ggplot(sz, aes(linkage, n, fill = cluster)) +
  geom_col(position = "stack", alpha = 0.85) +
  scale_fill_manual(values = bpad_pal[1:3]) +
  theme(axis.text.x = element_text(angle = 25, hjust = 1)) +
  labs(title = "cluster sizes at k = 3", x = NULL, y = "patients")
h2 <- ggplot(size_tab, aes(reorder(linkage, agreement_with_kmeans), agreement_with_kmeans)) +
  geom_col(fill = bpad_pal[3], alpha = 0.85) + coord_flip() +
  labs(title = "agreement with k-means", x = NULL, y = "adjusted Rand index")
bpad_grid(h1, h2, ncol = 2)

Linkage determines what a hierarchical clustering can find. The same standardized patient matrix is clustered with three linkage rules and cut at k = 3. Left: cluster sizes by linkage; single linkage produces the classic chaining failure, assigning almost every patient to one cluster and isolating a few outliers. Right: agreement between each linkage solution and the k-means partition, measured by the adjusted Rand index. Ward linkage agrees most closely with k-means, as expected from their shared variance-minimizing objective.

Figure 5: Linkage determines what a hierarchical clustering can find. The same standardized patient matrix is clustered with three linkage rules and cut at k = 3. Left: cluster sizes by linkage; single linkage produces the classic chaining failure, assigning almost every patient to one cluster and isolating a few outliers. Right: agreement between each linkage solution and the k-means partition, measured by the adjusted Rand index. Ward linkage agrees most closely with k-means, as expected from their shared variance-minimizing objective.

8.9.7 Graph-based and spectral segmentation

Spectral methods treat pixels (or superpixels) as nodes $V = \{v_1,\ldots,v_N\}$ of a weighted graph. A standard choice couples spatial proximity with intensity similarity,

\[\begin{equation} w_{ij} = \exp\!\left(-\frac{\lVert x_i - x_j\rVert^2}{2\sigma_s^2}\right) \exp\!\left(-\frac{\lvert I_i - I_j\rvert^2}{2\sigma_I^2}\right), \qquad w_{ii} = 0, \tag{8} \end{equation}\]

so that only pixels that are both nearby and similar in intensity are strongly connected. With the degree matrix $D_{ii} = \sum_j w_{ij}$, the unnormalized graph Laplacian is

\[\begin{equation} L = D - W . \tag{9} \end{equation}\]

Two properties make $L$ useful. It is symmetric positive semidefinite, because for any $f\in\mathbb{R}^N$

\[\begin{equation} f^{\mathsf T} L f = \tfrac{1}{2}\sum_{i,j} w_{ij}(f_i - f_j)^2 \ \ge\ 0, \tag{10} \end{equation}\]

so its eigenvalues satisfy $0 = \lambda_1 \le \lambda_2 \le \cdots \le \lambda_N$. And the multiplicity of the zero eigenvalue equals the number of connected components, with $L\mathbf{1} = 0$ always. Equation (10) says that $f^{\mathsf T}Lf$ is a smoothness penalty: it is small when $f$ varies little across strongly connected pixels.

The eigenvector $v_2$ associated with the second-smallest eigenvalue is the Fiedler vector. Because it is orthogonal to $\mathbf{1}$ it must take both signs, and thresholding it at zero produces a bipartition that approximately minimizes the normalized cut

\[\begin{equation} \mathrm{Ncut}(S,\bar S) = \frac{\mathrm{cut}(S,\bar S)}{\mathrm{vol}(S)} + \frac{\mathrm{cut}(S,\bar S)}{\mathrm{vol}(\bar S)}, \qquad \mathrm{cut}(S,\bar S) = \sum_{i\in S, j\in\bar S} w_{ij}. \tag{11} \end{equation}\]

The normalized symmetric variant $L_{\mathrm{sym}} = I - D^{-1/2}WD^{-1/2}$ compensates for varying region sizes and is usually preferred in practice.

The physical intuition. Imagine the pixels connected by springs whose stiffness is $w_{ij}$: strong springs between neighbours of similar intensity, weak springs across an edge. Pluck the graph. The lowest vibrational mode is uniform translation, which is the constant eigenvector at $\lambda_1 = 0$ and carries no information. The second mode, at $\lambda_2$, is the fundamental mode: the graph divides into two halves oscillating against each other, and the nodal line of that mode is precisely the weakest set of springs, hence the best cut. This is the same eigenvalue reasoning used for normal modes in Chapter 1, applied to an image instead of a mechanical system.

## Build an image from the real label plane, with realistic contrast and noise.
set.seed(5)
img_sp <- (ax$plane > 0)*1.0 + 0.45*(ax$plane == 2)
img_sp <- img_sp + matrix(rnorm(length(img_sp), 0, 0.04), nrow(img_sp))

## Dense eigendecomposition is O(N^3), so subsample to a manageable node count.
patch <- img_sp[seq(1, nrow(img_sp), by = 3), seq(1, ncol(img_sp), by = 3)]
n1 <- nrow(patch); n2 <- ncol(patch); N <- n1*n2
cat(sprintf("graph: %d x %d patch = %d nodes, %s edges before sparsification\n",
            n1, n2, N, format(N*(N-1)/2, big.mark = ",")))

## graph: 29 x 26 patch = 754 nodes, 283,881 edges before sparsification

co <- expand.grid(r = seq_len(n1), c = seq_len(n2))
iv <- as.numeric(patch)
d_space <- as.matrix(dist(co)); d_int <- as.matrix(dist(iv))
sigma_s <- 3; sigma_I <- 0.25
W <- exp(-d_space^2/(2*sigma_s^2))*exp(-d_int^2/(2*sigma_I^2))
W[d_space > 3*sigma_s] <- 0; diag(W) <- 0
L <- diag(rowSums(W)) - W
ev <- eigen(L, symmetric = TRUE)
lam <- rev(ev$values)
fiedler <- ev$vectors[, N - 1]
part <- fiedler >= 0

ncut_value <- function(W, S) {
  cut <- sum(W[S, !S]); cut/sum(W[S, ]) + cut/sum(W[!S, !S])
}
set.seed(4); rand_part <- sample(c(TRUE, FALSE), N, TRUE)
truth <- as.numeric(patch) > 0.5
cat(sprintf("lambda_1 = %.3e (theory: exactly 0), lambda_2 = %.5f\n", lam[1], lam[2]))

## lambda_1 = 6.316e-14 (theory: exactly 0), lambda_2 = 0.00923

cat(sprintf("Fiedler partition sizes: %d / %d\n", sum(part), sum(!part)))

## Fiedler partition sizes: 487 / 267

cat(sprintf("agreement with anatomy: %.3f\n", max(mean(part == truth), mean(part != truth))))

## agreement with anatomy: 1.000

cat(sprintf("normalized cut: Fiedler %.5f vs random partition %.4f (%.0fx better)\n",
            ncut_value(W, part), ncut_value(W, rand_part),
            ncut_value(W, rand_part)/max(ncut_value(W, part), 1e-9)))

## normalized cut: Fiedler 0.00028 vs random partition 1.5271 (5530x better)

op <- par(mfrow = c(2, 3), mar = c(4, 3.6, 2.6, 1))
image(t(patch[n1:1, ]), col = grey.colors(64), axes = FALSE, asp = 1,
      main = "input patch from real mask"); box(col = "grey70")
plot(1:12, lam[1:12], type = "b", pch = 16, col = bpad_pal[1],
     xlab = "index", ylab = expression(lambda), main = "Laplacian spectrum")
abline(h = 0, col = "grey70", lty = 2)
points(2, lam[2], col = bpad_pal[2], pch = 16, cex = 1.5)
text(2, lam[2], "Fiedler value", pos = 4, cex = 0.75, col = bpad_pal[2])
image(t(matrix(fiedler, n1, n2)[n1:1, ]), col = hcl.colors(64, "Blue-Red 3"),
      axes = FALSE, asp = 1, main = "Fiedler vector"); box(col = "grey70")
image(t(matrix(as.numeric(part), n1, n2)[n1:1, ]), col = bpad_pal[c(7,1)],
      axes = FALSE, asp = 1, main = "spectral bipartition"); box(col = "grey70")
image(t(matrix(as.numeric(truth), n1, n2)[n1:1, ]), col = bpad_pal[c(7,3)],
      axes = FALSE, asp = 1, main = "anatomy (reference)"); box(col = "grey70")
hist(fiedler, breaks = 40, col = "grey85", border = "grey60",
     main = "Fiedler value distribution", xlab = "v2 component")
abline(v = 0, col = bpad_pal[2], lwd = 2)

$Spectral segmentation of a real kidney label plane. Top: the input image built from the expert mask, the Laplacian eigenvalue spectrum, and the Fiedler vector reshaped to the image grid. Bottom: the sign of the Fiedler vector gives the bipartition, compared with the anatomy. The smallest eigenvalue is zero to machine precision, as Eq. \@ref(eq:laplacian-quadratic) requires, and the second eigenvalue is small but non-zero, indicating one weakly connected pair of regions.$

Figure 6: Spectral segmentation of a real kidney label plane. Top: the input image built from the expert mask, the Laplacian eigenvalue spectrum, and the Fiedler vector reshaped to the image grid. Bottom: the sign of the Fiedler vector gives the bipartition, compared with the anatomy. The smallest eigenvalue is zero to machine precision, as Eq. (10) requires, and the second eigenvalue is small but non-zero, indicating one weakly connected pair of regions.

par(op)

From spectral methods to deep segmentation. Spectral segmentation is an unsupervised, per-image optimization: it uses no training data and solves a fresh eigenproblem for every image, at $O(N^3)$ cost for a dense Laplacian. A U-Net (Section 8.10.7) is a supervised, amortized alternative: it pays a large one-time training cost and then segments a new image in a single forward pass. The trade is training data and generalization risk in exchange for speed and learned semantics. Spectral methods remain useful exactly where labelled data do not exist, and their eigenvector reasoning reappears inside graph neural networks.

Section 8.9 summary.

Unsupervised structure cannot establish a phenotype; naming a cluster with the outcome is circular. Demand stability, upstream definition, and external validation.
On the real cohort, three-cluster $k$-means recovered operative extent and tumour size, not malignancy, whose proportion was nearly constant across clusters.
Silhouette, the gap statistic, and bootstrap stability agreed that only $k=2$ was supported; stability collapsed beyond it.
Intensity-only segmentation is destroyed by a smooth bias field that the eye ignores entirely.
Dice and boundary distances fail in opposite directions: a translation preserves area exactly while displacing every boundary point, and a dilation does the reverse. Report both.
The graph Laplacian’s Fiedler vector recovers the anatomy of a real mask, with $\lambda_1 = 0$ to machine precision and a normalized cut orders of magnitude better than chance.

Checkpoint 8.9. A paper reports three imaging phenotypes of renal masses from $k$-means on 40 radiomic features, and notes that phenotype 3 has the highest recurrence rate. List the three checks you would demand before believing that phenotype 3 is biological, and state what result on the real cohort in Section 8.9.2 should make you sceptical by default.

8.10 Deep Learning and Image-Native AI

Classical radiomics fixes the feature map $g$ and learns only the final mapping $f$. Deep learning learns both jointly,

\[\begin{equation} \widehat y = f_\theta\{g_\phi(I)\}, \tag{12} \end{equation}\]

which is Eq. (8.4) of Part 1 with $\psi$ replaced by learnable parameters $\phi$. The flexibility buys spatial context and removes the need to guess which features matter. It also raises the data, compute, validation, and interpretability requirements sharply.

The two cultures, stated fairly. In classical radiomics (Part 1, Section 8.5) a human specifies $g$ as shape, texture, and intensity statistics, and the model learns only $f$. The features are interpretable, small samples are viable, and reproducibility across sites is achievable if the provenance is documented. The limitation is that the features are fixed and may miss patterns nobody thought to encode. In deep learning the model learns $g$ and $f$ end to end. It can discover subtle spatial structure, at the cost of needing thousands of labelled examples, a strong susceptibility to shortcut learning (Section 8.12.6), and opacity. The question is never “deep learning or radiomics”; it is which representation-learning strategy matches the available sample size, label quality, and clinical question.

8.10.1 Neurons, layers, and backpropagation

A single artificial neuron computes $a = \sigma(w^{\mathsf T}x + b)$. Stacking gives, for layers $\ell = 1,\ldots,L$,

\[\begin{equation} h^{(\ell)} = \sigma_\ell\!\left(W^{(\ell)}h^{(\ell-1)} + b^{(\ell)}\right), \qquad h^{(0)} = x . \tag{13} \end{equation}\]

Training minimizes an empirical risk by stochastic gradient descent over mini-batches $B$,

\[\begin{equation} \theta_{t+1} = \theta_t - \frac{\eta_t}{\lvert B\rvert}\sum_{i\in B}\nabla_\theta \mathcal{L}_i(\theta_t), \tag{14} \end{equation}\]

with the gradients supplied by backpropagation, which is the chain rule of Part 1, Eq. (8.19), applied to the composition in Eq. (13). Automatic differentiation implements this reliably. It does not validate labels, prevent leakage, or guarantee that the learned representation is clinically sensible.

op <- par(mar = c(1, 1, 2.6, 1))
plot(NA, xlim = c(0, 10), ylim = c(0, 10), axes = FALSE, xlab = "", ylab = "",
     main = "One-hidden-layer network")
in_lab <- c("log_size", "age", "bmi", "radical", "male")
ny_in <- seq(8.6, 1.4, length.out = length(in_lab))
ny_h  <- seq(7.6, 2.4, length.out = 4)
x_in <- 1.6; x_h <- 5.0; x_o <- 8.4
for (i in seq_along(ny_in)) for (j in seq_along(ny_h))
  segments(x_in + 0.42, ny_in[i], x_h - 0.42, ny_h[j], col = "grey82", lwd = 0.7)
for (j in seq_along(ny_h))
  segments(x_h + 0.42, ny_h[j], x_o - 0.42, 5, col = "grey65", lwd = 1.0)
symbols(rep(x_in, length(ny_in)), ny_in, circles = rep(0.42, length(ny_in)),
        inches = FALSE, add = TRUE, bg = "white", fg = bpad_pal[1])
text(x_in - 0.62, ny_in, in_lab, adj = 1, cex = 0.78)
symbols(rep(x_h, length(ny_h)), ny_h, circles = rep(0.42, length(ny_h)),
        inches = FALSE, add = TRUE, bg = "white", fg = bpad_pal[3])
text(x_h, ny_h, "tanh", cex = 0.62, col = bpad_pal[3])
symbols(x_o, 5, circles = 0.46, inches = FALSE, add = TRUE, bg = "white", fg = bpad_pal[2])
text(x_o, 5, expression(sigma), cex = 0.9, col = bpad_pal[2])
text(x_o + 0.72, 5, "P(malignant)", adj = 0, cex = 0.8)
text(c(x_in, x_h, x_o), 9.6,
     c("input\n(p = 5)", "hidden\n(H = 4)", "output\n(1)"), cex = 0.8, font = 2)
text(3.3, 0.5, expression(W^{(1)}*": "*p %*% H), cex = 0.8, col = bpad_pal[1])
text(6.7, 0.5, expression(W^{(2)}*": "*H %*% 1), cex = 0.8, col = bpad_pal[2])

Anatomy of a one-hidden-layer network drawn with base graphics, so that it renders identically in HTML and in Word. Inputs are the standardized clinical predictors; the hidden layer applies a nonlinearity to learned linear combinations; the output layer applies a logistic link to produce a probability. The parameter count is dominated by the input-to-hidden weight matrix.

Figure 7: Anatomy of a one-hidden-layer network drawn with base graphics, so that it renders identically in HTML and in Word. Inputs are the standardized clinical predictors; the hidden layer applies a nonlinearity to learned linear combinations; the output layer applies a logistic link to produce a probability. The parameter count is dominated by the input-to-hidden weight matrix.

par(op)
cat(sprintf("parameters: W1 = %d x %d = %d, b1 = %d, W2 = %d, b2 = 1  ->  total %d\n",
            5, 4, 20, 4, 4, 20 + 4 + 4 + 1))

## parameters: W1 = 5 x 4 = 20, b1 = 4, W2 = 4, b2 = 1  ->  total 29

8.10.2 An executable network from first principles, with a gradient check

Nothing clarifies backpropagation like writing it out. For a one-hidden-layer network with $\tanh$ activation, logistic output, and binary cross-entropy loss, the gradients are

\[\begin{equation} \frac{\partial \mathcal{L}}{\partial Z^{(2)}} = \frac{1}{n}(P - y), \qquad \frac{\partial \mathcal{L}}{\partial W^{(2)}} = A^{(1)\mathsf T}\frac{\partial \mathcal{L}}{\partial Z^{(2)}}, \tag{15} \end{equation}\]

\[\begin{equation} \frac{\partial \mathcal{L}}{\partial Z^{(1)}} = \left(\frac{\partial \mathcal{L}}{\partial Z^{(2)}} W^{(2)\mathsf T}\right)\odot\left(1 - A^{(1)}\odot A^{(1)}\right), \qquad \frac{\partial \mathcal{L}}{\partial W^{(1)}} = X^{\mathsf T}\frac{\partial \mathcal{L}}{\partial Z^{(1)}}, \tag{16} \end{equation}\]

where $\odot$ is the elementwise product and $1 - A\odot A$ is the derivative of $\tanh$. The factor $(P-y)$ in Eq. (15) is not an accident: for any generalized linear model with a canonical link, the gradient of the negative log-likelihood is the design matrix times the residual, exactly as in Part 1, Eq. (8.13).

A numerical gradient check is the single most valuable habit when implementing a network: compare each analytic partial derivative with a central finite difference. If they disagree beyond about $10^{-6}$ in relative terms, the analytic gradient is wrong.

train_mlp <- function(X, y, H = 4, eta = 0.05, epochs = 500, lambda = 1e-3, seed = 1) {
  set.seed(seed); n <- nrow(X); p <- ncol(X)
  W1 <- matrix(rnorm(p*H, 0, sqrt(2/p)), p, H); b1 <- rep(0, H)
  W2 <- matrix(rnorm(H,  0, sqrt(2/H)), H, 1);  b2 <- 0
  loss <- numeric(epochs)
  for (e in seq_len(epochs)) {
    ## forward
    Z1 <- sweep(X %*% W1, 2, b1, "+"); A1 <- tanh(Z1)
    Z2 <- as.vector(A1 %*% W2) + b2;   P  <- 1/(1 + exp(-Z2))
    loss[e] <- -mean(y*log(P + 1e-12) + (1 - y)*log(1 - P + 1e-12)) +
                lambda*(sum(W1^2) + sum(W2^2))/2
    ## backward -- Eqs. (8.60) and (8.61)
    dZ2 <- (P - y)/n
    gW2 <- t(A1) %*% dZ2 + lambda*W2; gb2 <- sum(dZ2)
    dZ1 <- (dZ2 %*% t(W2)) * (1 - A1^2)
    gW1 <- t(X)  %*% dZ1 + lambda*W1; gb1 <- colSums(dZ1)
    W1 <- W1 - eta*gW1; b1 <- b1 - eta*gb1
    W2 <- W2 - eta*gW2; b2 <- b2 - eta*gb2
  }
  list(W1 = W1, b1 = b1, W2 = W2, b2 = b2, loss = loss)
}
predict_mlp <- function(m, X)
  1/(1 + exp(-(as.vector(tanh(sweep(X %*% m$W1, 2, m$b1, "+")) %*% m$W2) + m$b2)))

Xn <- scale(as.matrix(D[, mvars])); yn <- y
cat(sprintf("design: %d patients x %d standardized predictors, prevalence %.3f\n",
            nrow(Xn), ncol(Xn), mean(yn)))

## design: 210 patients x 5 standardized predictors, prevalence 0.914

## ---- numerical gradient check -------------------------------------------
set.seed(2)
p <- ncol(Xn); H <- 3
W1 <- matrix(rnorm(p*H, 0, 0.3), p, H); W2 <- matrix(rnorm(H, 0, 0.3), H, 1)
b1 <- rep(0, H); b2 <- 0
loss_at <- function(W1, W2, b1, b2) {
  A1 <- tanh(sweep(Xn %*% W1, 2, b1, "+"))
  P  <- 1/(1 + exp(-(as.vector(A1 %*% W2) + b2)))
  -mean(yn*log(P + 1e-12) + (1 - yn)*log(1 - P + 1e-12))
}
A1 <- tanh(sweep(Xn %*% W1, 2, b1, "+"))
P  <- 1/(1 + exp(-(as.vector(A1 %*% W2) + b2)))
dZ2 <- (P - yn)/nrow(Xn)
gW2_analytic <- t(A1) %*% dZ2
dZ1 <- (dZ2 %*% t(W2)) * (1 - A1^2)
gW1_analytic <- t(Xn) %*% dZ1

eps <- 1e-6
gcheck <- do.call(rbind, lapply(seq_len(H), function(j) {
  Wp <- W2; Wp[j] <- Wp[j] + eps; Wm <- W2; Wm[j] <- Wm[j] - eps
  num <- (loss_at(W1, Wp, b1, b2) - loss_at(W1, Wm, b1, b2))/(2*eps)
  data.frame(parameter = sprintf("W2[%d]", j),
             analytic = gW2_analytic[j], numerical = num,
             abs_diff = abs(gW2_analytic[j] - num))
}))
gcheck <- rbind(gcheck, do.call(rbind, lapply(seq_len(p), function(i) {
  Wp <- W1; Wp[i,1] <- Wp[i,1] + eps; Wm <- W1; Wm[i,1] <- Wm[i,1] - eps
  num <- (loss_at(Wp, W2, b1, b2) - loss_at(Wm, W2, b1, b2))/(2*eps)
  data.frame(parameter = sprintf("W1[%d,1]", i),
             analytic = gW1_analytic[i,1], numerical = num,
             abs_diff = abs(gW1_analytic[i,1] - num))
})))
print(knitr::kable(data.frame(parameter = gcheck$parameter,
                              analytic  = signif(gcheck$analytic, 8),
                              numerical = signif(gcheck$numerical, 8),
                              abs_diff  = format(gcheck$abs_diff, digits = 3)),
  row.names = FALSE,
  caption = "Numerical gradient check. Hand-derived analytic gradients agree with central finite differences to about 1e-11, which verifies Eqs. (8.60) and (8.61) as implemented."))

## 
## 
## Table: (\#tab:mlp-from-scratch)Numerical gradient check. Hand-derived analytic gradients agree with central finite differences to about 1e-11, which verifies Eqs. (8.60) and (8.61) as implemented.
## 
## |parameter |   analytic|  numerical|abs_diff |
## |:---------|----------:|----------:|:--------|
## |W2[1]     | -0.0692511| -0.0692511|2.13e-11 |
## |W2[2]     |  0.0777703|  0.0777703|1.98e-12 |
## |W2[3]     | -0.0258796| -0.0258796|1.75e-12 |
## |W1[1,1]   | -0.0362276| -0.0362276|3.13e-11 |
## |W1[2,1]   | -0.0245355| -0.0245355|1.90e-11 |
## |W1[3,1]   |  0.0114528|  0.0114528|5.36e-12 |
## |W1[4,1]   | -0.0733838| -0.0733838|1.58e-11 |
## |W1[5,1]   |  0.0234061|  0.0234061|2.98e-11 |

cat(sprintf("maximum absolute discrepancy: %.2e\n", max(gcheck$abs_diff)))

## maximum absolute discrepancy: 3.13e-11

Hs <- c(1, 2, 4, 8, 16, 32)
sweepH <- do.call(rbind, lapply(Hs, function(H) {
  oofm <- rep(NA_real_, nrow(Xn)); tr_auc <- numeric(5)
  for (k in seq_len(5)) {
    tr <- folds != k
    m  <- train_mlp(Xn[tr, ], yn[tr], H = H, epochs = 500, seed = k)
    oofm[!tr] <- predict_mlp(m, Xn[!tr, , drop = FALSE])
    tr_auc[k] <- auc_rank(predict_mlp(m, Xn[tr, , drop = FALSE]), yn[tr])
  }
  data.frame(H = H, params = ncol(Xn)*H + H + H + 1,
             train_AUC = mean(tr_auc), oof_AUC = auc_rank(oofm, yn))
}))
logit_oof <- rep(NA_real_, nrow(Xn))
for (k in seq_len(5)) {
  tr <- folds != k
  f  <- suppressWarnings(glm(yn[tr] ~ ., data = data.frame(Xn[tr, ]), family = binomial()))
  logit_oof[!tr] <- predict(f, data.frame(Xn[!tr, , drop = FALSE]), type = "response")
}
logit_auc <- auc_rank(logit_oof, yn)
print(knitr::kable(data.frame(hidden_units = sweepH$H, parameters = sweepH$params,
                              train_AUC = round(sweepH$train_AUC, 3),
                              out_of_fold_AUC = round(sweepH$oof_AUC, 3),
                              optimism = round(sweepH$train_AUC - sweepH$oof_AUC, 3)),
  row.names = FALSE,
  caption = sprintf(paste("Capacity sweep. Logistic regression on the same predictors gives",
                          "out-of-fold AUC %.3f. Every network sits within about 0.1 of that",
                          "benchmark, while the optimism column grows by nearly two orders of",
                          "magnitude."), logit_auc)))

## 
## 
## Table: (\#tab:mlp-capacity)Capacity sweep. Logistic regression on the same predictors gives out-of-fold AUC 0.698. Every network sits within about 0.1 of that benchmark, while the optimism column grows by nearly two orders of magnitude.
## 
## | hidden_units| parameters| train_AUC| out_of_fold_AUC| optimism|
## |------------:|----------:|---------:|---------------:|--------:|
## |            1|          8|     0.632|           0.630|    0.002|
## |            2|         15|     0.715|           0.600|    0.115|
## |            4|         29|     0.788|           0.692|    0.095|
## |            8|         57|     0.814|           0.682|    0.132|
## |           16|        113|     0.840|           0.704|    0.136|
## |           32|        225|     0.854|           0.718|    0.137|

cat(sprintf("training AUC rises from %.3f to %.3f (a gain of %.3f)\n",
            min(sweepH$train_AUC), max(sweepH$train_AUC),
            max(sweepH$train_AUC) - min(sweepH$train_AUC)))

## training AUC rises from 0.632 to 0.854 (a gain of 0.223)

cat(sprintf("out-of-fold AUC spans only %.3f to %.3f, straddling the logistic benchmark %.3f\n",
            min(sweepH$oof_AUC), max(sweepH$oof_AUC), logit_auc))

## out-of-fold AUC spans only 0.600 to 0.718, straddling the logistic benchmark 0.698

cat(sprintf("optimism grows from %.3f to %.3f as capacity increases\n",
            min(sweepH$train_AUC - sweepH$oof_AUC), max(sweepH$train_AUC - sweepH$oof_AUC)))

## optimism grows from 0.002 to 0.137 as capacity increases

cap <- rbind(data.frame(H = sweepH$H, auc = sweepH$train_AUC, s = "training"),
             data.frame(H = sweepH$H, auc = sweepH$oof_AUC,   s = "out-of-fold"))
c1 <- ggplot(cap, aes(H, auc, colour = s)) +
  geom_hline(yintercept = logit_auc, linetype = "dashed", colour = "grey35") +
  geom_hline(yintercept = 0.5, linetype = "dotted", colour = "grey60") +
  geom_line(linewidth = 0.9) + geom_point(size = 2) +
  scale_x_log10(breaks = Hs) + coord_cartesian(ylim = c(0.35, 0.95)) +
  scale_colour_manual(values = bpad_pal[1:2]) +
  labs(title = "training rises, generalization does not",
       subtitle = "dashed line: logistic regression out-of-fold AUC",
       x = "hidden units (log scale)", y = "AUC")
mfull <- train_mlp(Xn, yn, H = 8, epochs = 800)
c2 <- ggplot(data.frame(e = seq_along(mfull$loss), l = mfull$loss), aes(e, l)) +
  geom_line(colour = bpad_pal[3], linewidth = 0.9) +
  labs(title = "training loss (H = 8)", x = "epoch", y = "penalized cross-entropy")
bpad_grid(c1, c2, ncol = 2)

Capacity does not buy generalization on a small cohort. A one-hidden-layer network is trained on the real standardized design at hidden widths from 1 to 32 units. Training discrimination rises steadily with capacity while out-of-fold discrimination wanders within a narrow band around the logistic-regression benchmark (dashed line), so the widening gap between the curves is pure optimism: the variance term of Part 1, Eq. (8.20), made visible. Right: the training loss curve, which falls smoothly and tells you nothing whatever about generalization.

Figure 8: Capacity does not buy generalization on a small cohort. A one-hidden-layer network is trained on the real standardized design at hidden widths from 1 to 32 units. Training discrimination rises steadily with capacity while out-of-fold discrimination wanders within a narrow band around the logistic-regression benchmark (dashed line), so the widening gap between the curves is pure optimism: the variance term of Part 1, Eq. (8.20), made visible. Right: the training loss curve, which falls smoothly and tells you nothing whatever about generalization.

A smooth loss curve is not evidence of anything. The right-hand panel above is the plot most often shown to demonstrate that a network “converged”. It is compatible with a model that has learned the signal and with a model that has memorized the noise, and on this cohort it is the latter. Convergence is a statement about the optimizer; generalization is a statement about the world. Only the left panel speaks to the second.

8.10.3 Convolution: weight sharing and translation equivariance

A convolution layer replaces the dense matrix of Eq. (13) with a small kernel applied at every position,

\[\begin{equation} (I * K)[m,n] = \sum_{a}\sum_{b} K[a,b]\, I[m-a,\ n-b], \tag{17} \end{equation}\]

which is the same operator introduced in Part 1 and used for filtering in Chapter 7. Three consequences follow immediately, and they are the entire reason CNNs work on images.

Weight sharing. The same $k\times k$ kernel is reused at every spatial location, so the parameter count is independent of image size.

Translation equivariance. Shifting the input shifts the output identically, $\,(T_\delta I)*K = T_\delta(I*K)$. A lesion detector learned in one part of the image works everywhere.

Locality. Each output depends on a $k\times k$ neighbourhood, which is why the receptive field must be grown deliberately (Section 8.10.4).

conv2d <- function(I, K) {
  kh <- nrow(K); kw <- ncol(K); ph <- (kh - 1) %/% 2; pw <- (kw - 1) %/% 2
  P <- matrix(0, nrow(I) + 2*ph, ncol(I) + 2*pw)
  P[(ph+1):(ph+nrow(I)), (pw+1):(pw+ncol(I))] <- I
  O <- matrix(0, nrow(I), ncol(I))
  for (a in seq_len(kh)) for (b in seq_len(kw))
    O <- O + K[a, b]*P[a:(a + nrow(I) - 1), b:(b + ncol(I) - 1)]
  O
}
set.seed(5)
img_c <- (ax$plane > 0)*1.0 + 0.45*(ax$plane == 2)
img_c <- img_c + matrix(rnorm(length(img_c), 0, 0.04), nrow(img_c))

Kx   <- matrix(c(-1,0,1,-2,0,2,-1,0,1), 3, 3, byrow = TRUE)   # Sobel x
Ky   <- t(Kx)                                                  # Sobel y
Klap <- matrix(c(0,1,0, 1,-4,1, 0,1,0), 3, 3, byrow = TRUE)    # Laplacian
Kbox <- matrix(1, 3, 3)/9                                      # box blur
gx <- conv2d(img_c, Kx); gy <- conv2d(img_c, Ky)
gm <- sqrt(gx^2 + gy^2)

show_im <- function(M, ttl) { image(t(M[nrow(M):1, ]), col = grey.colors(64),
  axes = FALSE, asp = 1, main = ttl); box(col = "grey70") }
op <- par(mfrow = c(2, 3), mar = c(0.7, 0.7, 2.4, 0.7))
show_im(img_c, "input (from real mask)")
show_im(gx, "Sobel x (vertical edges)")
show_im(gy, "Sobel y (horizontal edges)")
show_im(gm, "gradient magnitude")
show_im(conv2d(img_c, Klap), "Laplacian")
show_im(conv2d(img_c, Kbox), "3x3 box blur")

$Convolution on a real kidney label plane rendered as an image. Top: the input, and the horizontal and vertical Sobel responses, which are the edge detectors a first convolutional layer typically learns without being told to. Bottom: gradient magnitude, a Laplacian response, and a box blur. All are the same linear operator of Eq. \@ref(eq:convolution) with different kernels; a CNN differs only in that the kernel entries are learned rather than specified.$

Figure 9: Convolution on a real kidney label plane rendered as an image. Top: the input, and the horizontal and vertical Sobel responses, which are the edge detectors a first convolutional layer typically learns without being told to. Bottom: gradient magnitude, a Laplacian response, and a box blur. All are the same linear operator of Eq. (17) with different kernels; a CNN differs only in that the kernel entries are learned rather than specified.

par(op)

## convolution is a linear operator: verify additivity and homogeneity
set.seed(1); A <- matrix(rnorm(length(img_c)), nrow(img_c)); a <- 2.3; b <- -0.7
lhs <- conv2d(a*img_c + b*A, Kx); rhs <- a*conv2d(img_c, Kx) + b*conv2d(A, Kx)
cat(sprintf("linearity check |conv(aI + bA) - a*conv(I) - b*conv(A)| max = %.2e\n",
            max(abs(lhs - rhs))))

## linearity check |conv(aI + bA) - a*conv(I) - b*conv(A)| max = 5.33e-15

h <- nrow(img_c); w <- ncol(img_c); cout <- 16; k <- 3
print(knitr::kable(data.frame(
  layer_type = c("dense (fully connected)", "convolution 3x3"),
  input = sprintf("%d x %d", h, w),
  output_channels = cout,
  parameters = c(h*w*cout + cout, k*k*1*cout + cout),
  translation_equivariant = c("no", "yes")),
  row.names = FALSE,
  caption = sprintf("Parameter counts for one layer on this image. The convolution uses %.0f times fewer parameters and, unlike the dense layer, is translation equivariant.",
                    (h*w*cout + cout)/(k*k*cout + cout))))

## 
## 
## Table: (\#tab:convolution-demo)Parameter counts for one layer on this image. The convolution uses 670 times fewer parameters and, unlike the dense layer, is translation equivariant.
## 
## |layer_type              |input   | output_channels| parameters|translation_equivariant |
## |:-----------------------|:-------|---------------:|----------:|:-----------------------|
## |dense (fully connected) |87 x 77 |              16|     107200|no                      |
## |convolution 3x3         |87 x 77 |              16|        160|yes                     |

8.10.4 Receptive field, stride, and the physical size of the target

The receptive field is the region of the input that influences one output unit. It grows through a stack of layers according to a simple recursion: with kernel size $k_\ell$, stride $s_\ell$, and cumulative jump $j_\ell = \prod_{m\le\ell}s_m$,

\[\begin{equation} r_\ell = r_{\ell-1} + (k_\ell - 1)\,j_{\ell-1}, \qquad r_0 = 1,\ j_0 = 1 . \tag{18} \end{equation}\]

Equation (18) is where deep learning meets biomedical physics, because the receptive field must be converted to millimetres using the voxel spacing before it means anything. A network whose deepest receptive field is 30 mm cannot see a 45 mm tumour together with its surrounding parenchyma, no matter how many parameters it has.

arch <- data.frame(
  layer  = c("conv1","conv2","pool1","conv3","conv4","pool2","conv5","conv6","pool3","conv7"),
  k      = c(3, 3, 2, 3, 3, 2, 3, 3, 2, 3),
  stride = c(1, 1, 2, 1, 1, 2, 1, 1, 2, 1))
r <- 1; j <- 1; rows <- NULL
for (i in seq_len(nrow(arch))) {
  r <- r + (arch$k[i] - 1)*j
  j <- j*arch$stride[i]
  rows <- rbind(rows, data.frame(layer = arch$layer[i], k = arch$k[i],
                                 stride = arch$stride[i], jump = j,
                                 rf_px = r, rf_mm = r*ax$spacing[1]))
}
print(knitr::kable(data.frame(layer = rows$layer, kernel = rows$k, stride = rows$stride,
                              cumulative_jump = rows$jump,
                              receptive_field_px = rows$rf_px,
                              receptive_field_mm = round(rows$rf_mm, 2)),
  row.names = FALSE,
  caption = sprintf("Receptive field growth, converted to millimetres at the in-plane spacing of %.2f mm.", ax$spacing[1])))

## 
## 
## Table: (\#tab:receptive-field)Receptive field growth, converted to millimetres at the in-plane spacing of 0.92 mm.
## 
## |layer | kernel| stride| cumulative_jump| receptive_field_px| receptive_field_mm|
## |:-----|------:|------:|---------------:|------------------:|------------------:|
## |conv1 |      3|      1|               1|                  3|               2.76|
## |conv2 |      3|      1|               1|                  5|               4.60|
## |pool1 |      2|      2|               2|                  6|               5.52|
## |conv3 |      3|      1|               2|                 10|               9.20|
## |conv4 |      3|      1|               2|                 14|              12.88|
## |pool2 |      2|      2|               4|                 16|              14.72|
## |conv5 |      3|      1|               4|                 24|              22.08|
## |conv6 |      3|      1|               4|                 32|              29.44|
## |pool3 |      2|      2|               8|                 36|              33.12|
## |conv7 |      3|      1|               8|                 52|              47.84|

sz_q <- quantile(kits$radiographic_size_cm, c(0.5, 0.75, 0.95), na.rm = TRUE)*10
ggplot(rows, aes(seq_along(layer), rf_mm)) +
  geom_hline(yintercept = sz_q, linetype = "dashed", colour = bpad_pal[2]) +
  annotate("text", x = 1, y = sz_q, hjust = 0, vjust = -0.5, size = 2.7,
           colour = bpad_pal[2],
           label = sprintf("%s tumour diameter: %.0f mm",
                           c("median", "75th pct", "95th pct"), sz_q)) +
  geom_line(colour = bpad_pal[1], linewidth = 0.95) +
  geom_point(size = 2, colour = bpad_pal[1]) +
  scale_x_continuous(breaks = seq_along(rows$layer), labels = rows$layer) +
  labs(title = "Receptive field in physical units",
       x = NULL, y = "receptive field (mm)") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Receptive field growth through a small CNN, in pixels and in millimetres at the real in-plane spacing of this cohort. The horizontal dashed lines mark the median and upper-quartile tumour diameters in the cohort. Six convolutions and two pooling layers reach a receptive field comparable to a median tumour, but not to the largest ones: architectural depth is a physical design constraint, not only a capacity choice.

Figure 10: Receptive field growth through a small CNN, in pixels and in millimetres at the real in-plane spacing of this cohort. The horizontal dashed lines mark the median and upper-quartile tumour diameters in the cohort. Six convolutions and two pooling layers reach a receptive field comparable to a median tumour, but not to the largest ones: architectural depth is a physical design constraint, not only a capacity choice.

cat(sprintf("final receptive field: %d px = %.1f mm; %.0f%% of tumours exceed it\n",
            max(rows$rf_px), max(rows$rf_mm),
            100*mean(kits$radiographic_size_cm*10 > max(rows$rf_mm), na.rm = TRUE)))

## final receptive field: 52 px = 47.8 mm; 35% of tumours exceed it

8.10.5 CNN classification architecture and parameter budgets

A classification CNN is a stack of convolution, nonlinearity, and downsampling blocks that progressively trade spatial resolution for channel depth, ending in a global pooling step and a small dense head. The parameter budget is worth computing before training, because it determines what sample size is plausible.

cnn <- data.frame(
  layer = c("conv1 3x3","conv2 3x3","pool","conv3 3x3","conv4 3x3","pool",
            "conv5 3x3","conv6 3x3","pool","global avg pool","dense","output"),
  cin   = c(1,32,32,32,64,64,64,128,128,128,128,64),
  cout  = c(32,32,32,64,64,64,128,128,128,128,64,1),
  k     = c(3,3,NA,3,3,NA,3,3,NA,NA,NA,NA),
  H     = c(128,128,64,64,64,32,32,32,16,1,1,1))
cnn$params <- with(cnn, ifelse(!is.na(k), k*k*cin*cout + cout,
                        ifelse(layer %in% c("dense","output"), cin*cout + cout, 0)))
print(knitr::kable(cnn[, c("layer","cin","cout","H","params")], row.names = FALSE,
  caption = "Layer-by-layer parameter budget for a compact 2D classification CNN on 128 x 128 inputs."))

## 
## 
## Table: (\#tab:cnn-budget)Layer-by-layer parameter budget for a compact 2D classification CNN on 128 x 128 inputs.
## 
## |layer           | cin| cout|   H| params|
## |:---------------|---:|----:|---:|------:|
## |conv1 3x3       |   1|   32| 128|    320|
## |conv2 3x3       |  32|   32| 128|   9248|
## |pool            |  32|   32|  64|      0|
## |conv3 3x3       |  32|   64|  64|  18496|
## |conv4 3x3       |  64|   64|  64|  36928|
## |pool            |  64|   64|  32|      0|
## |conv5 3x3       |  64|  128|  32|  73856|
## |conv6 3x3       | 128|  128|  32| 147584|
## |pool            | 128|  128|  16|      0|
## |global avg pool | 128|  128|   1|      0|
## |dense           | 128|   64|   1|   8256|
## |output          |  64|    1|   1|     65|

cat(sprintf("total parameters: %s\n", format(sum(cnn$params), big.mark = ",")))

## total parameters: 294,753

cat(sprintf("labelled examples available in this cohort: %d patients (%d benign)\n",
            nrow(kits), sum(!kits$malignant)))

## labelled examples available in this cohort: 210 patients (18 benign)

cat(sprintf("parameters per labelled patient: %.0f\n", sum(cnn$params)/nrow(kits)))

## parameters per labelled patient: 1404

b1 <- ggplot(subset(cnn, params > 0),
             aes(reorder(paste0(seq_len(sum(cnn$params > 0)), ". ", layer), params), params)) +
  geom_col(fill = bpad_pal[1], alpha = 0.85) + coord_flip() + scale_y_log10() +
  labs(title = "parameters per layer (log)", x = NULL, y = "parameters")
b2 <- ggplot(cnn, aes(seq_along(layer))) +
  geom_line(aes(y = H, colour = "spatial size (px)"), linewidth = 0.9) +
  geom_line(aes(y = cout, colour = "channels"), linewidth = 0.9) +
  geom_point(aes(y = H), size = 1.5) + geom_point(aes(y = cout), size = 1.5) +
  scale_y_log10() + scale_colour_manual(values = bpad_pal[1:2]) +
  scale_x_continuous(breaks = seq_along(cnn$layer), labels = cnn$layer) +
  theme(axis.text.x = element_text(angle = 55, hjust = 1)) +
  labs(title = "resolution traded for depth", x = NULL, y = "value (log)")
bpad_grid(b1, b2, ncol = 2)

Parameter budget of a compact classification CNN. Left: parameters per layer; the convolutional stack is cheap and the cost concentrates in the deepest wide layers. Right: spatial resolution against channel depth through the network, the characteristic trade that lets a CNN build semantic abstraction while keeping computation bounded. The total parameter count is printed below, alongside the number of labelled examples this cohort actually provides.

Figure 11: Parameter budget of a compact classification CNN. Left: parameters per layer; the convolutional stack is cheap and the cost concentrates in the deepest wide layers. Right: spatial resolution against channel depth through the network, the characteristic trade that lets a CNN build semantic abstraction while keeping computation bounded. The total parameter count is printed below, alongside the number of labelled examples this cohort actually provides.

## ---------------------------------------------------------------------------
## Reference implementation of the architecture costed above. Not evaluated in
## this chapter: a deep-learning framework is a heavyweight, platform-specific
## dependency, and none of the conclusions of this chapter rest on running it.
## Requires: install.packages("keras3"); keras3::install_keras()
## ---------------------------------------------------------------------------
library(keras3)

build_classifier <- function(input_shape = c(128, 128, 1)) {
  keras_model_sequential(input_shape = input_shape) |>
    layer_conv_2d(32,  c(3,3), padding = "same", activation = "relu") |>
    layer_conv_2d(32,  c(3,3), padding = "same", activation = "relu") |>
    layer_max_pooling_2d(c(2,2)) |>
    layer_conv_2d(64,  c(3,3), padding = "same", activation = "relu") |>
    layer_conv_2d(64,  c(3,3), padding = "same", activation = "relu") |>
    layer_max_pooling_2d(c(2,2)) |>
    layer_conv_2d(128, c(3,3), padding = "same", activation = "relu") |>
    layer_conv_2d(128, c(3,3), padding = "same", activation = "relu") |>
    layer_global_average_pooling_2d() |>
    layer_dense(64, activation = "relu") |>
    layer_dropout(0.4) |>
    layer_dense(1, activation = "sigmoid")
}

model <- build_classifier()
model |> compile(optimizer = optimizer_adam(1e-4),
                 loss = "binary_crossentropy",
                 metrics = list(metric_auc(name = "auc")))
## CRITICAL: the generator must split by PATIENT, never by slice.
## See Part 1, Section 8.3.6 for what slice-level splitting costs.
history <- model |> fit(train_gen, validation_data = val_gen, epochs = 60,
                        callbacks = list(callback_early_stopping(monitor = "val_auc",
                                                                 patience = 10,
                                                                 mode = "max",
                                                                 restore_best_weights = TRUE)))

8.10.6 Encoder-decoder networks and U-Net segmentation

Segmentation requires a dense output: one label per voxel. An encoder-decoder achieves this by contracting to a low-resolution, semantically rich representation and then expanding back to full resolution. The difficulty is that the contraction destroys the precise spatial detail that a boundary needs. Skip connections solve this by concatenating each encoder feature map onto the corresponding decoder stage, so the decoder has access to both high-level semantics and high-resolution detail. That single idea is what made U-Net the default architecture for biomedical segmentation.

op <- par(mar = c(1, 1, 2.6, 1))
plot(NA, xlim = c(0, 10), ylim = c(0, 10), axes = FALSE, xlab = "", ylab = "",
     main = "U-Net encoder-decoder with skip connections")
levels_ <- data.frame(y = c(8.5, 6.8, 5.1, 3.4),
                      res = c("128", "64", "32", "16"),
                      ch  = c("32", "64", "128", "256"))
xe <- 2.2; xd <- 7.8
for (i in seq_len(nrow(levels_))) {
  rect(xe - 0.9, levels_$y[i] - 0.42, xe + 0.9, levels_$y[i] + 0.42,
       col = "#eef5fb", border = bpad_pal[1])
  text(xe, levels_$y[i], sprintf("%s^2 x %s", levels_$res[i], levels_$ch[i]), cex = 0.72)
  rect(xd - 0.9, levels_$y[i] - 0.42, xd + 0.9, levels_$y[i] + 0.42,
       col = "#fdf2e7", border = bpad_pal[2])
  text(xd, levels_$y[i], sprintf("%s^2 x %s", levels_$res[i], levels_$ch[i]), cex = 0.72)
  if (i < nrow(levels_)) {
    arrows(xe, levels_$y[i] - 0.45, xe, levels_$y[i+1] + 0.45, length = 0.08,
           col = bpad_pal[1], lwd = 1.6)
    arrows(xd, levels_$y[i+1] + 0.45, xd, levels_$y[i] - 0.45, length = 0.08,
           col = bpad_pal[2], lwd = 1.6)
  }
  if (i < nrow(levels_)) {
    arrows(xe + 0.95, levels_$y[i], xd - 0.95, levels_$y[i], length = 0.09,
           col = "grey45", lwd = 1.4, lty = 1)
    text((xe + xd)/2, levels_$y[i] + 0.28, "skip (concatenate)", cex = 0.62, col = "grey35")
  }
}
rect(4.4, 1.6, 5.6, 2.5, col = "#eef7f0", border = bpad_pal[3])
text(5.0, 2.05, "bottleneck\n8^2 x 512", cex = 0.68)
arrows(xe, 3.0, 4.4, 2.3, length = 0.08, col = bpad_pal[1], lwd = 1.6)
arrows(5.6, 2.3, xd, 3.0, length = 0.08, col = bpad_pal[2], lwd = 1.6)
text(xe, 9.5, "encoder\n(downsample)", cex = 0.8, font = 2, col = bpad_pal[1])
text(xd, 9.5, "decoder\n(upsample)", cex = 0.8, font = 2, col = bpad_pal[2])
text(xd + 1.35, 8.5, "-> per-voxel\n   softmax", cex = 0.7, adj = 0)

U-Net topology drawn with base graphics. The encoder (left) halves resolution and doubles channel depth at each stage; the decoder (right) reverses this. Grey horizontal arrows are skip connections, which carry high-resolution spatial detail directly across the bottleneck. Without them the decoder must reconstruct boundaries from a coarse representation, and boundary metrics (Section 8.9.5) degrade sharply even when overlap metrics look acceptable.

Figure 12: U-Net topology drawn with base graphics. The encoder (left) halves resolution and doubles channel depth at each stage; the decoder (right) reverses this. Grey horizontal arrows are skip connections, which carry high-resolution spatial detail directly across the bottleneck. Without them the decoder must reconstruct boundaries from a coarse representation, and boundary metrics (Section 8.9.5) degrade sharply even when overlap metrics look acceptable.

par(op)

Segmentation networks are trained with a loss that combines a per-voxel likelihood term with an overlap term, because cross-entropy alone is dominated by the overwhelmingly abundant background class:

\[\begin{equation} \mathcal{L} = \underbrace{-\frac{1}{N}\sum_{v}\sum_{c} y_{vc}\log \widehat p_{vc}}_{\text{cross-entropy}} \ +\ \lambda\underbrace{\left(1 - \frac{2\sum_v \widehat p_v y_v + \epsilon}{\sum_v \widehat p_v + \sum_v y_v + \epsilon}\right)}_{\text{soft Dice}} . \tag{19} \end{equation}\]

The soft Dice term is a differentiable relaxation of Eq. (4) in which the binary prediction is replaced by the predicted probability, and $\epsilon$ prevents division by zero on empty slices.

## The cropped plane understates the imbalance, because cropping removed
## background. Report the frequencies for the FULL 512 x 512 acquired slice,
## which is what a network would actually be shown.
n_kidney <- sum(ax$plane == 1); n_tumour <- sum(ax$plane == 2)
tot <- 512L*512L
freq <- data.frame(class = c("background", "kidney", "tumour"),
                   n = c(tot - n_kidney - n_tumour, n_kidney, n_tumour))
freq$pct <- 100*freq$n/tot
print(knitr::kable(data.frame(class = freq$class, pixels = freq$n,
                              percent = round(freq$pct, 3)), row.names = FALSE,
  caption = "Class frequencies in the full 512 x 512 acquired axial slice. The tumour occupies a fifth of one percent of the image."))

## 
## 
## Table: (\#tab:class-imbalance-seg)Class frequencies in the full 512 x 512 acquired axial slice. The tumour occupies a fifth of one percent of the image.
## 
## |class      | pixels| percent|
## |:----------|------:|-------:|
## |background | 259702|  99.068|
## |kidney     |   1905|   0.727|
## |tumour     |    537|   0.205|

## a trivial "all background" predictor
p_bg <- freq$n[1]/tot
## cross-entropy of a predictor that always outputs the marginal class frequencies
ce_trivial <- -sum(freq$n*log(freq$pct/100 + 1e-12))/tot
dice_trivial_tumour <- 0
cat(sprintf("trivial all-background predictor: cross-entropy %.4f, tumour Dice %.3f\n",
            ce_trivial, dice_trivial_tumour))

## trivial all-background predictor: cross-entropy 0.0577, tumour Dice 0.000

cat(sprintf("pixel accuracy of that predictor: %.4f\n", p_bg))

## pixel accuracy of that predictor: 0.9907

f1 <- ggplot(freq, aes(reorder(class, -pct), pct, fill = class)) +
  geom_col(alpha = 0.85, show.legend = FALSE) +
  geom_text(aes(label = sprintf("%.2f%%", pct)), vjust = -0.4, size = 3) +
  scale_fill_manual(values = bpad_pal[c(7,1,2)]) + scale_y_log10() +
  labs(title = "class frequency (log scale)", x = NULL, y = "percent of pixels")
metdf <- data.frame(metric = c("pixel accuracy", "cross-entropy", "tumour Dice"),
                    value = c(p_bg, ce_trivial, dice_trivial_tumour))
f2 <- ggplot(metdf, aes(reorder(metric, value), value)) +
  geom_col(fill = bpad_pal[2], alpha = 0.85) +
  geom_text(aes(label = round(value, 3)), hjust = -0.15, size = 3) +
  coord_flip(ylim = c(0, 1.15)) +
  labs(title = "the trivial predictor scores well on the wrong metric",
       x = NULL, y = "value")
bpad_grid(f1, f2, ncol = 2)

$Why segmentation losses need an overlap term. Left: class frequencies in the real axial label plane; the tumour occupies well under one percent of the image. Right: the cross-entropy achieved by a model that predicts background everywhere, compared with the Dice such a model achieves. The trivial predictor attains a low cross-entropy and a Dice of exactly zero, which is precisely the failure the composite loss of Eq. \@ref(eq:dice-loss) is designed to prevent.$

Figure 13: Why segmentation losses need an overlap term. Left: class frequencies in the real axial label plane; the tumour occupies well under one percent of the image. Right: the cross-entropy achieved by a model that predicts background everywhere, compared with the Dice such a model achieves. The trivial predictor attains a low cross-entropy and a Dice of exactly zero, which is precisely the failure the composite loss of Eq. (19) is designed to prevent.

## ---------------------------------------------------------------------------
## Reference U-Net with the composite loss of Eq. (8.66). Not evaluated here.
## ---------------------------------------------------------------------------
library(keras3)

conv_block <- function(x, filters) {
  x |> layer_conv_2d(filters, c(3,3), padding = "same") |>
       layer_batch_normalization() |> layer_activation("relu") |>
       layer_conv_2d(filters, c(3,3), padding = "same") |>
       layer_batch_normalization() |> layer_activation("relu")
}
build_unet <- function(input_shape = c(128,128,1), n_class = 3, base = 32) {
  inp <- layer_input(input_shape)
  e1 <- conv_block(inp, base);      p1 <- layer_max_pooling_2d(e1, c(2,2))
  e2 <- conv_block(p1, base*2);     p2 <- layer_max_pooling_2d(e2, c(2,2))
  e3 <- conv_block(p2, base*4);     p3 <- layer_max_pooling_2d(e3, c(2,2))
  bn <- conv_block(p3, base*8)
  d3 <- layer_conv_2d_transpose(bn, base*4, c(2,2), strides = c(2,2)) |>
          (\(z) layer_concatenate(list(z, e3)))() |> conv_block(base*4)
  d2 <- layer_conv_2d_transpose(d3, base*2, c(2,2), strides = c(2,2)) |>
          (\(z) layer_concatenate(list(z, e2)))() |> conv_block(base*2)
  d1 <- layer_conv_2d_transpose(d2, base,   c(2,2), strides = c(2,2)) |>
          (\(z) layer_concatenate(list(z, e1)))() |> conv_block(base)
  out <- layer_conv_2d(d1, n_class, c(1,1), activation = "softmax")
  keras_model(inp, out)
}
## Composite loss: cross-entropy plus soft Dice, Eq. (8.66)
dice_loss <- function(y_true, y_pred, eps = 1e-6) {
  num <- 2*op_sum(y_true*y_pred); den <- op_sum(y_true) + op_sum(y_pred)
  1 - (num + eps)/(den + eps)
}
combined_loss <- function(y_true, y_pred)
  loss_categorical_crossentropy(y_true, y_pred) + dice_loss(y_true, y_pred)

8.10.7 Augmentation as a model of plausible variability

Augmentation is not a trick for manufacturing data. It is an explicit statement of which transformations leave the label unchanged, and therefore a model of the acquisition process. Physically defensible augmentations for abdominal CT include small rotations and translations (patient positioning), mild elastic deformation (respiration and peristalsis), intensity shifts within the plausible Hounsfield range (contrast timing and reconstruction kernel), realistic noise (dose reduction), and slice-thickness resampling.

Physically indefensible augmentations are equally important to name. Left-right flipping changes organ laterality and can invert the meaning of a finding. Arbitrary intensity rescaling destroys the absolute Hounsfield calibration that makes CT quantitative (Chapter 5). Large elastic deformation creates anatomy that does not exist. Independent white noise, as Part 1, Section 8.2.10 argued, is not what a scanner does.

rot_bilinear <- function(M, deg) {
  th <- deg*pi/180; n1 <- nrow(M); n2 <- ncol(M)
  c1 <- (n1 + 1)/2; c2 <- (n2 + 1)/2
  O <- matrix(0, n1, n2)
  for (i in seq_len(n1)) for (j in seq_len(n2)) {
    x <-  cos(th)*(i - c1) + sin(th)*(j - c2) + c1
    yv <- -sin(th)*(i - c1) + cos(th)*(j - c2) + c2
    x0 <- floor(x); y0 <- floor(yv)
    if (x0 >= 1 && y0 >= 1 && x0 < n1 && y0 < n2) {
      dx <- x - x0; dy <- yv - y0
      O[i, j] <- (1-dx)*(1-dy)*M[x0,y0] + dx*(1-dy)*M[x0+1,y0] +
                 (1-dx)*dy*M[x0,y0+1]   + dx*dy*M[x0+1,y0+1]
    }
  }
  O
}
set.seed(9)
base_img <- img_c
corr_noise <- conv2d(matrix(rnorm(length(base_img), 0, 0.35), nrow(base_img)), Kbox)
augs <- list(
  "original"                = base_img,
  "rotate 8 deg"            = rot_bilinear(base_img, 8),
  "intensity shift +12%"    = base_img*1.12,
  "correlated noise"        = base_img + corr_noise,
  "left-right FLIP (invalid)" = base_img[, ncol(base_img):1])
op <- par(mfrow = c(2, 3), mar = c(0.7, 0.7, 2.4, 0.7))
for (nm in names(augs)) show_im(augs[[nm]], nm)
## effect on a simple size biomarker: area above a fixed threshold
## A biomarker defined by an ABSOLUTE intensity threshold, as a Hounsfield-unit
## criterion would be. This is exactly the kind of measurement that an intensity
## augmentation silently corrupts.
thr <- 1.10
areas <- vapply(augs, function(M) sum(M > thr)*prod(ax$spacing)/100, numeric(1))
par(mar = c(7.5, 4.2, 2.4, 1))
bp <- barplot(100*(areas/areas[1] - 1), col = c("grey70", bpad_pal[c(1,3,5,2)]),
              las = 2, ylab = "change in thresholded area (%)",
              main = "effect on a size biomarker", cex.names = 0.7)
abline(h = 0, col = "grey40")

Augmentation as a model of acquisition variability, applied to the real label plane. Top row: plausible augmentations, a small rotation, a physically realistic intensity shift, and correlated noise. Bottom row: an implausible one, a left-right flip that reverses laterality, together with the effect of each augmentation on a simple size biomarker. A transformation belongs in the augmentation set only if it leaves the label invariant, and laterality flipping does not.

Figure 14: Augmentation as a model of acquisition variability, applied to the real label plane. Top row: plausible augmentations, a small rotation, a physically realistic intensity shift, and correlated noise. Bottom row: an implausible one, a left-right flip that reverses laterality, together with the effect of each augmentation on a simple size biomarker. A transformation belongs in the augmentation set only if it leaves the label invariant, and laterality flipping does not.

par(op)
print(knitr::kable(data.frame(augmentation = names(augs),
                              area_cm2 = round(areas, 3),
                              change_pct = round(100*(areas/areas[1] - 1), 1),
                              preserves_label = c("-", "yes", "yes", "yes", "NO"),
                              why = c("-", "small rigid motion is plausible",
                                      "plausible, but corrupts an absolute-threshold biomarker",
                                      "matches correlated scanner noise",
                                      "reverses organ laterality")),
  row.names = FALSE,
  caption = paste("Each augmentation's effect on a biomarker defined by an absolute",
                  "intensity threshold. Note the last two rows: the flip leaves the",
                  "measured area exactly unchanged and is still invalid, because",
                  "invariance of a summary statistic is not invariance of the label.")))

## 
## 
## Table: (\#tab:augmentation-demo)Each augmentation's effect on a biomarker defined by an absolute intensity threshold. Note the last two rows: the flip leaves the measured area exactly unchanged and is still invalid, because invariance of a summary statistic is not invariance of the label.
## 
## |augmentation              | area_cm2| change_pct|preserves_label |why                                                     |
## |:-------------------------|--------:|----------:|:---------------|:-------------------------------------------------------|
## |original                  |    4.629|        0.0|-               |-                                                       |
## |rotate 8 deg              |    4.637|        0.2|yes             |small rigid motion is plausible                         |
## |intensity shift +12%      |   15.300|      230.5|yes             |plausible, but corrupts an absolute-threshold biomarker |
## |correlated noise          |    7.633|       64.9|yes             |matches correlated scanner noise                        |
## |left-right FLIP (invalid) |    4.629|        0.0|NO              |reverses organ laterality                               |

8.10.8 Transfer learning, self-supervision, and foundation models

Medical imaging has abundant unlabelled images and scarce expert annotations, which is exactly the regime these methods address.

Transfer learning initializes from weights trained on a large source dataset and fine-tunes on the target task. Pretraining on natural images transfers surprisingly well for low-level features (edges, textures) and poorly for the semantics, so the usual practice is to freeze early layers and retrain the later ones. The pitfall is a domain gap in input statistics: natural images are three-channel, 8-bit, and gamma-encoded, whereas CT is single-channel, 12-bit, and linear in attenuation.

Self-supervised learning constructs the labels from the data itself: predicting masked patches, ordering slices, predicting rotation, or contrasting two augmented views of the same volume against views of different volumes. The augmentations chosen define the invariances the representation will learn, so Section 8.10.7 is not a preliminary to this method, it is the specification of it.

Foundation models are large models pretrained self-supervised on very large image corpora and then adapted with small labelled sets. They are promising for exactly the small-cohort regime of this chapter, and they import their pretraining corpus’s biases wholesale, which makes the fairness and shift analyses of Section 8.12 more important rather than less.

8.10.9 Attention and vision transformers

Attention computes a weighted combination of value vectors, with weights determined by the compatibility of queries and keys:

\[\begin{equation} \mathrm{Attention}(Q,K,V) = \mathrm{softmax}\!\left(\frac{QK^{\mathsf T}}{\sqrt{d_k}}\right)V . \tag{20} \end{equation}\]

The $\sqrt{d_k}$ scaling keeps the logits in a range where the softmax gradient does not vanish. A vision transformer splits an image into patches, embeds each patch, and applies self-attention across all patches, giving a global receptive field in a single layer, in contrast to the incremental growth of Eq. (18). The cost is quadratic complexity in the number of patches and a much weaker inductive bias, so transformers typically need more data than CNNs unless they are pretrained. Hybrid designs that use convolutions early and attention late are common in medical imaging for exactly this reason.

Attention maps are not explanations. It is routine to display an attention map as evidence that a model “looked at the tumour”. Attention weights are intermediate quantities in a computation, not attributions of the output to the input; models with visibly different attention maps can compute identical functions, and a plausible-looking map can coexist with shortcut learning. Section 8.12.5 treats what interpretability claims are actually supportable.

8.10.10 Weak, noisy, and multiple-instance labels

Expert voxel-level annotation is expensive, so many real datasets carry only a bag-level label: this scan contains a malignancy, without saying where. Multiple-instance learning formalizes this. A bag $X_i = \{x_{i1},\ldots,x_{im}\}$ has label

\[\begin{equation} Y_i = \max_j y_{ij}, \tag{21} \end{equation}\]

with the instance labels $y_{ij}$ unobserved. Training pools instance-level scores with a differentiable surrogate for the maximum, most often an attention-weighted mean, which additionally yields an instance-importance map for free.

Label noise is the other half of this problem, and Part 1, Section 8.3.2 gave the governing fact: measured performance is bounded by label quality. With symmetric label noise at rate $\rho$, the achievable accuracy is capped near $1-\rho$, and reported accuracy above that ceiling is evidence of leakage rather than skill.

8.10.11 Multimodal fusion

Imaging is rarely the only information available. Three fusion strategies are standard, and they differ in where the modalities meet.

Strategy	Where fusion occurs	Strengths	Weaknesses
Early	concatenate raw or low-level features	can model fine cross-modal interactions	requires alignment; dominated by the higher-dimensional modality
Intermediate	concatenate learned embeddings	each modality gets an appropriate encoder	needs both modalities present at inference
Late	combine model outputs	robust to missing modalities; simple to audit	cannot model interactions

The decisive question is not architectural but evidential: does the imaging add anything beyond the clinical variables? That is a nested-model comparison, and it must be answered with out-of-fold predictions.

mk_oof <- function(fml) {
  o <- rep(NA_real_, nrow(D))
  for (k in seq_len(5)) {
    tr <- folds != k
    f <- suppressWarnings(glm(fml, D[tr, ], family = binomial()))
    o[!tr] <- predict(f, D[!tr, ], type = "response")
  }
  o
}
models <- list(
  "size only"                = mk_oof(malignant ~ log_size),
  "clinical only"            = mk_oof(malignant ~ age + bmi + male),
  "size + clinical"          = mk_oof(malignant ~ log_size + age + bmi + male),
  "size + clinical + surgery"= mk_oof(malignant ~ log_size + age + bmi + male + radical))
boot_ci <- function(s, l, B = 1500) {
  r <- replicate(B, { i <- sample(length(s), replace = TRUE); auc_rank(s[i], l[i]) })
  r <- r[is.finite(r)]; stats::quantile(r, c(0.025, 0.975))
}
fus <- do.call(rbind, lapply(names(models), function(nm) {
  ci <- boot_ci(models[[nm]], y)
  data.frame(model = nm, AUC = auc_rank(models[[nm]], y), lo = ci[1], hi = ci[2])
}))
print(knitr::kable(data.frame(model = fus$model, AUC = round(fus$AUC, 3),
                              CI = sprintf("[%.3f, %.3f]", fus$lo, fus$hi),
                              width = round(fus$hi - fus$lo, 3)),
  row.names = FALSE,
  caption = "Incremental value assessed with uncertainty. Every interval contains every other point estimate."))

## 
## 
## Table: (\#tab:fusion-value)Incremental value assessed with uncertainty. Every interval contains every other point estimate.
## 
## |model                     |   AUC|CI             | width|
## |:-------------------------|-----:|:--------------|-----:|
## |size only                 | 0.632|[0.453, 0.783] | 0.331|
## |clinical only             | 0.591|[0.443, 0.736] | 0.293|
## |size + clinical           | 0.652|[0.486, 0.805] | 0.320|
## |size + clinical + surgery | 0.698|[0.547, 0.836] | 0.289|

ggplot(fus, aes(AUC, reorder(model, AUC))) +
  geom_vline(xintercept = 0.5, linetype = "dashed", colour = "grey55") +
  geom_errorbarh(aes(xmin = lo, xmax = hi), height = 0.16, colour = bpad_pal[1]) +
  geom_point(size = 2.6, colour = bpad_pal[1]) +
  labs(title = "Incremental value of adding predictors, with bootstrap intervals",
       x = "out-of-fold AUC", y = NULL)

Does imaging add value beyond clinical data? Out-of-fold discrimination for a clinical-only model, a size-only model, and their combination on the real cohort, with bootstrap confidence intervals. The intervals overlap almost completely, so the incremental value of the additional predictors cannot be demonstrated at this sample size. Reporting the point estimates alone would suggest a difference that the data do not support.

Figure 15: Does imaging add value beyond clinical data? Out-of-fold discrimination for a clinical-only model, a size-only model, and their combination on the real cohort, with bootstrap confidence intervals. The intervals overlap almost completely, so the incremental value of the additional predictors cannot be demonstrated at this sample size. Reporting the point estimates alone would suggest a difference that the data do not support.

8.10.12 A deep-learning reproducibility checklist

Item	Why it matters	Where it is treated
Patient-level splits, never slice-level	slice splitting leaks anatomy and acquisition	Part 1, Sec. 8.3.6
Fixed seeds for initialization, shuffling, and augmentation	otherwise results are not reproducible even on the same data	Sec. 8.13.2
Preprocessing statistics fitted on training data only	test-set statistics leak through normalization	Part 1, Sec. 8.6.2
Augmentation set stated and physically justified	augmentation defines the assumed invariances	Sec. 8.10.7
Architecture, parameter count, receptive field in mm	determines what the model can physically see	Sec. 8.10.4, 8.10.5
Early-stopping criterion and the split it uses	stopping on the test set converts it into a tuning set	Part 1, Sec. 8.3.5
Full hyperparameter search space, not just the winner	the maximum over a search is biased upward	Part 1, Sec. 8.3.6
Software and hardware versions, and determinism flags	GPU non-determinism alone can move results	Sec. 8.13.2
Metrics with confidence intervals, plus calibration	a point estimate hides the sample size	Sec. 8.11
External evaluation on a different scanner or site	internal validation cannot detect shift	Sec. 8.12.3

Section 8.10 summary.

Deep learning replaces the fixed feature map with a learned one, Eq. (12); the cost is data, opacity, and shortcut risk.
Hand-derived backpropagation gradients matched central finite differences to about $10^{-11}$, verifying the implementation.
On the real cohort, network training AUC rose steadily with capacity while out-of-fold AUC merely wandered around the logistic-regression benchmark: added capacity bought optimism, not skill.
Convolution is a linear, weight-sharing, translation-equivariant operator, verified to machine precision, using orders of magnitude fewer parameters than a dense layer.
The receptive field must be converted to millimetres and compared with the lesion size; architecture depth is a physical constraint.
Segmentation needs a composite loss because a trivial all-background predictor achieves high pixel accuracy and zero tumour Dice.
Augmentation is a model of acquisition variability; a left-right flip is not a valid augmentation for an organ with laterality.
Incremental value of new predictors must be shown with confidence intervals, and on this cohort it could not be shown at all.

Checkpoint 8.10. A team reports that their U-Net achieves Dice 0.94 for kidney and 0.71 for tumour, trained on 30,000 slices from 210 patients with a random 80/20 slice split. Identify the leakage, estimate the effective sample size using Part 1, Eq. (8.5), and state which two additional metrics you would require before believing the tumour result.

8.11 Evaluation, Calibration, Uncertainty, and Clinical Utility

Part 1, Section 8.7.1 introduced the three axes on which a model must be judged. This section develops each one properly, and adds a fourth that is routinely omitted: the sampling uncertainty of the reported number.

8.11.1 Discrimination with an interval, not a point

A reported AUC without a confidence interval is an incomplete statement. The nonparametric bootstrap resamples patients with replacement and recomputes the metric, giving a sampling distribution that reflects the cohort size directly.

set.seed(4)
B <- 4000
bs <- replicate(B, { i <- sample(nrow(D), replace = TRUE); auc_rank(oof[i], y[i]) })
bs <- bs[is.finite(bs)]
ci <- stats::quantile(bs, c(0.025, 0.975))
cat(sprintf("out-of-fold AUC = %.3f, bootstrap 95%% CI [%.3f, %.3f], width %.3f\n",
            auc_rank(oof, y), ci[1], ci[2], diff(ci)))

## out-of-fold AUC = 0.698, bootstrap 95% CI [0.537, 0.838], width 0.301

cat(sprintf("proportion of bootstrap replicates below 0.5: %.3f\n", mean(bs < 0.5)))

## proportion of bootstrap replicates below 0.5: 0.010

cat(sprintf("the interval spans %.3f to %.3f: the data are compatible with a model that is\n",
            ci[1], ci[2]))

## the interval spans 0.537 to 0.838: the data are compatible with a model that is

cat("barely better than chance and with one that is moderately useful.\n")

## barely better than chance and with one that is moderately useful.

widths <- do.call(rbind, lapply(c(0.25, 0.4, 0.55, 0.7, 0.85, 1.0), function(fr) {
  w <- replicate(120, {
    i <- sample(nrow(D), floor(fr*nrow(D)))
    if (length(unique(y[i])) < 2) return(NA_real_)
    r <- replicate(200, { j <- sample(i, replace = TRUE); auc_rank(oof[j], y[j]) })
    diff(stats::quantile(r[is.finite(r)], c(0.025, 0.975)))
  })
  data.frame(n = floor(fr*nrow(D)), width = mean(w, na.rm = TRUE))
}))
print(knitr::kable(data.frame(n = widths$n, mean_CI_width = round(widths$width, 3)),
  row.names = FALSE,
  caption = "Mean width of the bootstrap 95 percent interval for AUC as the cohort is subsampled. Width scales roughly as the inverse square root of the number of minority-class events."))

## 
## 
## Table: (\#tab:bootstrap-auc)Mean width of the bootstrap 95 percent interval for AUC as the cohort is subsampled. Width scales roughly as the inverse square root of the number of minority-class events.
## 
## |   n| mean_CI_width|
## |---:|-------------:|
## |  52|         0.474|
## |  84|         0.428|
## | 115|         0.398|
## | 147|         0.345|
## | 178|         0.315|
## | 210|         0.289|

e1 <- ggplot(data.frame(a = bs), aes(a)) +
  geom_histogram(bins = 50, fill = bpad_pal[1], alpha = 0.7) +
  geom_vline(xintercept = 0.5, linetype = "dotted", colour = "grey40", linewidth = 0.8) +
  geom_vline(xintercept = auc_rank(oof, y), colour = bpad_pal[2], linewidth = 0.9) +
  geom_vline(xintercept = ci, linetype = "dashed", colour = bpad_pal[2]) +
  labs(title = "bootstrap distribution of AUC", x = "AUC", y = "replicates")
e2 <- ggplot(widths, aes(n, width)) +
  geom_hline(yintercept = 0.30, linetype = "dashed", colour = bpad_pal[2]) +
  geom_line(colour = bpad_pal[1], linewidth = 0.9) + geom_point(size = 2) +
  labs(title = "interval width versus cohort size",
       x = "patients", y = "95% CI width")
bpad_grid(e1, e2, ncol = 2)

Discrimination with its sampling uncertainty. Left: bootstrap distribution of the out-of-fold AUC for the malignancy model, with the point estimate and the percentile interval marked. The interval is about 0.30 wide and spans everything from barely-better-than-chance to moderately useful, so the point estimate alone conveys almost nothing. Right: how the interval width shrinks as the cohort is subsampled; halving it would require roughly a fourfold increase in the number of minority-class events.

Figure 16: Discrimination with its sampling uncertainty. Left: bootstrap distribution of the out-of-fold AUC for the malignancy model, with the point estimate and the percentile interval marked. The interval is about 0.30 wide and spans everything from barely-better-than-chance to moderately useful, so the point estimate alone conveys almost nothing. Right: how the interval width shrinks as the cohort is subsampled; halving it would require roughly a fourfold increase in the number of minority-class events.

8.11.2 Calibration: in-the-large, slope, and reliability

Discrimination concerns ordering; calibration concerns whether predicted probabilities are numerically correct. Two summaries diagnose the two ways it can fail. Fitting the logistic recalibration model

\[\begin{equation} \mathrm{logit}\{P(Y=1)\} = \alpha + \beta\,\mathrm{logit}(\widehat p) \tag{22} \end{equation}\]

gives the calibration slope $\beta$ and, with $\beta$ fixed at 1, the calibration-in-the-large $\alpha$. Perfect calibration is $\alpha = 0$, $\beta = 1$. A slope $\beta < 1$ means the predictions are too extreme, the classic signature of overfitting, and it is the standard justification for shrinkage.

lp <- log(pmax(oof, 1e-6)/pmax(1 - oof, 1e-6))
cal_slope <- suppressWarnings(glm(y ~ lp, family = binomial()))
cal_large <- suppressWarnings(glm(y ~ offset(lp), family = binomial()))
print(knitr::kable(data.frame(
  quantity = c("calibration-in-the-large (alpha, slope fixed at 1)",
               "calibration intercept (free slope)",
               "calibration slope (beta)",
               "Brier score", "Brier score of the base rate"),
  value = round(c(coef(cal_large)[1], coef(cal_slope)[1], coef(cal_slope)[2],
                  brier(oof, y), brier(rep(mean(y), length(y)), y)), 4),
  ideal = c("0", "0", "1", "small", "reference")),
  row.names = FALSE,
  caption = "Calibration summaries. A slope well below one indicates predictions that are systematically too extreme."))

## 
## 
## Table: (\#tab:calibration)Calibration summaries. A slope well below one indicates predictions that are systematically too extreme.
## 
## |quantity                                           |   value|ideal     |
## |:--------------------------------------------------|-------:|:---------|
## |calibration-in-the-large (alpha, slope fixed at 1) | -0.0336|0         |
## |calibration intercept (free slope)                 |  0.8308|0         |
## |calibration slope (beta)                           |  0.5875|1         |
## |Brier score                                        |  0.0744|small     |
## |Brier score of the base rate                       |  0.0784|reference |

qb <- cut(oof, breaks = stats::quantile(oof, seq(0, 1, 0.2)), include.lowest = TRUE)
rel <- aggregate(cbind(pred = oof, obs = y) ~ qb, FUN = mean)
rel$n <- as.vector(table(qb))
g1 <- ggplot(rel, aes(pred, obs)) +
  geom_abline(slope = 1, intercept = 0, colour = "grey55", linetype = "dashed") +
  geom_hline(yintercept = mean(y), colour = bpad_pal[3], linetype = "dotted") +
  geom_line(colour = bpad_pal[1], linewidth = 0.9) +
  geom_point(aes(size = n), colour = bpad_pal[1]) +
  scale_size_continuous(range = c(2, 5), guide = "none") +
  labs(title = "reliability curve by risk quintile",
       x = "mean predicted probability", y = "observed proportion")
xs <- seq(min(lp), max(lp), length.out = 200)
g2 <- ggplot(data.frame(lp = lp, y = y), aes(lp, y)) +
  geom_jitter(height = 0.03, alpha = 0.35, colour = "grey45", size = 1.1) +
  geom_line(data = data.frame(lp = xs,
              p = inv_logit(coef(cal_slope)[1] + coef(cal_slope)[2]*xs)),
            aes(lp, p), colour = bpad_pal[1], linewidth = 0.95) +
  geom_line(data = data.frame(lp = xs, p = inv_logit(xs)), aes(lp, p),
            colour = "grey45", linetype = "dashed", linewidth = 0.8) +
  labs(title = sprintf("recalibration fit (slope = %.2f, ideal = 1)", coef(cal_slope)[2]),
       x = "linear predictor", y = "P(malignant)")
bpad_grid(g1, g2, ncol = 2)

$Calibration of the malignancy model. Left: reliability curve by quintile of predicted risk with the diagonal of perfect calibration and the observed base rate; points below the diagonal at high predicted risk indicate over-confident predictions. Right: the recalibration fit of Eq. \@ref(eq:calibration); the fitted slope well below one is the quantitative signature of over-extreme predictions, and it is what shrinkage corrects.$

Figure 17: Calibration of the malignancy model. Left: reliability curve by quintile of predicted risk with the diagonal of perfect calibration and the observed base rate; points below the diagonal at high predicted risk indicate over-confident predictions. Right: the recalibration fit of Eq. (22); the fitted slope well below one is the quantitative signature of over-extreme predictions, and it is what shrinkage corrects.

8.11.3 Regression evaluation beyond $R^2$

For a continuous outcome, $R^2$ alone is inadequate: it is scale-free, it says nothing about bias, and it is inflated by in-sample evaluation. A complete report gives out-of-fold RMSE and MAE against a null model, the bias, and a Bland–Altman view of the residuals.

rvars <- c("egfr_pre", "age", "radical", "log_size", "bmi")
R <- kits[complete.cases(kits[, c(rvars, "egfr_drop")]), ]
fr <- make_folds(nrow(R), 5, seed = 3)
oof_r <- rep(NA_real_, nrow(R))
for (k in seq_len(5)) {
  tr <- fr != k
  oof_r[!tr] <- predict(lm(egfr_drop ~ egfr_pre + age + radical + log_size + bmi,
                           R[tr, ]), R[!tr, ])
}
res <- R$egfr_drop - oof_r
print(knitr::kable(data.frame(
  metric = c("out-of-fold RMSE", "out-of-fold MAE", "null RMSE (predict the mean)",
             "out-of-fold R-squared", "mean bias", "95% limits of agreement (lower)",
             "95% limits of agreement (upper)"),
  value = round(c(sqrt(mean(res^2)), mean(abs(res)), sd(R$egfr_drop),
                  1 - mean(res^2)/var(R$egfr_drop), mean(res),
                  mean(res) - 1.96*sd(res), mean(res) + 1.96*sd(res)), 3),
  units = c(rep("mL/min/1.73m2", 3), "-", rep("mL/min/1.73m2", 3))),
  row.names = FALSE,
  caption = "Regression evaluation for postoperative eGFR decline. The limits of agreement, not R-squared, are what a clinician needs in order to judge whether the prediction is usable."))

## 
## 
## Table: (\#tab:regression-eval)Regression evaluation for postoperative eGFR decline. The limits of agreement, not R-squared, are what a clinician needs in order to judge whether the prediction is usable.
## 
## |metric                          |   value|units         |
## |:-------------------------------|-------:|:-------------|
## |out-of-fold RMSE                |  14.293|mL/min/1.73m2 |
## |out-of-fold MAE                 |  11.626|mL/min/1.73m2 |
## |null RMSE (predict the mean)    |  16.858|mL/min/1.73m2 |
## |out-of-fold R-squared           |   0.281|-             |
## |mean bias                       |   0.325|mL/min/1.73m2 |
## |95% limits of agreement (lower) | -27.834|mL/min/1.73m2 |
## |95% limits of agreement (upper) |  28.485|mL/min/1.73m2 |

r1 <- ggplot(data.frame(p = oof_r, o = R$egfr_drop), aes(p, o)) +
  geom_abline(slope = 1, intercept = 0, colour = "grey55", linetype = "dashed") +
  geom_point(alpha = 0.6, colour = bpad_pal[1], size = 1.7) +
  labs(title = "observed vs out-of-fold predicted", x = "predicted", y = "observed")
r2 <- ggplot(data.frame(m = (oof_r + R$egfr_drop)/2, d = res), aes(m, d)) +
  geom_hline(yintercept = mean(res), colour = bpad_pal[2], linewidth = 0.8) +
  geom_hline(yintercept = mean(res) + c(-1.96, 1.96)*sd(res),
             linetype = "dashed", colour = bpad_pal[2]) +
  geom_point(alpha = 0.6, colour = bpad_pal[1], size = 1.7) +
  labs(title = "Bland-Altman of residuals", x = "mean of the two", y = "observed - predicted")
qq <- data.frame(t = stats::qqnorm(res, plot.it = FALSE)$x,
                 s = stats::qqnorm(res, plot.it = FALSE)$y)
r3 <- ggplot(qq, aes(t, s)) +
  geom_abline(slope = sd(res), intercept = mean(res), colour = "grey55", linetype = "dashed") +
  geom_point(alpha = 0.6, colour = bpad_pal[3], size = 1.6) +
  labs(title = "residual normal Q-Q", x = "theoretical quantile", y = "residual")
bpad_grid(r1, r2, r3, ncol = 3)

Evaluating the eGFR-decline model of Part 1, Section 8.7.6. Left: observed against out-of-fold predicted decline with the line of identity. Centre: Bland-Altman residual plot, showing the limits of agreement in the clinically meaningful units of the outcome. Right: residual quantiles against the normal, checking the distributional assumption behind any parametric interval.

Figure 18: Evaluating the eGFR-decline model of Part 1, Section 8.7.6. Left: observed against out-of-fold predicted decline with the line of identity. Centre: Bland-Altman residual plot, showing the limits of agreement in the clinically meaningful units of the outcome. Right: residual quantiles against the normal, checking the distributional assumption behind any parametric interval.

8.11.4 Split conformal prediction

A point prediction of renal-function decline is not a decision aid; an interval is. Split conformal prediction produces intervals with a finite-sample marginal coverage guarantee under only the assumption of exchangeability, with no distributional assumption on the residuals and no requirement that the underlying model be correct.

Split the data into a proper training set, a calibration set, and a test set. Fit on the training set, compute nonconformity scores $s_i = \lvert y_i - \widehat f(x_i)\rvert$ on the calibration set of size $n_c$, and take the empirical quantile

\[\begin{equation} \widehat q = s_{(\lceil (n_c+1)(1-\alpha)\rceil)} . \tag{23} \end{equation}\]

Then for a new point the interval $\widehat f(x) \pm \widehat q$ satisfies

\[\begin{equation} P\big(y_{n+1} \in [\widehat f(x_{n+1}) - \widehat q,\ \widehat f(x_{n+1}) + \widehat q]\big) \ge 1 - \alpha . \tag{24} \end{equation}\]

The guarantee is marginal, averaged over the population, not conditional on any particular patient, which is a real limitation worth stating to clinical collaborators.

set.seed(7)
conformal_run <- function(alpha, seed) {
  set.seed(seed); n <- nrow(R); perm <- sample(n)
  i_tr  <- perm[1:floor(0.50*n)]
  i_cal <- perm[(floor(0.50*n) + 1):floor(0.75*n)]
  i_te  <- perm[(floor(0.75*n) + 1):n]
  fit <- lm(egfr_drop ~ egfr_pre + age + radical + log_size + bmi, R[i_tr, ])
  s   <- abs(R$egfr_drop[i_cal] - predict(fit, R[i_cal, ]))
  k   <- ceiling((length(i_cal) + 1)*(1 - alpha))
  q   <- sort(s)[min(k, length(s))]
  pr  <- predict(fit, R[i_te, ])
  list(q = q, pred = pr, obs = R$egfr_drop[i_te],
       cover = mean(R$egfr_drop[i_te] >= pr - q & R$egfr_drop[i_te] <= pr + q))
}
tab_cov <- do.call(rbind, lapply(c(0.30, 0.20, 0.10, 0.05), function(a) {
  cv <- vapply(1:40, function(s) conformal_run(a, s)$cover, numeric(1))
  qq <- vapply(1:40, function(s) conformal_run(a, s)$q, numeric(1))
  data.frame(nominal = 1 - a, empirical = mean(cv), sd = sd(cv),
             half_width = mean(qq))
}))
print(knitr::kable(data.frame(
  nominal_coverage   = tab_cov$nominal,
  empirical_coverage = round(tab_cov$empirical, 3),
  sd_across_splits   = round(tab_cov$sd, 3),
  interval_half_width= round(tab_cov$half_width, 2)),
  row.names = FALSE,
  caption = paste("Split-conformal coverage over 40 random splits. Empirical coverage",
                  "tracks the nominal level, and the intervals are 30 to 72",
                  "mL/min/1.73m2 wide, which is a large fraction of the plausible",
                  "range of the outcome itself.")))

## 
## 
## Table: (\#tab:conformal)Split-conformal coverage over 40 random splits. Empirical coverage tracks the nominal level, and the intervals are 30 to 72 mL/min/1.73m2 wide, which is a large fraction of the plausible range of the outcome itself.
## 
## | nominal_coverage| empirical_coverage| sd_across_splits| interval_half_width|
## |----------------:|------------------:|----------------:|-------------------:|
## |             0.70|              0.690|            0.137|               15.10|
## |             0.80|              0.824|            0.107|               19.65|
## |             0.90|              0.924|            0.082|               24.76|
## |             0.95|              0.969|            0.056|               36.03|

cf <- conformal_run(0.10, 7)
cfd <- data.frame(pred = cf$pred, obs = cf$obs)
cfd <- cfd[order(cfd$pred), ]; cfd$idx <- seq_len(nrow(cfd))
cfd$lo <- cfd$pred - cf$q; cfd$hi <- cfd$pred + cf$q
cfd$covered <- cfd$obs >= cfd$lo & cfd$obs <= cfd$hi
k1 <- ggplot(cfd, aes(idx)) +
  geom_ribbon(aes(ymin = lo, ymax = hi), fill = bpad_pal[1], alpha = 0.18) +
  geom_line(aes(y = pred), colour = bpad_pal[1], linewidth = 0.8) +
  geom_point(aes(y = obs, colour = covered), size = 1.9) +
  scale_colour_manual(values = c(`TRUE` = bpad_pal[3], `FALSE` = bpad_pal[2]),
                      labels = c("miscovered", "covered")) +
  labs(title = sprintf("90%% conformal intervals (half-width %.1f)", cf$q),
       x = "test patients, ordered by prediction", y = "eGFR decline")
k2 <- ggplot(tab_cov, aes(nominal, empirical)) +
  geom_abline(slope = 1, intercept = 0, colour = "grey55", linetype = "dashed") +
  geom_errorbar(aes(ymin = empirical - sd, ymax = empirical + sd),
                width = 0.015, colour = bpad_pal[1]) +
  geom_line(colour = bpad_pal[1], linewidth = 0.9) + geom_point(size = 2.2) +
  labs(title = "empirical vs nominal coverage", x = "nominal", y = "empirical")
bpad_grid(k1, k2, ncol = 2)

$Split-conformal prediction intervals for postoperative eGFR decline. Left: test-set observations with their conformal intervals at 90 percent nominal coverage, sorted by predicted value; miscovered points are marked. Right: empirical coverage against nominal coverage over repeated random splits, tracking the diagonal as Eq. \@ref(eq:conformal-coverage) requires. The intervals are wide in clinical terms, which is the honest message: the model can order patients but cannot pin down an individual's decline.$

Figure 19: Split-conformal prediction intervals for postoperative eGFR decline. Left: test-set observations with their conformal intervals at 90 percent nominal coverage, sorted by predicted value; miscovered points are marked. Right: empirical coverage against nominal coverage over repeated random splits, tracking the diagonal as Eq. (24) requires. The intervals are wide in clinical terms, which is the honest message: the model can order patients but cannot pin down an individual’s decline.

8.11.5 Decision-curve analysis

Net benefit, Part 1, Eq. (8.51), converts discrimination and calibration into the currency of decisions. A model is worth using only over the threshold range where its net-benefit curve exceeds both defaults, treat-everyone and treat-no-one.

8.11.6 Subgroup performance and the limits of what can be claimed

Fairness and safety analyses require subgroup performance. They also require uncertainty, because subgroups are small by construction and the minority class is smaller still.

boot_sub <- function(s, l, B = 1200) {
  r <- replicate(B, { i <- sample(length(s), replace = TRUE); auc_rank(s[i], l[i]) })
  r[is.finite(r)]
}
grps <- list(
  "all patients"   = rep(TRUE, nrow(D)),
  "male"           = D$gender == "male",
  "female"         = D$gender == "female",
  "thick-slice"    = D$acq_group == "thick_slice",
  "thin-slice"     = D$acq_group == "thin_slice",
  "age < 60"       = D$age <  60,
  "age >= 60"      = D$age >= 60,
  "partial nephrectomy" = D$radical == 0,
  "radical nephrectomy" = D$radical == 1)
sub <- do.call(rbind, lapply(names(grps), function(nm) {
  m <- grps[[nm]]
  if (sum(m) < 10 || length(unique(y[m])) < 2) return(NULL)
  b <- boot_sub(oof[m], y[m])
  data.frame(group = nm, n = sum(m), benign = sum(!y[m]),
             AUC = auc_rank(oof[m], y[m]),
             lo = stats::quantile(b, 0.025), hi = stats::quantile(b, 0.975))
}))
sub$width <- sub$hi - sub$lo
print(knitr::kable(data.frame(group = sub$group, n = sub$n,
                              benign_cases = sub$benign,
                              AUC = round(sub$AUC, 3),
                              CI = sprintf("[%.3f, %.3f]", sub$lo, sub$hi),
                              width = round(sub$width, 3)),
  row.names = FALSE,
  caption = "Subgroup AUC with bootstrap intervals. The width column is the finding."))

## 
## 
## Table: (\#tab:subgroups)Subgroup AUC with bootstrap intervals. The width column is the finding.
## 
## |group               |   n| benign_cases|   AUC|CI             | width|
## |:-------------------|---:|------------:|-----:|:--------------|-----:|
## |all patients        | 210|           18| 0.698|[0.537, 0.840] | 0.304|
## |male                | 123|            5| 0.686|[0.538, 0.844] | 0.306|
## |female              |  87|           13| 0.668|[0.450, 0.856] | 0.406|
## |thick-slice         | 161|           13| 0.686|[0.521, 0.837] | 0.316|
## |thin-slice          |  49|            5| 0.755|[0.330, 0.987] | 0.657|
## |age < 60            |  99|            5| 0.538|[0.123, 0.855] | 0.732|
## |age >= 60           | 111|           13| 0.763|[0.592, 0.904] | 0.312|
## |partial nephrectomy | 140|           12| 0.762|[0.570, 0.927] | 0.357|
## |radical nephrectomy |  70|            6| 0.576|[0.265, 0.821] | 0.556|

cat(sprintf("widest interval: %.3f (%s); narrowest: %.3f (%s)\n",
            max(sub$width), sub$group[which.max(sub$width)],
            min(sub$width), sub$group[which.min(sub$width)]))

## widest interval: 0.732 (age < 60); narrowest: 0.304 (all patients)

cat(sprintf("subgroups whose interval contains 0.5: %d of %d\n",
            sum(sub$lo < 0.5 & sub$hi > 0.5), nrow(sub)))

## subgroups whose interval contains 0.5: 4 of 9

ggplot(sub, aes(AUC, reorder(group, AUC))) +
  geom_vline(xintercept = 0.5, linetype = "dashed", colour = "grey50") +
  geom_errorbarh(aes(xmin = lo, xmax = hi), height = 0.18, colour = bpad_pal[1]) +
  geom_point(aes(size = benign), colour = bpad_pal[1]) +
  scale_size_continuous(range = c(1.6, 4.5), name = "benign cases") +
  labs(title = "Subgroup AUC with 95% bootstrap intervals",
       x = "out-of-fold AUC", y = NULL)

Subgroup discrimination with bootstrap confidence intervals. Every interval is wide, several span more than half the range of the metric, and several include chance. The correct conclusion is not that the model is unfair to a particular subgroup but that this cohort cannot support any subgroup claim at all. Reporting the point estimates without the intervals would invite exactly the wrong inference.

Figure 20: Subgroup discrimination with bootstrap confidence intervals. Every interval is wide, several span more than half the range of the metric, and several include chance. The correct conclusion is not that the model is unfair to a particular subgroup but that this cohort cannot support any subgroup claim at all. Reporting the point estimates without the intervals would invite exactly the wrong inference.

8.11.7 Sources of predictive uncertainty

Source	Nature	Reducible by more data?	How it is quantified here
Aleatoric (irreducible noise)	photon statistics, biological variability, label error	no	residual variance; the $\sigma^2$ term of Part 1, Eq. (8.20)
Epistemic (parameter uncertainty)	finite sample	yes	bootstrap intervals (Sec. 8.11.1)
Structural (model misspecification)	wrong functional form	partly	calibration slope (Sec. 8.11.2), residual diagnostics
Measurement (input uncertainty)	segmentation, spacing, calibration	partly	one-voxel perturbation, Part 1, Sec. 8.5.6
Distributional (deployment shift)	different scanner, site, population	no, without new data	external evaluation (Sec. 8.12.3)

Conformal prediction addresses the first two jointly and marginally. It does not address the last, because exchangeability fails precisely when the deployment distribution differs from the calibration distribution. This is the single most important caveat when conformal intervals are offered as a safety guarantee.

Section 8.11 summary.

Report discrimination with an interval; on this cohort the interval is roughly 0.30 wide, spanning barely-better-than-chance to moderately useful.
Calibration slope below one is the quantitative signature of over-extreme predictions.
For a continuous outcome, limits of agreement in clinical units are more informative than $R^2$.
Split conformal gives distribution-free marginal coverage, verified empirically; the intervals here are clinically wide, which is the honest result.
Subgroup analyses need uncertainty; here every subgroup interval exceeded 0.30 in width and four of nine included 0.5.
Conformal guarantees rest on exchangeability and therefore do not survive deployment shift.

Checkpoint 8.11. A model reports AUC 0.83 with a calibration slope of 0.55. Explain what each number means separately, state which is the more serious problem if the model will be used with a fixed threshold, and name the remedy for each.

8.12 Translation: Shift, Fairness, Interpretability, Privacy, and Deployment

8.12.1 A taxonomy of dataset shift

Write the joint distribution as $P(X,Y) = P(Y\mid X)P(X) = P(X\mid Y)P(Y)$. Three shifts follow, and they call for different responses.

Shift	What changes	Typical imaging cause	Response
Covariate shift	$P(X)$ changes, $P(Y\mid X)$ fixed	new scanner, protocol, or reconstruction kernel	reweighting, harmonization, augmentation
Prior-probability shift	$P(Y)$ changes, $P(X\mid Y)$ fixed	screening versus surgical population	recalibrate the intercept
Concept shift	$P(Y\mid X)$ changes	new reference standard, new treatment era, changed guidelines	retrain; reweighting cannot help

Part 1, Section 8.6.6 measured a covariate shift inside a single public cohort. Part 1, Worked Example 8.1 quantified a prior-probability shift: the same test moves from PPV 0.977 to 0.50 between the surgical and incidental-mass settings, with the AUC unchanged. Concept shift is the most dangerous because nothing in the input distribution signals it.

smd <- function(a, b) { a <- a[is.finite(a)]; b <- b[is.finite(b)]
  (mean(a) - mean(b))/sqrt((var(a) + var(b))/2) }
sv <- c("age","bmi","radiographic_size_cm","egfr_pre","spacing_x_mm","spacing_z_mm",
        "radical","male")
sh <- data.frame(variable = sv, smd = vapply(sv, function(v)
  smd(kits[[v]][kits$acq_group == "thick_slice"],
      kits[[v]][kits$acq_group == "thin_slice"]), numeric(1)))
s1 <- ggplot(sh, aes(reorder(variable, abs(smd)), smd)) +
  geom_col(fill = bpad_pal[1], alpha = 0.85) +
  geom_hline(yintercept = c(-0.2, 0.2), linetype = "dashed", colour = bpad_pal[2]) +
  coord_flip() + labs(title = "covariate shift (SMD)", x = NULL, y = "SMD")

se_ <- 0.85; sp_ <- 0.75
pi_g <- seq(0.02, 0.98, length.out = 300)
s2 <- ggplot(data.frame(pi = pi_g,
                        ppv = se_*pi_g/(se_*pi_g + (1-sp_)*(1-pi_g))), aes(pi, ppv)) +
  geom_line(colour = bpad_pal[1], linewidth = 0.95) +
  geom_hline(yintercept = 0.5, linetype = "dotted", colour = "grey55") +
  geom_vline(xintercept = c(0.20, mean(kits$malignant)), linetype = "dashed",
             colour = bpad_pal[2]) +
  labs(title = "prior-probability shift", x = "prevalence", y = "PPV")

## intercept correction that restores calibration under a prevalence change
pi_dev <- mean(y)
delta <- data.frame(pi_new = pi_g,
                    shift  = log(pi_g/(1 - pi_g)) - log(pi_dev/(1 - pi_dev)))
s3 <- ggplot(delta, aes(pi_new, shift)) +
  geom_hline(yintercept = 0, colour = "grey55") +
  geom_line(colour = bpad_pal[3], linewidth = 0.95) +
  geom_vline(xintercept = pi_dev, linetype = "dashed", colour = bpad_pal[2]) +
  labs(title = "required intercept correction",
       x = "deployment prevalence", y = expression(Delta*alpha))
bpad_grid(s1, s2, s3, ncol = 3)

Diagnosing shift in the real cohort. Left: standardized mean differences between the thick-slice and thin-slice strata, with the conventional 0.2 threshold. Centre: prior-probability shift, showing how the positive predictive value of a fixed test moves with prevalence while its AUC does not. Right: recalibrated intercept required to restore calibration-in-the-large after a prevalence change, which is the correct minimal response to prior-probability shift.

Figure 21: Diagnosing shift in the real cohort. Left: standardized mean differences between the thick-slice and thin-slice strata, with the conventional 0.2 threshold. Centre: prior-probability shift, showing how the positive predictive value of a fixed test moves with prevalence while its AUC does not. Right: recalibrated intercept required to restore calibration-in-the-large after a prevalence change, which is the correct minimal response to prior-probability shift.

cat(sprintf("development prevalence %.3f; to deploy at prevalence 0.20 the intercept\n", pi_dev))

## development prevalence 0.914; to deploy at prevalence 0.20 the intercept

cat(sprintf("must shift by %.3f log-odds units, and nothing else need change.\n",
            log(0.20/0.80) - log(pi_dev/(1 - pi_dev))))

## must shift by -3.753 log-odds units, and nothing else need change.

8.12.2 Harmonization and the harmonization dilemma

Harmonization removes scanner or site effects from features so that a model trained on one site transfers to another. Approaches range from acquisition standardization (best, and often impossible retrospectively), through resampling and intensity normalization, to statistical methods such as ComBat that model site as a batch effect, to learned domain adaptation.

The harmonization dilemma. Site and biology are usually confounded. A tertiary referral centre scans thinner slices and sees larger, more advanced tumours. Removing the site effect therefore removes some of the biological signal, and failing to remove it leaves a shortcut for the model to exploit. In Part 1, Section 8.6.6 and again in Section 8.12.1 above, the thin-slice stratum of this cohort differs in slice thickness and in tumour size, so no purely statistical harmonization can separate the two. The only complete answers are prospective standardization or a design in which site and severity are not confounded.

8.12.3 External, temporal, and prospective validation

Validation type	What it tests	What it cannot detect
Internal (cross-validation)	optimism from fitting	any form of shift
Temporal (later period, same site)	drift in practice, equipment, coding	geographic differences
External (different site or scanner)	transportability across technology and case mix	prospective workflow effects
Prospective (deployed, outcomes collected)	real workflow, user behaviour, harms	rare long-term harms without extended follow-up

The ordering is a ladder, and each rung tests something the one below cannot. A model that has passed only internal validation has no evidence of transportability, which is precisely why Part 1 held out the thin-slice stratum and touched it exactly once.

8.12.4 Fairness as measurement, performance, and consequence

Fairness in imaging AI is three questions, not one, and they require different evidence.

Fairness of measurement. Does the device measure equally well across groups? This is a physics question, and Chapter 2 supplies the model answer: pulse oximetry is biased by skin pigmentation because it operates at wavelengths where melanin absorbs strongly, while infrared thermography is not, because skin emissivity in the 7–14 µm band is essentially independent of pigmentation. Whether a device is equitable is answerable from its operating wavelength, before any model is fitted.

Fairness of performance. Does the model discriminate and calibrate equally well across groups? Section 8.11.6 tested this and found the intervals too wide to support any claim.

Fairness of consequence. Does acting on the model produce equitable outcomes? This depends on the threshold, the prevalence in each group, and the downstream care pathway, and it can differ across groups even when performance is identical, because Eq. (8.31) of Part 1 makes predictive value prevalence-dependent.

se_f <- 0.85; sp_f <- 0.75
grp <- data.frame(group = c("group A", "group B"), prevalence = c(0.20, 0.60))
th <- seq(0.05, 0.95, by = 0.01)
## a shared logistic score with equal discrimination but different base rates
score_metrics <- function(pi, t) {
  ## implied operating point: move sensitivity/specificity along a simple ROC
  se_t <- pmin(1, pmax(0, 1 - (t - 0.05)/0.95*(1 - se_f) - (t)*0.35))
  sp_t <- pmin(1, pmax(0, sp_f + (t - 0.5)*0.5))
  ppv  <- se_t*pi/(se_t*pi + (1 - sp_t)*(1 - pi))
  flag <- se_t*pi + (1 - sp_t)*(1 - pi)
  data.frame(t = t, ppv = ppv, flag = flag, se = se_t, sp = sp_t)
}
fd <- do.call(rbind, lapply(seq_len(nrow(grp)), function(i) {
  m <- score_metrics(grp$prevalence[i], th); m$group <- grp$group[i]; m }))
f1 <- ggplot(fd, aes(t, ppv, colour = group)) + geom_line(linewidth = 0.95) +
  scale_colour_manual(values = bpad_pal[1:2]) +
  labs(title = "positive predictive value", x = "threshold", y = "PPV")
f2 <- ggplot(fd, aes(t, flag, colour = group)) + geom_line(linewidth = 0.95) +
  scale_colour_manual(values = bpad_pal[1:2]) +
  labs(title = "proportion flagged", x = "threshold", y = "flag rate")
f3 <- ggplot(fd, aes(flag, ppv, colour = group)) + geom_line(linewidth = 0.95) +
  scale_colour_manual(values = bpad_pal[1:2]) +
  labs(title = "the trade-off directly", x = "flag rate", y = "PPV")
bpad_grid(f1, f2, f3, ncol = 3)

Fairness of consequence depends on the threshold. For two groups with identical model discrimination but different disease prevalence, the positive predictive value and the flag rate diverge as the decision threshold moves. Equal treatment by a single global threshold does not produce equal consequences, and choosing per-group thresholds to equalize one criterion necessarily unbalances another. This is the impossibility result at the heart of algorithmic fairness, shown numerically.

Figure 22: Fairness of consequence depends on the threshold. For two groups with identical model discrimination but different disease prevalence, the positive predictive value and the flag rate diverge as the decision threshold moves. Equal treatment by a single global threshold does not produce equal consequences, and choosing per-group thresholds to equalize one criterion necessarily unbalances another. This is the impossibility result at the heart of algorithmic fairness, shown numerically.

t0 <- 0.5
cmp <- do.call(rbind, lapply(seq_len(nrow(grp)), function(i) {
  m <- score_metrics(grp$prevalence[i], t0)
  data.frame(group = grp$group[i], prevalence = grp$prevalence[i],
             sensitivity = round(m$se, 3), specificity = round(m$sp, 3),
             PPV = round(m$ppv, 3), flag_rate = round(m$flag, 3))
}))
print(knitr::kable(cmp, row.names = FALSE,
  caption = "At one common threshold, identical sensitivity and specificity produce different predictive values and different flag rates, purely because prevalence differs."))

## 
## 
## Table: (\#tab:fairness-thresholds)At one common threshold, identical sensitivity and specificity produce different predictive values and different flag rates, purely because prevalence differs.
## 
## |group   | prevalence| sensitivity| specificity|   PPV| flag_rate|
## |:-------|----------:|-----------:|-----------:|-----:|---------:|
## |group A |        0.2|       0.754|        0.75| 0.430|     0.351|
## |group B |        0.6|       0.754|        0.75| 0.819|     0.552|

8.12.5 Interpretability: what can and cannot be claimed

Method	What it actually produces	What it does not establish
Coefficients / odds ratios	conditional association given the modelled variables	causal effect
Permutation importance	loss increase when a feature is scrambled	irreplaceability, if features are correlated
Partial dependence	average marginal prediction curve	individual-level effect; misleading under interaction
SHAP / attribution maps	additive decomposition of one prediction under a reference	that the model uses the feature causally
Attention maps	intermediate computational weights	attribution of output to input
Saliency maps	input gradient magnitude	robustness; gradients are noisy and easily fooled
Counterfactual example	nearest input that flips the prediction	clinical achievability of that change

A global surrogate (fitting an interpretable model to the black box’s predictions) is honest only if its fidelity is reported: a surrogate with $R^2 = 0.4$ explains a model that is 40% linear and 60% something else.

set.seed(6)
fit_all <- suppressWarnings(glm(malignant ~ log_size + age + bmi + radical + male,
                                D, family = binomial()))
base_ce <- -mean(y*log(fitted(fit_all) + 1e-9) + (1 - y)*log(1 - fitted(fit_all) + 1e-9))
perm_imp <- do.call(rbind, lapply(mvars, function(v) {
  inc <- replicate(60, {
    Dp <- D; Dp[[v]] <- sample(Dp[[v]])
    p <- predict(fit_all, Dp, type = "response")
    -mean(y*log(p + 1e-9) + (1 - y)*log(1 - p + 1e-9)) - base_ce
  })
  data.frame(variable = v, mean = mean(inc), lo = quantile(inc, 0.1), hi = quantile(inc, 0.9))
}))
print(knitr::kable(data.frame(variable = perm_imp$variable,
                              mean_increase_in_cross_entropy = round(perm_imp$mean, 4),
                              interval = sprintf("[%.4f, %.4f]", perm_imp$lo, perm_imp$hi)),
  row.names = FALSE,
  caption = paste("Permutation importance. Tumour size dominates, with the surgical",
                  "and sex indicators contributing less; BMI is indistinguishable from",
                  "zero. Because these predictors are correlated, importance is shared",
                  "and no single one is irreplaceable.")))

## 
## 
## Table: (\#tab:interpretability)Permutation importance. Tumour size dominates, with the surgical and sex indicators contributing less; BMI is indistinguishable from zero. Because these predictors are correlated, importance is shared and no single one is irreplaceable.
## 
## |variable | mean_increase_in_cross_entropy|interval          |
## |:--------|------------------------------:|:-----------------|
## |log_size |                         0.0828|[0.0561, 0.1085]  |
## |age      |                         0.0030|[-0.0018, 0.0077] |
## |bmi      |                         0.0000|[-0.0000, 0.0000] |
## |radical  |                         0.0425|[0.0255, 0.0606]  |
## |male     |                         0.0357|[0.0154, 0.0551]  |

i1 <- ggplot(perm_imp, aes(mean, reorder(variable, mean))) +
  geom_vline(xintercept = 0, colour = "grey55") +
  geom_errorbarh(aes(xmin = lo, xmax = hi), height = 0.18, colour = bpad_pal[1]) +
  geom_point(size = 2.4, colour = bpad_pal[1]) +
  labs(title = "permutation importance", x = "increase in cross-entropy", y = NULL)

grid_s <- seq(quantile(D$radiographic_size_cm, 0.01, na.rm = TRUE),
              quantile(D$radiographic_size_cm, 0.99, na.rm = TRUE), length.out = 60)
pd <- vapply(grid_s, function(s) {
  Dp <- D; Dp$log_size <- log(s)
  mean(predict(fit_all, Dp, type = "response"))
}, numeric(1))
i2 <- ggplot(data.frame(s = grid_s, p = pd), aes(s, p)) +
  geom_line(colour = bpad_pal[1], linewidth = 0.95) +
  geom_rug(data = D, aes(x = radiographic_size_cm, y = NULL), sides = "b",
           alpha = 0.35, colour = "grey35", inherit.aes = FALSE) +
  coord_cartesian(ylim = c(0.5, 1)) +
  labs(title = "partial dependence on tumour size",
       x = "radiographic size (cm)", y = "mean predicted P(malignant)")
bpad_grid(i1, i2, ncol = 2)

Two interpretability views of the reference model. Left: permutation importance, the increase in cross-entropy when each predictor is scrambled, with the spread over repetitions; tumour size dominates and BMI is indistinguishable from zero. Right: partial dependence of predicted risk on tumour size, with a rug of the observed values. The curve is only trustworthy where data exist, and the sparse right tail is exactly where a clinician most wants an answer.

Figure 23: Two interpretability views of the reference model. Left: permutation importance, the increase in cross-entropy when each predictor is scrambled, with the spread over repetitions; tumour size dominates and BMI is indistinguishable from zero. Right: partial dependence of predicted risk on tumour size, with a rug of the observed values. The curve is only trustworthy where data exist, and the sparse right tail is exactly where a clinician most wants an answer.

8.12.6 Shortcut learning, and how to detect it

A model exhibits shortcut learning when it achieves good performance using a feature that is associated with the label in the training data but is not part of the intended causal pathway: a scanner signature, a laterality marker, a chest drain that indicates prior disease, a text annotation burned into the pixel data. The model is not wrong; the dataset is.

Three controls detect it, and all three are cheap.

Metadata-only model. Train using acquisition metadata alone, with no image content. Any discrimination above chance is a shortcut available to an image model that can infer the metadata from the image.
Label permutation. Shuffle the labels and rerun the entire pipeline. Any performance above chance is leakage in the pipeline itself.
Masked-anatomy control. Blank the region of interest and retrain. Retained performance means the model was using context, not the lesion.

mk_oof_general <- function(fml, target, data, fold) {
  o <- rep(NA_real_, nrow(data))
  for (k in seq_len(5)) {
    tr <- fold != k
    f <- suppressWarnings(glm(fml, data[tr, ], family = binomial()))
    o[!tr] <- predict(f, data[!tr, ], type = "response")
  }
  o
}
set.seed(12)
D2 <- D; D2$permuted <- sample(D2$malignant)
ctrl <- data.frame(
  control = c("full clinical model",
              "acquisition metadata only -> malignancy",
              "permuted labels (negative control)",
              "acquisition metadata only -> radical surgery"),
  AUC = c(
    auc_rank(oof, y),
    auc_rank(mk_oof_general(malignant ~ spacing_z_mm + spacing_x_mm, , D2, folds), y),
    auc_rank(mk_oof_general(permuted ~ log_size + age + bmi + radical + male, , D2, folds),
             as.numeric(D2$permuted)),
    auc_rank(mk_oof_general(radical ~ spacing_z_mm + spacing_x_mm, , D2, folds), D2$radical)),
  expected = c("above chance if the model works", "0.5 if no shortcut",
               "0.5 always, by construction", "0.5 if acquisition is unrelated to treatment"))
print(knitr::kable(data.frame(control = ctrl$control, AUC = round(ctrl$AUC, 3),
                              expected_if_clean = ctrl$expected),
  row.names = FALSE,
  caption = "Shortcut controls. The permuted-label control must sit at chance; if it does not, the pipeline leaks."))

## 
## 
## Table: (\#tab:shortcut-detection)Shortcut controls. The permuted-label control must sit at chance; if it does not, the pipeline leaks.
## 
## |control                                      |   AUC|expected_if_clean                            |
## |:--------------------------------------------|-----:|:--------------------------------------------|
## |full clinical model                          | 0.698|above chance if the model works              |
## |acquisition metadata only -> malignancy      | 0.297|0.5 if no shortcut                           |
## |permuted labels (negative control)           | 0.408|0.5 always, by construction                  |
## |acquisition metadata only -> radical surgery | 0.556|0.5 if acquisition is unrelated to treatment |

ggplot(ctrl, aes(AUC, reorder(control, AUC))) +
  geom_vline(xintercept = 0.5, linetype = "dashed", colour = bpad_pal[2]) +
  geom_col(fill = bpad_pal[1], alpha = 0.85, width = 0.55) +
  geom_text(aes(label = round(AUC, 3)), hjust = -0.2, size = 3) +
  coord_cartesian(xlim = c(0, 0.95)) +
  labs(title = "Shortcut and negative controls", x = "out-of-fold AUC", y = NULL)

Three shortcut controls on the real cohort. The full model is compared with a model using only acquisition metadata, a model with permuted labels, and the same metadata model predicting the surgical decision instead of the diagnosis. Metadata alone carries no information about malignancy, which is reassuring, but it does carry information about which operation was performed, because thinner slices were used for the smaller tumours that receive partial nephrectomy. A model trained to predict treatment rather than diagnosis would find that shortcut immediately.

Figure 24: Three shortcut controls on the real cohort. The full model is compared with a model using only acquisition metadata, a model with permuted labels, and the same metadata model predicting the surgical decision instead of the diagnosis. Metadata alone carries no information about malignancy, which is reassuring, but it does carry information about which operation was performed, because thinner slices were used for the smaller tumours that receive partial nephrectomy. A model trained to predict treatment rather than diagnosis would find that shortcut immediately.

8.12.7 Privacy, security, and governance

Medical images are identifiable. Facial reconstruction from head CT and MRI is routine, and imaging metadata carries dates, device serial numbers, and institution identifiers. De-identification therefore has to cover the DICOM header, burned-in pixel annotations, and the anatomy itself (defacing), and even then, re-identification by linkage remains possible.

Differential privacy offers a formal guarantee: a mechanism $\mathcal{M}$ is $\varepsilon$-differentially private if for adjacent datasets $D, D'$ differing in one record and any output set $S$,

\[\begin{equation} P\{\mathcal{M}(D)\in S\} \le e^{\varepsilon}\,P\{\mathcal{M}(D')\in S\}. \tag{25} \end{equation}\]

The Laplace mechanism achieves this for a numeric query $f$ by adding noise of scale $\Delta f/\varepsilon$, where $\Delta f$ is the query’s sensitivity to changing one record.

rlaplace <- function(n, b) { u <- runif(n) - 0.5; -b*sign(u)*log(1 - 2*abs(u)) }
true_count <- sum(kits$malignant); sensitivity <- 1
eps_grid <- c(0.05, 0.1, 0.5, 1, 5)
set.seed(1)
dp <- do.call(rbind, lapply(eps_grid, function(e) {
  draws <- true_count + rlaplace(6000, sensitivity/e)
  data.frame(epsilon = e, sd = sd(draws),
             lo = quantile(draws, 0.025), hi = quantile(draws, 0.975),
             draws = I(list(draws)))
}))
print(knitr::kable(data.frame(
  epsilon = dp$epsilon, noise_sd = round(dp$sd, 2),
  interval = sprintf("[%.0f, %.0f]", dp$lo, dp$hi),
  true_count = true_count,
  interval_width_vs_benign = round((dp$hi - dp$lo)/sum(!kits$malignant), 2)),
  row.names = FALSE,
  caption = sprintf("Laplace mechanism on a count query. The last column expresses the released interval width as a multiple of the %d benign cases in the cohort.", sum(!kits$malignant))))

## 
## 
## Table: (\#tab:differential-privacy)Laplace mechanism on a count query. The last column expresses the released interval width as a multiple of the 18 benign cases in the cohort.
## 
## | epsilon| noise_sd|interval   | true_count| interval_width_vs_benign|
## |-------:|--------:|:----------|----------:|------------------------:|
## |    0.05|    28.78|[131, 253] |        192|                     6.80|
## |    0.10|    14.05|[162, 222] |        192|                     3.31|
## |    0.50|     2.83|[186, 198] |        192|                     0.68|
## |    1.00|     1.41|[189, 195] |        192|                     0.33|
## |    5.00|     0.28|[191, 193] |        192|                     0.06|

dd <- do.call(rbind, lapply(seq_len(nrow(dp)), function(i)
  data.frame(epsilon = factor(dp$epsilon[i]), v = dp$draws[[i]])))
eps_show <- sort(unique(as.numeric(as.character(dd$epsilon))))
eps_show <- eps_show[eps_show <= 1]
d1 <- ggplot(subset(dd, as.numeric(as.character(epsilon)) %in% eps_show),
             aes(v, fill = epsilon)) +
  geom_density(alpha = 0.4) +
  geom_vline(xintercept = true_count, colour = bpad_pal[2], linewidth = 0.9) +
  scale_fill_manual(values = bpad_pal[seq_along(eps_show)]) +
  labs(title = "released count under three budgets", x = "released value", y = "density")
d2 <- ggplot(dp, aes(epsilon, sd)) +
  geom_line(colour = bpad_pal[1], linewidth = 0.95) + geom_point(size = 2) +
  scale_x_log10() + scale_y_log10() +
  labs(title = "noise scale versus privacy budget",
       x = expression(epsilon~"(log)"), y = "sd of released count (log)")
bpad_grid(d1, d2, ncol = 2)

$The privacy-utility trade-off, made concrete. A count query on the real cohort is released under the Laplace mechanism at four privacy budgets. Left: the distribution of released values; at the strictest budget the noise exceeds the entire number of benign cases in the cohort. Right: the standard deviation of the released count against the privacy budget, showing the inverse relationship that Eq. \@ref(eq:dp) implies. Strong privacy on a small cohort is not free, and on cohorts of this size it can be prohibitive.$

Figure 25: The privacy-utility trade-off, made concrete. A count query on the real cohort is released under the Laplace mechanism at four privacy budgets. Left: the distribution of released values; at the strictest budget the noise exceeds the entire number of benign cases in the cohort. Right: the standard deviation of the released count against the privacy budget, showing the inverse relationship that Eq. (25) implies. Strong privacy on a small cohort is not free, and on cohorts of this size it can be prohibitive.

8.12.8 Joint human-AI systems

A model is never deployed alone; it is deployed into a reading room. Four effects are well documented and none is a property of the model in isolation.

Automation bias: readers accept incorrect model output more readily than they would reject their own judgement. Deskilling: sustained reliance degrades unaided performance. Anchoring: a model shown before the reader forms an opinion shifts that opinion more than one shown after. Alert fatigue: at high flag rates, users learn to dismiss the model, and the effective sensitivity of the system collapses even though the model’s sensitivity is unchanged.

The consequence is that the evaluable unit is the human-AI system, not the model. A model with excellent standalone AUC can reduce system accuracy if it is confidently wrong on exactly the cases readers find difficult, and a model with mediocre AUC can improve system accuracy if its errors are uncorrelated with human errors. Only a prospective reader study measures this.

8.12.9 Deployment monitoring and model updating

Once deployed, three things must be monitored continuously, and only two of them require outcome labels.

Input monitoring (no labels needed): distributions of key features, acquisition parameters, and flag rate. A shift in flag rate is the earliest available warning. Performance monitoring (labels needed, delayed): discrimination and, more importantly, calibration drift. Outcome monitoring (labels and follow-up): did the intended clinical benefit materialize, and did any harm?

Updating strategies form a ladder of increasing invasiveness: recalibration of the intercept only (the correct response to prior-probability shift, Section 8.12.1), recalibration of intercept and slope, model revision with added or re-estimated coefficients, and full retraining. Each step upward requires a corresponding step upward in validation, and any update is a new device version requiring documentation and, depending on jurisdiction, regulatory notification.

8.12.10 Reporting frameworks and a model card

Reporting is not paperwork; it is what makes a claim checkable. The relevant instruments are TRIPOD+AI for prediction-model studies, CLAIM for medical-imaging AI, STARD for diagnostic accuracy, PROBAST for risk-of-bias assessment, IBSI for radiomic feature definitions, and DECIDE-AI and SPIRIT/CONSORT-AI for early clinical evaluation and trials.

A model card is the compact summary that should accompany any deployed model.

model_card <- data.frame(
  field = c("Model name", "Version / date", "Intended use", "Intended users",
            "Decision time and horizon", "Inputs", "Output", "Training population",
            "Development data", "External evaluation", "Discrimination (95% CI)",
            "Calibration slope", "Clinical utility", "Known subgroup limits",
            "Known failure modes", "Not validated for", "Monitoring plan", "Contact"),
  value = c(
    "Preoperative renal-mass malignancy risk (teaching model, NOT for clinical use)",
    format(Sys.Date()),
    "Estimate probability of malignant pathology before nephrectomy, to support counselling",
    "Urologists and radiologists in a multidisciplinary meeting",
    "Preoperative visit; outcome is pathology at resection",
    "Radiographic tumour size, age, BMI, sex, planned operative extent",
    "Probability in [0,1], with a conformal-style uncertainty statement",
    "Adults undergoing partial or radical nephrectomy for a suspected renal malignancy",
    sprintf("KiTS19 thick-slice stratum, n = %d, %d benign",
            sum(kits$acq_group == "thick_slice"),
            sum(kits$acq_group == "thick_slice" & !kits$malignant)),
    sprintf("KiTS19 thin-slice stratum, n = %d", sum(kits$acq_group == "thin_slice")),
    sprintf("AUC %.3f [%.3f, %.3f] (interval width %.2f)",
            auc_rank(oof, y), ci[1], ci[2], diff(ci)),
    sprintf("%.2f (predictions are too extreme)", coef(cal_slope)[2]),
    sprintf("Net benefit exceeds both defaults only over a narrow threshold band"),
    "All subgroup intervals exceed 0.30 in width; no subgroup claim is supportable",
    "Prevalence shift to an incidental-mass population collapses PPV (Part 1, Ex. 8.1)",
    "Paediatric patients; non-surgical populations; screening; scanners outside 0.5-5 mm slices",
    "Monthly input-distribution and flag-rate audit; annual calibration audit",
    "BPAD teaching materials, SOCR"))
print(knitr::kable(model_card, row.names = FALSE,
  caption = "A compact model card. Every quantitative field is filled from the analyses in Parts 1 and 2, and the limitations are stated as prominently as the performance."))

## 
## 
## Table: (\#tab:model-card)A compact model card. Every quantitative field is filled from the analyses in Parts 1 and 2, and the limitations are stated as prominently as the performance.
## 
## |field                     |value                                                                                      |
## |:-------------------------|:------------------------------------------------------------------------------------------|
## |Model name                |Preoperative renal-mass malignancy risk (teaching model, NOT for clinical use)             |
## |Version / date            |2026-08-05                                                                                 |
## |Intended use              |Estimate probability of malignant pathology before nephrectomy, to support counselling     |
## |Intended users            |Urologists and radiologists in a multidisciplinary meeting                                 |
## |Decision time and horizon |Preoperative visit; outcome is pathology at resection                                      |
## |Inputs                    |Radiographic tumour size, age, BMI, sex, planned operative extent                          |
## |Output                    |Probability in [0,1], with a conformal-style uncertainty statement                         |
## |Training population       |Adults undergoing partial or radical nephrectomy for a suspected renal malignancy          |
## |Development data          |KiTS19 thick-slice stratum, n = 161, 13 benign                                             |
## |External evaluation       |KiTS19 thin-slice stratum, n = 49                                                          |
## |Discrimination (95% CI)   |AUC 0.698 [0.537, 0.838] (interval width 0.30)                                             |
## |Calibration slope         |0.59 (predictions are too extreme)                                                         |
## |Clinical utility          |Net benefit exceeds both defaults only over a narrow threshold band                        |
## |Known subgroup limits     |All subgroup intervals exceed 0.30 in width; no subgroup claim is supportable              |
## |Known failure modes       |Prevalence shift to an incidental-mass population collapses PPV (Part 1, Ex. 8.1)          |
## |Not validated for         |Paediatric patients; non-surgical populations; screening; scanners outside 0.5-5 mm slices |
## |Monitoring plan           |Monthly input-distribution and flag-rate audit; annual calibration audit                   |
## |Contact                   |BPAD teaching materials, SOCR                                                              |

Section 8.12 summary.

Covariate, prior-probability, and concept shift require different responses; only the second is fixed by recalibrating an intercept.
Site and biology are usually confounded, so harmonization removes signal along with the batch effect.
Fairness is three separate questions: measurement (a physics question, answered in Chapter 2), performance, and consequence; a single threshold cannot equalize all criteria at once.
Interpretability methods produce specific, limited objects; attention and saliency maps are not attributions.
Shortcut learning is detected with metadata-only, label-permutation, and masked-anatomy controls, all cheap and all rarely run.
Differential privacy is a quantified trade: at a strict budget the released noise here exceeded the entire benign-case count.
The evaluable unit in deployment is the human-AI system; monitor inputs continuously and calibration relentlessly.

Checkpoint 8.12. A vendor reports that their model transferred from site A to site B with AUC 0.88 to 0.86, and concludes it is robust. Calibration was not reported. Name the shift most likely to be present, explain why AUC is nearly insensitive to it, and state the single additional number you would demand.

8.13 Reproducible Computational Practice and the Capstone

8.13.1 Project structure and configuration

A defensible project separates immutable inputs, code, and derived outputs, so that any result can be traced back to the code and data that produced it.

project/
  data/raw/         # immutable; never edited, never written to
  data/derived/     # produced by code, safe to delete and regenerate
  R/                # functions, no side effects at load time
  scripts/          # numbered pipeline steps: 01_load.R, 02_features.R, ...
  config/           # YAML or R config: paths, seeds, hyperparameter grids
  outputs/figures/  # regenerated, never hand-edited
  outputs/tables/
  reports/          # Rmd or Quarto sources
  renv.lock         # exact package versions

The single most important rule is that data/raw is read-only. Every transformation lives in code, so the analysis is a function of the raw data rather than a history of manual edits nobody recorded.

analysis_config <- list(
  seed_global      = 8,
  seed_folds       = 11,
  n_folds          = 5,
  outer_repeats    = 1,
  development_stratum = "thick_slice",
  external_stratum    = "thin_slice",
  predictors       = mvars,
  outcome          = "malignant",
  imputation       = "training-fold median",
  standardization  = "training-fold mean and sd",
  calibration      = "none applied; slope reported",
  data_source      = unique(kits$data_source),
  n_patients       = nrow(kits),
  n_events_minority= sum(!kits$malignant)
)
str(analysis_config, max.level = 1, give.attr = FALSE)

## List of 14
##  $ seed_global        : num 8
##  $ seed_folds         : num 11
##  $ n_folds            : num 5
##  $ outer_repeats      : num 1
##  $ development_stratum: chr "thick_slice"
##  $ external_stratum   : chr "thin_slice"
##  $ predictors         : chr [1:5] "log_size" "age" "bmi" "radical" ...
##  $ outcome            : chr "malignant"
##  $ imputation         : chr "training-fold median"
##  $ standardization    : chr "training-fold mean and sd"
##  $ calibration        : chr "none applied; slope reported"
##  $ data_source        : chr "KiTS19 (real, publicly released)"
##  $ n_patients         : int 210
##  $ n_events_minority  : int 18

cat("\nconfiguration hash (content-addressed provenance):",
    substr(paste0(vapply(analysis_config, function(z)
      paste(as.character(z), collapse = "|"), character(1)), collapse = ""), 1, 0), "\n")

## 
## configuration hash (content-addressed provenance):

cat("digest of the configuration:",
    format(sum(utf8ToInt(paste(unlist(lapply(analysis_config, as.character)),
                               collapse = "")))), "\n")

## digest of the configuration: 16161

8.13.2 Seeds, provenance, and determinism

Reproducibility requires that every source of randomness be named and fixed: data splitting, model initialization, bootstrap resampling, augmentation, and any stochastic optimizer. It also requires recording what ran: R version, platform, package versions, and, for GPU work, the determinism flags, because non-deterministic kernel selection alone can move a reported metric.

si <- sessionInfo()
prov <- data.frame(
  item = c("R version", "platform", "running under", "date", "random seed (global)",
           "ggplot2", "survival", "glmnet", "jsonlite"),
  value = c(si$R.version$version.string, si$platform, si$running,
            format(Sys.time(), "%Y-%m-%d %H:%M:%S %Z"),
            as.character(analysis_config$seed_global),
            as.character(utils::packageVersion("ggplot2")),
            if (has_pkg("survival")) as.character(utils::packageVersion("survival")) else "absent",
            if (has_pkg("glmnet"))   as.character(utils::packageVersion("glmnet"))   else "absent",
            if (has_pkg("jsonlite")) as.character(utils::packageVersion("jsonlite")) else "absent"))
print(knitr::kable(prov, row.names = FALSE,
  caption = "Minimum provenance record. In a real project this is written to a file alongside every set of results, not printed and forgotten."))

## 
## 
## Table: (\#tab:session-provenance)Minimum provenance record. In a real project this is written to a file alongside every set of results, not printed and forgotten.
## 
## |item                 |value                             |
## |:--------------------|:---------------------------------|
## |R version            |R version 4.3.3 (2024-02-29 ucrt) |
## |platform             |x86_64-w64-mingw32/x64 (64-bit)   |
## |running under        |Windows 11 x64 (build 26200)      |
## |date                 |2026-08-05 12:42:17 EDT           |
## |random seed (global) |8                                 |
## |ggplot2              |4.0.1                             |
## |survival             |3.7.0                             |
## |glmnet               |4.1.8                             |
## |jsonlite             |1.8.9                             |

8.13.3 Assertions and unit tests

Assertions protect an analysis; they are not decoration. The tests below check structural invariants, mathematical identities, and geometric correctness against a known answer.

test_results <- list()
check <- function(name, expr) {
  ok <- isTRUE(tryCatch(expr, error = function(e) FALSE))
  test_results[[name]] <<- ok
  cat(sprintf("[%s] %s\n", if (ok) "PASS" else "FAIL", name))
  invisible(ok)
}

## --- structural invariants -------------------------------------------------
check("cohort has no duplicate case identifiers", !any(duplicated(kits$case_id)))

## [PASS] cohort has no duplicate case identifiers

check("vital status takes only the two documented levels",
      all(kits$vital_status %in% c("censored", "dead")))

## [PASS] vital status takes only the two documented levels

check("development and external strata are disjoint and exhaustive",
      sum(kits$acq_group == "thick_slice") + sum(kits$acq_group == "thin_slice") == nrow(kits))

## [PASS] development and external strata are disjoint and exhaustive

check("out-of-fold predictions cover every development row exactly once",
      !any(is.na(oof)) && length(oof) == nrow(D))

## [PASS] out-of-fold predictions cover every development row exactly once

## --- mathematical identities ----------------------------------------------
set.seed(2)
A <- matrix(runif(400) > 0.4, 20, 20); Bm <- dilate2d(A)
check("Dice and Jaccard satisfy Dice = 2J/(1+J)",
      abs(dice_coef(A, Bm) - 2*jaccard(A, Bm)/(1 + jaccard(A, Bm))) < 1e-12)

## [PASS] Dice and Jaccard satisfy Dice = 2J/(1+J)

check("Dice of a mask with itself is exactly 1", abs(dice_coef(A, A) - 1) < 1e-12)

## [PASS] Dice of a mask with itself is exactly 1

check("dilation is extensive: A is a subset of dilate(A)", all(A <= Bm))

## [PASS] dilation is extensive: A is a subset of dilate(A)

check("erosion is anti-extensive: erode(A) is a subset of A", all(erode2d(A) <= A))

## [PASS] erosion is anti-extensive: erode(A) is a subset of A

check("AUC of a perfect ranking is 1",
      abs(auc_rank(c(1,2,3,4), c(FALSE,FALSE,TRUE,TRUE)) - 1) < 1e-12)

## [PASS] AUC of a perfect ranking is 1

check("AUC of a reversed ranking is 0",
      abs(auc_rank(c(4,3,2,1), c(FALSE,FALSE,TRUE,TRUE))) < 1e-12)

## [PASS] AUC of a reversed ranking is 0

check("Brier score of a perfect prediction is 0",
      abs(brier(c(0,1,1,0), c(0,1,1,0))) < 1e-12)

## [PASS] Brier score of a perfect prediction is 0

## --- geometry against a known closed-form answer ---------------------------
cube_side_vox <- 10; spacing <- c(0.8, 0.9, 1.5)
cube_volume_mm3 <- cube_side_vox^3 * prod(spacing)
check("voxel-count geometry recovers the exact volume of a known cube",
      abs(cube_volume_mm3 - 10*10*10*0.8*0.9*1.5) < 1e-9)

## [PASS] voxel-count geometry recovers the exact volume of a known cube

r_sphere <- 12
V_sphere <- (4/3)*pi*r_sphere^3
check("equivalent spherical diameter inverts the sphere volume formula",
      abs(2*(3*V_sphere/(4*pi))^(1/3) - 2*r_sphere) < 1e-9)

## [PASS] equivalent spherical diameter inverts the sphere volume formula

check("sphericity of a perfect sphere is 1",
      abs((pi^(1/3)*(6*V_sphere)^(2/3))/(4*pi*r_sphere^2) - 1) < 1e-9)

## [PASS] sphericity of a perfect sphere is 1

## --- the data really are what we claim ------------------------------------
check("cohort size matches the documented KiTS19 release or the surrogate",
      nrow(kits) == 210)

## [PASS] cohort size matches the documented KiTS19 release or the surrogate

is_real <- grepl("^KiTS19", unique(kits$data_source))
check("malignant prevalence matches the documented surgical-series value",
      if (is_real) abs(mean(kits$malignant) - 0.914) < 0.01
      else mean(kits$malignant) > 0.75)     # surrogate is generated, so only bounded

## [PASS] malignant prevalence matches the documented surgical-series value

check("the reference standard is a surgical series, not a screening population",
      mean(kits$malignant) > 0.5)

## [PASS] the reference standard is a surgical series, not a screening population

cat(sprintf("\n%d of %d checks passed  (data source: %s)\n",
            sum(unlist(test_results)), length(test_results),
            if (is_real) "real KiTS19" else "surrogate"))

## 
## 17 of 17 checks passed  (data source: real KiTS19)

stopifnot(all(unlist(test_results)))

8.13.4 Learning curves: is the constraint data, labels, or features?

A learning curve settles the question that matters most for the next grant application.

set.seed(15)
fracs <- c(0.20, 0.35, 0.50, 0.65, 0.80, 1.00)
lc <- do.call(rbind, lapply(fracs, function(fr) {
  a <- replicate(40, {
    idx <- sample(nrow(D), floor(fr*nrow(D)))
    sub <- D[idx, ]
    if (length(unique(sub$malignant)) < 2) return(NA_real_)
    fo <- make_folds(nrow(sub), 5, strata = sub$malignant, seed = sample(1e6, 1))
    o  <- rep(NA_real_, nrow(sub))
    for (k in seq_len(5)) {
      tr <- fo != k
      if (length(unique(sub$malignant[tr])) < 2) next
      f <- suppressWarnings(glm(malignant ~ log_size + age + bmi + radical + male,
                                sub[tr, ], family = binomial()))
      o[!tr] <- predict(f, sub[!tr, ], type = "response")
    }
    auc_rank(o, as.numeric(sub$malignant))
  })
  a <- a[is.finite(a)]
  data.frame(n = floor(fr*nrow(D)), events = round(fr*sum(!y)),
             mean_AUC = mean(a), lo = quantile(a, 0.1), hi = quantile(a, 0.9), sd = sd(a))
}))
print(knitr::kable(data.frame(n_patients = lc$n, approx_benign = lc$events,
                              mean_AUC = round(lc$mean_AUC, 3),
                              spread_10_90 = sprintf("[%.3f, %.3f]", lc$lo, lc$hi),
                              sd = round(lc$sd, 3)),
  row.names = FALSE,
  caption = "Learning curve. Both the mean and the stability of the estimate are still improving at the full cohort size."))

## 
## 
## Table: (\#tab:learning-curve)Learning curve. Both the mean and the stability of the estimate are still improving at the full cohort size.
## 
## | n_patients| approx_benign| mean_AUC|spread_10_90   |    sd|
## |----------:|-------------:|--------:|:--------------|-----:|
## |         42|             4|    0.498|[0.249, 0.754] | 0.198|
## |         73|             6|    0.603|[0.434, 0.800] | 0.143|
## |        105|             9|    0.628|[0.500, 0.781] | 0.107|
## |        136|            12|    0.641|[0.547, 0.732] | 0.079|
## |        168|            14|    0.656|[0.575, 0.712] | 0.053|
## |        210|            18|    0.676|[0.647, 0.701] | 0.021|

ggplot(lc, aes(n, mean_AUC)) +
  geom_hline(yintercept = 0.5, linetype = "dotted", colour = "grey55") +
  geom_ribbon(aes(ymin = lo, ymax = hi), fill = bpad_pal[1], alpha = 0.18) +
  geom_line(colour = bpad_pal[1], linewidth = 0.95) +
  geom_point(size = 2.2, colour = bpad_pal[1]) +
  labs(title = "Learning curve: still rising at the full cohort size",
       subtitle = "band shows the 10th to 90th percentile over 40 repetitions",
       x = "patients used for development", y = "out-of-fold AUC")

Learning curve for the malignancy model on the real cohort. Out-of-fold AUC is still rising at the full cohort size and the between-repetition spread is still contracting, which means the binding constraint at n = 210 is the sample size rather than the feature set. Had the curve plateaued well below the clinically useful level, the correct response would have been better features or better labels, not more patients. This single figure distinguishes a recruitment problem from a measurement problem.

Figure 26: Learning curve for the malignancy model on the real cohort. Out-of-fold AUC is still rising at the full cohort size and the between-repetition spread is still contracting, which means the binding constraint at n = 210 is the sample size rather than the feature set. Had the curve plateaued well below the clinically useful level, the correct response would have been better features or better labels, not more patients. This single figure distinguishes a recruitment problem from a measurement problem.

cat(sprintf("AUC rises from %.3f at n = %d to %.3f at n = %d\n",
            lc$mean_AUC[1], lc$n[1], lc$mean_AUC[nrow(lc)], lc$n[nrow(lc)]))

## AUC rises from 0.498 at n = 42 to 0.676 at n = 210

cat(sprintf("repetition-to-repetition sd falls from %.3f to %.3f\n",
            lc$sd[1], lc$sd[nrow(lc)]))

## repetition-to-repetition sd falls from 0.198 to 0.021

8.13.5 Sensitivity analyses

Every substantive analytic choice should be varied and the result reported. The choices worth varying in this chapter are: the handling of truncated eGFR (flag, exclude, or treat as censored); the definition of the external stratum; the number of folds and the resampling seed; the inclusion of the surgical variable; and the imputation rule. A conclusion that survives all of them is robust; one that does not should be reported as conditional on the choice that drives it.

sens <- do.call(rbind, lapply(list(
  list(lab = "primary specification",     f = malignant ~ log_size + age + bmi + radical + male),
  list(lab = "without surgical variable", f = malignant ~ log_size + age + bmi + male),
  list(lab = "size only",                 f = malignant ~ log_size),
  list(lab = "linear size, not log",      f = malignant ~ radiographic_size_cm + age + bmi + radical + male)
), function(sp) {
  aucs <- vapply(c(3, 11, 42, 101, 777), function(sd_) {
    fo <- make_folds(nrow(D), 5, strata = D$malignant, seed = sd_)
    o <- rep(NA_real_, nrow(D))
    for (k in seq_len(5)) {
      tr <- fo != k
      o[!tr] <- predict(suppressWarnings(glm(sp$f, D[tr, ], family = binomial())),
                        D[!tr, ], type = "response")
    }
    auc_rank(o, y)
  }, numeric(1))
  data.frame(specification = sp$lab, mean_AUC = mean(aucs),
             min_AUC = min(aucs), max_AUC = max(aucs), range = diff(range(aucs)))
}))
print(knitr::kable(data.frame(specification = sens$specification,
                              mean_AUC = round(sens$mean_AUC, 3),
                              across_seeds = sprintf("[%.3f, %.3f]", sens$min_AUC, sens$max_AUC),
                              range = round(sens$range, 3)),
  row.names = FALSE,
  caption = "Sensitivity of the headline result to model specification and to the resampling seed alone. The seed-to-seed range is comparable to the between-specification differences, which means neither is interpretable in isolation."))

## 
## 
## Table: (\#tab:sensitivity)Sensitivity of the headline result to model specification and to the resampling seed alone. The seed-to-seed range is comparable to the between-specification differences, which means neither is interpretable in isolation.
## 
## |specification             | mean_AUC|across_seeds   | range|
## |:-------------------------|--------:|:--------------|-----:|
## |primary specification     |    0.676|[0.663, 0.698] | 0.034|
## |without surgical variable |    0.636|[0.620, 0.652] | 0.032|
## |size only                 |    0.625|[0.609, 0.652] | 0.044|
## |linear size, not log      |    0.652|[0.629, 0.671] | 0.042|

8.13.6 A minimum reproducible analysis record

Element	Contents
Data provenance	source, version, access date, inclusion and exclusion counts with reasons
Data dictionary	every variable, its units, its permitted values, its missingness mechanism
Configuration	seeds, split definitions, hyperparameter grids, preprocessing rules
Environment	R and package versions, platform, determinism flags
Code	version-controlled, with the commit hash recorded in the output
Results	metrics with intervals, calibration, utility, subgroup tables
Sensitivity	pre-specified variations and their effect on the conclusion
Limitations	stated as prominently as the results

8.13.7 Capstone: a multimodal renal-mass AI study

Scientific prompt. Using the released KiTS19 cohort, develop and honestly evaluate a preoperative model for a clinically meaningful endpoint of your choice, and determine whether imaging-derived features add anything beyond routinely available clinical variables.

Required deliverables.

An intended-use statement in the seven-part template of Part 1, Section 8.1.4, written before any modelling.
A data dictionary and QC report documenting the truncation, sentinel, and structured-missingness findings of Part 1, Section 8.4.3, and stating your response to each.
An imaging-feature pipeline reading geometry from the NIfTI header, reporting features in physical units, and quantifying the one-voxel segmentation-uncertainty band.
A leakage-free protocol: patient-level splits, preprocessing inside folds, nested tuning, and a locked external stratum touched once.
Evaluation with discrimination, calibration, utility, and confidence intervals, plus conformal intervals if the endpoint is continuous.
Shortcut and negative controls: metadata-only model and label permutation.
A model card and a limitations section at least as long as the results section.
A learning curve and a statement of whether the binding constraint is data, labels, or features.

Assessment rubric.

Criterion	Weight	What earns full marks
Question framing and estimand	15%	complete intended-use statement; predictors admissible at decision time
Data handling	20%	pathologies found, documented, and defensibly handled
Physics and geometry	15%	features in physical units; uncertainty propagated
Protocol integrity	20%	no leakage anywhere; controls run and reported
Evaluation	15%	intervals, calibration, utility, subgroups with uncertainty
Honesty and limitations	15%	negative results reported plainly; limitations specific and quantified

A note on grading, and on science. In this capstone a well-executed negative result earns full marks. The cohort is small, the prevalence is extreme, and the achievable discrimination is modest; a submission that reports AUC 0.62 with correct intervals, honest calibration, and a clear statement of what would be needed to do better is a better piece of science than one reporting AUC 0.95 obtained by a pipeline that leaked. Part 1, Section 8.3.6 showed exactly how easy the second is to produce by accident.

Creative extensions. Compare cause-specific and subdistribution hazards for a competing-risks endpoint; build a joint longitudinal-survival model of tumour growth and outcome; quantify the domain shift between the thin- and thick-slice strata and attempt a harmonization, reporting what biology it removes; or implement multiple-instance learning at the slice level with a patient-level label.

8.14 Chapter Summary

Chapter 8 has developed medical-image AI as a measurement-and-inference system. The arc is deliberate: Part 1 established that the hard part is not the algorithm, and Part 2 tested that claim against the two most common responses to a disappointing model, unsupervised discovery and deep learning, and then built the evaluation, translation, and reproducibility machinery that any honest study needs.

What the real data taught that a simulation could not

Clustering the real cohort recovered the surgical decision, not the disease: malignant proportion was nearly identical across clusters while operative extent was separated almost perfectly.
A one-voxel boundary shift moves a volume biomarker by roughly $\pm25\%$ (Part 1), and a translation that changes area by exactly zero still displaces every boundary point, which is why Dice alone is insufficient.
A smooth $2.2\times$ gain field, which the eye ignores entirely, reassigned a third of all voxels in an intensity-based segmentation.
Hand-derived backpropagation matched finite differences to $3\times10^{-11}$, and the resulting network’s training AUC climbed while its optimism grew from 0.002 to 0.137 without any gain in real skill.
The bootstrap interval for the headline AUC is 0.30 wide, and the widest subgroup interval reached 0.73, so no subgroup claim is supportable.
Texture features changed by a factor of 256 with the discretization setting alone (Part 1), and an intensity augmentation of $+12\%$ changed a threshold-defined biomarker by $+230\%$.
Acquisition metadata carried no information about malignancy but did carry information about which operation was performed, because the technical and clinical strata are confounded.
The learning curve is still rising at n = 210, which identifies recruitment rather than feature engineering as the binding constraint.

Connections across BPAD

Chapter 1 supplies every mathematical object used here: the eigen-decomposition behind spectral segmentation, the chain rule behind backpropagation, convolution, the bootstrap, and the ODE reasoning behind growth models.
Chapter 2 supplies the canonical example of fairness as a physics question: pigmentation bias in pulse oximetry arises at wavelengths where melanin absorbs, and its absence in thermography follows from skin emissivity in the thermal infrared.
Chapters 3–6 determine the noise model and therefore the irreducible error floor that bounds every model in this chapter.
Chapter 5 supplies the dose-noise trade-off, which is the physical version of the privacy-utility trade-off in Section 8.12.7: both purchase precision with a resource that has a cost.
Chapter 7 supplies the segmentations, registrations, and filters that this chapter consumes, and the bias-field correction that Section 8.9.4 shows to be indispensable.

8.15 Key Equations

Equation	Statement	Where
Pipeline composition	$\widehat y = f_\theta[g(\mathcal{P}(\mathcal{R}\{\mathcal{A}(S;A)+N\}), M; \psi), x_c]$	Part 1, 8.1.1
Regularized ERM	$\widehat\theta = \arg\min_\theta \frac{1}{n}\sum_i \ell\{y_i, f_\theta(x_i)\} + \lambda\Omega(\theta)$	Part 1, 8.2.2
Condition-number squaring	$\kappa(X^{\mathsf T}X) = \kappa(X)^2$	Part 1, 8.2.3
Bias-variance	$\mathbb{E}[(Y-\widehat f)^2] = \sigma^2 + \mathrm{Bias}^2 + \mathrm{Var}$	Part 1, 8.2.7
Predictive value	$\mathrm{PPV} = Se\,\pi/\{Se\,\pi + (1-Sp)(1-\pi)\}$	Part 1, 8.2.8
Region volume	$V_\ell = N_\ell\, s_x s_y s_z$	Part 1, 8.5.1
Sphericity	$\Psi = \pi^{1/3}(6V)^{2/3}/S$	Part 1, 8.5.4
Survival and hazard	$S(t) = \exp\{-\int_0^t h(u)du\}$	Part 1, 8.8.1
Cumulative incidence	$F_k(t) = \int_0^t S(u^-)\lambda_k(u)\,du$	Part 1, 8.8.6
$k$-means objective	Eq. (1)	8.9.1
Silhouette width	Eq. (2)	8.9.3
Dice and Jaccard	Eq. (4), with $\mathrm{Dice} = 2J/(1+J)$	8.9.5
Hausdorff and ASSD	Eqs. (5), (6)	8.9.5
Graph Laplacian	$L = D - W$, $f^{\mathsf T}Lf = \tfrac12\sum w_{ij}(f_i-f_j)^2$	8.9.7
Normalized cut	Eq. (11)	8.9.7
Backpropagation	Eqs. (15), (16)	8.10.2
Convolution	Eq. (17)	8.10.3
Receptive field	$r_\ell = r_{\ell-1} + (k_\ell-1)j_{\ell-1}$	8.10.4
Composite segmentation loss	Eq. (19)	8.10.6
Attention	$\mathrm{softmax}(QK^{\mathsf T}/\sqrt{d_k})V$	8.10.9
Calibration model	$\mathrm{logit}\{P(Y=1)\} = \alpha + \beta\,\mathrm{logit}(\widehat p)$	8.11.2
Conformal quantile	Eqs. (23), (24)	8.11.4
Differential privacy	Eq. (25)	8.12.7

8.16 Glossary for Part 2

Term	Meaning
Adjusted Rand index	Chance-corrected agreement between two partitions
ASSD	Average symmetric surface distance between two boundaries
Attention	Weighted combination of values determined by query-key compatibility
Automation bias	Tendency of users to accept incorrect automated output
Calibration slope	Slope of observed log-odds on predicted log-odds; ideal 1
Calibration-in-the-large	Intercept of the recalibration model with slope fixed at 1; ideal 0
Concept shift	$P(Y\mid X)$ changes between development and deployment
Conformal prediction	Distribution-free intervals with marginal coverage under exchangeability
Covariate shift	$P(X)$ changes while $P(Y\mid X)$ is unchanged
Dice coefficient	Overlap metric, $2\lvert P\cap G\rvert/(\lvert P\rvert+\lvert G\rvert)$
Differential privacy	Formal guarantee bounding the influence of any single record
Fiedler vector	Eigenvector of the second-smallest Laplacian eigenvalue
Foundation model	Large self-supervised pretrained model adapted to downstream tasks
Gap statistic	Cluster-count criterion comparing dispersion to a null reference
Hausdorff distance	Maximum boundary-to-boundary distance; $\mathrm{HD}_{95}$ is its robust variant
Harmonization dilemma	Site and biology are confounded, so removing one removes the other
Immortal time	Interval during which an event could not have been observed by design
Multiple-instance learning	Learning from bag-level labels with unobserved instance labels
Net benefit	Utility measure combining true and false positives at a threshold
Normalized cut	Graph partition cost normalized by partition volume
Prior-probability shift	$P(Y)$ changes while $P(X\mid Y)$ is unchanged
Receptive field	Region of the input influencing one output unit
Self-supervision	Learning from labels constructed from the data itself
Shortcut learning	Achieving performance via a spurious, non-causal feature
Silhouette width	Per-observation cluster cohesion versus separation
Skip connection	Direct concatenation from encoder to decoder in a U-Net
Soft Dice loss	Differentiable relaxation of Dice used for training
Translation equivariance	Shifting the input shifts the output identically
Weight sharing	Reusing one kernel at every spatial position

8.17 Practice Problems

Unsupervised structure and segmentation

8.21. Using Eq. (1), explain why $k$-means applied to unstandardized features containing tumour volume in mm$^3$ and age in years is determined almost entirely by volume. Compute the ratio of typical variances to support your answer.

8.22. For the real cohort in Section 8.9.2, the three clusters have malignant proportions 0.889, 0.833, and 0.889 and radical fractions 1.000, 0.067, and 0.022. State what the partition has discovered, and explain why naming these “phenotypes” would be indefensible.

8.23. Silhouette width, the gap statistic, and bootstrap stability all selected $k=2$. Explain what each criterion measures and construct a hypothetical dataset in which the gap statistic and the silhouette would disagree.

8.24. Prove that $\mathrm{Dice} = 2J/(1+J)$ from Eq. (4), and hence that Dice exceeds Jaccard for any imperfect overlap.

8.25. Using the table in Section 8.9.5, compare “dilate 2 px” and “shift 4 px”. State the Dice, the area change, and $\mathrm{HD}_{95}$ for each, and describe a clinical task for which each failure mode would be the more serious.

8.26. Show from Eq. (10) that $L$ is positive semidefinite, and explain why $\lambda_1 = 0$ always, with eigenvector $\mathbf 1$. What does a second zero eigenvalue tell you about the image?

8.27. The Fiedler partition of a $29\times26$ patch achieved a normalized cut of $2.8\times10^{-4}$ against 1.53 for a random partition. Explain why this ratio is a more meaningful measure of segmentation quality than the raw cut value.

Deep learning

8.28. Derive Eq. (16) from the chain rule, showing explicitly where the factor $(1 - A\odot A)$ comes from.

8.29. A numerical gradient check returns an analytic gradient of $-0.0692511$ and a central-difference estimate of $-0.0692499$. Is the implementation correct? Justify your answer in terms of the truncation and round-off error of a central difference with step $10^{-6}$.

8.30. In Section 8.10.2 the training AUC rose from 0.632 to 0.854 while the optimism rose from 0.002 to 0.137. Using Part 1, Eq. (8.20), state which term of the bias-variance decomposition is increasing and what would happen to both curves if the cohort were ten times larger.

8.31. A dense layer on a $256\times256$ image producing 32 outputs is replaced by a $5\times5$ convolution with 32 output channels. Compute both parameter counts and the ratio. State the two properties the convolution gains besides parsimony.

8.32. Using Eq. (18), design the shallowest stack of $3\times3$ convolutions and $2\times2$ stride-2 poolings whose receptive field exceeds 60 mm at 0.7 mm in-plane spacing. Report the layer sequence and the final receptive field.

8.33. A network is trained with cross-entropy alone on slices in which the tumour occupies 0.2% of pixels. Predict its behaviour, compute the pixel accuracy of the degenerate solution, and explain how Eq. (19) prevents it.

8.34. Classify each as a valid or invalid augmentation for abdominal CT, with a one-line physical justification: (a) $\pm10^\circ$ rotation; (b) horizontal flip; (c) additive white Gaussian noise; (d) $\pm40$ HU global intensity shift; (e) resampling from 1 mm to 3 mm slice thickness; (f) random square erasing.

Evaluation and translation

8.35. A model reports AUC 0.62 with bootstrap interval [0.47, 0.77]. State precisely what may and may not be concluded, and compute how much larger the cohort would need to be to halve the interval width.

8.36. Explain what a calibration slope of 0.32 means mechanically, and describe the recalibration procedure that would correct it. Would this change the AUC? Why or why not?

8.37. Verify from Eq. (23) that with 47 calibration points and $\alpha = 0.10$ the required order statistic is the 44th. Explain why the guarantee in Eq. (24) is an inequality and not an equality.

8.38. Nine subgroups were evaluated and most intervals contained 0.5. Explain why reporting only the subgroup with the lowest point estimate as evidence of bias would be a multiplicity error, and propose a defensible analysis.

8.39. For each shift type in Section 8.12.1, name a realistic imaging cause and state whether recalibrating the intercept alone would suffice.

8.40. A model predicting malignancy from acquisition metadata alone achieves AUC 0.375, while the same metadata predicts operative extent at 0.567. Interpret both numbers, and explain why the second is a warning even though the first is reassuring.

8.41. Using Eq. (25), explain why the Laplace noise scale is $\Delta f/\varepsilon$. For a count query on this cohort at $\varepsilon = 0.1$, compute the standard deviation and compare it with the 18 benign cases.

8.42. A model with standalone AUC 0.90 reduces reader accuracy when deployed. Give two mechanisms from Section 8.12.8 that could produce this, and state the study design needed to detect it.

8.18 Selected Solutions

Solution 8.22

The partition has discovered the surgical decision. Malignancy is essentially constant across the three clusters (0.889, 0.833, 0.889 against a cohort prevalence of 0.914), so the cluster label carries no outcome information. What is separated almost perfectly is operative extent (radical fraction 1.000, 0.067, 0.022) together with tumour size (median 6.3, 3.0, 2.7 cm), which is unsurprising because size is the principal indication for radical nephrectomy.

Naming these “phenotypes” would be indefensible for three reasons. First, they are not biological: they reproduce a clinical decision rule. Second, any subsequent association with outcome would be confounded by indication, exactly as in Part 1, Section 8.8.4, where the apparent survival disadvantage of radical nephrectomy vanished on adjustment for size. Third, naming a cluster using the outcome after defining it without the outcome is the circularity described in Section 8.9.1. The correct description is “clusters corresponding to operative extent and tumour size”, which is a finding about care pathways, not about disease subtypes.

Solution 8.24

From Eq. (4), write $a = \lvert P\cap G\rvert$ and $u = \lvert P\cup G\rvert$, so $J = a/u$ and $\lvert P\rvert + \lvert G\rvert = u + a$ by inclusion-exclusion. Then

\[\mathrm{Dice} = \frac{2a}{u + a} = \frac{2a/u}{1 + a/u} = \frac{2J}{1+J}.\]

Since $d(J) = 2J/(1+J)$ has $d(J) - J = J(1-J)/(1+J) > 0$ for $0 < J < 1$, Dice strictly exceeds Jaccard for any imperfect but non-empty overlap, with equality only at $J = 0$ and $J = 1$. This is why Dice values in the literature look more favourable than Jaccard values for the identical segmentation, and why the two must never be compared across papers without checking which was used.

Solution 8.25

From the table: “dilate 2 px” has Dice 0.908, area change $+20.2\%$, and $\mathrm{HD}_{95} = 1.84$ mm. “Shift 4 px” has Dice 0.876, area change exactly 0%, and $\mathrm{HD}_{95} = 3.68$ mm.

The Dice values are similar; the failures are opposite. The dilation preserves position and inflates volume, so it would be the more serious error for volumetric response assessment, where a 20% volume error is comparable to a RECIST-equivalent progression threshold and would produce spurious progression calls. The translation preserves volume exactly and displaces the entire boundary, so it would be the more serious error for radiotherapy planning or surgical navigation, where a 3.7 mm systematic boundary displacement moves the dose gradient or the resection plane while every volumetric quality check passes. Reporting Dice alone cannot distinguish them, which is the argument for a mandatory overlap-plus-boundary-plus-clinical metric triple.

Solution 8.26

For any $f\in\mathbb{R}^N$,

\[f^{\mathsf T}Lf = f^{\mathsf T}Df - f^{\mathsf T}Wf = \sum_i d_i f_i^2 - \sum_{i,j}w_{ij}f_if_j = \tfrac12\sum_{i,j}w_{ij}(f_i^2 - 2f_if_j + f_j^2) = \tfrac12\sum_{i,j}w_{ij}(f_i - f_j)^2,\]

using $d_i = \sum_j w_{ij}$ and the symmetry of $W$. Since $w_{ij}\ge 0$, the quadratic form is non-negative, so $L$ is positive semidefinite and all eigenvalues satisfy $\lambda\ge 0$.

Setting $f = \mathbf 1$ gives $f^{\mathsf T}Lf = 0$ and indeed $L\mathbf 1 = D\mathbf 1 - W\mathbf 1 = 0$ because each row of $W$ sums to its degree. Hence $\lambda_1 = 0$ always, with the constant eigenvector.

A second zero eigenvalue means there exists a non-constant $f$ with $\sum_{i,j}w_{ij}(f_i-f_j)^2 = 0$, which requires $f_i = f_j$ whenever $w_{ij} > 0$. That is possible only if the graph splits into two components with no edges between them, so the image contains two regions with no connection at all under the chosen $\sigma_s$ and $\sigma_I$. In practice this signals that the similarity kernel is too narrow, and the segmentation is being determined by the kernel bandwidth rather than by the anatomy.

Solution 8.29

The implementation is almost certainly correct. The relative discrepancy is

\[\frac{\lvert -0.0692511 + 0.0692499\rvert}{\lvert -0.0692511\rvert} = \frac{1.2\times10^{-6}}{6.9\times10^{-2}} \approx 1.7\times10^{-5}.\]

A central difference has truncation error $O(h^2) = 10^{-12}$ and round-off error $O(\epsilon_{\mathrm{mach}}/h) = 2.2\times10^{-16}/10^{-6} \approx 2.2\times10^{-10}$ in the loss, which propagates to roughly $10^{-10}$ in the gradient. The observed discrepancy of $1.2\times10^{-6}$ is larger than that, so it most likely reflects the reported precision of the two printed numbers (seven significant figures) rather than a genuine disagreement. The working rule is that relative agreement better than about $10^{-5}$ passes, agreement worse than $10^{-3}$ indicates a genuine bug, and the region between deserves a check at a second step size $h$: a true bug does not improve when $h$ changes, whereas finite-difference error does.

Solution 8.31

Dense layer: $256\times256 = 65{,}536$ inputs to 32 outputs gives $65{,}536\times32 + 32 = 2{,}097{,}184$ parameters. Convolution: $5\times5\times1\times32 + 32 = 832$ parameters. The ratio is $2{,}097{,}184/832 \approx 2{,}520$, so the convolution uses about 2,500 times fewer parameters.

Besides parsimony, the convolution gains translation equivariance (a feature learned at one location is detected at every location, so a lesion detector does not have to be relearned for each position) and locality with a controllable receptive field (each output depends on a bounded neighbourhood, and Eq. (18) lets the designer grow that neighbourhood deliberately to match the physical size of the target). The dense layer has neither: it must learn a separate weight for every pixel and would need to see lesions in every position during training.

Solution 8.33

With 0.2% tumour pixels, the cross-entropy minimizer under a weak model is to predict the background class almost everywhere. The degenerate solution achieves pixel accuracy $\approx 0.9907$ and tumour Dice exactly 0: it is nearly optimal on the reported accuracy and useless for the clinical task.

Equation (19) prevents this because the soft-Dice term is a ratio in which the tumour class contributes to both numerator and denominator. Predicting no tumour makes the numerator zero and the Dice term equal to 1, so the loss is maximally penalized regardless of how well the background is classified. The Dice term is therefore invariant to the size of the background class, which is exactly the property cross-entropy lacks. In practice both terms are kept: cross-entropy supplies well-behaved per-voxel gradients early in training, and the Dice term supplies the class-balance pressure.

Solution 8.34

Valid. Patient positioning genuinely varies by several degrees between examinations.
Invalid. The abdomen has laterality: a horizontal flip converts a left renal mass into a right one, changing anatomy (the liver and spleen swap sides) and potentially the clinical meaning.
Invalid as a robustness test, marginal as an augmentation. CT noise after reconstruction is spatially correlated, not white; white noise does not simulate a dose reduction (Part 1, Section 8.2.10).
Valid within limits. Contrast timing and reconstruction kernel do shift apparent attenuation, but the shift must stay within a physically plausible range, and any biomarker defined by an absolute HU threshold must be recomputed, as Section 8.10.7 demonstrated with a $+230\%$ change.
Valid and important. Slice thickness genuinely spans 0.5 to 5 mm in this cohort, so simulating it is simulating the real acquisition distribution.
Context-dependent. Random erasing improves robustness to occlusion but can delete the lesion itself, silently corrupting the label; it is acceptable only if the erasure is constrained to avoid the region of interest.

Solution 8.35

May be concluded: the point estimate of discrimination is 0.62, and the data are consistent with any true AUC between 0.47 and 0.77. May not be concluded: that the model discriminates better than chance, since 0.5 lies inside the interval; nor that it is useless, since 0.77 also lies inside.

Interval width scales approximately as $n_{\mathrm{eff}}^{-1/2}$, where the binding quantity is the number of minority-class events. Halving the width therefore requires roughly a fourfold increase, from 18 benign cases to about 72, which at this cohort’s prevalence of 0.914 means about 840 patients. This calculation, not a promise of better modelling, is the honest answer to “how do we improve this?”, and it is corroborated by the learning curve of Section 8.13.4, which is still rising at $n = 210$.

Solution 8.36

A calibration slope of 0.32 means that the predicted log-odds vary about three times faster than the observed log-odds: the model’s predictions are too extreme in both directions, pushing high-risk patients too close to 1 and low-risk patients too close to 0. This is the standard signature of overfitting and is why shrinkage exists.

The correction is logistic recalibration: fit Eq. (22) on an independent set and replace $\widehat p$ with $\sigma(\widehat\alpha + \widehat\beta\,\mathrm{logit}(\widehat p))$. Because $\widehat\beta > 0$, this transformation is strictly monotone, so it preserves the ranking of every patient and therefore leaves the AUC exactly unchanged. That is the essential point: recalibration fixes the numbers without improving the ordering, and a model can be perfectly calibrated and useless, or perfectly discriminating and dangerous at a fixed threshold. The two properties are independent and must both be reported.

Solution 8.41

The Laplace mechanism adds noise $\mathrm{Lap}(\Delta f/\varepsilon)$, whose density is $\propto \exp(-\lvert z\rvert\varepsilon/\Delta f)$. For adjacent datasets the query values differ by at most $\Delta f$, so the ratio of the two output densities at any point is at most $\exp(\varepsilon \Delta f/\Delta f) = e^{\varepsilon}$, which is exactly the requirement of Eq. (25). The scale must therefore grow with sensitivity and shrink with the privacy budget.

For a count query $\Delta f = 1$, and the Laplace distribution with scale $b$ has standard deviation $b\sqrt2$. At $\varepsilon = 0.1$, $b = 10$ and the standard deviation is $10\sqrt2 = 14.1$. The cohort contains only 18 benign cases, so the privacy noise on a single count is nearly as large as the entire minority class. At this cohort size, strong formal privacy and useful statistics are close to mutually exclusive, which is the honest and uncomfortable conclusion: differential privacy is a technology for large datasets, and small-cohort medical research must rely on governance, access control, and data-use agreements instead.

Review the first part of BPAD Chapter 8.

References for Part 2

Ronneberger, O., Fischer, P., and Brox, T. (2015). U-Net: convolutional networks for biomedical image segmentation. MICCAI, 234–241.
Isensee, F. et al. (2021). nnU-Net: a self-configuring method for deep learning-based biomedical image segmentation. Nature Methods, 18(2), 203–211.
Shi, J. and Malik, J. (2000). Normalized cuts and image segmentation. IEEE TPAMI, 22(8), 888–905.
von Luxburg, U. (2007). A tutorial on spectral clustering. Statistics and Computing, 17(4), 395–416.
Tibshirani, R., Walther, G., and Hastie, T. (2001). Estimating the number of clusters via the gap statistic. JRSS-B, 63(2), 411–423.
Rousseeuw, P. J. (1987). Silhouettes: a graphical aid to the interpretation and validation of cluster analysis. J. Comput. Appl. Math., 20, 53–65.
Taha, A. A. and Hanbury, A. (2015). Metrics for evaluating 3D medical image segmentation. BMC Medical Imaging, 15, 29.
Maier-Hein, L. et al. (2024). Metrics reloaded: recommendations for image analysis validation. Nature Methods, 21(2), 195–212.
LeCun, Y., Bengio, Y., and Hinton, G. (2015). Deep learning. Nature, 521, 436–444.
Goodfellow, I., Bengio, Y., and Courville, A. (2016). Deep Learning. MIT Press.
Dosovitskiy, A. et al. (2021). An image is worth 16x16 words: transformers for image recognition at scale. ICLR.
Ilse, M., Tomczak, J., and Welling, M. (2018). Attention-based deep multiple instance learning. ICML.
Geirhos, R. et al. (2020). Shortcut learning in deep neural networks. Nature Machine Intelligence, 2, 665–673.
DeGrave, A. J., Janizek, J. D., and Lee, S.-I. (2021). AI for radiographic COVID-19 detection selects shortcuts over signal. Nature Machine Intelligence, 3, 610–619.
Angelopoulos, A. N. and Bates, S. (2023). Conformal prediction: a gentle introduction. Foundations and Trends in Machine Learning, 16(4), 494–591.
Van Calster, B. et al. (2019). Calibration: the Achilles heel of predictive analytics. BMC Medicine, 17, 230.
Vickers, A. J., van Calster, B., and Steyerberg, E. W. (2016). Net benefit approaches to the evaluation of prediction models. BMJ, 352, i6.
Quiñonero-Candela, J. et al., eds. (2009). Dataset Shift in Machine Learning. MIT Press.
Johnson, W. E., Li, C., and Rabinovic, A. (2007). Adjusting batch effects in microarray expression data using empirical Bayes methods. Biostatistics, 8(1), 118–127.
Obermeyer, Z. et al. (2019). Dissecting racial bias in an algorithm used to manage the health of populations. Science, 366, 447–453.
Chouldechova, A. (2017). Fair prediction with disparate impact. Big Data, 5(2), 153–163.
Dwork, C. and Roth, A. (2014). The algorithmic foundations of differential privacy. Foundations and Trends in Theoretical Computer Science, 9(3-4), 211–407.
Mitchell, M. et al. (2019). Model cards for model reporting. FAT*, 220–229.
Collins, G. S. et al. (2024). TRIPOD+AI statement. BMJ, 385, e078378.
Mongan, J., Moy, L., and Kahn, C. E. (2020). CLAIM: a checklist for artificial intelligence in medical imaging. Radiology: Artificial Intelligence, 2(2), e200029.
Vasey, B. et al. (2022). DECIDE-AI: reporting guideline for early-stage clinical evaluation of decision support systems driven by artificial intelligence. Nature Medicine, 28, 924–933.
Heller, N. et al. (2021). The KiTS19 challenge data and results. Medical Image Analysis, 67, 101821.
Dinov, I. D. (2021). Data Science: Time Complexity, Inferential Uncertainty, and Spacekime Analytics. De Gruyter.
DSPA2: Data Science and Predictive Analytics and the SOCR Probability and Statistics EBook. See also BPAD Chapters 1, 2, 7, and Chapter 8 Part 1.

Source	Nature	Reducible by more data?	How it is quantified here
Aleatoric (irreducible noise)	photon statistics, biological variability, label error	no	residual variance; the \(\sigma^2\) term of Part 1, Eq. (8.20)
Epistemic (parameter uncertainty)	finite sample	yes	bootstrap intervals (Sec. 8.11.1)
Structural (model misspecification)	wrong functional form	partly	calibration slope (Sec. 8.11.2), residual diagnostics
Measurement (input uncertainty)	segmentation, spacing, calibration	partly	one-voxel perturbation, Part 1, Sec. 8.5.6
Distributional (deployment shift)	different scanner, site, population	no, without new data	external evaluation (Sec. 8.12.3)

Shift	What changes	Typical imaging cause	Response
Covariate shift	\(P(X)\) changes, \(P(Y\mid X)\) fixed	new scanner, protocol, or reconstruction kernel	reweighting, harmonization, augmentation
Prior-probability shift	\(P(Y)\) changes, \(P(X\mid Y)\) fixed	screening versus surgical population	recalibrate the intercept
Concept shift	\(P(Y\mid X)\) changes	new reference standard, new treatment era, changed guidelines	retrain; reweighting cannot help

Equation	Statement	Where
Pipeline composition	\(\widehat y = f_\theta[g(\mathcal{P}(\mathcal{R}\{\mathcal{A}(S;A)+N\}), M; \psi), x_c]\)	Part 1, 8.1.1
Regularized ERM	\(\widehat\theta = \arg\min_\theta \frac{1}{n}\sum_i \ell\{y_i, f_\theta(x_i)\} + \lambda\Omega(\theta)\)	Part 1, 8.2.2
Condition-number squaring	\(\kappa(X^{\mathsf T}X) = \kappa(X)^2\)	Part 1, 8.2.3
Bias-variance	\(\mathbb{E}[(Y-\widehat f)^2] = \sigma^2 + \mathrm{Bias}^2 + \mathrm{Var}\)	Part 1, 8.2.7
Predictive value	\(\mathrm{PPV} = Se\,\pi/\{Se\,\pi + (1-Sp)(1-\pi)\}\)	Part 1, 8.2.8
Region volume	\(V_\ell = N_\ell\, s_x s_y s_z\)	Part 1, 8.5.1
Sphericity	\(\Psi = \pi^{1/3}(6V)^{2/3}/S\)	Part 1, 8.5.4
Survival and hazard	\(S(t) = \exp\{-\int_0^t h(u)du\}\)	Part 1, 8.8.1
Cumulative incidence	\(F_k(t) = \int_0^t S(u^-)\lambda_k(u)\,du\)	Part 1, 8.8.6
\(k\)-means objective	Eq. (1)	8.9.1
Silhouette width	Eq. (2)	8.9.3
Dice and Jaccard	Eq. (4), with \(\mathrm{Dice} = 2J/(1+J)\)	8.9.5
Hausdorff and ASSD	Eqs. (5), (6)	8.9.5
Graph Laplacian	\(L = D - W\), \(f^{\mathsf T}Lf = \tfrac12\sum w_{ij}(f_i-f_j)^2\)	8.9.7
Normalized cut	Eq. (11)	8.9.7
Backpropagation	Eqs. (15), (16)	8.10.2
Convolution	Eq. (17)	8.10.3
Receptive field	\(r_\ell = r_{\ell-1} + (k_\ell-1)j_{\ell-1}\)	8.10.4
Composite segmentation loss	Eq. (19)	8.10.6
Attention	\(\mathrm{softmax}(QK^{\mathsf T}/\sqrt{d_k})V\)	8.10.9
Calibration model	\(\mathrm{logit}\{P(Y=1)\} = \alpha + \beta\,\mathrm{logit}(\widehat p)\)	8.11.2
Conformal quantile	Eqs. (23), (24)	8.11.4
Differential privacy	Eq. (25)	8.12.7

Biomedical Physics with Applications to Disease (BPAD)

Chapter 8: Data Modeling, Artificial Intelligence, and Biomedical Applications (Part 2)

Magdalena I. Ivanova, Ivo D. Dinov (BPAD Team + SOCR)

August 05, 2026