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

A predictive model does not begin at the first line of machine-learning code. It begins at the photon, the spin, or the pressure wave, and it inherits every assumption made between there and the number it is handed. Medical-image AI is a measurement-and-inference system, not a classifier bolted onto an image.

How to Use This Chapter

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

Chapters 2 through 6 explained how physical interactions become medical images. Chapter 7 explained how images become measurements. Table 1 maps every earlier result this chapter depends on to the section that uses it. This chapter takes the final scientific step: deciding what those measurements imply about disease, treatment, physiology, and future outcome.

That step is not separable from biomedical physics. The chain runs

\[\text{patient} \rightarrow \text{physical signal} \rightarrow \text{sampled data} \rightarrow \text{reconstructed image} \rightarrow \text{region or representation} \rightarrow \text{prediction or decision},\]

and a change anywhere along it, in dose, field strength, pulse sequence, reconstruction kernel, voxel spacing, contrast timing, segmentation boundary, or reference standard, changes the numbers a model receives. A model that is excellent on one scanner can be useless on another for reasons that are entirely physical.

Chapter-level clinical puzzle. A patient has a renal mass on contrast-enhanced CT. Five different questions are now on the table, and they are not interchangeable. Segmentation: which voxels are kidney and which are tumour? Characterization: is this mass malignant? Regression: how much kidney function will be lost after surgery? Prognosis: what is the risk of death, and over what horizon? Decision support: does any of this change what the surgeon should do? Each has a different target, a different loss function, a different sampling unit, a different performance measure, and a different evidentiary standard. Confusing them is the most common failure in the imaging-AI literature.

Sections covered in Part 1

Part 1 of Chapter 8 covers Sections 8.1 through 8.8: the conceptual framing, the mathematical foundations, study design, the real kidney case study, image-to-feature extraction, preprocessing and representation, supervised learning, and time-to-event and longitudinal modelling. Sections 8.9 onward (unsupervised learning and segmentation, deep learning, evaluation and clinical utility, translation and deployment, and reproducible practice) follow separately and build on the same cohort.

A two-track reading plan

The core track is Sections 8.1–8.5 plus 8.7: the framing, the mathematics, the design rules, the case study, feature extraction, and classical prediction. The advanced track adds Sections 8.6 and 8.8: copulas and rank dependence, acquisition shift, survival analysis, competing risks, and longitudinal growth modelling. Every advanced topic stays attached to the same kidney problem, so the chapter deepens rather than becoming a catalogue of algorithms.

Each section closes with a Section summary suitable for a lecture slide and a Checkpoint question. Try it yourself boxes ask you to change code and observe what moves. The most valuable exercises are the ones that make a result get worse in a predictable way, because that is how you learn which assumption was carrying the result.

Expectations, stated plainly. Real cohorts do not produce tidy figures. In this surgical series roughly nine of every ten resected masses are malignant, only 21 deaths are observed in 210 patients, laboratory values are truncated at the top of their reportable range, slice thickness varies tenfold, and the outcome most amenable to modelling is available in fewer than half the cases. Discrimination statistics are correspondingly modest, and several models in this chapter barely beat their own base rate. Resist the urge to tune until the numbers look publishable. Diagnosing why a real dataset resists a method is the actual skill being taught.

Learning Objectives

#	Objective	Section
1	Express the imaging-AI pipeline as a composition of physical, computational, and statistical operators, and decompose prediction error across its stages.	8.1
2	Identify patient, examination, lesion, slice, and voxel levels in a hierarchical dataset and choose the correct splitting unit.	8.1, 8.3
3	Write a complete intended-use statement and distinguish prediction, estimation, and causal targets.	8.1, 8.7
4	Formulate supervised learning as regularized empirical risk minimization and derive squared-error and cross-entropy losses from their likelihoods.	8.2
5	Explain why $\kappa(X^{\mathsf T}X) = \kappa(X)^2$ makes normal-equation inversion unsafe for radiomic features, and use QR or SVD instead.	8.2
6	Decompose expected prediction error into bias, variance, and irreducible noise, and demonstrate the trade-off computationally.	8.2
7	Convert sensitivity and specificity into a post-test probability and explain why AUC is prevalence-invariant while predictive value is not.	8.2, 8.7
8	Propagate measurement uncertainty through a feature pipeline by the delta method and by Monte Carlo perturbation.	8.2, 8.5
9	Match a likelihood to an imaging noise model (Gaussian, Poisson, Rician) and explain why additive Gaussian noise is not a universal robustness test.	8.2
10	Define an estimand, target population, and unit of analysis, and design a reference standard whose error is documented.	8.3
11	Classify missingness as MCAR, MAR, or MNAR, and recognize structurally informative missingness in real records.	8.3, 8.4
12	Distinguish censoring, truncation, and competing events, and state when a binary endpoint is inadmissible.	8.3, 8.8
13	Construct a leakage-free train/tune/test workflow with nested resampling, and quantify the optimism that leakage produces.	8.3
14	Load real clinical records while preserving units, truncation, semantic sentinels, and provenance.	8.4
15	Read a NIfTI volume from its byte-level header and convert voxel counts into physical volumes.	8.5
16	Compute shape, first-order, and texture features in physical units, and state which are scale-invariant.	8.5
17	Quantify the segmentation-uncertainty floor of a radiomic biomarker by boundary perturbation.	8.5
18	Fit preprocessing on training data only, and diagnose collinearity using correlation, VIF, and the singular-value spectrum.	8.6
19	Connect PCA to the SVD and distinguish dimensionality reduction from feature selection.	8.6
20	Separate marginal behaviour from dependence structure using the empirical copula.	8.6
21	Detect and characterize acquisition-linked domain shift in a real cohort.	8.6
22	Fit and honestly evaluate logistic, penalized, nearest-neighbour, tree-based, and kernel models on a small imbalanced cohort.	8.7
23	Report discrimination, calibration, and clinical utility, and explain when precision-recall is preferable to ROC.	8.7
24	Recognize when a censored outcome requires survival methods, and interpret Kaplan–Meier, log-rank, and Cox output including the proportional-hazards check.	8.8
25	Explain immortal-time bias and construct a landmark analysis.	8.8
26	Show why $1 - \mathrm{KM}$ overestimates risk under competing events and use a cumulative incidence function instead.	8.8
27	Fit exponential, logistic, and Gompertz growth models to longitudinal tumour volumes and explain why they are indistinguishable at short follow-up.	8.8

Notation

Symbol	Meaning
$X \in \mathbb{R}^{n\times p}$	design matrix, $n$ observations by $p$ predictors
$y,\ \widehat y$	outcome and prediction
$\theta,\ \beta$	model parameters; $\beta$ for linear predictors
$\ell(y,\widehat y),\ \Omega(\theta),\ \lambda$	loss, penalty, penalty strength
$\sigma(\eta) = (1+e^{-\eta})^{-1}$	logistic (sigmoid) link
$\kappa(A)$	condition number of $A$
$I \in \mathbb{R}^{n_x\times n_y\times n_z}$	image volume (tensor)
$M$	segmentation mask; $M_v \in \{0,1,2\}$ for background, kidney, tumour
$z = g(I,M;\psi)$	extracted features, with extraction choices $\psi$
$\Delta_{\mathrm{acq}},\Delta_{\mathrm{recon}},\Delta_{\mathrm{seg}},\Delta_{\mathrm{model}},\Delta_{\mathrm{shift}}$	error contributions by pipeline stage
$D,\ T$	development and deployment distributions
$R_D(f),\ R_T(f)$	risk under development and deployment
$\widetilde T = \min(T,C),\ \delta = \mathbb{I}(T\le C)$	observed time and event indicator
$S(t),\ h(t),\ H(t)$	survival, hazard, cumulative hazard
$F_k(t)$	cause-specific cumulative incidence function
$\pi,\ Se,\ Sp$	prevalence, sensitivity, specificity
$V(t),\ K,\ k$	tumour volume, carrying capacity, growth rate

Bridges from Earlier Chapters

Table 1: Where this chapter draws on the rest of BPAD.
Earlier result	Chapter	Role here	Section
Vectors, matrices, tensors	1	design matrices and voxel arrays	8.2, 8.5
Eigenvalues, SVD, condition number	1	PCA, low-rank structure, numerical stability of least squares	8.2, 8.6
Probability axioms and Bayes’ theorem	1	likelihoods, posterior risk, prevalence dependence	8.2, 8.7
Derivatives, gradients, Hessians, optimization	1	gradient descent and backpropagation	8.2
Convolution and the point spread function	1, 7	learned spatial filters; instrument blur as a feature confound	8.5
Fourier analysis, sampling, aliasing	1, 7	why voxel spacing changes texture features	8.5, 8.6
Noise models: Gaussian, Poisson, Rician	1, 2, 4, 5, 6	matching a likelihood to a modality	8.2
SNR, CNR, and the $\sqrt{N}$ law	1, 5, 6	detectability limits on any downstream model	8.2
Error propagation, delta method, Monte Carlo	1, 2	feature robustness and predictive uncertainty	8.2, 8.5
Maximum likelihood and confidence intervals	1	estimation, inference, calibration	8.2, 8.7
ODEs: exponential, logistic, Gompertz growth	1	longitudinal tumour-volume modelling	8.8
Compartmental models and AUC	1, 6	exposure summaries as predictors	8.8
Registration, segmentation, quantification	7	the upstream source of every feature used here	8.5
Optical and thermal quantification, calibration bias	2	pigmentation bias as a physically grounded fairness problem	8.7

Computational Environment, Data Provenance, and Reproducibility

Requirements. A recent R installation. The analytic spine uses base R plus jsonlite to read the released metadata. Optional sections use survival, glmnet, randomForest, e1071, pROC, and cmprsk; each is guarded by has_pkg(), so a missing package degrades a demonstration gracefully rather than stopping the knit.

Data. The clinical cohort is the 210 publicly released KiTS19 training cases, read from kits.json (about 0.5 MB, no authentication). Imaging features are computed from the corresponding expert segmentation.nii.gz label volumes (about 0.8 MB each). Everything is cached under bpad.cache_dir.

Two cohorts that must not be merged. KiTS19 and KiTS21 are related but distinct challenge releases with different case counts, annotation protocols, and metadata fields. This chapter uses KiTS19 throughout, and the released volumes are contrast-enhanced abdominal CT, not MRI. Descriptions of one release must not be silently attached to data from the other, and a chapter that mixes them will report morphometry that cannot be reproduced.

Option	Default	Controls
`bpad.cache_dir`	session temp directory	where downloads are cached; set to a persistent folder to download once
`bpad.allow_network`	`TRUE`	set `FALSE` to run purely from an existing cache or the surrogate cohort
`bpad.n_mask_cases`	`6`	how many real segmentation volumes are streamed into imaging features
`bpad.run_heavy`	`FALSE`	enable longer-running demonstrations

## Cache persistently and process more segmentation volumes.
options(
  bpad.cache_dir    = "~/bpad_cache",
  bpad.allow_network = TRUE,
  bpad.n_mask_cases = 30,
  bpad.run_heavy    = TRUE
)
rmarkdown::render("BPAD_Chap08_DataModeling_AI_ML.Rmd",
                  output_format = "bookdown::html_document2")

On interactive demonstrations. Earlier drafts of this chapter embedded Shiny applications directly in executable chunks. That is not compatible with a static knit: a shinyApp() object at the top level of a chunk cannot render in a Word document and will stall or error in a plain HTML build. In this edition every interactive demonstration appears as a non-evaluated code listing that a reader can copy into app.R and run, and each is paired with a static figure that carries the same teaching point in both HTML and DOCX. Readers who want the apps embedded in a self-contained page should use the shinylive route, which packages a WebAssembly build of R into the HTML; note that this substantially increases page size and first-load time.

8.1 From Physical Measurement to a Learning Problem

8.1.1 A compositional model of the imaging-AI chain

Let $S$ denote an underlying biological state, $A$ the acquisition settings, $N$ measurement noise, $\mathcal{R}$ reconstruction, $\mathcal{P}$ preprocessing, $g$ feature extraction or representation learning, and $f_\theta$ a prediction rule. The observed image is

\[\begin{equation} I = \mathcal{R}\{\mathcal{A}(S;A) + N\}, \tag{1} \end{equation}\]

features are extracted from the image and a mask $M$ by

\[\begin{equation} z = g(I, M; \psi), \tag{2} \end{equation}\]

where $\psi$ collects every extraction choice (resampling grid, intensity discretization, bin count, mask definition), and the prediction combines imaging with clinical information $x_c$:

\[\begin{equation} \widehat y = f_\theta(z, x_c). \tag{3}. \end{equation}\]

Composing these gives the whole pipeline in one line,

\[\begin{equation} \widehat y = f_\theta\Big[\,g\big(\mathcal{P}(\mathcal{R}\{\mathcal{A}(S;A)+N\}),\,M;\,\psi\big),\ x_c\Big]. \tag{4} \end{equation}\]

Equation (4) is the organizing statement of this chapter. The fitted model depends not only on $\theta$ but on $A$, $\mathcal{R}$, $\mathcal{P}$, $M$, and $\psi$, none of which appear in the usual machine-learning formulation and all of which are physics or engineering choices made before any learning occurs.

8.1.2 Error decomposition across the pipeline

A useful accounting device separates the sources of discrepancy between prediction and truth:

\[\begin{equation} \widehat y - y = \underbrace{\Delta_{\mathrm{acq}}}_{\text{acquisition}} + \underbrace{\Delta_{\mathrm{recon}}}_{\text{reconstruction}} + \underbrace{\Delta_{\mathrm{seg}}}_{\text{segmentation}} + \underbrace{\Delta_{\mathrm{model}}}_{\text{statistical model}} + \underbrace{\Delta_{\mathrm{shift}}}_{\text{deployment shift}} + \underbrace{\epsilon}_{\text{residual noise}}. \tag{5} \end{equation}\]

The terms are neither additive nor independent in general, so Eq. (5) is a reasoning device, not an exact identity. Its purpose is to prevent the default habit of attributing all uncertainty to the final algorithm. When an imaging-AI model transfers poorly, the first hypothesis should usually be $\Delta_{\mathrm{acq}}$ or $\Delta_{\mathrm{shift}}$, not $\Delta_{\mathrm{model}}$. Section 8.5.6 measures $\Delta_{\mathrm{seg}}$ directly for a real biomarker, and Section 8.6.6 measures $\Delta_{\mathrm{shift}}$ directly in the real cohort.

Key idea. Every quantity a model consumes is a measurement with a provenance and an uncertainty. Chapters 2 through 7 established what those uncertainties are; this chapter is about not throwing that knowledge away the moment the data enters a spreadsheet.

8.1.3 Units of observation and hierarchical dependence

Index patients by $i$, examinations by $j$, lesions by $k$, and voxels by $v$, so an observation is $I_{ijkv}$. Voxels within a lesion are far more alike than voxels from different patients; slices within one examination share anatomy, contrast phase, and reconstruction kernel. Treating these as independent inflates the nominal sample size without adding information.

The correct notion is the effective sample size. For $m$ correlated observations per patient with intraclass correlation $\rho$,

\[\begin{equation} n_{\mathrm{eff}} \approx \frac{n\,m}{1 + (m-1)\rho}, \tag{6} \end{equation}\]

so at $\rho = 0.8$, one hundred slices from each of 200 patients carry an effective sample size close to $200 \times 100 / 80 = 250$, not 20,000. Every confidence interval, every $p$-value, and every claim of statistical power must be computed against $n_{\mathrm{eff}}$.

Leakage warning. If slices from one patient are divided between training and testing sets, the model is tested on anatomy and acquisition characteristics it has already seen. Split at the highest relevant independent unit: usually the patient, and sometimes the site or the time period. Section 8.3.6 shows what this costs when it is ignored.

8.1.4 Prediction target, intended use, and action

A target must be attached to a time point. A preoperative model cannot use postoperative variables. A risk model must state its horizon. A segmentation model must state its label ontology. A decision-support model must state the action it is meant to inform.

A complete intended-use statement names seven things: population, decision time, permissible inputs, output, horizon, user, and action. A convenient template is

In [target population], at [decision time], use [available predictors] to estimate [target] over [time horizon] in order to support [clinical or scientific action].

For this chapter’s running problem:

In adults undergoing evaluation of a renal mass, at the preoperative visit, use contrast-enhanced CT and routinely available clinical variables to estimate the probability of malignant pathology, in order to support surgical planning and patient counselling.

That single sentence fixes the estimand and immediately disqualifies pathologic size, ISUP grade, T stage, ischaemia time, and operative blood loss as predictors, because none of them exists at the decision time. Section 8.4 shows that the real cohort contains all five, which is exactly why the statement must be written before the data are opened.

8.1.5 Learning paradigms

In supervised learning, paired inputs and labels $\{(x_i,y_i)\}_{i=1}^{n}$ estimate a mapping from $X$ to $Y$; classification, regression, survival prediction, detection, and segmentation are all supervised. In unsupervised learning the model sees $X$ alone: PCA, clustering, density estimation, and most anomaly detection. In semi-supervised learning a small labelled set is combined with a much larger unlabelled one. In self-supervised learning the targets are constructed from the data itself, for example by reconstructing masked patches or contrasting augmented views, before task-specific fine-tuning. Medical imaging is unusually well suited to the last of these, because unlabelled images are abundant while expert annotations are scarce and expensive.

8.1.6 Mechanistic, statistical, predictive, and causal models

These four words are routinely used interchangeably and mean four different things.

Model type	Goal	Success criterion	Example in this chapter
Mechanistic	encode physical or biological law	consistency with known dynamics	logistic tumour growth (Sec. 8.8.7)
Statistical	describe association and quantify uncertainty	valid coverage and calibration	eGFR decline regression (Sec. 8.7.6)
Predictive	forecast future observations	out-of-sample accuracy and utility	malignancy classifier (Sec. 8.7.2)
Causal	estimate the effect of an intervention	unbiasedness under explicit assumptions	effect of radical versus partial nephrectomy (Sec. 8.7.10)

A feature may predict well precisely because it records a downstream consequence of treatment. That makes it a good predictor and an invalid baseline confounder. Conversely a causal estimate can be scientifically decisive while adding little to individual-level prediction accuracy.

Prediction is not causation. Predicting who receives radical nephrectomy is not the same as estimating whether assigning radical nephrectomy changes outcome. Section 8.8.4 shows both computations on the real cohort and they point in opposite directions: radical nephrectomy is associated with worse survival on a log-rank test and with no detectable effect once tumour size is adjusted for, because larger tumours are the reason the radical operation was chosen.

Section 8.1 summary.

The pipeline is a composition, Eq. (4); the fitted model depends on acquisition, reconstruction, mask, and extraction choices, not only on $\theta$.
Prediction error decomposes across pipeline stages; failures to transfer are usually acquisition or shift terms, not model terms.
Effective sample size, not row count, governs inference; split at the highest independent unit.
An intended-use statement fixes the estimand and disqualifies predictors that do not exist at decision time.
Mechanistic, statistical, predictive, and causal models have different success criteria and are not substitutes.

Checkpoint 8.1. A study reports a lesion-classification AUC of 0.94 from 40,000 CT slices drawn from 180 patients, split at the slice level. Using Eqs. (6) and (5), give two independent reasons the reported number is likely to be optimistic, and state which one you would investigate first.

8.2 Mathematical Foundations of Statistical Learning

8.2.1 Design matrices and parameterized prediction rules

For $n$ patients and $p$ predictors,

\[\begin{equation} X = \begin{bmatrix} 1 & x_{11} & \cdots & x_{1p}\\ 1 & x_{21} & \cdots & x_{2p}\\ \vdots & \vdots & \ddots & \vdots\\ 1 & x_{n1} & \cdots & x_{np} \end{bmatrix} \in \mathbb{R}^{n\times(p+1)}, \qquad y = \begin{bmatrix} y_1 & \cdots & y_n\end{bmatrix}^{\mathsf T}. \tag{7} \end{equation}\]

A linear predictor is $\eta = X\beta$. Regression uses $\widehat y = X\widehat\beta$ directly; logistic regression maps the linear predictor to a probability through

\[\begin{equation} \sigma(\eta_i) = \frac{1}{1 + e^{-\eta_i}}. \tag{8} \end{equation}\]

An image is itself a tensor, so a neural network is a parameterized sequence of tensor transformations ending in a score, a probability, or a mask. The distinction between “classical” and “deep” models is one of parameterization and capacity, not of a different underlying framework.

8.2.2 Empirical risk minimization

Most fitting procedures solve

\[\begin{equation} \widehat\theta = \arg\min_{\theta} \left[\underbrace{\frac{1}{n}\sum_{i=1}^{n}\ell\{y_i, f_\theta(x_i)\}}_{\text{fidelity}} + \lambda\underbrace{\Omega(\theta)}_{\text{regularization}}\right]. \tag{9} \end{equation}\]

The first term measures misfit through a loss $\ell$; the second penalizes complexity; the tuning parameter $\lambda \ge 0$ sets the trade-off. Almost every method in this chapter is a choice of $\ell$, $\Omega$, and the function class $\{f_\theta\}$.

Task	Loss $\ell$	Typical penalty $\Omega$
Linear regression	squared error $(y-\widehat y)^2$	none, ridge, lasso
Logistic regression	binary cross-entropy	none, ridge, lasso
Support vector machine	hinge $\max(0, 1-y\widehat y)$	squared norm
Cox model	negative partial log-likelihood	ridge, lasso
Neural classifier	cross-entropy	weight decay, dropout, augmentation
Segmentation network	cross-entropy, Dice, or a sum	weight decay, augmentation

Choose the loss to match the noise model. Squared error is the maximum-likelihood loss under additive Gaussian noise; cross-entropy is the maximum-likelihood loss for Bernoulli outcomes; the Poisson deviance is correct for photon counts. Section 8.2.10 makes the modality-specific version of this argument. A loss is not a stylistic preference, it is a distributional assumption.

8.2.3 Least squares, likelihood, and numerical conditioning

For the Gaussian linear model $y = X\beta + \varepsilon$ with $\varepsilon \sim \mathcal{N}(0,\sigma^2 I)$, least squares minimizes $Q(\beta) = \lVert y - X\beta\rVert_2^2$. Differentiating and setting the gradient to zero gives the normal equations

\[\begin{equation} X^{\mathsf T}X\,\widehat\beta = X^{\mathsf T}y, \tag{10} \end{equation}\]

with solution $\widehat\beta = (X^{\mathsf T}X)^{-1}X^{\mathsf T}y$ when $X^{\mathsf T}X$ is invertible. Geometrically, $X\widehat\beta$ is the orthogonal projection of $y$ onto the column space of $X$; probabilistically, it is the maximum-likelihood estimate.

Numerically, however, Eq. (10) is a trap for radiomic data. Forming the cross-product squares the condition number:

\[\begin{equation} \kappa(X^{\mathsf T}X) = \kappa(X)^2. \tag{11} \end{equation}\]

Radiomic features are heavily collinear (volume, surface area, bounding-box diagonal, and equivalent diameter all measure size), so $\kappa(X)$ is already large, and Eq. (11) turns a merely awkward problem into an unstable one. Two factorizations avoid it.

QR decomposition. Write $X = QR$ with $Q^{\mathsf T}Q = I$ and $R$ upper triangular. Then Eq. (10) reduces to $R\widehat\beta = Q^{\mathsf T}y$, solved by back-substitution. Because orthogonal transformations preserve lengths, the condition number is never squared. This is what lm() actually does.

Singular value decomposition. Write $X = U\Sigma V^{\mathsf T}$. Collinear features appear as small singular values $\sigma_j$, and the pseudo-inverse

\[\begin{equation} X^{+} = V\Sigma^{+}U^{\mathsf T}, \qquad \Sigma^{+}_{jj} = \begin{cases}1/\sigma_j, & \sigma_j > \tau\\ 0, & \sigma_j \le \tau\end{cases} \tag{12} \end{equation}\]

truncates the unstable directions explicitly rather than inverting them. Truncated SVD is PCA regression in disguise (Section 8.6.4), and the choice of $\tau$ is a regularization decision that must be made inside resampling.

set.seed(3)
n <- 120
make_design <- function(eps) {
  ## A size-like feature and three near-duplicates of it, plus two independent ones
  size   <- rnorm(n, 4, 1.2)
  cbind(volume   = size^3 + rnorm(n, 0, eps),
        surface  = size^2 + rnorm(n, 0, eps),
        diameter = size    + rnorm(n, 0, eps),
        bbox_diag= 1.7*size + rnorm(n, 0, eps),
        age      = rnorm(n, 60, 11),
        bmi      = rnorm(n, 29, 5))
}
X0 <- make_design(1e-3)
sv <- svd(scale(X0))$d
cat("singular values:", signif(sv, 4), "\n")

## singular values: 21.72 11.91 9.714 2.457 0.1654 0.008743

cat(sprintf("kappa(X) = %.3g   kappa(X'X) = %.3g   ratio = %.3g\n",
            max(sv)/min(sv), (max(sv)/min(sv))^2, (max(sv)/min(sv))))

## kappa(X) = 2.48e+03   kappa(X'X) = 6.17e+06   ratio = 2.48e+03

## Coefficient recovery error for three solvers as collinearity worsens
beta_true <- c(0.02, -0.05, 0.8, 0.1, 0.03, -0.02)
eps_grid  <- 10^seq(0, -6, length.out = 25)
err <- t(sapply(eps_grid, function(e) {
  X <- make_design(e)
  y <- as.vector(X %*% beta_true) + rnorm(n, 0, 0.5)
  b_norm <- tryCatch(as.vector(solve(t(X) %*% X, t(X) %*% y)),
                     error = function(z) rep(NA, ncol(X)))
  b_qr   <- as.vector(qr.solve(X, y))
  s      <- svd(X); keep <- s$d > max(s$d)*1e-8
  b_svd  <- as.vector(s$v[, keep, drop = FALSE] %*%
                      ((t(s$u[, keep, drop = FALSE]) %*% y)/s$d[keep]))
  c(normal = sqrt(mean((b_norm - beta_true)^2)),
    QR     = sqrt(mean((b_qr   - beta_true)^2)),
    SVD    = sqrt(mean((b_svd  - beta_true)^2)))
}))
kap <- sapply(eps_grid, function(e) { d <- svd(scale(make_design(e)))$d; max(d)/min(d) })

p1 <- ggplot(data.frame(k = seq_along(sv), sv = sv), aes(k, sv)) +
  geom_col(fill = bpad_pal[1], width = 0.65) +
  scale_y_log10() +
  labs(title = "singular-value spectrum (log scale)",
       x = "index", y = expression(sigma[j]))

ed <- do.call(rbind, lapply(colnames(err), function(cn)
  data.frame(kappa = kap, err = err[, cn], solver = cn)))
p2 <- ggplot(ed, aes(kappa, err, colour = solver)) +
  geom_line(linewidth = 0.9, na.rm = TRUE) +
  scale_x_log10() + scale_y_log10() +
  scale_colour_manual(values = bpad_pal[1:3]) +
  labs(title = "coefficient RMSE versus conditioning",
       x = expression(kappa(X)~"(log)"), y = "RMSE (log)")
bpad_grid(p1, p2, ncol = 2)

$Numerical conditioning of least squares with collinear radiomic-style features. Left: the singular-value spectrum of a design matrix in which three size descriptors are near-duplicates; the small singular values are the collinear directions. Right: coefficient error versus the collinearity level for three solvers. Normal-equation inversion degrades fastest because it squares the condition number, Eq. \@ref(eq:condition-squared).$

Figure 1: Numerical conditioning of least squares with collinear radiomic-style features. Left: the singular-value spectrum of a design matrix in which three size descriptors are near-duplicates; the small singular values are the collinear directions. Right: coefficient error versus the collinearity level for three solvers. Normal-equation inversion degrades fastest because it squares the condition number, Eq. (11).

8.2.4 Logistic regression and cross-entropy

For a binary outcome $Y_i \in \{0,1\}$ with $p_i = \sigma(x_i^{\mathsf T}\beta)$, the Bernoulli likelihood is $L(\beta) = \prod_i p_i^{y_i}(1-p_i)^{1-y_i}$, and the negative log-likelihood is the binary cross-entropy

\[\begin{equation} -\log L(\beta) = -\sum_{i=1}^{n}\big[y_i\log p_i + (1-y_i)\log(1-p_i)\big]. \tag{13} \end{equation}\]

A coefficient $\beta_j$ is the change in conditional log-odds per unit increase in $x_j$, holding the other modelled variables fixed. It is not a causal effect without design assumptions, and it is not an absolute risk without the intercept and the covariate values.

Logistic growth is not logistic regression. Both involve a sigmoid. Logistic growth (Section 8.8.7) is a dynamical model for a bounded quantity evolving in time, derived from an ODE. Logistic regression is a probability model for a binary outcome at a fixed time. They share a functional form and nothing else.

8.2.5 Regularization

Ridge uses $\Omega(\beta) = \lVert\beta\rVert_2^2$ and shrinks correlated coefficients toward each other, stabilizing collinear designs without removing variables. Lasso uses $\Omega(\beta) = \lVert\beta\rVert_1$ and can set coefficients exactly to zero, performing selection. Elastic net combines both:

\[\begin{equation} \Omega_{\alpha}(\beta) = \alpha\lVert\beta\rVert_1 + \tfrac{1}{2}(1-\alpha)\lVert\beta\rVert_2^2 . \tag{14} \end{equation}\]

Regularization is part of fitting, not a post-processing step. The penalty strength must be selected inside resampling. Choosing $\lambda$ on the complete dataset and then reporting cross-validated performance is leakage, and Section 8.3.6 quantifies how much optimism it buys.

Ridge or lasso for radiomics? With many strongly correlated size descriptors, lasso will pick one essentially at random and drop the rest, so the selected feature set is unstable across resamples even when predictive performance is stable. Ridge keeps all of them with shrunken coefficients, which is worse for interpretation and often better for prediction. If a stable, interpretable subset is required, prefer explicit grouping or clustering of correlated features (Section 8.6.3) over relying on lasso to make the choice.

8.2.6 Optimization, gradient descent, and backpropagation

For differentiable risk $J(\theta)$, gradient descent iterates

\[\begin{equation} \theta^{(m+1)} = \theta^{(m)} - \alpha_m \nabla_\theta J\{\theta^{(m)}\}, \tag{15} \end{equation}\]

with step size $\alpha_m$. For a composition $f = f_L\circ\cdots\circ f_1$ with intermediate activations $h_\ell$, the chain rule gives backpropagation

\[\begin{equation} \frac{\partial f}{\partial\theta_\ell} = \frac{\partial f_L}{\partial h_{L-1}}\cdots \frac{\partial h_{\ell+1}}{\partial h_\ell}\frac{\partial h_\ell}{\partial\theta_\ell}. \tag{16} \end{equation}\]

Equation (16) is nothing more than the multivariable chain rule of Chapter 1 applied to a deep composition, and its product structure is why very deep networks suffer vanishing or exploding gradients: a product of $L$ factors each slightly below or above one behaves exponentially in $L$.

8.2.7 Bias, variance, and irreducible noise

For squared-error prediction at a point $x$,

\[\begin{equation} \mathbb{E}\big[(Y - \widehat f(x))^2\big] = \underbrace{\sigma^2}_{\text{irreducible}} + \underbrace{\big(\mathbb{E}[\widehat f(x)] - f(x)\big)^2}_{\text{bias}^2} + \underbrace{\mathbb{E}\big[(\widehat f(x) - \mathbb{E}[\widehat f(x)])^2\big]}_{\text{variance}} . \tag{17} \end{equation}\]

Bias is systematic error from approximating a complex truth with a simpler model (underfitting). Variance is sensitivity to which particular training sample was drawn (overfitting). The noise term $\sigma^2$ is a floor no model can cross, and in imaging it is set by the physics: photon statistics, thermal noise, and reference-standard error (Section 8.3.2).

The demonstration below both shows the trade-off curve and estimates the three terms separately by repeated sampling, which is the only way to see that bias and variance really are moving in opposite directions.

set.seed(8)
true_f <- function(x) sin(2*pi*x)
sigma_noise <- 0.30
n_tr <- 45
x_te <- seq(0.02, 0.98, length.out = 200); f_te <- true_f(x_te)
degrees <- 1:14

## (a) single training sample: apparent versus generalization error
x_tr <- sort(runif(n_tr)); y_tr <- true_f(x_tr) + rnorm(n_tr, 0, sigma_noise)
tr_mse <- te_mse <- numeric(length(degrees))
for (j in seq_along(degrees)) {
  fit <- lm(y_tr ~ poly(x_tr, degrees[j], raw = TRUE))
  tr_mse[j] <- mean(residuals(fit)^2)
  te_mse[j] <- mean((f_te - predict(fit, data.frame(x_tr = x_te)))^2)
}

## (b) repeated sampling: separate bias^2 from variance
R <- 200
bias2 <- vari <- numeric(length(degrees))
for (j in seq_along(degrees)) {
  P <- matrix(NA_real_, R, length(x_te))
  for (r in seq_len(R)) {
    xs <- sort(runif(n_tr)); ys <- true_f(xs) + rnorm(n_tr, 0, sigma_noise)
    fr <- lm(ys ~ poly(xs, degrees[j], raw = TRUE))
    P[r, ] <- predict(fr, data.frame(xs = x_te))
  }
  mbar <- colMeans(P)
  bias2[j] <- mean((mbar - f_te)^2)
  vari[j]  <- mean(apply(P, 2, var))
}

print(knitr::kable(data.frame(
  degree = degrees, train_MSE = round(tr_mse, 4), test_MSE = round(te_mse, 4),
  bias2 = round(bias2, 4), variance = round(vari, 4),
  bias2_plus_var = round(bias2 + vari, 4))[c(1,2,3,4,6,9,12,14), ],
  row.names = FALSE,
  caption = "Bias-squared and variance estimated from 200 independent training samples. Their sum tracks the test error; the irreducible noise variance is 0.09."))

## 
## 
## Table: (\#tab:bias-variance-demo)Bias-squared and variance estimated from 200 independent training samples. Their sum tracks the test error; the irreducible noise variance is 0.09.
## 
## | degree| train_MSE| test_MSE|  bias2| variance| bias2_plus_var|
## |------:|---------:|--------:|------:|--------:|--------------:|
## |      1|    0.2933|   0.1866| 0.1741|   0.0134|         0.1875|
## |      2|    0.2931|   0.1864| 0.1723|   0.0233|         0.1956|
## |      3|    0.0901|   0.0083| 0.0036|   0.0086|         0.0122|
## |      4|    0.0901|   0.0083| 0.0035|   0.0113|         0.0148|
## |      6|    0.0754|   0.0240| 0.0001|   0.0178|         0.0179|
## |      9|    0.0696|   0.0405| 0.0008|   0.2715|         0.2723|
## |     12|    0.0565|   0.6199| 0.0211|   2.9783|         2.9994|
## |     14|    0.0552|   0.1563| 0.0562|   5.4426|         5.4988|

p1 <- ggplot(rbind(data.frame(d = degrees, e = tr_mse, s = "training (apparent)"),
                   data.frame(d = degrees, e = te_mse, s = "test (generalization)")),
             aes(d, e, colour = s)) +
  geom_line(linewidth = 0.9) + geom_point(size = 1.6) +
  scale_y_log10() + scale_colour_manual(values = bpad_pal[1:2]) +
  labs(title = "apparent error is not generalization error",
       x = "polynomial degree", y = "MSE (log)")

p2 <- ggplot(rbind(data.frame(d = degrees, v = bias2,        s = "bias^2"),
                   data.frame(d = degrees, v = vari,         s = "variance"),
                   data.frame(d = degrees, v = bias2 + vari, s = "sum")),
             aes(d, v, colour = s)) +
  geom_line(linewidth = 0.9) + geom_point(size = 1.4) +
  scale_y_log10() + scale_colour_manual(values = bpad_pal[c(1,2,3)]) +
  labs(title = "the two components move oppositely",
       x = "polynomial degree", y = "contribution to MSE (log)")
bpad_grid(p1, p2, ncol = 2)

$The bias-variance trade-off, Eq. \@ref(eq:bias-variance). Left: training error falls monotonically with polynomial degree while test error is U-shaped. Right: bias-squared and variance estimated by 200 independent training samples, with their sum plus the irreducible noise reproducing the test error. The optimum is where the two curves cross in influence, not where either is smallest.$

Figure 2: The bias-variance trade-off, Eq. (17). Left: training error falls monotonically with polynomial degree while test error is U-shaped. Right: bias-squared and variance estimated by 200 independent training samples, with their sum plus the irreducible noise reproducing the test error. The optimum is where the two curves cross in influence, not where either is smallest.

cat(sprintf("optimal degree by test MSE: %d;  by bias^2 + variance: %d\n",
            degrees[which.min(te_mse)], degrees[which.min(bias2 + vari)]))

## optimal degree by test MSE: 3;  by bias^2 + variance: 3

8.2.8 Bayes’ theorem, prevalence, and clinical interpretation

With prevalence $\pi$, sensitivity $Se$, and specificity $Sp$, the positive and negative predictive values are

\[\begin{equation} \mathrm{PPV} = P(D\mid +) = \frac{Se\,\pi}{Se\,\pi + (1-Sp)(1-\pi)}, \qquad \mathrm{NPV} = P(\bar D\mid -) = \frac{Sp\,(1-\pi)}{Sp\,(1-\pi) + (1-Se)\pi}. \tag{18} \end{equation}\]

Sensitivity and specificity are properties of the test. Predictive values are properties of the test and the population. AUC, being a function of $Se$ and $Sp$ alone, is prevalence-invariant; PPV is not. This is why a paper that reports only AUC can be simultaneously accurate and clinically misleading.

The point is sharpest in this chapter’s own cohort, where the malignancy prevalence in a surgical series is about 0.91, while the prevalence among all incidentally detected small renal masses is closer to 0.20.

ppv <- function(pi, se, sp) se*pi/(se*pi + (1 - sp)*(1 - pi))
npv <- function(pi, se, sp) sp*(1 - pi)/(sp*(1 - pi) + (1 - se)*pi)
pi_g <- seq(0.005, 0.99, length.out = 400)
se <- 0.80; sp <- 0.80

pv <- rbind(data.frame(pi = pi_g, v = ppv(pi_g, se, sp), q = "positive predictive value"),
            data.frame(pi = pi_g, v = npv(pi_g, se, sp), q = "negative predictive value"))
ggplot(pv, aes(pi, v, colour = q)) +
  geom_line(linewidth = 0.95) +
  geom_vline(xintercept = c(0.20, 0.914), linetype = "dashed", colour = "grey45") +
  annotate("text", x = 0.20, y = 0.05, label = "incidental mass\n0.20", size = 2.7,
           hjust = -0.08, colour = "grey30") +
  annotate("text", x = 0.914, y = 0.05, label = "surgical series\n0.91", size = 2.7,
           hjust = 1.08, colour = "grey30") +
  scale_colour_manual(values = bpad_pal[1:2]) +
  labs(title = "Predictive value depends on where the test is deployed",
       x = "disease prevalence", y = "probability")

$Predictive values as functions of prevalence for a test with 80 percent sensitivity and 80 percent specificity, Eq. \@ref(eq:ppv-npv). Vertical lines mark the surgical-series prevalence of this chapter's cohort (0.91) and a plausible incidental-mass prevalence (0.20). The same test is nearly conclusive in one setting and close to useless in the other, while its AUC is identical in both.$

Figure 3: Predictive values as functions of prevalence for a test with 80 percent sensitivity and 80 percent specificity, Eq. (18). Vertical lines mark the surgical-series prevalence of this chapter’s cohort (0.91) and a plausible incidental-mass prevalence (0.20). The same test is nearly conclusive in one setting and close to useless in the other, while its AUC is identical in both.

print(knitr::kable(data.frame(
  prevalence = c(0.05, 0.10, 0.20, 0.50, 0.914),
  PPV = round(ppv(c(0.05,0.10,0.20,0.50,0.914), se, sp), 3),
  NPV = round(npv(c(0.05,0.10,0.20,0.50,0.914), se, sp), 3),
  AUC = "unchanged"),
  caption = "A fixed test at five prevalences. PPV moves by a factor of eight; AUC does not move at all."))

## 
## 
## Table: (\#tab:prevalence-demo)A fixed test at five prevalences. PPV moves by a factor of eight; AUC does not move at all.
## 
## | prevalence|   PPV|   NPV|AUC       |
## |----------:|-----:|-----:|:---------|
## |      0.050| 0.174| 0.987|unchanged |
## |      0.100| 0.308| 0.973|unchanged |
## |      0.200| 0.500| 0.941|unchanged |
## |      0.500| 0.800| 0.800|unchanged |
## |      0.914| 0.977| 0.273|unchanged |

Worked Example 8.1 (why the surgical AUC misleads). A malignancy classifier with $Se = Sp = 0.80$ is developed in this chapter’s surgical cohort, where $\pi = 0.914$. Its PPV there is

\[\mathrm{PPV} = \frac{0.80\times0.914}{0.80\times0.914 + 0.20\times0.086} = \frac{0.731}{0.749} = 0.977,\]

so a positive result is almost conclusive. Deployed for incidentally detected small masses with $\pi = 0.20$,

\[\mathrm{PPV} = \frac{0.80\times0.20}{0.80\times0.20 + 0.20\times0.80} = \frac{0.16}{0.32} = 0.50,\]

a coin flip. The AUC is unchanged. The clinical consequence is that a model validated only in resected specimens cannot be used to decide whether to resect, because the population that would receive the test is precisely the one from which the surgical cohort was selected. This is selection on the outcome, and no amount of algorithmic sophistication repairs it.

8.2.9 Measurement error and uncertainty propagation

Suppose the recorded predictor is $X^{\ast} = X + \delta$ with $\mathbb{E}[\delta] = 0$ and $\mathrm{Var}(\delta) = \sigma_\delta^2$. Two consequences follow. Attenuation: classical independent measurement error in a predictor shrinks its estimated coefficient toward zero, by a factor approximately $\sigma_X^2/(\sigma_X^2 + \sigma_\delta^2)$ in simple regression, so effect sizes are systematically understated. Selection instability: near a decision boundary, small perturbations cause lasso or stepwise procedures to select different members of a collinear group on repeated scans of the same patient.

For a differentiable pipeline $g(U)$ with input covariance $\Sigma_U$, the delta method of Chapter 1 gives

\[\begin{equation} \mathrm{Var}\{g(U)\} \approx \nabla g(\mu_U)^{\mathsf T}\,\Sigma_U\,\nabla g(\mu_U), \tag{19} \end{equation}\]

with the off-diagonal terms of $\Sigma_U$ mandatory rather than optional. When $g$ is non-differentiable, which is the normal case in imaging because masks are discrete and thresholds are hard, Monte Carlo perturbation is both more honest and easier: perturb plausible acquisition or segmentation inputs, recompute features and predictions, and summarize the resulting distribution. Section 8.5.6 does exactly this for a real tumour volume.

What a model confidence interval usually omits. A coefficient interval from glm() conditions on the features as if they were exact. It says nothing about voxel spacing, intensity calibration, contrast timing, or where the annotator drew the boundary. In imaging AI these are frequently the dominant uncertainty terms, and a reported interval that excludes them is not conservative, it is wrong in a known direction.

8.2.10 Imaging noise models, likelihoods, SNR, and CNR

The physics of image formation determines the noise distribution, and the noise distribution determines the correct likelihood. Chapter 1 developed these distributions; Chapters 2, 4, 5, and 6 attached them to modalities.

Model	Statement	Signal dependence	Modality
Gaussian	$X = \mu + \epsilon$, $\epsilon\sim\mathcal{N}(0,\sigma^2)$	none	thermal and amplifier noise (CT, optical)
Poisson	$N\sim\mathrm{Poisson}(\lambda)$, $\mathbb{E}=\mathrm{Var}=\lambda$	$\sigma = \sqrt{\lambda}$	photon counting (PET, SPECT, low-dose CT)
Rician	$M = \sqrt{S_{\mathrm{re}}^2 + S_{\mathrm{im}}^2}$	biased high at low SNR	magnitude MRI

The standard quality metrics are

\[\begin{equation} \mathrm{SNR} = \frac{\mu_{\mathrm{signal}}}{\sigma_{\mathrm{noise}}}, \qquad \mathrm{CNR} = \frac{\lvert\mu_1 - \mu_2\rvert}{\sigma_{\mathrm{noise}}}, \tag{20} \end{equation}\]

and for photon-limited imaging $\mathrm{SNR} = \sqrt{\lambda}$, so doubling SNR costs a fourfold increase in dose or time. This is a hard ceiling on any downstream model: no classifier can recover information that the counting statistics never delivered, and a model that appears to do so is reading something else, usually a scanner signature.

set.seed(18)
N <- 20000; mu <- 25; sg <- 5
gauss <- mu + rnorm(N, 0, sg)
pois  <- rpois(N, mu)
rice  <- sqrt((mu + rnorm(N, 0, sg))^2 + rnorm(N, 0, sg)^2)
rice0 <- sqrt(rnorm(N, 0, sg)^2 + rnorm(N, 0, sg)^2)   # true signal = 0

nd <- rbind(data.frame(v = gauss, m = "Gaussian"),
            data.frame(v = pois,  m = "Poisson"),
            data.frame(v = rice,  m = "Rician (magnitude MRI)"))
ggplot(nd, aes(v, fill = m)) +
  geom_histogram(bins = 60, alpha = 0.55, position = "identity") +
  scale_fill_manual(values = bpad_pal[1:3]) +
  labs(title = "Noise distribution follows the physics of detection",
       x = "measured value", y = "count")

Three imaging noise models at a common mean signal of 25 counts or units. Gaussian noise is symmetric and signal-independent; Poisson noise is discrete with variance equal to its mean; the magnitude-MRI (Rician) distribution is right-skewed and, at low signal, has a positive mean even when the true signal is zero. Fitting a Gaussian likelihood to the right-hand panel biases every derived quantity.

Figure 4: Three imaging noise models at a common mean signal of 25 counts or units. Gaussian noise is symmetric and signal-independent; Poisson noise is discrete with variance equal to its mean; the magnitude-MRI (Rician) distribution is right-skewed and, at low signal, has a positive mean even when the true signal is zero. Fitting a Gaussian likelihood to the right-hand panel biases every derived quantity.

print(knitr::kable(data.frame(
  model = c("Gaussian", "Poisson", "Rician (A = 25)", "Rician (A = 0)"),
  mean  = round(c(mean(gauss), mean(pois), mean(rice), mean(rice0)), 3),
  sd    = round(c(sd(gauss),   sd(pois),   sd(rice),   sd(rice0)), 3),
  theory = c("mean 25, sd 5", sprintf("sd = sqrt(25) = %.2f", sqrt(25)),
             "approx Gaussian at high SNR",
             sprintf("mean = sigma*sqrt(pi/2) = %.2f", sg*sqrt(pi/2)))),
  caption = paste("Empirical moments against theory. The last row is the",
                  "low-SNR magnitude bias: the mean is positive although the",
                  "true signal is zero.")))

## 
## 
## Table: (\#tab:noise-models-demo)Empirical moments against theory. The last row is the low-SNR magnitude bias: the mean is positive although the true signal is zero.
## 
## |model           |   mean|    sd|theory                         |
## |:---------------|------:|-----:|:------------------------------|
## |Gaussian        | 24.964| 4.991|mean 25, sd 5                  |
## |Poisson         | 25.056| 5.005|sd = sqrt(25) = 5.00           |
## |Rician (A = 25) | 25.490| 4.946|approx Gaussian at high SNR    |
## |Rician (A = 0)  |  6.284| 3.272|mean = sigma*sqrt(pi/2) = 6.27 |

cat(sprintf("SNR (Gaussian) = %.1f;  SNR (Poisson) = %.1f = sqrt(lambda)\n",
            mu/sg, mean(pois)/sd(pois)))

## SNR (Gaussian) = 5.0;  SNR (Poisson) = 5.0 = sqrt(lambda)

Augmentation must match the physics. Adding independent Gaussian noise to images is not a universal robustness test. Reconstruction, iterative denoising, and spatial filtering all introduce spatially correlated noise, so real acquisition variability is coloured, not white. Testing a model with white noise and declaring it robust to low-dose CT or low-field MRI is a category error: the perturbation you applied is not the perturbation the scanner applies.

8.2.11 Basis expansions and transformations

A linear model can represent nonlinear dependence through a basis expansion,

\[\begin{equation} f(x) = \beta_0 + \sum_{m=1}^{M}\beta_m\,\phi_m(x), \tag{21} \end{equation}\]

where $\phi_m$ may be polynomials, splines, interactions, or physically motivated transformations. Three cautions apply.

Global polynomials misbehave. High-degree polynomials exhibit Runge’s phenomenon (violent oscillation near the domain boundary) and non-local influence: moving one data point changes predictions everywhere. Natural cubic splines and B-splines use localized piecewise bases and avoid both.

Log transformation changes the error model. If $\log Y = X\beta + \epsilon$ with $\epsilon\sim\mathcal{N}(0,\sigma^2)$, then $Y = e^{X\beta}e^{\epsilon}$, which is multiplicative error. Because $\mathbb{E}[e^{\epsilon}] = e^{\sigma^2/2} \ne 1$, naive back-transformation estimates the median, not the mean; the mean requires the correction

\[\begin{equation} \widehat{\mathbb{E}}[Y\mid X] = e^{X\widehat\beta}\,e^{\widehat\sigma^2/2}. \tag{22} \end{equation}\]

This matters directly in Section 8.8.7, where tumour volume is modelled on the log scale precisely because growth is multiplicative.

Complexity is a hyperparameter. Spline knot placement, polynomial degree, and Box–Cox exponents must be selected inside each training fold. Choosing them on the full dataset is leakage of exactly the kind Section 8.3.6 quantifies.

Section 8.2 summary.

Fitting is regularized empirical risk minimization, Eq. (9); the loss encodes a distributional assumption.
Never invert $X^{\mathsf T}X$ for collinear radiomic features: $\kappa(X^{\mathsf T}X) = \kappa(X)^2$. Use QR or a truncated SVD.
Expected error decomposes into bias$^2$, variance, and an irreducible floor set by the imaging physics.
Sensitivity and specificity belong to the test; predictive values belong to the test and the population; AUC hides the difference.
Measurement error attenuates coefficients and destabilizes selection; propagate uncertainty by the delta method or, for discrete pipelines, Monte Carlo.
Match the likelihood to the modality’s noise model, and make augmentation physically plausible rather than merely Gaussian.

Checkpoint 8.2. A radiomic study reports 1200 features from 90 patients and fits ordinary least squares by inverting $X^{\mathsf T}X$. Identify three distinct failures this guarantees, name the section of this chapter that addresses each, and state the minimum change that would make the analysis defensible.

) # 8.3 Study Design Before Model Fitting

A biomedical AI study should be designed backward from its intended use, not forward from an available dataset. Everything in this section happens before a single model is fitted, and most published failures in imaging AI are failures here rather than failures of algorithm choice.

8.3.1 Estimand, target population, and unit of analysis

Let $D$ be the development distribution and $T$ the deployment distribution. The quantity of practical interest is the deployment risk

\[\begin{equation} R_T(f) = \mathbb{E}_{(X,Y)\sim T}\big[L\{Y, f(X)\}\big], \tag{23} \end{equation}\]

whereas training minimizes an estimate of $R_D(f)$. A model can perform excellently under $D$ and fail under $T$ whenever scanners, sites, prevalence, patient mix, annotation practice, or clinical workflow differ. Section 8.6.6 measures a real instance of $D \ne T$ inside a single public cohort.

The unit of analysis must match the independent sampling unit.

Data structure	Splitting unit	Reason
One scan per patient	patient	prevents a person entering two partitions
Multiple scans per patient	patient	longitudinal images are strongly dependent
Multiple lesions per patient	patient, with lesion-level modelling	lesions share biology and acquisition
Multiple slices or patches	patient or volume	slice-level splitting is severe leakage
Multiple sites	site-aware patient split	needed to assess transportability
Repeated readers	patient plus reader-aware analysis	labels and predictions are reader-correlated

8.3.2 Outcome and reference-standard design

The reference standard must be defined before inspecting model results. If the observed label $Y$ is a noisy version of a latent truth $Y^{\ast}$, then $P(Y \ne Y^{\ast}) > 0$ and measured performance is bounded above by label quality. A model cannot be more accurate than its labels, and reported accuracy above the label-agreement ceiling is evidence of leakage, not of skill.

Differential label error is the dangerous case. Pathology is preferentially available for suspicious lesions; long-term follow-up is preferentially incomplete in low-risk patients. Both make the label depend on the predictors through the referral process rather than through biology.

A rigorous protocol records: who established the label and whether they were blinded to the model inputs; the interval between imaging and reference measurement; how indeterminate cases were handled; inter-reader reliability when labels are subjective; and whether any information appears in both the predictors and the label definition. For segmentation the reference is a spatial object, so consensus masks with uncertainty bands are more honest than a single nominally exact boundary.

8.3.3 Missing data and informative observation

Let $R_j = 1$ indicate that $X_j$ is observed. The classical mechanisms are missing completely at random (MCAR) ($R_j \perp (X,Y)$), missing at random (MAR) ($R_j \perp X_j \mid X_{-j}, Y$), and missing not at random (MNAR) (missingness still depends on the unobserved value after conditioning). Complete-case analysis is generally valid only under restrictive conditions, and mean imputation performed before splitting leaks the test distribution into training. A defensible workflow estimates imputation parameters on each training resample and applies them unchanged to the corresponding validation fold (Section 8.6.2).

Missingness may itself be informative and predictive. An advanced sequence may be omitted because the patient could not tolerate a long examination; a laboratory value may be absent because the clinician judged it unnecessary. Missingness indicators can improve prediction while encoding local workflow rather than biology, so they must be tested for transportability.

Absence is not a biological zero. A missing laboratory value, an unperformed sequence, an out-of-range truncated value, and a measured value of zero are four different states. Section 8.4.3 finds all four in the real cohort, and shows that a single call to as.numeric() silently destroys the distinction in every one of them.

8.3.4 Censoring, truncation, and competing events

For prognosis let $T$ be the event time and $C$ the censoring time. We observe

\[\begin{equation} \widetilde T = \min(T, C), \qquad \delta = \mathbb{I}(T \le C). \tag{24} \end{equation}\]

A binary endpoint such as “dead at last contact” discards follow-up duration and misclassifies everyone censored before the horizon. Standard survival methods handle incomplete follow-up but require conditionally independent censoring,

\[\begin{equation} T \perp C \mid X, \tag{25} \end{equation}\]

which fails when patients are lost to follow-up because their disease is progressing. Inverse probability of censoring weighting or joint longitudinal-survival modelling is then required.

Left truncation arises when a patient enters the risk set at time $V_i > 0$ after the biological origin. Resetting the origin to zero produces immortal-time bias, because patients who would have had the event before entry are never observed. The risk set must be

\[\begin{equation} Y_i(t) = \mathbb{I}(V_i < t \le \widetilde T_i). \tag{26} \end{equation}\]

Competing events occur when one event prevents another, for example non-cancer death preventing observed recurrence. Two hazards must then be distinguished. The cause-specific hazard

\[\begin{equation} \lambda_k(t) = \lim_{\Delta t\to0}\frac{P(t \le T < t+\Delta t,\ D = k \mid T \ge t)}{\Delta t} \tag{27} \end{equation}\]

addresses aetiology, whereas the subdistribution hazard underlying the Fine–Gray model keeps patients who experienced a competing event in the risk set and addresses absolute prognosis. Section 8.8.6 shows numerically that using $1 - \mathrm{KM}$ under competing events overestimates absolute risk substantially, and that the cumulative incidence function must be used instead.

8.3.5 Partition roles and nested resampling

Each partition gets exactly one role, and roles are not recyclable.

Partition	Permitted	Forbidden
Training folds	fit model and all preprocessing	any final performance claim
Tuning folds	select hyperparameters and thresholds	repeated redesign until favourable
Locked internal test	one-time assessment in the development setting	feature selection, calibration fitting
External test	assess a genuinely different site, scanner, time, or population	model revision followed by calling it untouched

For small cohorts, nested cross-validation replaces a fixed tuning split: an outer loop estimates generalization while an inner loop chooses hyperparameters. Repeated cross-validation measures resampling variability; it does not manufacture independent patients.

8.3.6 Leakage is a property of the pipeline, not of a variable

The most consequential idea in this section is that leakage is not a list of forbidden columns. It is any information path from held-out data into a fitting decision. Feature screening, imputation, standardization, threshold selection, hyperparameter choice, and even the decision to try another model after seeing test performance are all fitting decisions.

The demonstration below makes the cost concrete. It constructs data in which the outcome is independent of every predictor by construction, so the only honest answer is $\mathrm{AUC} = 0.5$. Then it computes cross-validated AUC two ways: screening the 500 features on all the data before cross-validating (the common shortcut), and screening strictly inside each training fold.

set.seed(42)

leak_experiment <- function(n = 100, p = 500, K = 5, n_keep = 10, reps = 40) {
  out <- matrix(NA_real_, reps, 2, dimnames = list(NULL, c("leaky", "nested")))
  for (b in seq_len(reps)) {
    X <- matrix(rnorm(n*p), n, p)
    y <- rbinom(n, 1, 0.5) == 1                       # independent of X by construction
    fold <- sample(rep_len(seq_len(K), n))

    ## (1) LEAKY: screen on all n rows, then cross-validate the chosen features
    sc  <- abs(apply(X, 2, function(z) cor(z, as.numeric(y))))
    top <- order(sc, decreasing = TRUE)[seq_len(n_keep)]
    o1  <- rep(NA_real_, n)
    for (k in seq_len(K)) {
      tr <- fold != k
      f  <- suppressWarnings(glm(y[tr] ~ ., data = data.frame(X[tr, top, drop = FALSE]),
                                 family = binomial()))
      o1[!tr] <- predict(f, newdata = data.frame(X[!tr, top, drop = FALSE]),
                         type = "response")
    }

    ## (2) NESTED: screen inside each training fold only
    o2 <- rep(NA_real_, n)
    for (k in seq_len(K)) {
      tr  <- fold != k
      s   <- abs(apply(X[tr, ], 2, function(z) cor(z, as.numeric(y[tr]))))
      tp  <- order(s, decreasing = TRUE)[seq_len(n_keep)]
      f   <- suppressWarnings(glm(y[tr] ~ ., data = data.frame(X[tr, tp, drop = FALSE]),
                                  family = binomial()))
      o2[!tr] <- predict(f, newdata = data.frame(X[!tr, tp, drop = FALSE]),
                         type = "response")
    }
    out[b, ] <- c(auc_rank(o1, y), auc_rank(o2, y))
  }
  out
}

L <- leak_experiment()
print(knitr::kable(data.frame(
  procedure = c("screen on all data, then cross-validate", "screen inside training folds only"),
  mean_AUC  = round(colMeans(L), 3),
  sd_AUC    = round(apply(L, 2, sd), 3),
  min_AUC   = round(apply(L, 2, min), 3),
  max_AUC   = round(apply(L, 2, max), 3),
  truth     = c(0.5, 0.5)),
  row.names = FALSE,
  caption = "Cross-validated AUC on pure noise, 40 replicates. The leaky pipeline never once reports the truth."))

## 
## 
## Table: (\#tab:leakage-simulation)Cross-validated AUC on pure noise, 40 replicates. The leaky pipeline never once reports the truth.
## 
## |procedure                               | mean_AUC| sd_AUC| min_AUC| max_AUC| truth|
## |:---------------------------------------|--------:|------:|-------:|-------:|-----:|
## |screen on all data, then cross-validate |    0.820|  0.031|   0.738|   0.889|   0.5|
## |screen inside training folds only       |    0.484|  0.079|   0.325|   0.646|   0.5|

ld <- rbind(data.frame(auc = L[, "leaky"],  s = "screen on all data"),
            data.frame(auc = L[, "nested"], s = "screen inside folds"))
ggplot(ld, aes(auc, fill = s)) +
  geom_histogram(bins = 26, alpha = 0.6, position = "identity") +
  geom_vline(xintercept = 0.5, linetype = "dashed", colour = "grey30", linewidth = 0.8) +
  scale_fill_manual(values = bpad_pal[2:1]) +
  labs(title = "Cross-validated AUC on data with no signal whatsoever",
       x = "out-of-fold AUC", y = "replicates")

What leakage buys. Both procedures cross-validate on data in which the label is independent of all 500 features, so the true AUC is exactly 0.5 (dashed line). Screening features on the full dataset before cross-validation manufactures an apparent AUC above 0.8 across every replicate. Screening inside each training fold recovers the truth. No model, metric, or resampling scheme can detect this after the fact.

Figure 5: What leakage buys. Both procedures cross-validate on data in which the label is independent of all 500 features, so the true AUC is exactly 0.5 (dashed line). Screening features on the full dataset before cross-validation manufactures an apparent AUC above 0.8 across every replicate. Screening inside each training fold recovers the truth. No model, metric, or resampling scheme can detect this after the fact.

Read the numbers above carefully. The leaky pipeline is not slightly optimistic; it converts pure noise into an AUC around 0.81 with a standard deviation of 0.03, so it looks both strong and stable. Stability is not validity. This is why the ordering of operations in a pipeline is a scientific claim and must be reported as one.

8.3.7 Sample size, events, dimensionality, and learning curves

Three quantities constrain what a cohort can support.

Events, not rows. For binary and survival outcomes the binding constraint is the number of events in the minority class. A useful rule of thumb is at least 10 to 20 events per candidate predictor for stable coefficient estimation, and modern shrinkage-based sample-size formulae are stricter still. This chapter’s cohort has 210 patients, 18 benign masses, and 21 deaths: a handful of predictors, not a hundred.

Dimensionality. As $p$ grows at fixed $n$, points become nearly equidistant and neighbourhood-based methods lose their meaning. The demonstration below shows the concentration of pairwise distances directly.

set.seed(1)
dims <- c(2, 3, 5, 10, 20, 50, 100, 200, 400)
cd <- t(sapply(dims, function(p) {
  X <- matrix(runif(400*p), 400, p)
  d <- as.matrix(dist(X)); d <- d[upper.tri(d)]
  c(ratio = min(d)/max(d), cv = sd(d)/mean(d))
}))
print(knitr::kable(data.frame(dimension = dims,
                        min_over_max = round(cd[, "ratio"], 3),
                        coef_variation = round(cd[, "cv"], 3)),
             caption = "Pairwise-distance concentration for 400 uniform points."))

## 
## 
## Table: (\#tab:curse-dimensionality)Pairwise-distance concentration for 400 uniform points.
## 
## | dimension| min_over_max| coef_variation|
## |---------:|------------:|--------------:|
## |         2|        0.001|          0.479|
## |         3|        0.010|          0.369|
## |         5|        0.047|          0.277|
## |        10|        0.128|          0.195|
## |        20|        0.256|          0.135|
## |        50|        0.472|          0.084|
## |       100|        0.590|          0.059|
## |       200|        0.714|          0.041|
## |       400|        0.780|          0.030|

p1 <- ggplot(data.frame(p = dims, r = cd[, "ratio"]), aes(p, r)) +
  geom_line(colour = bpad_pal[1], linewidth = 0.95) + geom_point(size = 1.8) +
  scale_x_log10() +
  labs(title = "min / max pairwise distance", x = "dimension (log)", y = "ratio")
p2 <- ggplot(data.frame(p = dims, cv = cd[, "cv"]), aes(p, cv)) +
  geom_line(colour = bpad_pal[2], linewidth = 0.95) + geom_point(size = 1.8) +
  scale_x_log10() + scale_y_log10() +
  labs(title = "coefficient of variation of distances",
       x = "dimension (log)", y = "sd / mean (log)")
bpad_grid(p1, p2, ncol = 2)

Distance concentration in high dimensions. Left: for uniformly distributed points, the ratio of the minimum to the maximum pairwise distance rises toward one as dimension increases, so 'nearest' and 'farthest' become nearly indistinguishable. Right: the same effect as the coefficient of variation of pairwise distances collapsing. This is why k-nearest-neighbour methods degrade on wide radiomic feature sets unless the dimension is first reduced.

Figure 6: Distance concentration in high dimensions. Left: for uniformly distributed points, the ratio of the minimum to the maximum pairwise distance rises toward one as dimension increases, so ‘nearest’ and ‘farthest’ become nearly indistinguishable. Right: the same effect as the coefficient of variation of pairwise distances collapsing. This is why k-nearest-neighbour methods degrade on wide radiomic feature sets unless the dimension is first reduced.

Learning curves. Plotting performance against training-set size distinguishes two very different situations: a curve still rising at the full sample size means more data would help, whereas a curve that has plateaued below the clinical requirement means the features or labels, not the sample size, are the limiting factor. This distinction determines whether the next investment should be recruitment or measurement.

Section 8.3 summary.

Design backward from intended use; the estimand is a deployment risk, Eq. (23), not a training risk.
Split at the independent unit; label quality caps achievable performance.
Distinguish MCAR, MAR, and MNAR, and never impute before splitting.
Censoring, truncation, and competing events each require a specific estimator; a binary endpoint is usually inadmissible.
Leakage is a pipeline property: screening on all data turned pure noise into AUC 0.81.
Events, not rows, and dimension, not row count, set what a cohort can support.

Checkpoint 8.3. You are given 210 patients, 21 deaths, and 900 radiomic features, and asked to build a survival model. State the maximum number of candidate predictors you would consider, the resampling scheme, and the single design change that would most improve the credibility of the result.

8.4 Running Case Study: the Real KiTS19 Kidney Cohort

8.4.1 Provenance

The executable spine of this chapter uses the 210 publicly released KiTS19 training cases: adults who underwent partial or radical nephrectomy for a suspected renal malignancy, with contrast-enhanced preoperative CT, expert voxel-level reference segmentations (background, kidney, tumour), clinical records, and surgical outcomes. The clinical file is kits.json; the reference masks are segmentation.nii.gz per case.

Nothing in the analytic spine is simulated. Where a demonstration requires something the cohort does not record, such as cause of death (Section 8.8.6), the text says so explicitly and labels the computation as a designed simulation.

8.4.2 Loading real records while preserving their pathologies

Real clinical files contain values that are not numbers even though they occupy numeric fields. This loader keeps every one of them, splitting each such field into a numeric value and a flag that records why the value is absent. This is the single most important habit in this section: as.numeric() applied blindly converts information into NA and the analyst never learns what was lost.

## ---------------------------------------------------------------------------
## Parse a field that is "usually numeric" into (value, flag). Anything that is
## not a plain number is preserved verbatim as a flag level, never coerced away.
## ---------------------------------------------------------------------------
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(...) {
    path_ <- c(...)
    vapply(j, function(z) {
      v <- z
      for (p in path_) { v <- v[[p]]; if (is.null(v)) return(NA_character_) }
      if (length(v) != 1) return(NA_character_)
      as.character(v)
    }, character(1))
  }
  egfr_pre  <- num_with_flag(pull("last_preop_egfr",  "value"))
  egfr_post <- num_with_flag(pull("last_postop_egfr", "value"))
  isch      <- num_with_flag(pull("ischemia_time"))
  hosp      <- num_with_flag(pull("hospitalization"))

  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"))),
    smoking              = pull("smoking_history"),
    radiographic_size_cm = suppressWarnings(as.numeric(pull("radiographic_size"))),
    pathologic_size_cm   = suppressWarnings(as.numeric(pull("pathologic_size"))),
    malignant            = pull("malignant") == "TRUE",
    t_stage              = pull("pathology_t_stage"),
    isup                 = suppressWarnings(as.numeric(pull("tumor_isup_grade"))),
    histology            = pull("tumor_histologic_subtype"),
    ## THREE DISTINCT SURGICAL FIELDS -- see Section 8.4.3
    procedure            = pull("surgical_procedure"),   # partial vs radical (extent)
    access               = pull("surgery_type"),         # open/laparoscopic/robotic
    route                = pull("surgical_approach"),    # transperitoneal vs retroperitoneal
    egfr_pre             = egfr_pre$value,  egfr_pre_flag  = egfr_pre$flag,
    egfr_post            = egfr_post$value, egfr_post_flag = egfr_post$flag,
    ischemia_min         = isch$value,      ischemia_flag  = isch$flag,
    hosp_days            = hosp$value,      hosp_flag      = hosp$flag,
    vital_status         = pull("vital_status"),
    vital_days           = suppressWarnings(as.numeric(pull("vital_days_after_surgery"))),
    ckd                  = pull("comorbidities", "chronic_kidney_disease") == "TRUE",
    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)
}

## ---------------------------------------------------------------------------
## Documented surrogate cohort, used ONLY if the real file cannot be reached.
## It reproduces the released schema, the marginal distributions, AND the four
## data pathologies of Section 8.4.3, so every lesson below still lands. It is
## labelled in `data_source` so no number in this chapter is ever ambiguous.
## ---------------------------------------------------------------------------
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, shape = 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)))
  egfr_post_raw <- pmax(5, egfr - drop)
  malig <- runif(n) < inv_logit(-0.4 + 0.55*size)
  lp    <- 0.02*(age - 60) + 0.16*size + 0.35*(procedure == "radical_nephrectomy")
  tt    <- rexp(n, rate = 0.02*exp(lp)); cc <- runif(n, 0.05, 9)
  time  <- pmin(tt, cc); dead <- tt <= cc
  sz    <- sample(c(0.5, 1.0, 1.25, 2.0, 3.0, 5.0), n, TRUE, c(20, 29, 60, 55, 30, 16)/n)
  ## inject the real pathologies
  egfr_chr      <- ifelse(egfr >= 90, ">=90", as.character(egfr))
  egfr_chr[sample(n, 57)] <- NA
  egfr_post_chr <- ifelse(egfr_post_raw >= 90, ">=90", as.character(egfr_post_raw))
  egfr_post_chr[age < 16] <- "age<16"
  egfr_post_chr[sample(n, 53)] <- NA
  isch_chr <- ifelse(procedure == "radical_nephrectomy", "not_applicable",
                     as.character(round(pmax(4, rnorm(n, 22, 8)))))
  isch_chr[procedure == "partial_nephrectomy"][1:10] <- NA
  hosp_chr <- as.character(pmax(1, round(rlnorm(n, 1.0, 0.6)))); hosp_chr[1] <- "died_before_discharge"
  ep <- num_with_flag(egfr_chr); eo <- num_with_flag(egfr_post_chr)
  ic <- num_with_flag(isch_chr); hp <- num_with_flag(hosp_chr)
  data.frame(
    case_id = sprintf("case_%05d", seq_len(n) - 1L), age = age,
    gender = sample(c("male","female"), n, TRUE, c(0.62, 0.38)),
    bmi = round(pmax(16, rnorm(n, 29.5, 6.4)), 1),
    smoking = sample(c("never_smoked","previous_smoker","current_smoker"),
                     n, TRUE, c(100, 76, 34)/n),
    radiographic_size_cm = round(size, 1),
    pathologic_size_cm   = round(pmax(0.3, size + rnorm(n, -0.05, 1.17)), 1),
    malignant = malig, t_stage = NA_character_, isup = NA_real_,
    histology = NA_character_, procedure = procedure,
    access = sample(c("robotic","open","laparoscopic"), n, TRUE, c(122,60,28)/n),
    route  = sample(c("Transperitoneal","Retroperitoneal"), n, TRUE, c(171,39)/n),
    egfr_pre = ep$value, egfr_pre_flag = ep$flag,
    egfr_post = eo$value, egfr_post_flag = eo$flag,
    ischemia_min = ic$value, ischemia_flag = ic$flag,
    hosp_days = hp$value, hosp_flag = hp$flag,
    vital_status = ifelse(dead, "dead", "censored"),
    vital_days = round(time*365.25),
    ckd = runif(n) < 0.11,
    spacing_x_mm = round(runif(n, 0.6, 1.0), 4), spacing_z_mm = sz,
    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()

## --- derived analysis variables: explicit, documented, reversible ----------
kits$event      <- as.integer(kits$vital_status == "dead")   # right-censored outcome
kits$time_yr    <- kits$vital_days/365.25
kits$radical    <- as.integer(kits$procedure == "radical_nephrectomy")
kits$open       <- as.integer(kits$access    == "open")
kits$egfr_drop  <- kits$egfr_pre - kits$egfr_post
## A REAL technical stratum: KiTS19 slice thickness spans 0.5-5.0 mm.
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 cases, %d fields\n", nrow(kits), ncol(kits)))

## n = 210 cases, 34 fields

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

## malignant prevalence = 0.914 (18 benign)

cat(sprintf("deaths = %d, censored = %d, follow-up to %.1f years\n",
            sum(kits$event), sum(!kits$event), max(kits$time_yr, na.rm = TRUE)))

## deaths = 21, censored = 189, follow-up to 8.4 years

8.4.3 Data-quality control: five real findings

Assertions should protect an analysis, not decorate it. The checks below are written so that the structural requirements must hold, while the substantive findings are reported rather than suppressed.

stopifnot(nrow(kits) > 0, !any(duplicated(kits$case_id)))
stopifnot(all(kits$vital_status %in% c("censored", "dead")))
stopifnot(all(kits$procedure %in% c("partial_nephrectomy", "radical_nephrectomy")))

cat("=== FINDING 1: the age range includes paediatric nephrectomy ===\n")

## === FINDING 1: the age range includes paediatric nephrectomy ===

print(summary(kits$age))

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1.00   50.00   61.00   58.35   68.00   90.00

cat(sprintf("cases under 18: %d\n\n", sum(kits$age < 18, na.rm = TRUE)))

## cases under 18: 4

cat("=== FINDING 2: laboratory truncation at the top of the reportable range ===\n")

## === FINDING 2: laboratory truncation at the top of the reportable range ===

print(table(preoperative = kits$egfr_pre_flag))

## preoperative
##    >=90 missing numeric 
##      33      57     120

print(table(postoperative = kits$egfr_post_flag))

## postoperative
##    >=90  age<16 missing numeric 
##      22       2      53     133

cat("'>=90' is a right-censored laboratory value, not a missing one.\n\n")

## '>=90' is a right-censored laboratory value, not a missing one.

cat("=== FINDING 3: missingness is structural, not random ===\n")

## === FINDING 3: missingness is structural, not random ===

print(table(procedure = kits$procedure, ischemia_time = kits$ischemia_flag))

##                      ischemia_time
## procedure             missing not_applicable numeric
##   partial_nephrectomy      10              0     130
##   radical_nephrectomy       0             70       0

cat("\n=== FINDING 4: three distinct surgical fields are easily confused ===\n")

## 
## === FINDING 4: three distinct surgical fields are easily confused ===

print(table(procedure_extent = kits$procedure))

## procedure_extent
## partial_nephrectomy radical_nephrectomy 
##                 140                  70

print(table(operative_access = kits$access))

## operative_access
## laparoscopic         open      robotic 
##           28           60          122

print(table(anatomic_route   = kits$route))

## anatomic_route
## Retroperitoneal Transperitoneal 
##              39             171

cat("\n=== FINDING 5: string sentinels hiding inside a numeric field ===\n")

## 
## === FINDING 5: string sentinels hiding inside a numeric field ===

print(table(hospitalization = kits$hosp_flag))

## hospitalization
## died_before_discharge               numeric 
##                     1                   209

Each finding changes what may legitimately be done next.

1. Paediatric cases. The cohort spans a very wide age range and includes at least one Wilms tumour. A model trained on adults and applied to children is extrapolating across a different disease. Either exclude them with a documented rule or model age flexibly and report the subgroup separately.

2. Truncated laboratory values. An eGFR reported as ">=90" is right-censored at 90, not missing. Imputing 90 understates the true value; imputing the mean of the numeric cases understates it further and induces a spurious negative correlation with the truncation indicator. The defensible options are to analyse eGFR as a censored outcome, to dichotomize at a clinically meaningful threshold, or to carry the truncation flag as a covariate. Section 8.7.6 does the last of these and reports the consequence.

3. Structured missingness. Ischaemia time is "not_applicable" for exactly the radical nephrectomies and for none of the partials, because no vascular clamping is performed when the whole kidney is removed. This is not MCAR, MAR, or MNAR in the usual sense: it is a deterministic function of the procedure. Imputing it would fabricate a clamping time for an operation in which no clamp was applied, and including it as a predictor would smuggle the procedure into the model under another name.

4. Three surgical fields. This is the subtlest trap in the file. surgical_procedure records the extent of parenchyma removed (partial or radical); surgery_type records the operative access (open, laparoscopic, robotic); surgical_approach records the anatomic route (transperitoneal or retroperitoneal). They are three different clinical variables with three different meanings.

A bug worth naming. A natural-looking line such as radical <- as.numeric(surgery_type == "radical") is syntactically valid, silently wrong, and produces an indicator that is identically zero, because surgery_type never takes the value "radical" in this file. The model then fits a constant column, the coefficient is unidentifiable or dropped, and nothing in the output announces the error. The same trap applies to surgical_approach == "open". Always tabulate a categorical variable before encoding it, and prefer %in% against an explicit, asserted level set.

5. Semantic sentinels. Length of stay is recorded as "died_before_discharge" for one patient. That is not a missing number, it is the most severe possible outcome, and coercing it to NA deletes the worst case in the cohort from any length-of-stay analysis.

8.4.4 Exploratory analysis and the radiographic-pathologic gap

k <- kits
p1 <- ggplot(k, aes(procedure, radiographic_size_cm, fill = procedure)) +
  geom_boxplot(alpha = 0.55, outlier.size = 0.8, show.legend = FALSE) +
  scale_fill_manual(values = bpad_pal[1:2]) +
  labs(title = "(a) size by operative extent", x = NULL, y = "radiographic size (cm)")

brk <- cut(k$radiographic_size_cm, breaks = c(0, 2, 4, 7, 10, Inf),
           labels = c("<2", "2-4", "4-7", "7-10", ">10"))
agg <- aggregate(malignant ~ brk, data = data.frame(brk, malignant = k$malignant), FUN = mean)
agg$n <- as.vector(table(brk)[as.character(agg$brk)])
p2 <- ggplot(agg, aes(brk, malignant)) +
  geom_col(fill = bpad_pal[3], alpha = 0.8) +
  geom_text(aes(label = paste0("n=", n)), vjust = -0.4, size = 2.9) +
  scale_y_continuous(limits = c(0, 1.08)) +
  labs(title = "(b) malignant proportion by size", x = "size (cm)", y = "proportion malignant")

p3 <- ggplot(k, aes(time_yr, fill = factor(event, labels = c("censored", "died")))) +
  geom_histogram(bins = 30, alpha = 0.7, position = "stack") +
  scale_fill_manual(values = c(bpad_pal[7], bpad_pal[2])) +
  labs(title = "(c) follow-up time and events", x = "years after surgery", y = "patients")

ba <- k[complete.cases(k[, c("radiographic_size_cm", "pathologic_size_cm")]), ]
ba$avg <- (ba$radiographic_size_cm + ba$pathologic_size_cm)/2
ba$dif <- ba$radiographic_size_cm - ba$pathologic_size_cm
mu_d <- mean(ba$dif); sd_d <- sd(ba$dif)
p4 <- ggplot(ba, aes(avg, dif)) +
  geom_point(alpha = 0.5, colour = bpad_pal[1], size = 1.5) +
  geom_hline(yintercept = mu_d, colour = bpad_pal[2], linewidth = 0.8) +
  geom_hline(yintercept = mu_d + c(-1.96, 1.96)*sd_d, linetype = "dashed",
             colour = bpad_pal[2]) +
  labs(title = "(d) Bland-Altman: radiographic vs pathologic size",
       x = "mean of the two measurements (cm)", y = "radiographic - pathologic (cm)")
bpad_grid(p1, p2, p3, p4, ncol = 2)

Exploratory view of the real cohort. (a) Tumour size distribution by procedure: radical nephrectomy is performed for markedly larger tumours, which will confound every unadjusted comparison in Section 8.8.4. (b) Malignancy proportion rises with size, but even the smallest masses in this surgical series are mostly malignant. (c) Follow-up is short relative to the event rate. (d) Bland-Altman agreement between preoperative radiographic size and postoperative pathologic size: correlation is high yet the limits of agreement span more than four centimetres.

Figure 7: Exploratory view of the real cohort. (a) Tumour size distribution by procedure: radical nephrectomy is performed for markedly larger tumours, which will confound every unadjusted comparison in Section 8.8.4. (b) Malignancy proportion rises with size, but even the smallest masses in this surgical series are mostly malignant. (c) Follow-up is short relative to the event rate. (d) Bland-Altman agreement between preoperative radiographic size and postoperative pathologic size: correlation is high yet the limits of agreement span more than four centimetres.

cat(sprintf("Pearson r (radiographic, pathologic) = %.3f\n",
            cor(ba$radiographic_size_cm, ba$pathologic_size_cm)))

## Pearson r (radiographic, pathologic) = 0.931

cat(sprintf("mean bias = %+.3f cm; limits of agreement = [%+.2f, %+.2f] cm\n",
            mu_d, mu_d - 1.96*sd_d, mu_d + 1.96*sd_d))

## mean bias = +0.055 cm; limits of agreement = [-2.25, +2.36] cm

cat(sprintf("median tumour size = %.1f cm: the agreement band (%.2f cm wide) is wider\n",
            median(ba$radiographic_size_cm), 3.92*sd_d))

## median tumour size = 4.0 cm: the agreement band (4.61 cm wide) is wider

cat("than the median tumour itself, so the two measurements are not interchangeable\n")

## than the median tumour itself, so the two measurements are not interchangeable

cat("for any threshold-based decision.\n")

## for any threshold-based decision.

High correlation is not agreement. The two size measurements correlate at about $r = 0.93$, which many papers would report as validation. The Bland–Altman limits of agreement, however, span more than four centimetres, comparable to the median tumour size itself. Correlation measures whether two variables move together; agreement measures whether they can be substituted for one another. For a threshold-based clinical decision at 4 cm, only agreement matters. Chapter 7’s discussion of measurement uncertainty is the right lens here, not Chapter 1’s correlation coefficient.

8.4.5 Locking development and external-test roles

The external test set is defined by a real, documented technical property rather than an invented “site”: thin-slice ($\le 1$ mm) versus thick-slice ($> 1$ mm) reconstructions. Slice thickness changes partial-volume averaging, and therefore changes every spatial and texture feature computed downstream. This is precisely the kind of shift a deployed model meets.

dev_set <- kits[kits$acq_group == "thick_slice", ]
ext_set <- kits[kits$acq_group == "thin_slice",  ]

cmp <- function(v, label) {
  d <- dev_set[[v]]; e <- ext_set[[v]]
  data.frame(variable = label,
             development = sprintf("%.2f", median(d, na.rm = TRUE)),
             external    = sprintf("%.2f", median(e, na.rm = TRUE)),
             wilcoxon_p  = signif(suppressWarnings(
               stats::wilcox.test(d, e)$p.value), 3))
}
print(knitr::kable(rbind(
  cmp("spacing_z_mm",        "slice thickness (mm)"),
  cmp("spacing_x_mm",        "in-plane spacing (mm)"),
  cmp("radiographic_size_cm","tumour size (cm)"),
  cmp("age",                 "age (years)"),
  cmp("bmi",                 "BMI")),
  row.names = FALSE,
  caption = "The held-out group is genuinely technically different, and the difference is not confined to geometry: tumour size also differs, so the shift is simultaneously technical and clinical."))

## 
## 
## Table: (\#tab:lock-partitions)The held-out group is genuinely technically different, and the difference is not confined to geometry: tumour size also differs, so the shift is simultaneously technical and clinical.
## 
## |variable              |development |external | wilcoxon_p|
## |:---------------------|:-----------|:--------|----------:|
## |slice thickness (mm)  |5.00        |0.50     |     0.0000|
## |in-plane spacing (mm) |0.78        |0.81     |     0.2780|
## |tumour size (cm)      |4.30        |2.80     |     0.0036|
## |age (years)           |61.00       |60.00    |     0.4260|
## |BMI                   |29.47       |32.09    |     0.1560|

cat(sprintf("development n = %d (malignant %.3f, deaths %d)\n",
            nrow(dev_set), mean(dev_set$malignant), sum(dev_set$event)))

## development n = 161 (malignant 0.919, deaths 20)

cat(sprintf("external    n = %d (malignant %.3f, deaths %d)\n",
            nrow(ext_set), mean(ext_set$malignant), sum(ext_set$event)))

## external    n = 49 (malignant 0.898, deaths 1)

cat("\nAll model development from here on uses `dev_set` only.\n")

## 
## All model development from here on uses `dev_set` only.

cat("`ext_set` is touched exactly once, in Section 8.7.2.\n")

## `ext_set` is touched exactly once, in Section 8.7.2.

Section 8.4 summary.

Real records contain truncated values, semantic sentinels, deterministic missingness, and near-identically-named fields with different meanings; a blind as.numeric() destroys all four.
Ischaemia time is not_applicable for exactly the radical nephrectomies: missingness can be a deterministic function of another variable, which makes both imputation and inclusion indefensible.
Three surgical fields (extent, access, route) are routinely confused, and confusing them produces a constant predictor with no error message.
Radiographic and pathologic tumour size correlate at 0.93 with limits of agreement wider than four centimetres: correlation is not agreement.
The development/external split uses a real acquisition property, and that property is entangled with tumour size, so technical and clinical shift arrive together.

Checkpoint 8.4. The external (thin-slice) subset has systematically smaller tumours than the development subset. Give one plausible clinical reason for this association, and explain why it means that a drop in external performance cannot be attributed to slice thickness alone.

8.5 From CT Volumes and Masks to Quantitative Features

8.5.1 Image tensors, coordinates, and label volumes

A CT volume is a third-order array $I \in \mathbb{R}^{n_x\times n_y\times n_z}$ (Chapter 1, Section 1.3), and a segmentation is a co-registered label volume $M$ with $M_v \in \{0,1,2\}$ for background, kidney, and tumour. The bridge from array indices to physical space is the voxel spacing $(s_x, s_y, s_z)$ in millimetres, so a single voxel occupies

\[\begin{equation} \Delta V = s_x s_y s_z \quad [\text{mm}^3], \tag{28} \end{equation}\]

and the volume of a labelled region containing $N_\ell$ voxels is

\[\begin{equation} V_\ell = N_\ell\, s_x s_y s_z . \tag{29} \end{equation}\]

Equation (29) is the entire reason a biomedical physicist rather than a programmer should compute imaging biomarkers: a voxel count is a dimensionless integer, and it becomes a clinical quantity only through geometry that lives in the file header.

Never take geometry from a spreadsheet. The KiTS19 metadata records voxel spacing in the order $(x, y, z)$, while the NIfTI header of the corresponding volume stores pixdim in a different axis order for these reconstructions. Reading spacing from one source and array dimensions from the other silently produces volumes that are wrong by the ratio of two spacings. The code below reads both from the same header, and Section 8.5.3 verifies the result against an independently recorded clinical measurement.

8.5.2 Reading a NIfTI volume from its byte-level header

Reading the header yourself is a biomedical-physics exercise, not an inconvenience. The bytes that turn a voxel count into a physical volume are these:

Bytes	Field	Meaning
0–3	`sizeof_hdr`	must be 348; also reveals byte order
40–55	`dim[0..7]`	`dim[0]` is the number of used dimensions
70–71	`datatype`	numeric type code of the voxel data
72–73	`bitpix`	bits per voxel
76–107	`pixdim[0..7]`	physical voxel size, millimetres
108–111	`vox_offset`	byte offset at which voxel data begin
112–119	`scl_slope`, `scl_inter`	affine intensity rescaling

## ---------------------------------------------------------------------------
## A NIfTI-1 reader written directly against the file specification.
## No imaging package required. Byte order is detected from sizeof_hdr.
## ---------------------------------------------------------------------------
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: sizeof_hdr = ", sz)
  gi <- function(off, n, size)
    readBin(hdr[(off + 1):(off + n*size)], "integer", n = n, size = size, endian = endian)
  gd <- function(off, n)
    readBin(hdr[(off + 1):(off + 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),
       bitpix     = gi(72, 1, 2),
       pixdim     = pd[2:(1 + d[1])],
       vox_offset = gd(108, 1),
       scl_slope  = gd(112, 1),
       scl_inter  = gd(116, 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 NIfTI datatype code ", dt))

8.5.3 Streaming morphometrics in one pass

A KiTS19 label volume holds of order $10^8$ voxels. Loading it entirely into memory is unnecessary: every quantity we need can be accumulated plane by plane in $O(n_x n_y)$ memory. The function below makes one streaming pass and returns volume, surface area, sphericity, equivalent diameter, cross-sectional-area profile, and slice extent, all in physical units.

Two design choices deserve comment. Surface area is accumulated as the total area of exposed voxel faces, counting in-plane boundaries with the appropriate in-plane spacings and through-plane boundaries with the through-plane spacing. This is a consistent, reproducible definition; it is not the only one, and it systematically overestimates a smooth surface because a voxelized boundary is staircase-shaped. The slice axis is identified from the header as the anisotropic axis rather than assumed to be the third, because reconstruction pipelines permute axes.

mask_morphometrics <- function(path, labels = c(kidney = 1, tumor = 2)) {
  h  <- read_nifti_header(path)
  tp <- nifti_type(h$datatype)
  d  <- h$dim[1:3]; nx <- d[1]; ny <- d[2]; nz <- d[3]
  sp <- h$pixdim[1:3]

  con <- if (grepl("\\.gz$", path)) gzfile(path, "rb") else file(path, "rb")
  on.exit(close(con))
  invisible(readBin(con, "raw", n = h$vox_offset))    # skip to the voxel data

  st   <- lapply(labels, function(L) list(n = 0, area = 0,
                 p1 = numeric(nx), p2 = numeric(ny), p3 = numeric(nz)))
  prev <- vector("list", length(labels))

  for (k in seq_len(nz)) {
    plane <- matrix(readBin(con, tp$what, n = nx*ny,
                            size = tp$size, signed = tp$signed), nx, ny)
    for (li in seq_along(labels)) {
      b <- plane == labels[li]; nv <- sum(b); o <- st[[li]]
      if (nv > 0) {
        o$n  <- o$n + nv
        o$p1 <- o$p1 + rowSums(b); o$p2 <- o$p2 + colSums(b); o$p3[k] <- nv
        ## in-plane exposed faces (edges of the array plus internal boundaries)
        o$area <- o$area +
          (sum(b[1, ]) + sum(b[nx, ]))*sp[2]*sp[3] +
          (sum(b[, 1]) + sum(b[, ny]))*sp[1]*sp[3] +
          (sum(b[-1, ]  & !b[-nx, ]) + sum(b[-nx, ] & !b[-1, ]))*sp[2]*sp[3] +
          (sum(b[, -1]  & !b[, -ny]) + sum(b[, -ny] & !b[, -1]))*sp[1]*sp[3]
      }
      ## through-plane exposed faces, comparing with the previous plane
      pb <- prev[[li]]
      if (is.null(pb)) { if (nv > 0) o$area <- o$area + nv*sp[1]*sp[2]
      } else            o$area <- o$area + (sum(b & !pb) + sum(pb & !b))*sp[1]*sp[2]
      prev[[li]] <- b; st[[li]] <- o
    }
  }
  for (li in seq_along(labels)) {
    pb <- prev[[li]]
    if (!is.null(pb)) st[[li]]$area <- st[[li]]$area + sum(pb)*sp[1]*sp[2]
  }

  ## the through-plane (slice) axis is the anisotropic one; do not assume axis 3
  slice_axis <- if (length(unique(round(sp, 3))) == 1) 3L
                else which.max(abs(sp - stats::median(sp)))
  inplane_mm2 <- prod(sp[-slice_axis])

  res <- list(dim = d, spacing_mm = sp, slice_axis = slice_axis,
              slice_thickness_mm = sp[slice_axis])
  for (li in seq_along(labels)) {
    nm <- names(labels)[li]; o <- st[[li]]
    V  <- o$n*prod(sp)              # voxel count times voxel volume
    prof <- list(o$p1, o$p2, o$p3)[[slice_axis]]*inplane_mm2
    idx  <- which(prof > 0)
    res[[paste0(nm, "_volume_cm3")]]        <- V/1000
    res[[paste0(nm, "_surface_cm2")]]       <- o$area/100
    res[[paste0(nm, "_sphericity")]]        <- if (o$area > 0)
                                                (pi^(1/3)*(6*V)^(2/3))/o$area else NA_real_
    res[[paste0(nm, "_eq_diameter_cm")]]    <- 2*(3*V/(4*pi))^(1/3)/10
    res[[paste0(nm, "_max_slice_area_cm2")]]<- if (length(idx)) max(prof)/100 else 0
    res[[paste0(nm, "_n_slices")]]          <- length(idx)
    res[[paste0(nm, "_extent_mm")]]         <- if (length(idx))
      (diff(range(idx)) + 1)*sp[slice_axis] else 0
    res[[paste0(nm, "_profile_cm2")]]       <- prof/1000*10   # cm^2 per slice
  }
  res
}

8.5.4 Shape features in physical units

Three dimensionless descriptors are used throughout. The equivalent spherical diameter

\[\begin{equation} d_{\mathrm{eq}} = 2\left(\frac{3V}{4\pi}\right)^{1/3} \tag{30} \end{equation}\]

converts a volume into the diameter a sphere of that volume would have, which is what makes a mask-derived measurement comparable with a radiologist’s calliper reading. Sphericity

\[\begin{equation} \Psi = \frac{\pi^{1/3}\,(6V)^{2/3}}{S} \in (0, 1] \tag{31} \end{equation}\]

equals one exactly for a sphere and decreases with surface irregularity. The surface-to-volume index $S/V^{2/3}$ is scale-invariant but unnormalized. All three are invariant to the choice of length unit provided area and volume use the same base unit, and all three are affected by voxel size through the staircase surface, which is why Section 8.5.8 matters.

mask_urls <- sprintf(
  "https://raw.githubusercontent.com/neheller/kits19/master/data/case_%05d/segmentation.nii.gz",
  0:(max(bpad_n_mask, 1) - 1))

morph_cache <- file.path(bpad_cache_dir, sprintf("morph_%d.rds", bpad_n_mask))
if (file.exists(morph_cache)) {
  morph_list <- readRDS(morph_cache)
} else {
  morph_list <- list()
  for (i in seq_along(mask_urls)) {
    pth <- bpad_fetch(mask_urls[i], sprintf("seg_%05d.nii.gz", i - 1))
    if (is.na(pth)) next
    m <- try(mask_morphometrics(pth), silent = TRUE)
    if (!inherits(m, "try-error")) {
      m$case_id <- sprintf("case_%05d", i - 1)
      morph_list[[length(morph_list) + 1]] <- m
    }
  }
  if (length(morph_list)) saveRDS(morph_list, morph_cache)
}
have_imaging <- length(morph_list) > 0
cat(sprintf("real segmentation volumes processed: %d\n", length(morph_list)))

## real segmentation volumes processed: 6

if (have_imaging) {
  img <- do.call(rbind, lapply(morph_list, function(m) data.frame(
    case_id             = m$case_id,
    slice_thickness_mm  = m$slice_thickness_mm,
    inplane_mm          = m$spacing_mm[-m$slice_axis][1],
    kidney_volume_cm3   = m$kidney_volume_cm3,
    tumor_volume_cm3    = m$tumor_volume_cm3,
    tumor_surface_cm2   = m$tumor_surface_cm2,
    tumor_sphericity    = m$tumor_sphericity,
    tumor_eq_diam_cm    = m$tumor_eq_diameter_cm,
    tumor_max_area_cm2  = m$tumor_max_slice_area_cm2,
    tumor_extent_mm     = m$tumor_extent_mm,
    stringsAsFactors    = FALSE)))
  img$surface_volume_index <- img$tumor_surface_cm2/(pmax(img$tumor_volume_cm3, 1e-6)^(2/3))
  print(knitr::kable(img[, c("case_id","slice_thickness_mm","kidney_volume_cm3",
                       "tumor_volume_cm3","tumor_eq_diam_cm","tumor_sphericity")],
               digits = 3, row.names = FALSE,
               caption = "Morphometry computed from real expert reference masks, in physical units."))

  ## external validation against the recorded radiographic size
  chk <- merge(img, kits[, c("case_id", "radiographic_size_cm", "pathologic_size_cm")],
               by = "case_id")
  cat(sprintf("\nmask-derived equivalent diameter vs recorded radiographic size:\n"))
  print(round(chk[, c("tumor_eq_diam_cm", "radiographic_size_cm", "pathologic_size_cm")], 2))
  if (nrow(chk) >= 3)
    cat(sprintf("mean absolute difference = %.2f cm over %d cases\n",
                mean(abs(chk$tumor_eq_diam_cm - chk$radiographic_size_cm)), nrow(chk)))

  ## visualize one real case
  m1 <- morph_list[[1]]
  prof <- m1$tumor_profile_cm2; kprof <- m1$kidney_profile_cm2
  pf <- rbind(data.frame(s = seq_along(prof),  a = prof,  lab = "tumour"),
              data.frame(s = seq_along(kprof), a = kprof, lab = "kidney"))
  g1 <- ggplot(pf, aes(s, a, colour = lab)) +
    geom_line(linewidth = 0.8) +
    scale_colour_manual(values = bpad_pal[c(2,1)]) +
    labs(title = sprintf("cross-sectional area, %s", m1$case_id),
         x = sprintf("index along slice axis %d", m1$slice_axis),
         y = expression("area (cm"^2*")"))
  g2 <- ggplot(chk, aes(radiographic_size_cm, tumor_eq_diam_cm)) +
    geom_abline(slope = 1, intercept = 0, colour = "grey55", linetype = "dashed") +
    geom_point(colour = bpad_pal[3], size = 2.6) +
    labs(title = "mask geometry vs clinical calliper",
         x = "recorded radiographic size (cm)",
         y = expression("mask "*d[eq]*" (cm)"))
  bpad_grid(g1, g2, ncol = 2)
} else {
  cat("No segmentation volume was reachable; the geometry pipeline is exercised\n")
  cat("below on a synthetic label volume instead (Section 8.5.6).\n")
}

## 
## 
## Table: (\#tab:real-morphometry)Morphometry computed from real expert reference masks, in physical units.
## 
## |case_id    | slice_thickness_mm| kidney_volume_cm3| tumor_volume_cm3| tumor_eq_diam_cm| tumor_sphericity|
## |:----------|------------------:|-----------------:|----------------:|----------------:|----------------:|
## |case_00000 |                0.5|           387.450|            7.824|            2.463|            0.557|
## |case_00001 |                0.5|           474.834|            7.265|            2.403|            0.429|
## |case_00002 |                1.0|           542.161|           34.976|            4.058|            0.612|
## |case_00003 |                1.0|           393.211|           10.225|            2.693|            0.616|
## |case_00004 |                4.0|           439.857|           18.848|            3.302|            0.644|
## |case_00005 |                0.5|           430.187|           60.293|            4.865|            0.571|
## 
## mask-derived equivalent diameter vs recorded radiographic size:
##   tumor_eq_diam_cm radiographic_size_cm pathologic_size_cm
## 1             2.46                  2.4                2.1
## 2             2.40                  2.2                2.4
## 3             4.06                  4.5                5.4
## 4             2.69                  3.0                2.9
## 5             3.30                  3.0                3.4
## 6             4.87                  4.9                5.4
## mean absolute difference = 0.23 cm over 6 cases

Figure 8: Morphometry streamed from real KiTS19 expert segmentations. Left: an axial view of the label volume, with kidney and tumour classes. Centre: the tumour cross-sectional-area profile along the identified slice axis, whose integral is the tumour volume. Right: mask-derived equivalent spherical diameter against the independently recorded radiographic size, with the line of identity. The two agree closely, which validates the geometry pipeline against an external measurement.

8.5.5 First-order and texture features, and the discretization trap

First-order features summarize the intensity distribution inside the mask: mean, standard deviation, skewness, kurtosis, entropy, and percentiles. In CT these are in Hounsfield units and therefore physically interpretable and comparable across scanners, which is a considerable advantage over MRI, where signal intensity has no absolute scale.

Texture features summarize spatial co-occurrence. For a grey-level co-occurrence matrix $P(i,j)$ built from neighbouring voxel pairs, the standard descriptors are

\[\begin{equation} \text{contrast} = \sum_{i,j}(i-j)^2 P(i,j), \qquad \text{homogeneity} = \sum_{i,j}\frac{P(i,j)}{1 + (i-j)^2}, \qquad \text{entropy} = -\sum_{i,j}P(i,j)\log P(i,j). \tag{32} \end{equation}\]

Every one of these depends on the discretization: the number of grey levels, the intensity range used for binning, and the neighbourhood offset. Texture values computed with 16 levels are not comparable with values computed with 64 levels, and neither is comparable across a different intensity window.

glcm_features <- function(img, mask = NULL, levels = 32, offset = c(1, 0)) {
  v <- if (is.null(mask)) img else ifelse(mask, img, NA)
  rng <- range(v, na.rm = TRUE)
  q <- floor((v - rng[1])/(diff(rng) + 1e-12)*(levels - 1)) + 1
  n1 <- nrow(q); n2 <- ncol(q)
  i1 <- 1:(n1 - offset[1]); j1 <- 1:(n2 - offset[2])
  a  <- q[i1, j1, drop = FALSE]
  b  <- q[i1 + offset[1], j1 + offset[2], drop = FALSE]
  ok <- !is.na(a) & !is.na(b)
  P  <- table(factor(a[ok], levels = 1:levels), factor(b[ok], levels = 1:levels))
  P  <- (P + t(P))/sum(P + t(P))                    # symmetric, normalized
  ii <- matrix(1:levels, levels, levels)
  jj <- t(ii); dd <- (ii - jj)
  Pn <- as.numeric(P); nz <- Pn > 0
  c(contrast    = sum(dd^2*P),
    homogeneity = sum(P/(1 + dd^2)),
    energy      = sum(P^2),
    entropy     = -sum(Pn[nz]*log(Pn[nz])))
}

set.seed(5)
N <- 96
xx <- outer(1:N, 1:N, function(i, j) i)
yy <- outer(1:N, 1:N, function(i, j) j)
patch <- 60 + 25*sin(2*pi*xx/13) + 15*cos(2*pi*yy/21) + matrix(rnorm(N*N, 0, 6), N)
lev <- c(8, 16, 32, 64, 128)
tx <- t(sapply(lev, function(L) glcm_features(patch, levels = L)))
print(knitr::kable(cbind(grey_levels = lev, round(tx, 4)), row.names = FALSE,
             caption = "The same image, five discretizations. Contrast and entropy are not comparable across rows."))

## 
## 
## Table: (\#tab:texture-discretization)The same image, five discretizations. Contrast and entropy are not comparable across rows.
## 
## | grey_levels| contrast| homogeneity| energy| entropy|
## |-----------:|--------:|-----------:|------:|-------:|
## |           8|   0.7760|      0.6963| 0.0618|  2.9672|
## |          16|   2.9346|      0.4765| 0.0149|  4.3745|
## |          32|  12.0683|      0.2788| 0.0036|  5.7904|
## |          64|  49.2011|      0.1528| 0.0009|  7.1508|
## |         128| 199.0298|      0.0803| 0.0003|  8.3771|

td <- rbind(data.frame(L = lev, v = tx[, "contrast"]/tx[1, "contrast"], f = "contrast"),
            data.frame(L = lev, v = tx[, "entropy"] /tx[1, "entropy"],  f = "entropy"),
            data.frame(L = lev, v = tx[, "homogeneity"]/tx[1, "homogeneity"], f = "homogeneity"))
ggplot(td, aes(L, v, colour = f)) +
  geom_line(linewidth = 0.9) + geom_point(size = 2) +
  scale_x_log10(breaks = lev) + scale_y_log10() +
  scale_colour_manual(values = bpad_pal[1:3]) +
  labs(title = "Texture features versus grey-level discretization",
       subtitle = "values relative to the 8-level setting; the image never changed",
       x = "number of grey levels (log)", y = "relative value (log)")

Texture features are not absolute numbers. The same synthetic region is quantized at 8 to 128 grey levels; GLCM contrast and entropy change monotonically and by large factors, while the underlying image is unchanged. Any radiomic report that omits the discretization settings is unreproducible in principle, not merely in practice.

Figure 9: Texture features are not absolute numbers. The same synthetic region is quantized at 8 to 128 grey levels; GLCM contrast and entropy change monotonically and by large factors, while the underlying image is unchanged. Any radiomic report that omits the discretization settings is unreproducible in principle, not merely in practice.

8.5.6 The segmentation-uncertainty floor

Every radiomic biomarker inherits the uncertainty of the boundary that defined it. The cleanest way to quantify this, following the Monte Carlo argument of Section 8.2.9, is to perturb the mask by one voxel of morphological dilation and erosion and recompute the biomarker. One voxel is not a pessimistic perturbation: it is smaller than the disagreement between two expert annotators.

## 4-connected morphology on a binary array, base R only
dilate3d <- function(B) {
  D <- B
  D[-1, , ]  <- D[-1, , ]  | B[-dim(B)[1], , ]
  D[-dim(B)[1], , ] <- D[-dim(B)[1], , ] | B[-1, , ]
  D[, -1, ]  <- D[, -1, ]  | B[, -dim(B)[2], ]
  D[, -dim(B)[2], ] <- D[, -dim(B)[2], ] | B[, -1, ]
  D[, , -1]  <- D[, , -1]  | B[, , -dim(B)[3]]
  D[, , -dim(B)[3]] <- D[, , -dim(B)[3]] | B[, , -1]
  D
}
erode3d <- function(B) !dilate3d(!B)

geom_of <- function(B, sp) {
  n <- sum(B); V <- n*prod(sp); d <- dim(B)
  ## exposed-face surface area, same definition as in mask_morphometrics()
  S <- 0
  S <- S + sum(B[1, , ])*sp[2]*sp[3] + sum(B[d[1], , ])*sp[2]*sp[3]
  S <- S + sum(B[, 1, ])*sp[1]*sp[3] + sum(B[, d[2], ])*sp[1]*sp[3]
  S <- S + sum(B[, , 1])*sp[1]*sp[2] + sum(B[, , d[3]])*sp[1]*sp[2]
  S <- S + (sum(B[-1, , ] & !B[-d[1], , ]) + sum(B[-d[1], , ] & !B[-1, , ]))*sp[2]*sp[3]
  S <- S + (sum(B[, -1, ] & !B[, -d[2], ]) + sum(B[, -d[2], ] & !B[, -1, ]))*sp[1]*sp[3]
  S <- S + (sum(B[, , -1] & !B[, , -d[3]]) + sum(B[, , -d[3]] & !B[, , -1]))*sp[1]*sp[2]
  c(volume_cm3   = V/1000,
    surface_cm2  = S/100,
    sphericity   = (pi^(1/3)*(6*V)^(2/3))/S,
    eq_diam_cm   = 2*(3*V/(4*pi))^(1/3)/10,
    sv_index     = (S/100)/((V/1000)^(2/3)))
}

## A synthetic but realistically irregular tumour, at real KiTS-like spacing
set.seed(11)
n <- 60; sp <- c(0.92, 0.92, 1.5)
g  <- expand.grid(i = 1:n, j = 1:n, k = 1:n)
r  <- sqrt(((g$i - 30)*sp[1])^2 + ((g$j - 30)*sp[2])^2 + ((g$k - 30)*sp[3])^2)
lob <- 1 + 0.18*sin(3*atan2(g$j - 30, g$i - 30)) + 0.12*cos(4*g$k/n*pi)
B <- array(r < 11*lob, c(n, n, n))

variants <- list(eroded = erode3d(B), reference = B, dilated = dilate3d(B))
G <- t(sapply(variants, geom_of, sp = sp))
rel <- 100*(G/matrix(G["reference", ], nrow(G), ncol(G), byrow = TRUE) - 1)
pct <- round(rel[, c("volume_cm3","surface_cm2","sphericity")], 1)
colnames(pct) <- c("volume_pct_change","surface_pct_change","sphericity_pct_change")
print(knitr::kable(cbind(round(G, 3), pct),
             caption = paste("One-voxel boundary perturbation of a single mask.",
                             "The last three columns are percentage changes relative",
                             "to the reference mask.")))

## 
## 
## Table: (\#tab:segmentation-uncertainty)One-voxel boundary perturbation of a single mask. The last three columns are percentage changes relative to the reference mask.
## 
## |          | volume_cm3| surface_cm2| sphericity| eq_diam_cm| sv_index| volume_pct_change| surface_pct_change| sphericity_pct_change|
## |:---------|----------:|-----------:|----------:|----------:|--------:|-----------------:|------------------:|---------------------:|
## |eroded    |      5.521|      26.022|      0.581|      2.193|    8.330|             -23.6|              -16.0|                  -0.5|
## |reference |      7.228|      30.986|      0.583|      2.399|    8.289|               0.0|                0.0|                   0.0|
## |dilated   |      9.158|      35.906|      0.590|      2.596|    8.203|              26.7|               15.9|                   1.0|

pd <- data.frame(
  feature = rep(c("volume", "surface area", "sphericity", "equivalent diameter",
                  "surface/volume index"), each = 2),
  variant = rep(c("eroded by 1 voxel", "dilated by 1 voxel"), 5),
  pct = c(rel["eroded","volume_cm3"],  rel["dilated","volume_cm3"],
          rel["eroded","surface_cm2"], rel["dilated","surface_cm2"],
          rel["eroded","sphericity"],  rel["dilated","sphericity"],
          rel["eroded","eq_diam_cm"],  rel["dilated","eq_diam_cm"],
          rel["eroded","sv_index"],    rel["dilated","sv_index"]))
ggplot(pd, aes(reorder(feature, abs(pct)), pct, fill = variant)) +
  geom_col(position = "dodge", alpha = 0.85) +
  geom_hline(yintercept = 0, colour = "grey40") +
  coord_flip() +
  scale_fill_manual(values = bpad_pal[1:2]) +
  labs(title = "Percentage change in each biomarker from a one-voxel boundary shift",
       x = NULL, y = "change (%)")

The segmentation-uncertainty floor. A synthetic tumour mask is dilated and eroded by a single voxel, and each derived biomarker is recomputed. Volume changes by tens of percent, surface area more, and sphericity least. No downstream model can be more reproducible than these bands, so a study reporting radiomic differences smaller than the one-voxel band has measured its annotator, not its patients.

Figure 10: The segmentation-uncertainty floor. A synthetic tumour mask is dilated and eroded by a single voxel, and each derived biomarker is recomputed. Volume changes by tens of percent, surface area more, and sphericity least. No downstream model can be more reproducible than these bands, so a study reporting radiomic differences smaller than the one-voxel band has measured its annotator, not its patients.

cat(sprintf("\nvolume is the most sensitive scalar: %+.1f%% to %+.1f%%\n",
            rel["eroded","volume_cm3"], rel["dilated","volume_cm3"]))

## 
## volume is the most sensitive scalar: -23.6% to +26.7%

cat(sprintf("equivalent diameter is a cube root of volume and so moves ~1/3 as much: %+.1f%% to %+.1f%%\n",
            rel["eroded","eq_diam_cm"], rel["dilated","eq_diam_cm"]))

## equivalent diameter is a cube root of volume and so moves ~1/3 as much: -8.6% to +8.2%

From equation to reporting standard. Because $d_{\mathrm{eq}} \propto V^{1/3}$, the relative uncertainty in diameter is one third of the relative uncertainty in volume: $\sigma_d/d = \tfrac{1}{3}\sigma_V/V$. Conversely, a $\pm5\%$ diameter error becomes a $\pm15\%$ volume error. This is the delta method of Chapter 1 applied to a power law, and it explains why volumetric response criteria are more sensitive and noisier than diameter-based ones. Sphericity is the most stable of the descriptors here precisely because a dilation inflates $V^{2/3}$ and $S$ in partially compensating ways.

8.5.7 Joining imaging and clinical tables safely

A join is a scientific operation, not a clerical one. Three rules prevent the usual failures. Join on an explicit key with asserted uniqueness on both sides. Use an inner join and report how many rows were lost, because a silent left join fills imaging columns with NA and any model that ignores them then silently changes its analysis population. State the analysis population explicitly: imaging-derived features exist only for the processed subset, so imaging models are fitted on a documented subsample and their performance is not comparable with clinical-only models fitted on all 210 cases.

if (have_imaging) {
  stopifnot(!any(duplicated(img$case_id)), !any(duplicated(kits$case_id)))
  joined <- merge(kits, img, by = "case_id", all = FALSE)     # explicit inner join
  cat(sprintf("clinical rows %d, imaging rows %d, joined rows %d\n",
              nrow(kits), nrow(img), nrow(joined)))
  cat(sprintf("cases dropped from the clinical table: %d\n", nrow(kits) - nrow(joined)))
  cat("\nAnalysis population for imaging models is the joined subset only.\n")
} else {
  joined <- kits[0, ]
  cat("no imaging features available; imaging models are skipped\n")
}

## clinical rows 210, imaging rows 6, joined rows 6
## cases dropped from the clinical table: 204
## 
## Analysis population for imaging models is the joined subset only.

8.5.8 Sampling, resampling, and aliasing

Chapter 1 established the Nyquist criterion; Chapter 7 applied it to interpolation. Its consequence here is specific and often overlooked: resampling a volume changes its texture features. A texture computed at 5 mm slice thickness and one computed at 0.5 mm are measurements of different physical quantities, because partial-volume averaging has already low-pass filtered the thick-slice data. Resampling the thick data to a fine grid does not restore the lost high frequencies; it interpolates them, and interpolation is not information.

f0 <- 5; fs <- 8
tc <- seq(0, 1, length.out = 1000)
ts <- seq(0, 1, by = 1/fs)
alias_f <- abs(f0 - fs*round(f0/fs))
ad <- rbind(data.frame(t = tc, y = sin(2*pi*f0*tc), s = "true 5 Hz structure"),
            data.frame(t = tc, y = sin(2*pi*alias_f*tc), s = sprintf("apparent %.0f Hz", alias_f)))
p1 <- ggplot(ad, aes(t, y, colour = s)) +
  geom_line(linewidth = 0.8) +
  geom_point(data = data.frame(t = ts, y = sin(2*pi*f0*ts), s = "samples"),
             colour = bpad_pal[3], size = 2.4) +
  scale_colour_manual(values = c(bpad_pal[7], bpad_pal[2], bpad_pal[3])) +
  labs(title = sprintf("undersampling at %d Hz folds 5 Hz to %.0f Hz", fs, alias_f),
       x = "position (arbitrary units)", y = "intensity")

xs <- seq(0, 10, by = 1.25)
ys <- 1/(1 + exp(-(xs - 5)*2.2))
xf <- seq(0, 10, length.out = 400)
ip <- rbind(
  data.frame(x = xf, y = approx(xs, ys, xf, method = "constant", rule = 2)$y, k = "nearest"),
  data.frame(x = xf, y = approx(xs, ys, xf, method = "linear",   rule = 2)$y, k = "linear"),
  data.frame(x = xf, y = spline(xs, ys, xout = xf)$y,                          k = "spline"))
p2 <- ggplot(ip, aes(x, y, colour = k)) +
  geom_line(linewidth = 0.85) +
  geom_point(data = data.frame(x = xs, y = ys, k = "samples"),
             colour = "black", size = 2) +
  scale_colour_manual(values = c(bpad_pal[1:3], "black")) +
  labs(title = "interpolation kernels on the same coarse samples",
       x = "position (mm)", y = "intensity")
bpad_grid(p1, p2, ncol = 2)

Two sampling failures that corrupt imaging features. Left: a 5 Hz structure sampled at 8 Hz is reconstructed as a 3 Hz structure, the folding predicted by the Nyquist criterion; no downstream processing can undo this. Right: three interpolation kernels applied to the same coarse samples of a smooth boundary. Nearest-neighbour introduces staircase edges that inflate surface-area and texture features; linear interpolation smooths them; spline interpolation overshoots at sharp transitions, creating intensity values that were never measured.

Figure 11: Two sampling failures that corrupt imaging features. Left: a 5 Hz structure sampled at 8 Hz is reconstructed as a 3 Hz structure, the folding predicted by the Nyquist criterion; no downstream processing can undo this. Right: three interpolation kernels applied to the same coarse samples of a smooth boundary. Nearest-neighbour introduces staircase edges that inflate surface-area and texture features; linear interpolation smooths them; spline interpolation overshoots at sharp transitions, creating intensity values that were never measured.

cat(sprintf("spline interpolation range: [%.3f, %.3f]; original sample range: [%.3f, %.3f]\n",
            min(subset(ip, k == "spline")$y), max(subset(ip, k == "spline")$y),
            min(ys), max(ys)))

## spline interpolation range: [-0.002, 1.002]; original sample range: [0.000, 1.000]

cat("Values outside the measured range are interpolation artefacts, not measurements.\n")

## Values outside the measured range are interpolation artefacts, not measurements.

8.5.9 Feature provenance

Every feature in a defensible study carries a record of: the image series and reconstruction kernel; the mask and who produced it; the resampling grid and interpolation kernel; the intensity range and discretization; the exact feature definition and software version; and the physical units. Two features with the same name and different provenance are different variables, and the community conventions (IBSI for radiomic definitions, and the reporting checklists of Section 8.11) exist precisely because this bookkeeping is not optional.

Section 8.5 summary.

A voxel count becomes a clinical quantity only through the geometry in the file header, Eq. (29); read spacing and dimensions from the same source.
Morphometry can be streamed in one pass in $O(n_x n_y)$ memory; the slice axis should be identified from pixdim, not assumed.
Mask-derived equivalent diameter agrees with the recorded clinical size, which validates the geometry pipeline against an independent measurement.
Texture features are meaningless without their discretization settings.
A one-voxel boundary shift moves volume by tens of percent; that is the reproducibility floor of every downstream model, and $\sigma_d/d = \tfrac{1}{3}\sigma_V/V$.
Resampling cannot restore frequencies that partial-volume averaging removed.

Checkpoint 8.5. A multi-centre study harmonizes all volumes to 1 mm isotropic spacing before extracting texture features, and reports that this removed the site effect. Using Sections 8.5.5 and 8.5.8, explain why resampling alone cannot achieve this, and name one measurement that would test whether the residual site effect has actually gone.

8.6 Exploratory Structure, Preprocessing, and Representation

8.6.1 Feature engineering with physical meaning

Every derived column should be a documented, reversible transformation of a recorded quantity, with a stated reason. Four kinds of reason recur.

Multiplicative processes go on the log scale. Tumour volume grows multiplicatively (Section 8.8.7), so $\log V$ is the natural modelling scale and a coefficient becomes a proportional effect.

Categorical variables need asserted level sets. As Section 8.4.3 showed, encoding against a level that does not exist produces a constant column and no error. Every indicator below is built with %in% against levels that have been tabulated first.

Scale-invariant shape indices are preferable to raw sizes when comparing across scanners, because they cancel the length unit.

Truncation and sentinel flags are carried, never imputed away.

engineer <- function(dat) {
  out <- dat
  ## asserted level sets -- see Section 8.4.3
  stopifnot(all(out$procedure %in% c("partial_nephrectomy", "radical_nephrectomy")))
  stopifnot(all(out$access    %in% c("open", "laparoscopic", "robotic")))
  out$radical      <- as.integer(out$procedure == "radical_nephrectomy")
  out$open_access  <- as.integer(out$access    == "open")
  out$male         <- ifelse(is.na(out$gender), NA_integer_,
                             as.integer(out$gender == "male"))
  out$smoke_ever   <- ifelse(is.na(out$smoking), NA_integer_,
                             as.integer(out$smoking %in% c("previous_smoker",
                                                           "current_smoker")))
  out$ckd_yes      <- as.integer(out$ckd)
  ## multiplicative size on the log scale
  out$log_size     <- log(pmax(out$radiographic_size_cm, 0.1))
  ## acquisition covariates: genuine technical properties of each scan
  out$thin_slice   <- as.integer(out$acq_group == "thin_slice")
  out$inplane_mm   <- out$spacing_x_mm
  ## truncation and sentinel flags carried explicitly
  out$egfr_pre_truncated  <- as.integer(out$egfr_pre_flag  == ">=90")
  out$egfr_post_truncated <- as.integer(out$egfr_post_flag == ">=90")
  out
}
kits_m <- engineer(kits)

cat("engineered indicators (checking none is constant, the Section 8.4.3 failure):\n")

## engineered indicators (checking none is constant, the Section 8.4.3 failure):

ind <- c("radical","open_access","male","smoke_ever","ckd_yes","thin_slice",
         "egfr_pre_truncated","egfr_post_truncated")
print(vapply(kits_m[ind], function(z) length(unique(na.omit(z))), integer(1)))

##             radical         open_access                male          smoke_ever 
##                   2                   2                   2                   2 
##             ckd_yes          thin_slice  egfr_pre_truncated egfr_post_truncated 
##                   2                   2                   2                   2

stopifnot(all(vapply(kits_m[ind], function(z) length(unique(na.omit(z))), integer(1)) > 1))
cat("all indicators vary; none is a silently constant column\n")

## all indicators vary; none is a silently constant column

if (have_imaging) {
  joined_m <- engineer(joined)
  joined_m$log_tumor_volume <- log1p(pmax(joined_m$tumor_volume_cm3, 0))
  cat(sprintf("\nimaging-model rows: %d\n", nrow(joined_m)))
}

## 
## imaging-model rows: 6

8.6.2 Training-only imputation and standardization

Preprocessing is part of the model. The estimator below is fitted on training rows and then applied elsewhere without refitting, which is the concrete implementation of the rule stated in Section 8.3.5.

fit_preprocessor <- function(dat, vars) {
  x <- dat[, vars, drop = FALSE]
  med <- vapply(x, function(z) stats::median(z, na.rm = TRUE), numeric(1))
  med[!is.finite(med)] <- 0
  xi <- x
  for (j in seq_along(xi)) xi[[j]][is.na(xi[[j]])] <- med[j]
  ctr <- vapply(xi, mean, numeric(1))
  scl <- vapply(xi, safe_sd, numeric(1))
  structure(list(vars = vars, median = med, center = ctr, scale = scl),
            class = "bpad_prep")
}
apply_preprocessor <- function(prep, dat) {
  x <- dat[, prep$vars, drop = FALSE]
  for (j in seq_along(x)) {
    x[[j]][is.na(x[[j]])] <- prep$median[j]
    x[[j]] <- (x[[j]] - prep$center[j])/prep$scale[j]
  }
  as.matrix(x)
}

dev_m <- engineer(dev_set); ext_m <- engineer(ext_set)

## Candidate predictors. Note that `thin_slice` is CONSTANT within the
## development set by construction, since that set is defined as the thick-slice
## stratum. Including it would create exactly the silently-constant column of
## Section 8.4.3, so the guard below detects and drops zero-variance columns
## rather than trusting the analyst to remember.
pred_vars <- c("age", "bmi", "log_size", "radical", "open_access", "male",
               "smoke_ever", "ckd_yes", "thin_slice", "inplane_mm")
nuniq <- vapply(dev_m[pred_vars], function(z) length(unique(na.omit(z))), integer(1))
if (any(nuniq < 2)) {
  cat("dropped as constant within the development set:",
      paste(names(nuniq)[nuniq < 2], collapse = ", "), "\n")
  pred_vars <- pred_vars[nuniq >= 2]
}

## dropped as constant within the development set: thin_slice

prep  <- fit_preprocessor(dev_m, pred_vars)      # fitted on development data ONLY
Xdev  <- apply_preprocessor(prep, dev_m)
Xext  <- apply_preprocessor(prep, ext_m)

print(knitr::kable(data.frame(variable = prep$vars,
                        train_median = round(prep$median, 3),
                        train_center = round(prep$center, 3),
                        train_scale  = round(prep$scale, 3)),
             row.names = FALSE,
             caption = "Preprocessing parameters estimated on the development set. These exact numbers, not recomputed ones, are applied to the external set."))

## 
## 
## Table: (\#tab:preprocessor)Preprocessing parameters estimated on the development set. These exact numbers, not recomputed ones, are applied to the external set.
## 
## |variable    | train_median| train_center| train_scale|
## |:-----------|------------:|------------:|-----------:|
## |age         |       61.000|       58.826|      14.098|
## |bmi         |       29.470|       30.949|       6.564|
## |log_size    |        1.459|        1.438|       0.564|
## |radical     |        0.000|        0.373|       0.485|
## |open_access |        0.000|        0.311|       0.464|
## |male        |        1.000|        0.609|       0.490|
## |smoke_ever  |        1.000|        0.559|       0.498|
## |ckd_yes     |        0.000|        0.075|       0.263|
## |inplane_mm  |        0.779|        0.793|       0.110|

cat(sprintf("\nafter transfer, external columns have mean %.3f and sd %.3f\n",
            mean(colMeans(Xext)), mean(apply(Xext, 2, sd))))

## 
## after transfer, external columns have mean -0.148 and sd 0.985

cat("Departures from 0 and 1 are the domain shift itself, and must not be 'corrected'\n")

## Departures from 0 and 1 are the domain shift itself, and must not be 'corrected'

cat("by re-standardizing the external set, which would erase the shift being measured.\n")

## by re-standardizing the external set, which would erase the shift being measured.

8.6.3 Collinearity: diagnosis and remedies

Radiomic and clinical features are correlated by construction, because several of them measure size. Three diagnostics should be reported together: the correlation matrix, the variance inflation factor

\[\begin{equation} \mathrm{VIF}_j = \frac{1}{1 - R_j^2}, \tag{33} \end{equation}\]

where $R_j^2$ is from regressing $x_j$ on the remaining predictors, and the singular-value spectrum of the standardized design, whose condition number is the quantity that Eq. (11) squares.

Cm <- suppressWarnings(cor(Xdev, use = "pairwise.complete.obs"))
cdf <- expand.grid(a = colnames(Cm), b = colnames(Cm))
cdf$r <- as.vector(abs(Cm))
p1 <- ggplot(cdf, aes(a, b, fill = r)) +
  geom_tile(colour = "white", linewidth = 0.3) +
  scale_fill_gradient(low = "white", high = bpad_pal[1], limits = c(0, 1)) +
  theme(axis.text.x = element_text(angle = 60, hjust = 1),
        legend.position = "right", legend.title = element_text()) +
  labs(title = "absolute correlation", x = NULL, y = NULL, fill = "|r|")

svv <- svd(Xdev)$d
p2 <- ggplot(data.frame(k = seq_along(svv), s = svv), aes(k, s)) +
  geom_col(fill = bpad_pal[2], width = 0.65) +
  labs(title = sprintf("singular values (kappa = %.1f)", max(svv)/min(svv)),
       x = "index", y = expression(sigma[j]))
bpad_grid(p1, p2, ncol = 2)

Collinearity diagnostics on the real development design matrix. Left: absolute correlations among predictors. Right: the singular-value spectrum; the condition number reported below is what normal-equation inversion would square. Three remedies are available and they answer different questions: drop a representative from each correlated group (interpretability), combine them into an index (parsimony), or shrink them jointly with ridge (prediction).

Figure 12: Collinearity diagnostics on the real development design matrix. Left: absolute correlations among predictors. Right: the singular-value spectrum; the condition number reported below is what normal-equation inversion would square. Three remedies are available and they answer different questions: drop a representative from each correlated group (interpretability), combine them into an index (parsimony), or shrink them jointly with ridge (prediction).

vif <- vapply(seq_len(ncol(Xdev)), function(j) {
  r2 <- summary(lm(Xdev[, j] ~ Xdev[, -j]))$r.squared
  1/(1 - r2)
}, numeric(1))
print(knitr::kable(data.frame(variable = colnames(Xdev), VIF = round(vif, 2))[order(-vif), ],
             row.names = FALSE,
             caption = "Variance inflation factors. Values above about 5 to 10 indicate that a coefficient is poorly identified even when the model predicts well."))

## 
## 
## Table: (\#tab:collinearity)Variance inflation factors. Values above about 5 to 10 indicate that a coefficient is poorly identified even when the model predicts well.
## 
## |variable    |  VIF|
## |:-----------|----:|
## |inplane_mm  | 2.04|
## |radical     | 1.85|
## |log_size    | 1.79|
## |bmi         | 1.77|
## |male        | 1.26|
## |age         | 1.23|
## |smoke_ever  | 1.14|
## |open_access | 1.13|
## |ckd_yes     | 1.07|

cat(sprintf("kappa(X) = %.1f, so kappa(X'X) = %.1f\n",
            max(svv)/min(svv), (max(svv)/min(svv))^2))

## kappa(X) = 2.6, so kappa(X'X) = 6.8

Three remedies, three different questions. Drop a representative from each correlated group when you need an interpretable coefficient for a named variable. Combine the group into a composite index (a mean of standardized members, or its first principal component) when you need parsimony and are willing to interpret a construct rather than a measurement. Shrink jointly with ridge when prediction is the only goal and interpretation is not claimed. Choosing lasso and reporting the surviving feature as “selected” is the one option that is defensible for prediction and indefensible for interpretation, because the survivor is arbitrary within its group.

8.6.4 PCA as the SVD of the centred data

For a centred design $X_c$, the sample covariance is $\Sigma = X_c^{\mathsf T}X_c/(n-1)$. Principal components are its eigenvectors, and equivalently the right singular vectors of $X_c = U\Sigma_{\mathrm{sv}}V^{\mathsf T}$, with

\[\begin{equation} \lambda_j = \frac{\sigma_j^2}{n-1}, \tag{34} \end{equation}\]

connecting the two computations exactly. This is the same identity used in Chapter 1 for low-rank denoising, applied now to a patients-by-features matrix instead of an image.

Xc <- scale(Xdev, center = TRUE, scale = FALSE)
pc <- prcomp(Xc)
sv <- svd(Xc)
stopifnot(isTRUE(all.equal(sv$d^2/(nrow(Xc) - 1), pc$sdev^2)))
cat("PCA via eigen(cov) and via SVD agree exactly (Eq. 8.34 verified)\n")

## PCA via eigen(cov) and via SVD agree exactly (Eq. 8.34 verified)

vex <- pc$sdev^2/sum(pc$sdev^2)
p1 <- ggplot(data.frame(k = seq_along(vex), v = 100*vex, c = 100*cumsum(vex)), aes(k)) +
  geom_col(aes(y = v), fill = bpad_pal[1], alpha = 0.8) +
  geom_line(aes(y = c), colour = bpad_pal[2], linewidth = 0.9) +
  geom_point(aes(y = c), colour = bpad_pal[2], size = 1.6) +
  labs(title = "scree and cumulative variance", x = "component", y = "percent")

ld <- rbind(data.frame(v = rownames(pc$rotation), l = pc$rotation[, 1], pc = "PC1"),
            data.frame(v = rownames(pc$rotation), l = pc$rotation[, 2], pc = "PC2"))
p2 <- ggplot(ld, aes(reorder(v, abs(l)), l, fill = pc)) +
  geom_col(position = "dodge") + coord_flip() +
  scale_fill_manual(values = bpad_pal[1:2]) +
  labs(title = "loadings", x = NULL, y = "loading")

sc <- data.frame(PC1 = pc$x[, 1], PC2 = pc$x[, 2],
                 malignant = ifelse(dev_m$malignant, "malignant", "benign"))
p3 <- ggplot(sc, aes(PC1, PC2, colour = malignant)) +
  geom_point(alpha = 0.7, size = 1.9) +
  scale_colour_manual(values = c(malignant = bpad_pal[2], benign = bpad_pal[1])) +
  labs(title = "patients in PC space")
bpad_grid(p1, p2, p3, ncol = 3)

Principal component analysis of the real development design. Left: scree plot with cumulative variance explained. Centre: loadings of the first two components, showing that PC1 is dominated by the size-related and acquisition-related variables. Right: patients projected onto the first two components, coloured by malignancy. The absence of visible separation is informative: the dominant directions of variance in this cohort are not the directions that discriminate the outcome.

Figure 13: Principal component analysis of the real development design. Left: scree plot with cumulative variance explained. Centre: loadings of the first two components, showing that PC1 is dominated by the size-related and acquisition-related variables. Right: patients projected onto the first two components, coloured by malignancy. The absence of visible separation is informative: the dominant directions of variance in this cohort are not the directions that discriminate the outcome.

print(knitr::kable(data.frame(component = 1:5,
                        variance_pct = round(100*vex[1:5], 1),
                        cumulative_pct = round(100*cumsum(vex)[1:5], 1)),
             caption = "Variance explained by the leading components."))

## 
## 
## Table: (\#tab:pca-demo)Variance explained by the leading components.
## 
## | component| variance_pct| cumulative_pct|
## |---------:|------------:|--------------:|
## |         1|         20.7|           20.7|
## |         2|         18.5|           39.2|
## |         3|         14.6|           53.7|
## |         4|         13.0|           66.8|
## |         5|          9.8|           76.5|

cat(sprintf("\ncomponents needed for 90%% of variance: %d of %d\n",
            which(cumsum(vex) >= 0.90)[1], length(vex)))

## 
## components needed for 90% of variance: 7 of 9

PCA is unsupervised, and that is a limitation as well as a feature. PCA maximizes variance, not discrimination. A scanner-intensity drift or a slice-thickness effect can dominate PC1 while carrying no biological signal at all, and a genuinely predictive but low-variance direction can land in PC12 and be discarded by a “keep 90% of variance” rule. The figure above shows exactly this: the leading components are size and acquisition, and malignancy is not visibly separated. If prediction is the goal, prefer supervised dimension reduction or penalized regression to unsupervised truncation.

8.6.5 Separating margins from dependence: the empirical copula

Correlation conflates two distinct things: the shape of each variable’s marginal distribution and the structure of their dependence. Sklar’s theorem separates them. For continuous marginals $F$ and $G$, the joint distribution factors as

\[\begin{equation} H(x,y) = C\{F(x), G(y)\}, \tag{35} \end{equation}\]

where the copula $C$ carries all the dependence and is invariant under any strictly monotone transformation of either variable. The empirical copula replaces each variable by its scaled rank,

\[\begin{equation} u_i = \frac{\mathrm{rank}(x_i)}{n+1}, \qquad v_i = \frac{\mathrm{rank}(y_i)}{n+1}, \tag{36} \end{equation}\]

which is why Spearman’s and Kendall’s coefficients are unchanged by log or square-root transformation while Pearson’s is not.

pair <- kits[complete.cases(kits[, c("radiographic_size_cm","pathologic_size_cm")]), ]
x0 <- pair$radiographic_size_cm; y0 <- pair$pathologic_size_cm
tf <- list(identity = function(z) z, log = function(z) log(z), sqrt = function(z) sqrt(z))

raw <- cop <- list(); stats_tab <- NULL
for (nm in names(tf)) {
  x <- tf[[nm]](x0); y <- tf[[nm]](y0)
  u <- rank(x, ties.method = "average")/(length(x) + 1)
  v <- rank(y, ties.method = "average")/(length(y) + 1)
  raw[[nm]] <- ggplot(data.frame(x, y), aes(x, y)) +
    geom_point(alpha = 0.45, colour = bpad_pal[1], size = 1.3) +
    labs(title = sprintf("%s scale", nm), x = "radiographic", y = "pathologic")
  cop[[nm]] <- ggplot(data.frame(u, v), aes(u, v)) +
    geom_point(alpha = 0.45, colour = bpad_pal[2], size = 1.3) +
    labs(title = sprintf("copula, %s", nm), x = "u = F(x)", y = "v = G(y)")
  stats_tab <- rbind(stats_tab, data.frame(
    scale = nm,
    Pearson  = round(cor(x, y), 4),
    Spearman = round(cor(x, y, method = "spearman"), 4),
    Kendall  = round(cor(x, y, method = "kendall"), 4)))
}
bpad_grid(raw$identity, raw$log, raw$sqrt,
          cop$identity, cop$log, cop$sqrt, ncol = 3)

Marginal shape versus dependence structure, using the real radiographic and pathologic tumour sizes. Top row: the raw scatter and its log and square-root transforms look substantially different. Bottom row: the corresponding empirical copulas are identical, because a monotone transformation permutes nothing. Pearson correlation moves with the transformation; Spearman and Kendall do not.

Figure 14: Marginal shape versus dependence structure, using the real radiographic and pathologic tumour sizes. Top row: the raw scatter and its log and square-root transforms look substantially different. Bottom row: the corresponding empirical copulas are identical, because a monotone transformation permutes nothing. Pearson correlation moves with the transformation; Spearman and Kendall do not.

print(knitr::kable(stats_tab, row.names = FALSE,
             caption = "Pearson correlation depends on the scale; rank-based measures are invariant, because they are functions of the copula alone."))

## 
## 
## Table: (\#tab:copula-demo)Pearson correlation depends on the scale; rank-based measures are invariant, because they are functions of the copula alone.
## 
## |scale    | Pearson| Spearman| Kendall|
## |:--------|-------:|--------:|-------:|
## |identity |  0.9311|     0.93|   0.797|
## |log      |  0.9318|     0.93|   0.797|
## |sqrt     |  0.9382|     0.93|   0.797|

## ---------------------------------------------------------------------------
## Interactive version of the figure above. Save as app.R and run with
## shiny::runApp(). Not evaluated here so that the chapter knits statically.
## ---------------------------------------------------------------------------
library(shiny); library(ggplot2)
ui <- fluidPage(
  titlePanel("Copula invariance explorer"),
  sidebarLayout(
    sidebarPanel(
      selectInput("v1", "Variable X:", names(numeric_features), names(numeric_features)[1]),
      selectInput("v2", "Variable Y:", names(numeric_features), names(numeric_features)[2]),
      radioButtons("tr", "Monotone transform:", c("none", "log", "sqrt")),
      helpText("The scatter changes shape; the copula and the rank correlations do not.")),
    mainPanel(plotOutput("scatter", height = "290px"),
              plotOutput("copula",  height = "290px"),
              verbatimTextOutput("stats"))))
server <- function(input, output) {
  dat <- reactive({
    x <- numeric_features[[input$v1]]; y <- numeric_features[[input$v2]]
    ok <- is.finite(x) & is.finite(y); x <- x[ok]; y <- y[ok]
    if (input$tr == "log")  { x <- log(x); y <- log(y) }
    if (input$tr == "sqrt") { x <- sqrt(x); y <- sqrt(y) }
    list(x = x, y = y,
         u = rank(x)/(length(x) + 1), v = rank(y)/(length(y) + 1))
  })
  output$scatter <- renderPlot({ d <- dat()
    ggplot(data.frame(x = d$x, y = d$y), aes(x, y)) + geom_point(alpha = 0.6) +
      theme_minimal() })
  output$copula <- renderPlot({ d <- dat()
    ggplot(data.frame(u = d$u, v = d$v), aes(u, v)) + geom_point(alpha = 0.6) +
      theme_minimal() })
  output$stats <- renderPrint({ d <- dat()
    c(Pearson  = cor(d$x, d$y),
      Spearman = cor(d$x, d$y, method = "spearman"),
      Kendall  = cor(d$x, d$y, method = "kendall")) })
}
shinyApp(ui, server)

8.6.6 Acquisition shift, measured

Section 8.4.5 defined the external set by slice thickness. Here we quantify the shift, which is the empirical content of the $\Delta_{\mathrm{shift}}$ term in Eq. (5).

smd <- function(a, b) {
  a <- a[is.finite(a)]; b <- b[is.finite(b)]
  (mean(a) - mean(b))/sqrt((var(a) + var(b))/2)
}
shift_vars <- c("age","bmi","radiographic_size_cm","spacing_x_mm","spacing_z_mm",
                "egfr_pre","radical","open_access","male","smoke_ever")
sd_tab <- data.frame(
  variable = shift_vars,
  smd = vapply(shift_vars, function(v) smd(dev_m[[v]], ext_m[[v]]), numeric(1)))
sd_tab$abs <- abs(sd_tab$smd)
print(knitr::kable(sd_tab[order(-sd_tab$abs), c("variable","smd")], row.names = FALSE,
             digits = 3,
             caption = "Standardized mean differences between development (thick-slice) and external (thin-slice) strata."))

## 
## 
## Table: (\#tab:acquisition-shift)Standardized mean differences between development (thick-slice) and external (thin-slice) strata.
## 
## |variable             |    smd|
## |:--------------------|------:|
## |spacing_z_mm         |  3.654|
## |radical              |  0.376|
## |smoke_ever           |  0.303|
## |open_access          |  0.244|
## |radiographic_size_cm |  0.204|
## |male                 |  0.198|
## |bmi                  | -0.150|
## |egfr_pre             |  0.139|
## |age                  |  0.138|
## |spacing_x_mm         | -0.135|

p1 <- ggplot(sd_tab, aes(reorder(variable, abs), 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 = "standardized mean difference", x = NULL, y = "SMD (dev - ext)")
p2 <- ggplot(rbind(data.frame(v = dev_m$radiographic_size_cm, g = "development (thick)"),
                   data.frame(v = ext_m$radiographic_size_cm, g = "external (thin)")),
             aes(v, fill = g)) +
  geom_density(alpha = 0.45) +
  scale_fill_manual(values = bpad_pal[1:2]) +
  labs(title = "tumour size by acquisition stratum",
       x = "radiographic size (cm)", y = "density")
bpad_grid(p1, p2, ncol = 2)

Acquisition shift in a single public cohort. Left: standardized mean differences between the thick-slice development set and the thin-slice external set; bars beyond the dashed lines at 0.2 exceed the conventional threshold for a non-negligible imbalance. Right: tumour size distributions in the two strata. The shift is not confined to geometry, so a performance drop on the external set cannot be attributed to slice thickness alone.

Figure 15: Acquisition shift in a single public cohort. Left: standardized mean differences between the thick-slice development set and the thin-slice external set; bars beyond the dashed lines at 0.2 exceed the conventional threshold for a non-negligible imbalance. Right: tumour size distributions in the two strata. The shift is not confined to geometry, so a performance drop on the external set cannot be attributed to slice thickness alone.

cat(sprintf("\nvariables with |SMD| > 0.2: %s\n",
            paste(sd_tab$variable[sd_tab$abs > 0.2], collapse = ", ")))

## 
## variables with |SMD| > 0.2: radiographic_size_cm, spacing_z_mm, radical, open_access, smoke_ever

Section 8.6 summary.

Engineer features for a stated physical reason, assert categorical level sets, and carry truncation flags rather than imputing them.
Fit imputation and standardization on training data only, and do not re-standardize the external set: that would erase the shift you are trying to measure.
Report correlation, VIF, and the singular-value spectrum together; the three remedies for collinearity answer three different questions.
PCA is the SVD of the centred data, Eq. (34), and it is unsupervised: high variance is not the same as high information.
The copula separates marginal shape from dependence, which is why rank correlations are transformation-invariant.
Acquisition shift in this cohort is technical and clinical, so external performance drops cannot be attributed to geometry alone.

Checkpoint 8.6. A colleague standardizes the combined development and external data together “so the scales match”, then splits and reports external performance. Identify the specific quantity that has leaked, predict the direction of the bias, and state the one-line fix.

8.7 Supervised Learning: Classification and Regression

8.7.1 What to report, and why three numbers are needed

A prediction model must be assessed on three logically independent axes, and reporting only one is the most common reporting failure in the field.

Discrimination asks whether the model ranks patients correctly. For a binary outcome the AUC equals the probability that a randomly chosen case outranks a randomly chosen control,

\[\begin{equation} \mathrm{AUC} = P\{\widehat p(X_1) > \widehat p(X_0) \mid Y_1 = 1, Y_0 = 0\}, \tag{37} \end{equation}\]

which is exactly the Mann–Whitney statistic implemented in auc_rank().

Calibration asks whether the predicted probabilities are numerically right. The Brier score

\[\begin{equation} \mathrm{BS} = \frac{1}{n}\sum_{i=1}^{n}\big(\widehat p_i - y_i\big)^2 \tag{38} \end{equation}\]

decomposes into calibration, refinement, and irreducible components, and it can be compared with the score of a model that always predicts the base rate, $\mathrm{BS}_0 = \bar y(1-\bar y)$. A model with excellent AUC and poor calibration ranks well and lies about absolute risk, which is worse than useless when a threshold is applied.

Clinical utility asks whether acting on the model helps. At a decision threshold $t$, net benefit is

\[\begin{equation} \mathrm{NB}(t) = \frac{\mathrm{TP}(t)}{n} - \frac{\mathrm{FP}(t)}{n}\cdot\frac{t}{1-t}, \tag{39} \end{equation}\]

where $t/(1-t)$ is the exchange rate between a false positive and a true positive implied by choosing that threshold. Section 8.7.8 computes it.

8.7.2 Malignancy classification on the real cohort

The intended-use statement of Section 8.1.4 restricts predictors to information available preoperatively. That excludes pathologic size, ISUP grade, T stage, ischaemia time, and blood loss, all of which are present in the file and all of which would improve apparent performance while making the model unusable for its stated purpose.

model_vars <- c("log_size", "age", "bmi", "radical", "open_access", "male", "smoke_ever")
d <- dev_m[complete.cases(dev_m[, c(model_vars, "malignant")]), ]
cat(sprintf("development rows = %d; benign = %d; malignant = %d\n",
            nrow(d), sum(!d$malignant), sum(d$malignant)))

## development rows = 161; benign = 13; malignant = 148

cat(sprintf("events per candidate predictor = %.1f (Section 8.3.7 rule of thumb: 10-20)\n",
            sum(!d$malignant)/length(model_vars)))

## events per candidate predictor = 1.9 (Section 8.3.7 rule of thumb: 10-20)

fml  <- as.formula(paste("malignant ~", paste(model_vars, collapse = " + ")))
fml1 <- malignant ~ log_size

folds <- make_folds(nrow(d), K = 5, strata = d$malignant, seed = 11)
oof_full <- oof_size <- rep(NA_real_, nrow(d))
for (k in seq_len(5)) {
  tr <- d[folds != k, ]; te <- d[folds == k, ]
  oof_full[folds == k] <- predict(suppressWarnings(glm(fml,  tr, family = binomial())),
                                  te, type = "response")
  oof_size[folds == k] <- predict(suppressWarnings(glm(fml1, tr, family = binomial())),
                                  te, type = "response")
}

base_rate <- mean(d$malignant)
perf <- data.frame(
  model = c("clinical model (7 predictors)", "tumour size alone", "always predict base rate"),
  AUC   = round(c(auc_rank(oof_full, d$malignant),
                  auc_rank(oof_size, d$malignant), 0.500), 3),
  Brier = round(c(brier(oof_full, d$malignant), brier(oof_size, d$malignant),
                  brier(rep(base_rate, nrow(d)), d$malignant)), 4))
perf$skill_vs_base <- round(1 - perf$Brier/perf$Brier[3], 3)
print(knitr::kable(perf, row.names = FALSE,
             caption = "Out-of-fold discrimination and calibration. The Brier skill score relative to the base rate is close to zero: in absolute-probability terms this model barely improves on quoting the prevalence."))

## 
## 
## Table: (\#tab:malignancy-model)Out-of-fold discrimination and calibration. The Brier skill score relative to the base rate is close to zero: in absolute-probability terms this model barely improves on quoting the prevalence.
## 
## |model                         |   AUC|  Brier| skill_vs_base|
## |:-----------------------------|-----:|------:|-------------:|
## |clinical model (7 predictors) | 0.621| 0.0750|        -0.011|
## |tumour size alone             | 0.684| 0.0717|         0.034|
## |always predict base rate      | 0.500| 0.0742|         0.000|

## external test, touched exactly once
e <- ext_m[complete.cases(ext_m[, c(model_vars, "malignant")]), ]
fit_final <- suppressWarnings(glm(fml, d, family = binomial()))
p_ext <- predict(fit_final, e, type = "response")
cat(sprintf("\nEXTERNAL (thin-slice) n = %d, benign = %d\n", nrow(e), sum(!e$malignant)))

## 
## EXTERNAL (thin-slice) n = 49, benign = 5

cat(sprintf("external AUC = %.3f   (development out-of-fold AUC = %.3f)\n",
            auc_rank(p_ext, e$malignant), auc_rank(oof_full, d$malignant)))

## external AUC = 0.727   (development out-of-fold AUC = 0.621)

cat(sprintf("external Brier = %.4f  (external base-rate Brier = %.4f)\n",
            brier(p_ext, e$malignant),
            brier(rep(mean(e$malignant), nrow(e)), e$malignant)))

## external Brier = 0.0800  (external base-rate Brier = 0.0916)

## ---- figures ----
roc_pts <- function(p, y) {
  o <- order(-p); y <- as.numeric(y)[o]
  data.frame(fpr = cumsum(1 - y)/sum(1 - y), tpr = cumsum(y)/sum(y))
}
r1 <- roc_pts(oof_full, d$malignant); r1$m <- "clinical model"
r2 <- roc_pts(oof_size, d$malignant); r2$m <- "size alone"
g1 <- ggplot(rbind(r1, r2), aes(fpr, tpr, colour = m)) +
  geom_abline(slope = 1, intercept = 0, colour = "grey60", linetype = "dashed") +
  geom_line(linewidth = 0.9) +
  scale_colour_manual(values = bpad_pal[1:2]) +
  labs(title = "ROC (out-of-fold)", x = "false positive rate", y = "true positive rate")

bins <- cut(oof_full, breaks = quantile(oof_full, seq(0, 1, 0.2)), include.lowest = TRUE)
cal <- aggregate(cbind(pred = oof_full, obs = as.numeric(d$malignant)) ~ bins, FUN = mean)
g2 <- ggplot(cal, aes(pred, obs)) +
  geom_abline(slope = 1, intercept = 0, colour = "grey60", linetype = "dashed") +
  geom_hline(yintercept = base_rate, colour = bpad_pal[3], linetype = "dotted") +
  geom_point(colour = bpad_pal[1], size = 2.6) + geom_line(colour = bpad_pal[1]) +
  coord_cartesian(xlim = c(0.6, 1), ylim = c(0.6, 1)) +
  labs(title = "calibration by quintile of risk",
       x = "mean predicted probability", y = "observed proportion")

g3 <- ggplot(data.frame(p = oof_full,
                        y = ifelse(d$malignant, "malignant", "benign")),
             aes(p, fill = y)) +
  geom_density(alpha = 0.5) +
  scale_fill_manual(values = c(malignant = bpad_pal[2], benign = bpad_pal[1])) +
  labs(title = "predicted risk by true class", x = "predicted probability", y = "density")
bpad_grid(g1, g2, g3, ncol = 3)

Out-of-fold malignancy prediction in the real cohort. Left: ROC curves for the clinical model and for tumour size alone; the two are close, indicating that most of the signal is size. Centre: calibration plot with the base rate marked; the model is close to the diagonal but spans a narrow range of predicted risk. Right: predicted-probability distributions by true class, showing extensive overlap. With 13 benign masses in the development set, this is close to the ceiling the labels can support.

Figure 16: Out-of-fold malignancy prediction in the real cohort. Left: ROC curves for the clinical model and for tumour size alone; the two are close, indicating that most of the signal is size. Centre: calibration plot with the base rate marked; the model is close to the diagonal but spans a narrow range of predicted risk. Right: predicted-probability distributions by true class, showing extensive overlap. With 13 benign masses in the development set, this is close to the ceiling the labels can support.

Read this result honestly. The out-of-fold AUC is modest and the Brier skill score relative to the base rate is essentially zero, meaning the model adds almost nothing to simply quoting the prevalence. The external AUC happens to come out higher than the development estimate, which is not evidence of good transfer: with five benign masses in the external set, that estimate has a confidence interval spanning most of the unit interval, and the sensible reading is that development and external performance are indistinguishable from each other and both close to uninformative. Note also that tumour size alone outperforms the seven-predictor model out of fold, which is what over-parameterization looks like at 1.9 events per predictor. This is not a failure of the code. With 13 benign masses among 161 development cases, a prevalence of 0.92, and preoperative predictors restricted to size and demographics, there is very little discriminable signal available. The correct conclusion is that this cohort cannot support a preoperative malignancy classifier of clinical value, and the correct next step is a better cohort (an unselected series of incidental masses) or better features (enhancement kinetics, texture on real intensity volumes), not a better optimizer.

8.7.3 Coefficients and their interpretation

ci <- suppressWarnings(confint.default(fit_final))
co <- data.frame(term = names(coef(fit_final)), est = coef(fit_final),
                 lo = ci[, 1], hi = ci[, 2])
co <- co[co$term != "(Intercept)", ]
co[, c("or","or_lo","or_hi")] <- exp(co[, c("est","lo","hi")])
print(knitr::kable(data.frame(term = co$term,
                        odds_ratio = round(co$or, 3),
                        CI_low = round(co$or_lo, 3),
                        CI_high = round(co$or_hi, 3),
                        crosses_1 = ifelse(co$or_lo < 1 & co$or_hi > 1, "yes", "no")),
             row.names = FALSE,
             caption = "Adjusted odds ratios. Note that `log_size` is on the log scale, so its odds ratio is a multiplicative effect per unit change in log size, that is per e-fold change in diameter."))

## 
## 
## Table: (\#tab:coefficients)Adjusted odds ratios. Note that `log_size` is on the log scale, so its odds ratio is a multiplicative effect per unit change in log size, that is per e-fold change in diameter.
## 
## |term        | odds_ratio| CI_low| CI_high|crosses_1 |
## |:-----------|----------:|------:|-------:|:---------|
## |log_size    |     10.723|  1.806|  63.671|no        |
## |age         |      0.987|  0.940|   1.036|yes       |
## |bmi         |      1.003|  0.917|   1.098|yes       |
## |radical     |      0.225|  0.035|   1.443|yes       |
## |open_access |      1.703|  0.332|   8.740|yes       |
## |male        |      2.390|  0.707|   8.078|yes       |
## |smoke_ever  |      0.980|  0.284|   3.375|yes       |

ggplot(co, aes(or, reorder(term, or))) +
  geom_vline(xintercept = 1, linetype = "dashed", colour = "grey50") +
  geom_errorbarh(aes(xmin = or_lo, xmax = or_hi), height = 0.18, colour = bpad_pal[1]) +
  geom_point(size = 2.4, colour = bpad_pal[1]) +
  scale_x_log10() +
  labs(title = "Adjusted odds ratios (log scale)", x = "odds ratio", y = NULL)

Odds ratios with 95 percent confidence intervals from the final logistic model, on a logarithmic axis. Only tumour size is clearly separated from the null value of one; every other interval is wide and crosses it, which is what 13 events across 7 predictors buys.

Figure 17: Odds ratios with 95 percent confidence intervals from the final logistic model, on a logarithmic axis. Only tumour size is clearly separated from the null value of one; every other interval is wide and crosses it, which is what 13 events across 7 predictors buys.

8.7.4 k-nearest neighbours and the bias-variance dial

kNN makes the bias-variance trade-off of Eq. (17) visible in a single tuning parameter. Small $k$ gives a flexible, low-bias, high-variance rule; large $k$ approaches the base rate, which is maximum bias and minimum variance.

knn_prob <- function(Xtr, ytr, Xte, k) {
  vapply(seq_len(nrow(Xte)), function(i) {
    d2 <- colSums((t(Xtr) - Xte[i, ])^2)
    mean(ytr[order(d2)[seq_len(k)]])
  }, numeric(1))
}
Xd <- apply_preprocessor(prep, d); yd <- as.numeric(d$malignant)
ks <- c(1, 3, 5, 7, 11, 15, 21, 31, 45, 65)
knn_auc <- vapply(ks, function(k) {
  oof <- rep(NA_real_, nrow(Xd))
  for (f in seq_len(5)) {
    tr <- folds != f
    oof[!tr] <- knn_prob(Xd[tr, , drop = FALSE], yd[tr], Xd[!tr, , drop = FALSE],
                         min(k, sum(tr)))
  }
  auc_rank(oof, d$malignant)
}, numeric(1))
print(knitr::kable(data.frame(k = ks, out_of_fold_AUC = round(knn_auc, 3)), row.names = FALSE,
             caption = "kNN out-of-fold AUC against the neighbourhood size."))

## 
## 
## Table: (\#tab:knn-demo)kNN out-of-fold AUC against the neighbourhood size.
## 
## |  k| out_of_fold_AUC|
## |--:|---------------:|
## |  1|           0.498|
## |  3|           0.509|
## |  5|           0.590|
## |  7|           0.646|
## | 11|           0.662|
## | 15|           0.693|
## | 21|           0.683|
## | 31|           0.558|
## | 45|           0.547|
## | 65|           0.478|

ggplot(data.frame(k = ks, auc = knn_auc), aes(k, auc)) +
  geom_hline(yintercept = 0.5, linetype = "dashed", colour = "grey55") +
  geom_line(colour = bpad_pal[1], linewidth = 0.9) + geom_point(size = 2) +
  scale_x_log10(breaks = ks) +
  labs(title = "kNN: the bias-variance trade-off as one tuning parameter",
       subtitle = "small k: high variance; large k: high bias, converging to the base rate",
       x = "k (log scale)", y = "out-of-fold AUC")

k-nearest neighbours on the standardized real development features. Out-of-fold AUC rises and then falls as k increases, tracing the bias-variance trade-off directly. At the largest k the classifier converges toward predicting the base rate for everyone, which yields chance discrimination. The optimum must be chosen inside resampling, never by reading this curve and reporting its peak.

Figure 18: k-nearest neighbours on the standardized real development features. Out-of-fold AUC rises and then falls as k increases, tracing the bias-variance trade-off directly. At the largest k the classifier converges toward predicting the base rate for everyone, which yields chance discrimination. The optimum must be chosen inside resampling, never by reading this curve and reporting its peak.

cat(sprintf("best k = %d at AUC %.3f; k = 1 gives %.3f; k = %d gives %.3f\n",
            ks[which.max(knn_auc)], max(knn_auc), knn_auc[1],
            ks[length(ks)], knn_auc[length(knn_auc)]))

## best k = 15 at AUC 0.693; k = 1 gives 0.498; k = 65 gives 0.478

8.7.5 Comparing model families

Different function classes are compared under an identical resampling scheme, identical preprocessing fitted inside folds, and identical predictors. Anything else compares pipelines rather than models.

oof_of <- function(fit_fun, pred_fun) {
  o <- rep(NA_real_, nrow(d))
  for (k in seq_len(5)) {
    tr <- folds != k
    pp <- fit_preprocessor(d[tr, ], pred_vars)         # preprocessing INSIDE the fold
    Xtr <- apply_preprocessor(pp, d[tr, ]); Xte <- apply_preprocessor(pp, d[!tr, ])
    m <- fit_fun(Xtr, yd[tr])
    o[!tr] <- pred_fun(m, Xte)
  }
  o
}
res <- list()
res[["logistic"]] <- oof_of(
  function(X, y) suppressWarnings(glm(y ~ ., data = data.frame(X, y = y), family = binomial())),
  function(m, X) predict(m, data.frame(X), type = "response"))
res[["kNN (k=11)"]] <- oof_of(function(X, y) list(X = X, y = y),
                              function(m, X) knn_prob(m$X, m$y, X, 11))
if (has_pkg("glmnet")) {
  ## glmnet warns that a class has fewer than 8 observations. That warning is a
  ## substantive fact about this cohort, not a nuisance, so it is reported here
  ## rather than silently suppressed.
  cat(sprintf("note: the minority (benign) class has %d cases in total, so an\n",
              sum(!d$malignant)))
  cat("inner training fold contains roughly 10. Penalized fitting is being asked\n")
  cat("to estimate a penalty from about ten events.\n")
  res[["ridge logistic"]] <- suppressWarnings(oof_of(
    function(X, y) glmnet::cv.glmnet(X, y, family = "binomial", alpha = 0, nfolds = 5),
    function(m, X) as.numeric(predict(m, X, s = "lambda.min", type = "response"))))
}

## note: the minority (benign) class has 13 cases in total, so an
## inner training fold contains roughly 10. Penalized fitting is being asked
## to estimate a penalty from about ten events.

if (has_pkg("randomForest")) res[["random forest"]] <- oof_of(
  function(X, y) randomForest::randomForest(X, factor(y), ntree = 500),
  function(m, X) predict(m, X, type = "prob")[, 2])
if (has_pkg("e1071")) res[["SVM (RBF)"]] <- oof_of(
  function(X, y) e1071::svm(X, factor(y), probability = TRUE, kernel = "radial"),
  function(m, X) attr(predict(m, X, probability = TRUE), "probabilities")[, "1"])

cmp <- data.frame(
  model = names(res),
  AUC   = round(vapply(res, function(p) auc_rank(p, d$malignant), numeric(1)), 3),
  Brier = round(vapply(res, function(p) brier(p, d$malignant), numeric(1)), 4))
cmp$skill <- round(1 - cmp$Brier/brier(rep(base_rate, nrow(d)), d$malignant), 3)
print(knitr::kable(cmp[order(-cmp$AUC), ], row.names = FALSE,
             caption = "Out-of-fold performance with preprocessing refitted inside every fold."))

## 
## 
## Table: (\#tab:classifier-comparison)Out-of-fold performance with preprocessing refitted inside every fold.
## 
## |model          |   AUC|  Brier|  skill|
## |:--------------|-----:|------:|------:|
## |kNN (k=11)     | 0.665| 0.0737|  0.007|
## |logistic       | 0.535| 0.0806| -0.086|
## |random forest  | 0.485| 0.0835| -0.125|
## |ridge logistic | 0.432| 0.0757| -0.020|
## |SVM (RBF)      | 0.350| 0.0771| -0.039|

## fold-to-fold variability, which is what makes the ranking above unreliable
per_fold <- do.call(rbind, lapply(names(res), function(nm)
  data.frame(model = nm,
             fold = seq_len(5),
             auc = vapply(seq_len(5), function(k)
               auc_rank(res[[nm]][folds == k], d$malignant[folds == k]), numeric(1)))))
ggplot(per_fold, aes(reorder(model, auc, FUN = median), auc)) +
  geom_hline(yintercept = 0.5, linetype = "dashed", colour = "grey55") +
  geom_boxplot(fill = bpad_pal[1], alpha = 0.35, outlier.shape = NA) +
  geom_jitter(width = 0.12, size = 1.8, colour = bpad_pal[2]) +
  coord_flip() +
  labs(title = "Per-fold AUC: between-model differences are smaller than between-fold spread",
       x = NULL, y = "AUC within fold")

Five model families under one resampling scheme on the real development set. Per-fold AUCs overlap heavily and straddle chance, so the between-family ordering is not interpretable: with 13 benign cases the fold-to-fold spread swamps any difference between function classes. Reporting the maximum over families without accounting for that selection is the multiplicity problem of Section 8.3.6 in miniature.

Figure 19: Five model families under one resampling scheme on the real development set. Per-fold AUCs overlap heavily and straddle chance, so the between-family ordering is not interpretable: with 13 benign cases the fold-to-fold spread swamps any difference between function classes. Reporting the maximum over families without accounting for that selection is the multiplicity problem of Section 8.3.6 in miniature.

8.7.6 Regression: predicting postoperative renal-function decline

This is the outcome the cohort actually supports. The target is the fall in estimated glomerular filtration rate,

\[\Delta\mathrm{eGFR} = \mathrm{eGFR}_{\mathrm{pre}} - \mathrm{eGFR}_{\mathrm{post}},\]

which is continuous, clinically meaningful, and mechanistically linked to the operation performed: removing a whole kidney removes roughly half the filtering parenchyma.

rvars <- c("egfr_pre", "age", "radical", "log_size", "bmi", "egfr_pre_truncated")
r <- kits_m[complete.cases(kits_m[, c(rvars, "egfr_drop")]), ]
cat(sprintf("eGFR analysis population: %d of %d cases (%.0f%%)\n",
            nrow(r), nrow(kits_m), 100*nrow(r)/nrow(kits_m)))

## eGFR analysis population: 93 of 210 cases (44%)

cat("The restriction is itself a finding: postoperative eGFR is recorded for fewer\n")

## The restriction is itself a finding: postoperative eGFR is recorded for fewer

cat("than half the cohort, so this analysis population is not the whole cohort.\n\n")

## than half the cohort, so this analysis population is not the whole cohort.

fit_r <- lm(egfr_drop ~ egfr_pre + age + radical + log_size + bmi + egfr_pre_truncated,
            data = r)
print(knitr::kable(round(summary(fit_r)$coefficients, 4),
             caption = "Linear model for postoperative eGFR decline (mL/min/1.73m2)."))

## 
## 
## Table: (\#tab:egfr-regression)Linear model for postoperative eGFR decline (mL/min/1.73m2).
## 
## |            | Estimate| Std. Error| t value| Pr(>&#124;t&#124;)|
## |:-----------|--------:|----------:|-------:|------------------:|
## |(Intercept) | -43.1083|    17.0443| -2.5292|             0.0132|
## |egfr_pre    |   0.5507|     0.1146|  4.8061|             0.0000|
## |age         |   0.1170|     0.1415|  0.8270|             0.4105|
## |radical     |  18.9697|     3.9832|  4.7624|             0.0000|
## |log_size    |  -2.9349|     3.0765| -0.9540|             0.3427|
## |bmi         |   0.1959|     0.2449|  0.7998|             0.4260|

cat(sprintf("\nR-squared = %.3f, adjusted = %.3f, residual SD = %.2f\n",
            summary(fit_r)$r.squared, summary(fit_r)$adj.r.squared,
            summary(fit_r)$sigma))

## 
## R-squared = 0.406, adjusted = 0.372, residual SD = 13.36

f2 <- make_folds(nrow(r), 5, seed = 3)
oof_r <- rep(NA_real_, nrow(r))
for (k in seq_len(5)) {
  tr <- f2 != k
  oof_r[!tr] <- predict(lm(egfr_drop ~ egfr_pre + age + radical + log_size + bmi +
                             egfr_pre_truncated, data = r[tr, ]), r[!tr, ])
}
cat(sprintf("out-of-fold RMSE = %.2f;  null RMSE (predict the mean) = %.2f\n",
            sqrt(mean((oof_r - r$egfr_drop)^2)), sd(r$egfr_drop)))

## out-of-fold RMSE = 14.29;  null RMSE (predict the mean) = 16.86

cat(sprintf("out-of-fold R-squared = %.3f\n",
            1 - mean((oof_r - r$egfr_drop)^2)/var(r$egfr_drop)))

## out-of-fold R-squared = 0.281

q1 <- ggplot(data.frame(obs = r$egfr_drop, pred = oof_r), aes(pred, obs)) +
  geom_abline(slope = 1, intercept = 0, colour = "grey55", linetype = "dashed") +
  geom_point(alpha = 0.6, colour = bpad_pal[1], size = 1.8) +
  labs(title = "out-of-fold prediction", x = "predicted decline", y = "observed decline")
q2 <- ggplot(data.frame(f = fitted(fit_r), e = residuals(fit_r)), aes(f, e)) +
  geom_hline(yintercept = 0, colour = "grey55") +
  geom_point(alpha = 0.6, colour = bpad_pal[2], size = 1.8) +
  geom_smooth(method = "loess", se = FALSE, colour = bpad_pal[3], linewidth = 0.8) +
  labs(title = "residuals versus fitted", x = "fitted", y = "residual")
q3 <- ggplot(r, aes(factor(radical, labels = c("partial", "radical")), egfr_drop,
                    fill = factor(radical))) +
  geom_boxplot(alpha = 0.5, show.legend = FALSE) +
  scale_fill_manual(values = bpad_pal[1:2]) +
  labs(title = "decline by operative extent", x = NULL, y = "eGFR decline")
bpad_grid(q1, q2, q3, ncol = 3)

Predicting postoperative renal-function decline on the real cohort. Left: observed against out-of-fold predicted decline, with the line of identity. Centre: residuals against fitted values, checking for heteroscedasticity. Right: decline by operative extent, showing the effect that dominates the model. Unlike the malignancy classifier, this model has real predictive skill, because the outcome is mechanistically determined by a recorded intervention.

Figure 20: Predicting postoperative renal-function decline on the real cohort. Left: observed against out-of-fold predicted decline, with the line of identity. Centre: residuals against fitted values, checking for heteroscedasticity. Right: decline by operative extent, showing the effect that dominates the model. Unlike the malignancy classifier, this model has real predictive skill, because the outcome is mechanistically determined by a recorded intervention.

From coefficient to physiology. The radical coefficient is the additional eGFR decline attributable to removing the whole kidney rather than part of it, adjusted for baseline function, age, tumour size, and BMI. Its magnitude is close to what a mechanistic argument predicts: if two kidneys contribute roughly equally to filtration and the contralateral kidney hypertrophies only partially, a nephrectomy should cost on the order of a third to a half of baseline eGFR. The egfr_pre coefficient, being well below one, reflects both regression to the mean and the ceiling imposed by the ">=90" truncation of Section 8.4.3; the truncation indicator is carried in the model precisely so that this artefact is visible rather than absorbed into the other coefficients.

8.7.7 Nested cross-validation

When a hyperparameter is tuned, the tuning must sit inside the performance estimate. The comparison below applies the Section 8.3.6 lesson to a real tuning decision rather than to simulated noise.

if (has_pkg("glmnet")) {
  set.seed(4)
  Xall <- apply_preprocessor(fit_preprocessor(d, pred_vars), d)
  ## (a) OPTIMISTIC: choose lambda on all development data, then cross-validate
  cv_all  <- suppressWarnings(
    glmnet::cv.glmnet(Xall, yd, family = "binomial", alpha = 0, nfolds = 5))
  lam_all <- cv_all$lambda.min
  oof_opt <- rep(NA_real_, nrow(d))
  for (k in seq_len(5)) {
    tr <- folds != k
    m  <- suppressWarnings(glmnet::glmnet(Xall[tr, ], yd[tr], family = "binomial",
                                          alpha = 0, lambda = lam_all))
    oof_opt[!tr] <- as.numeric(predict(m, Xall[!tr, ], type = "response"))
  }
  ## (b) HONEST: choose lambda inside each outer training fold
  oof_nest <- rep(NA_real_, nrow(d)); lam_used <- numeric(5)
  for (k in seq_len(5)) {
    tr <- folds != k
    inner <- suppressWarnings(glmnet::cv.glmnet(Xall[tr, ], yd[tr],
                                 family = "binomial", alpha = 0, nfolds = 5))
    lam_used[k] <- inner$lambda.min
    oof_nest[!tr] <- as.numeric(predict(inner, Xall[!tr, ], s = "lambda.min",
                                        type = "response"))
  }
  print(knitr::kable(data.frame(
    procedure = c("lambda chosen on all development data", "lambda chosen inside each fold"),
    AUC   = round(c(auc_rank(oof_opt, d$malignant), auc_rank(oof_nest, d$malignant)), 3),
    Brier = round(c(brier(oof_opt, d$malignant), brier(oof_nest, d$malignant)), 4)),
    row.names = FALSE,
    caption = "Nested versus non-nested tuning of the ridge penalty on the real cohort."))
  cat(sprintf("lambda selected on all data: %.4f\n", lam_all))
  cat(sprintf("lambda selected per fold  : %s\n",
              paste(sprintf("%.4f", lam_used), collapse = ", ")))
  cat("\nThe optimism here is small because one scalar hyperparameter leaks little.\n")
  cat("It grows with the number of tuned choices, and Section 8.3.6 shows the extreme.\n")
} else {
  cat("glmnet unavailable; nested-tuning comparison skipped\n")
}

## 
## 
## Table: (\#tab:nested-cv)Nested versus non-nested tuning of the ridge penalty on the real cohort.
## 
## |procedure                             |   AUC|  Brier|
## |:-------------------------------------|-----:|------:|
## |lambda chosen on all development data | 0.582| 0.0739|
## |lambda chosen inside each fold        | 0.439| 0.0762|
## lambda selected on all data: 0.3150
## lambda selected per fold  : 0.7622, 48.5689, 0.0517, 52.3252, 53.4274
## 
## The optimism here is small because one scalar hyperparameter leaks little.
## It grows with the number of tuned choices, and Section 8.3.6 shows the extreme.

8.7.8 Class imbalance, precision-recall, and decision thresholds

With a prevalence of 0.92 the benign class is the rare one. ROC treats the two errors symmetrically; precision-recall does not, and is the more informative curve when the quantity of interest concerns the rare class.

p_benign <- 1 - oof_full; y_benign <- !d$malignant
pr_pts <- function(p, y) {
  o <- order(-p); y <- as.numeric(y)[o]
  tp <- cumsum(y); fp <- cumsum(1 - y)
  data.frame(recall = tp/sum(y), precision = tp/(tp + fp))
}
rp <- roc_pts(p_benign, y_benign); pp <- pr_pts(p_benign, y_benign)
cat(sprintf("identifying the RARE benign class: AUC = %.3f, prevalence = %.3f\n",
            auc_rank(p_benign, y_benign), mean(y_benign)))

## identifying the RARE benign class: AUC = 0.621, prevalence = 0.081

cat(sprintf("average precision (area under PR) = %.3f; a no-skill classifier gives %.3f\n",
            mean(pp$precision, na.rm = TRUE), mean(y_benign)))

## average precision (area under PR) = 0.116; a no-skill classifier gives 0.081

g1 <- ggplot(rp, aes(fpr, tpr)) +
  geom_abline(slope = 1, intercept = 0, colour = "grey60", linetype = "dashed") +
  geom_line(colour = bpad_pal[1], linewidth = 0.9) +
  labs(title = "ROC for the rare class", x = "false positive rate", y = "true positive rate")
g2 <- ggplot(pp, aes(recall, precision)) +
  geom_hline(yintercept = mean(y_benign), linetype = "dashed", colour = "grey55") +
  geom_line(colour = bpad_pal[2], linewidth = 0.9) +
  coord_cartesian(ylim = c(0, 1)) +
  labs(title = "precision-recall for the rare class", x = "recall", y = "precision")

## net benefit for the primary task (identifying malignancy)
th <- seq(0.5, 0.99, by = 0.005)
nb <- vapply(th, function(t) {
  pos <- oof_full >= t
  tp <- sum(pos & d$malignant); fp <- sum(pos & !d$malignant)
  tp/nrow(d) - (fp/nrow(d))*(t/(1 - t))
}, numeric(1))
nb_all  <- base_rate - (1 - base_rate)*(th/(1 - th))
nbd <- rbind(data.frame(t = th, nb = nb,     s = "model"),
             data.frame(t = th, nb = nb_all, s = "treat everyone"),
             data.frame(t = th, nb = 0,      s = "treat no one"))
g3 <- ggplot(nbd, aes(t, nb, colour = s)) +
  geom_line(linewidth = 0.9) +
  coord_cartesian(ylim = c(-0.1, max(nb, na.rm = TRUE) + 0.05)) +
  scale_colour_manual(values = bpad_pal[c(1,2,7)]) +
  labs(title = "decision curve", x = "threshold probability", y = "net benefit")
bpad_grid(g1, g2, g3, ncol = 3)

$Two views of the same predictions, and their decision consequences. Left: ROC and precision-recall curves for identifying the rare benign class. The PR curve exposes how quickly precision collapses, information the ROC curve compresses into a small region near the origin. Right: net benefit, Eq. \@ref(eq:net-benefit), against threshold probability, compared with the two default strategies of treating everyone and treating no one. A model is clinically useful only over the threshold range where its curve is the highest of the three.$

Figure 21: Two views of the same predictions, and their decision consequences. Left: ROC and precision-recall curves for identifying the rare benign class. The PR curve exposes how quickly precision collapses, information the ROC curve compresses into a small region near the origin. Right: net benefit, Eq. (39), against threshold probability, compared with the two default strategies of treating everyone and treating no one. A model is clinically useful only over the threshold range where its curve is the highest of the three.

cat(sprintf("\nthreshold range where the model beats both defaults: %s\n",
            if (any(nb > pmax(nb_all, 0)))
              sprintf("%.2f to %.2f", min(th[nb > pmax(nb_all, 0)]),
                      max(th[nb > pmax(nb_all, 0)]))
            else "none"))

## 
## threshold range where the model beats both defaults: 0.50 to 0.95

8.7.9 Prediction versus causation, on real data

## Is `radical` a predictor of malignancy, a consequence of it, or a confounded proxy?
print(table(procedure = kits_m$procedure, malignant = kits_m$malignant))

##                      malignant
## procedure             FALSE TRUE
##   partial_nephrectomy    12  128
##   radical_nephrectomy     6   64

cat(sprintf("\nmalignancy: partial %.3f, radical %.3f\n",
            mean(kits_m$malignant[kits_m$radical == 0]),
            mean(kits_m$malignant[kits_m$radical == 1])))

## 
## malignancy: partial 0.914, radical 0.914

cat(sprintf("median tumour size: partial %.1f cm, radical %.1f cm\n",
            median(kits_m$radiographic_size_cm[kits_m$radical == 0], na.rm = TRUE),
            median(kits_m$radiographic_size_cm[kits_m$radical == 1], na.rm = TRUE)))

## median tumour size: partial 3.0 cm, radical 6.6 cm

m_un  <- suppressWarnings(glm(malignant ~ radical, kits_m, family = binomial()))
m_adj <- suppressWarnings(glm(malignant ~ radical + log_size + age, kits_m,
                              family = binomial()))
print(knitr::kable(data.frame(
  model = c("unadjusted", "adjusted for size and age"),
  OR_radical = round(exp(c(coef(m_un)["radical"], coef(m_adj)["radical"])), 3),
  p_value    = signif(c(summary(m_un)$coefficients["radical", 4],
                        summary(m_adj)$coefficients["radical", 4]), 3)),
  row.names = FALSE,
  caption = "The apparent association between operative extent and malignancy shrinks once tumour size is adjusted for, because size is the reason the operation was chosen."))

## 
## 
## Table: (\#tab:causal-vs-predictive)The apparent association between operative extent and malignancy shrinks once tumour size is adjusted for, because size is the reason the operation was chosen.
## 
## |model                     | OR_radical| p_value|
## |:-------------------------|----------:|-------:|
## |unadjusted                |      1.000|  1.0000|
## |adjusted for size and age |      0.261|  0.0599|

The operative extent is chosen after the imaging and partly because of it. It is therefore a legitimate predictor in a model whose inputs are permitted to include the surgical plan, and an illegitimate one in a model intended to inform whether to operate. Worse, it is a collider-adjacent variable: conditioning on a downstream consequence of the exposure can induce associations that are not causal. The only defence is the intended-use statement of Section 8.1.4, written before the data were opened.

Section 8.7 summary.

Report discrimination, calibration, and utility; a good AUC with poor calibration lies about absolute risk.
On this real cohort the malignancy model has modest discrimination and near-zero Brier skill, and the honest conclusion is about the cohort, not the algorithm.
kNN exposes the bias-variance trade-off in one tuning parameter; the optimum must be chosen inside resampling.
Between-model differences are smaller than between-fold spread when events, not capacity, are the constraint.
eGFR decline is the outcome this cohort supports, and the radical coefficient matches a mechanistic prediction.
Under strong imbalance, precision-recall and net benefit are more informative than ROC.
A variable chosen because of the imaging cannot be a predictor in a model that informs the same choice.

Checkpoint 8.7. The malignancy model achieves an out-of-fold AUC around 0.62 with a Brier skill score near zero. A colleague proposes adding 400 texture features to improve it. Using Sections 8.3.6, 8.3.7, and 8.5.6, give three quantitative reasons this is likely to make the reported performance better and the actual performance worse.

8.8 Time-to-Event and Longitudinal Modeling

8.8.1 Why a binary endpoint is usually inadmissible

The cohort records vital_status and vital_days_after_surgery. Reducing these to “dead at last contact” throws away follow-up duration and misclassifies every patient censored before the horizon of interest. With 21 deaths among 210 patients and follow-up ranging from days to more than eight years, that reduction would be indefensible.

The three functions that describe a time-to-event outcome are linked by identities worth stating once. The survival function is $S(t) = P(T > t)$; the hazard is the instantaneous event rate among survivors,

\[\begin{equation} h(t) = \lim_{\Delta t\to 0}\frac{P(t \le T < t+\Delta t \mid T \ge t)}{\Delta t} = -\frac{d}{dt}\log S(t), \tag{40} \end{equation}\]

and integrating gives the cumulative hazard and hence

\[\begin{equation} S(t) = \exp\{-H(t)\}, \qquad H(t) = \int_0^t h(u)\,du . \tag{41} \end{equation}\]

Equation (41) is a first-order ODE solution of exactly the kind developed in Chapter 1: a constant hazard $h$ gives $S(t) = e^{-ht}$, the same exponential law that governs radioactive decay in Chapter 6 and drug clearance in Chapter 1. A hazard that rises with time gives a Weibull; one that falls gives a heavy-tailed survival curve.

8.8.2 The Kaplan-Meier (KM) estimator

With distinct event times $t_1 < t_2 < \cdots$, $d_j$ events and $n_j$ at risk at $t_j$,

\[\begin{equation} \widehat S_{\mathrm{KM}}(t) = \prod_{t_j \le t}\left(1 - \frac{d_j}{n_j}\right), \tag{42} \end{equation}\]

a product over the conditional probabilities of surviving each event time. Censored observations leave the risk set without contributing an event, which is how incomplete follow-up is handled without discarding the patient.

if (has_pkg("survival")) {
  library(survival)
  surv_obj <- Surv(kits_m$time_yr, kits_m$event)
  km0 <- survfit(surv_obj ~ 1)
  cat(sprintf("n = %d, events = %d, median follow-up = %.2f years, maximum = %.2f\n",
              nrow(kits_m), sum(kits_m$event),
              median(kits_m$time_yr), max(kits_m$time_yr)))
  st <- summary(km0, times = c(1, 2, 3, 5))
  print(knitr::kable(data.frame(years = st$time, at_risk = st$n.risk, events = st$n.event,
                          survival = round(st$surv, 3),
                          lower = round(st$lower, 3), upper = round(st$upper, 3)),
               row.names = FALSE,
               caption = "Kaplan-Meier survival estimates. The number at risk falls quickly, which is why the late estimates are imprecise."))

  kmd <- data.frame(t = km0$time, s = km0$surv, lo = km0$lower, hi = km0$upper)
  kmd <- rbind(data.frame(t = 0, s = 1, lo = 1, hi = 1), kmd)
  s1 <- ggplot(kmd, aes(t, s)) +
    geom_ribbon(aes(ymin = lo, ymax = hi), fill = bpad_pal[1], alpha = 0.18) +
    geom_step(colour = bpad_pal[1], linewidth = 0.9) +
    coord_cartesian(ylim = c(0.5, 1)) +
    labs(title = "overall survival", x = "years after surgery", y = "S(t)")

  km1 <- survfit(surv_obj ~ radical, data = kits_m)
  gl  <- rep(names(km1$strata), km1$strata)
  kmg <- data.frame(t = km1$time, s = km1$surv,
                    g = ifelse(grepl("=0", gl), "partial", "radical"))
  kmg <- rbind(data.frame(t = 0, s = 1, g = c("partial","radical")), kmg)
  s2 <- ggplot(kmg, aes(t, s, colour = g)) +
    geom_step(linewidth = 0.9) +
    coord_cartesian(ylim = c(0.5, 1)) +
    scale_colour_manual(values = bpad_pal[1:2]) +
    labs(title = "by operative extent", x = "years after surgery", y = "S(t)")
  bpad_grid(s1, s2, ncol = 2)

  lr <- survdiff(surv_obj ~ radical, data = kits_m)
  cat(sprintf("\nlog-rank test: chi-square = %.2f on %d df, p = %.4f\n",
              lr$chisq, length(lr$n) - 1,
              stats::pchisq(lr$chisq, length(lr$n) - 1, lower.tail = FALSE)))
  print(data.frame(group = c("partial","radical"), n = as.vector(lr$n),
                   observed = as.vector(lr$obs), expected = round(as.vector(lr$exp), 2)))
}

## n = 210, events = 21, median follow-up = 2.16 years, maximum = 8.41
## 
## 
## Table: (\#tab:kaplan-meier)Kaplan-Meier survival estimates. The number at risk falls quickly, which is why the late estimates are imprecise.
## 
## | years| at_risk| events| survival| lower| upper|
## |-----:|-------:|------:|--------:|-----:|-----:|
## |     1|     152|      9|    0.949| 0.917| 0.982|
## |     2|     118|      7|    0.900| 0.854| 0.948|
## |     3|      71|      3|    0.870| 0.816| 0.928|
## |     5|      37|      0|    0.870| 0.816| 0.928|

Overall and stratified survival in the real cohort. Left: the Kaplan-Meier curve with 95 percent confidence band; the band widens markedly beyond four years, where few patients remain at risk. Right: survival stratified by operative extent, with the number at risk annotated. Radical nephrectomy is associated with worse survival on an unadjusted comparison, which Section 8.8.4 shows is largely explained by tumour size.

Figure 22: Overall and stratified survival in the real cohort. Left: the Kaplan-Meier curve with 95 percent confidence band; the band widens markedly beyond four years, where few patients remain at risk. Right: survival stratified by operative extent, with the number at risk annotated. Radical nephrectomy is associated with worse survival on an unadjusted comparison, which Section 8.8.4 shows is largely explained by tumour size.

## 
## log-rank test: chi-square = 6.53 on 1 df, p = 0.0106
##     group   n observed expected
## 1 partial 140        8    13.56
## 2 radical  70       13     7.44

8.8.3 The Cox proportional-hazards model

The Cox model specifies

\[\begin{equation} h(t \mid X) = h_0(t)\exp(X^{\mathsf T}\beta), \tag{43} \end{equation}\]

leaving the baseline hazard $h_0$ unspecified and estimating $\beta$ from the partial likelihood

\[\begin{equation} L(\beta) = \prod_{i: \delta_i = 1} \frac{\exp(X_i^{\mathsf T}\beta)}{\sum_{j \in \mathcal{R}(t_i)}\exp(X_j^{\mathsf T}\beta)}, \tag{44} \end{equation}\]

where $\mathcal{R}(t_i)$ is the risk set at $t_i$. Two consequences follow. Because $h_0$ cancels, $e^{\beta_j}$ is a hazard ratio, a relative quantity, and absolute risk requires an estimate of $h_0$ as well. And because only the ordering of event times enters Eq. (44), the model is invariant to any monotone transformation of the time axis.

The proportional-hazards assumption is testable, and it must be tested rather than assumed.

if (has_pkg("survival")) {
  cx <- coxph(Surv(time_yr, event) ~ age + log_size + radical + male, data = kits_m)
  sm <- summary(cx)
  print(knitr::kable(data.frame(
    term = rownames(sm$coefficients),
    hazard_ratio = round(sm$coefficients[, "exp(coef)"], 3),
    CI_low  = round(sm$conf.int[, "lower .95"], 3),
    CI_high = round(sm$conf.int[, "upper .95"], 3),
    p_value = signif(sm$coefficients[, "Pr(>|z|)"], 3)),
    row.names = FALSE,
    caption = "Cox proportional-hazards model for overall survival."))
  cat(sprintf("events = %d, concordance = %.3f (se %.3f)\n",
              sm$nevent, sm$concordance[1], sm$concordance[2]))
  cat(sprintf("events per parameter = %.1f\n", sm$nevent/length(coef(cx))))

  zph <- cox.zph(cx)
  cat("\nproportional-hazards test (Schoenfeld residuals):\n")
  print(round(zph$table, 4))
  cat(sprintf("\nglobal test p = %.3f: %s\n", zph$table["GLOBAL", "p"],
              ifelse(zph$table["GLOBAL", "p"] > 0.05,
                     "no evidence against proportional hazards",
                     "evidence of a time-varying effect")))

  cd <- data.frame(term = rownames(sm$conf.int),
                   hr = sm$conf.int[, "exp(coef)"],
                   lo = sm$conf.int[, "lower .95"],
                   hi = sm$conf.int[, "upper .95"])
  c1 <- ggplot(cd, aes(hr, reorder(term, hr))) +
    geom_vline(xintercept = 1, linetype = "dashed", colour = "grey50") +
    geom_errorbarh(aes(xmin = lo, xmax = hi), height = 0.18, colour = bpad_pal[1]) +
    geom_point(size = 2.4, colour = bpad_pal[1]) +
    scale_x_log10() +
    labs(title = "hazard ratios (log scale)", x = "hazard ratio", y = NULL)

  idx <- which(rownames(zph$table) == "log_size")
  sr  <- data.frame(t = zph$x, r = zph$y[, idx])
  c2 <- ggplot(sr, aes(t, r)) +
    geom_hline(yintercept = 0, colour = "grey55") +
    geom_point(alpha = 0.7, colour = bpad_pal[2], size = 1.8) +
    geom_smooth(method = "loess", se = TRUE, colour = bpad_pal[1], linewidth = 0.8) +
    labs(title = "scaled Schoenfeld residuals, log_size",
         x = "transformed time", y = "residual (beta(t))")
  bpad_grid(c1, c2, ncol = 2)
}

## 
## 
## Table: (\#tab:cox-model)Cox proportional-hazards model for overall survival.
## 
## |term     | hazard_ratio| CI_low| CI_high| p_value|
## |:--------|------------:|------:|-------:|-------:|
## |age      |        1.022|  0.986|   1.058|   0.233|
## |log_size |        2.789|  1.023|   7.603|   0.045|
## |radical  |        1.133|  0.336|   3.818|   0.840|
## |male     |        1.787|  0.690|   4.629|   0.232|
## events = 21, concordance = 0.741 (se 0.056)
## events per parameter = 5.2
## 
## proportional-hazards test (Schoenfeld residuals):
##           chisq df      p
## age      2.0237  1 0.1549
## log_size 1.2186  1 0.2696
## radical  1.4680  1 0.2257
## male     0.2016  1 0.6534
## GLOBAL   4.5408  4 0.3377
## 
## global test p = 0.338: no evidence against proportional hazards

Cox model on the real cohort. Left: hazard ratios with 95 percent confidence intervals on a logarithmic axis. Right: scaled Schoenfeld residuals against time for tumour size; a non-zero slope would indicate a time-varying effect and hence a violation of proportional hazards. The formal test is reported below.

Figure 23: Cox model on the real cohort. Left: hazard ratios with 95 percent confidence intervals on a logarithmic axis. Right: scaled Schoenfeld residuals against time for tumour size; a non-zero slope would indicate a time-varying effect and hence a violation of proportional hazards. The formal test is reported below.

8.8.4 Variable Confounding

Section 8.8.2 found a significant unadjusted survival difference by operative extent. Section 8.8.3 found no detectable adjusted effect. Both computations are correct, and the discrepancy is the whole lesson.

if (has_pkg("survival")) {
  fits <- list(
    "radical alone"                   = coxph(Surv(time_yr, event) ~ radical, kits_m),
    "+ tumour size"                   = coxph(Surv(time_yr, event) ~ radical + log_size, kits_m),
    "+ tumour size, age, sex"         = coxph(Surv(time_yr, event) ~ radical + log_size +
                                                age + male, kits_m))
  ct <- do.call(rbind, lapply(names(fits), function(nm) {
    s <- summary(fits[[nm]])
    data.frame(model = nm,
               HR    = s$conf.int["radical", "exp(coef)"],
               lo    = s$conf.int["radical", "lower .95"],
               hi    = s$conf.int["radical", "upper .95"],
               p     = s$coefficients["radical", "Pr(>|z|)"])
  }))
  print(knitr::kable(data.frame(model = ct$model, HR_radical = round(ct$HR, 3),
                          CI = sprintf("[%.2f, %.2f]", ct$lo, ct$hi),
                          p_value = signif(ct$p, 3)),
               row.names = FALSE,
               caption = "Hazard ratio for radical nephrectomy under progressive adjustment. The effect attenuates toward the null as the indication for surgery is controlled."))

  f1 <- ggplot(ct, aes(HR, reorder(model, HR))) +
    geom_vline(xintercept = 1, linetype = "dashed", colour = "grey50") +
    geom_errorbarh(aes(xmin = lo, xmax = hi), height = 0.16, colour = bpad_pal[1]) +
    geom_point(size = 2.6, colour = bpad_pal[1]) +
    scale_x_log10() +
    labs(title = "confounding by indication", x = "hazard ratio for radical (log scale)",
         y = NULL)
  f2 <- ggplot(kits_m, aes(factor(radical, labels = c("partial", "radical")),
                           radiographic_size_cm, fill = factor(radical))) +
    geom_boxplot(alpha = 0.5, show.legend = FALSE) +
    scale_fill_manual(values = bpad_pal[1:2]) +
    labs(title = "the mechanism: size drives the choice",
         x = NULL, y = "radiographic size (cm)")
  bpad_grid(f1, f2, ncol = 2)
}

## 
## 
## Table: (\#tab:confounding-demo)Hazard ratio for radical nephrectomy under progressive adjustment. The effect attenuates toward the null as the indication for surgery is controlled.
## 
## |model                   | HR_radical|CI           | p_value|
## |:-----------------------|----------:|:------------|-------:|
## |radical alone           |      3.003|[1.24, 7.28] |  0.0149|
## |+ tumour size           |      1.325|[0.40, 4.40] |  0.6450|
## |+ tumour size, age, sex |      1.133|[0.34, 3.82] |  0.8400|

The same variable, two conclusions. Left: unadjusted and adjusted hazard ratios for radical nephrectomy; adjustment for tumour size moves the estimate toward the null and widens the interval. Right: the mechanism, tumour size by operative extent, showing why the two analyses disagree. Radical nephrectomy is performed for larger tumours, and larger tumours carry worse prognosis, so the unadjusted comparison attributes the size effect to the operation.

Figure 24: The same variable, two conclusions. Left: unadjusted and adjusted hazard ratios for radical nephrectomy; adjustment for tumour size moves the estimate toward the null and widens the interval. Right: the mechanism, tumour size by operative extent, showing why the two analyses disagree. Radical nephrectomy is performed for larger tumours, and larger tumours carry worse prognosis, so the unadjusted comparison attributes the size effect to the operation.

Confounding by indication. The treatment was chosen because of a prognostic factor. No amount of data volume repairs this; only design does, through randomization, an instrument, or explicit adjustment for the full indication with a documented assumption of no unmeasured confounding. A machine-learning model with 400 features and a validated pipeline reproduces this bias faithfully and reports it with a tight confidence interval.

8.8.5 Immortal-time bias and landmark analysis

A frequent error is to classify patients by something that happens after the time origin, for example “received adjuvant therapy” versus “did not”. Patients in the treated group must survive long enough to be treated, so they are guaranteed event-free over that interval. That guarantee, not the treatment, produces the apparent benefit.

The landmark method fixes this: choose a landmark time $t_L$, exclude patients with an event or censoring before $t_L$, classify exposure using information available at $t_L$ only, and start the clock at $t_L$.

if (has_pkg("survival")) {
  set.seed(31)
  n <- 900
  Ttrue <- rexp(n, 0.25); Cc <- runif(n, 0, 8)
  obs <- pmin(Ttrue, Cc); ev <- Ttrue <= Cc
  ## exposure can only occur after 0.5 years, and has NO effect on survival
  t_expose <- 0.5
  exposed  <- (obs > t_expose) & (runif(n) < 0.55)
  naive <- survfit(Surv(obs, ev) ~ exposed)
  gl <- rep(names(naive$strata), naive$strata)
  nd <- rbind(data.frame(t = 0, s = 1, g = c("unexposed","exposed")),
              data.frame(t = naive$time, s = naive$surv,
                         g = ifelse(grepl("FALSE", gl), "unexposed", "exposed")))
  n1 <- ggplot(nd, aes(t, s, colour = g)) + geom_step(linewidth = 0.9) +
    scale_colour_manual(values = bpad_pal[1:2]) +
    labs(title = "naive: apparent benefit from immortal time",
         x = "years", y = "S(t)")

  keep <- obs > t_expose
  lm_fit <- survfit(Surv(obs[keep] - t_expose, ev[keep]) ~ exposed[keep])
  gl2 <- rep(names(lm_fit$strata), lm_fit$strata)
  ld <- rbind(data.frame(t = 0, s = 1, g = c("unexposed","exposed")),
              data.frame(t = lm_fit$time, s = lm_fit$surv,
                         g = ifelse(grepl("FALSE", gl2), "unexposed", "exposed")))
  n2 <- ggplot(ld, aes(t, s, colour = g)) + geom_step(linewidth = 0.9) +
    scale_colour_manual(values = bpad_pal[1:2]) +
    labs(title = "landmark at 0.5 years: no difference",
         x = "years after the landmark", y = "S(t)")
  bpad_grid(n1, n2, ncol = 2)

  hr_naive <- summary(coxph(Surv(obs, ev) ~ exposed))$conf.int
  hr_lm    <- summary(coxph(Surv(obs[keep] - t_expose, ev[keep]) ~ exposed[keep]))$conf.int
  print(knitr::kable(data.frame(
    analysis = c("naive (immortal time included)", "landmark at 0.5 years"),
    HR = round(c(hr_naive[1, 1], hr_lm[1, 1]), 3),
    CI = c(sprintf("[%.2f, %.2f]", hr_naive[1, 3], hr_naive[1, 4]),
           sprintf("[%.2f, %.2f]", hr_lm[1, 3], hr_lm[1, 4])),
    truth = c("HR = 1", "HR = 1")),
    row.names = FALSE,
    caption = "The exposure was simulated to have no effect whatsoever. The naive analysis nonetheless reports a large protective hazard ratio; the landmark analysis recovers the truth."))
}

Immortal-time bias and its repair, in a designed simulation in which the exposure has exactly no effect. Left: naive classification by an exposure that can only occur after six months creates a large apparent survival benefit, because exposed patients had to survive to be exposed. Right: the landmark analysis, starting the clock at six months and using only exposure status known by then, correctly recovers no difference.

Figure 25: Immortal-time bias and its repair, in a designed simulation in which the exposure has exactly no effect. Left: naive classification by an exposure that can only occur after six months creates a large apparent survival benefit, because exposed patients had to survive to be exposed. Right: the landmark analysis, starting the clock at six months and using only exposure status known by then, correctly recovers no difference.

## 
## 
## Table: (\#tab:landmark-bias)The exposure was simulated to have no effect whatsoever. The naive analysis nonetheless reports a large protective hazard ratio; the landmark analysis recovers the truth.
## 
## |analysis                       |    HR|CI           |truth  |
## |:------------------------------|-----:|:------------|:------|
## |naive (immortal time included) | 0.652|[0.55, 0.78] |HR = 1 |
## |landmark at 0.5 years          | 0.948|[0.78, 1.15] |HR = 1 |

8.8.6 Competing events and the cumulative incidence function

When several event types compete, the probability of experiencing cause $k$ by time $t$ is the cause-specific cumulative incidence function

\[\begin{equation} F_k(t) = \int_0^t S(u^-)\,\lambda_k(u)\,du, \tag{45} \end{equation}\]

where $S$ is the overall event-free survival across all causes. The Aalen–Johansen estimator is the empirical version. Treating competing events as ordinary censoring and reporting $1 - \widehat S_{\mathrm{KM}}$ for cause $k$ overestimates $F_k$, because Kaplan–Meier implicitly assumes that a patient censored by a competing death could still later experience cause $k$, which is impossible.

KiTS19 records vital status but not cause of death, so the demonstration below is a designed simulation with known truth, labelled as such.

if (has_pkg("survival") && has_pkg("cmprsk")) {
  set.seed(20)
  n <- 4000; h1 <- 0.06; h2 <- 0.20        # cause 1 rate, competing cause rate
  t1 <- rexp(n, h1); t2 <- rexp(n, h2); cc <- runif(n, 0, 12)
  tt <- pmin(t1, t2, cc)
  cause <- ifelse(tt == t1, 1L, ifelse(tt == t2, 2L, 0L))
  cat(sprintf("cause 1 events = %d, competing events = %d, censored = %d\n",
              sum(cause == 1), sum(cause == 2), sum(cause == 0)))

  grid <- c(1, 2, 3, 5, 8)
  naive <- 1 - summary(survfit(Surv(tt, cause == 1) ~ 1), times = grid)$surv
  ci    <- cmprsk::cuminc(tt, cause, cencode = 0)
  tp    <- cmprsk::timepoints(ci, grid)
  truth <- (h1/(h1 + h2))*(1 - exp(-(h1 + h2)*grid))   # closed form for this design
  print(knitr::kable(data.frame(
    years = grid, truth = round(truth, 4),
    Aalen_Johansen = round(tp$est[1, ], 4),
    naive_1_minus_KM = round(naive, 4),
    overestimate_pct = round(100*(naive/truth - 1), 0)),
    row.names = FALSE,
    caption = "Cause-1 cumulative incidence: the truth, the correct estimator, and the naive one."))

  tg <- seq(0.05, 10, by = 0.05)
  aj <- cmprsk::timepoints(ci, tg)
  nk <- 1 - summary(survfit(Surv(tt, cause == 1) ~ 1), times = tg, extend = TRUE)$surv
  cd <- rbind(
    data.frame(t = tg, v = (h1/(h1+h2))*(1 - exp(-(h1+h2)*tg)), s = "truth"),
    data.frame(t = tg, v = as.numeric(aj$est[1, ]),             s = "Aalen-Johansen CIF"),
    data.frame(t = tg, v = nk,                                  s = "naive 1 - KM"))
  w1 <- ggplot(cd, aes(t, v, colour = s, linetype = s)) +
    geom_line(linewidth = 0.9, na.rm = TRUE) +
    scale_colour_manual(values = c("grey30", bpad_pal[1], bpad_pal[2])) +
    scale_linetype_manual(values = c("dashed", "solid", "solid")) +
    labs(title = "cause-1 cumulative incidence", x = "years", y = "F_1(t)")
  w2 <- ggplot(rbind(
      data.frame(t = tg, v = as.numeric(aj$est[1, ]), s = "cause 1 (of interest)"),
      data.frame(t = tg, v = as.numeric(aj$est[2, ]), s = "cause 2 (competing)")),
      aes(t, v, fill = s)) +
    geom_area(alpha = 0.6, position = "stack", na.rm = TRUE) +
    scale_fill_manual(values = bpad_pal[c(1,5)]) +
    labs(title = "stacked cumulative incidence", x = "years", y = "probability")
  bpad_grid(w1, w2, ncol = 2)
} else {
  cat("survival and/or cmprsk unavailable; competing-risk demonstration skipped\n")
}

## cause 1 events = 668, competing events = 2114, censored = 1218
## 
## 
## Table: (\#tab:competing-risks)Cause-1 cumulative incidence: the truth, the correct estimator, and the naive one.
## 
## | years|  truth| Aalen_Johansen| naive_1_minus_KM| overestimate_pct|
## |-----:|------:|--------------:|----------------:|----------------:|
## |     1| 0.0528|         0.0632|           0.0698|               32|
## |     2| 0.0936|         0.1006|           0.1209|               29|
## |     3| 0.1250|         0.1345|           0.1769|               42|
## |     5| 0.1679|         0.1797|           0.2779|               66|
## |     8| 0.2019|         0.2062|           0.3700|               83|

Why the complement of Kaplan-Meier is wrong under competing events. A simulation with a known cause-1 cumulative incidence (dashed line) and a competing cause three times as frequent. The Aalen-Johansen estimator tracks the truth; the naive complement of Kaplan-Meier diverges upward and overstates cause-1 risk by more than 80 percent at eight years. The two stacked cumulative incidence curves on the right sum to the overall event probability.

Figure 26: Why the complement of Kaplan-Meier is wrong under competing events. A simulation with a known cause-1 cumulative incidence (dashed line) and a competing cause three times as frequent. The Aalen-Johansen estimator tracks the truth; the naive complement of Kaplan-Meier diverges upward and overstates cause-1 risk by more than 80 percent at eight years. The two stacked cumulative incidence curves on the right sum to the overall event probability.

Which hazard answers which question. Use the cause-specific hazard, Eq. (27), for aetiological questions: does this exposure act on the biology of this event? Use the cumulative incidence function, Eq. (45), or a Fine–Gray subdistribution model for prognostic questions: what is this patient’s absolute probability of this event by five years? In an older surgical cohort, where competing non-cancer death is common, the two can point in different directions, and reporting only the first while writing clinical conclusions about the second is a substantive error rather than a technical one.

8.8.7 Longitudinal tumour growth: the ODEs of Chapter 1

Serial imaging turns a static biomarker into a trajectory. The three canonical growth laws of Chapter 1 reappear here, now fitted to volumes rather than cell counts. Unconstrained exponential growth,

\[\begin{equation} \frac{dV}{dt} = kV \quad\Longrightarrow\quad V(t) = V_0 e^{kt}, \tag{46} \end{equation}\]

has the specific growth rate $k$ and the volume doubling time

\[\begin{equation} \mathrm{VDT} = \frac{\ln 2}{k}, \qquad k = \frac{\ln(V_2/V_1)}{t_2 - t_1}, \tag{47} \end{equation}\]

which is the standard clinical summary of a two-timepoint measurement. Logistic and Gompertz growth add saturation at a carrying capacity $K$:

\[\begin{equation} \frac{dV}{dt} = kV\left(1 - \frac{V}{K}\right), \qquad \frac{dV}{dt} = k_G V\ln\frac{K}{V}. \tag{48} \end{equation}\]

V0 <- 8; K <- 500; VDT <- 500
k  <- log(2)/VDT
kG <- k/log(K/V0)                              # matched initial specific growth rate
tg <- seq(0, 3650, by = 5)
Ve <- V0*exp(k*tg)
Vl <- K/(1 + (K/V0 - 1)*exp(-k*tg))
Vg <- K*exp(log(V0/K)*exp(-kG*tg))
d_of <- function(V) 2*(3*V/(4*pi))^(1/3)

gm <- rbind(data.frame(t = tg/365.25, V = Ve, m = "exponential"),
            data.frame(t = tg/365.25, V = Vl, m = "logistic"),
            data.frame(t = tg/365.25, V = Vg, m = "Gompertz"))
print(knitr::kable(data.frame(
  years = c(1, 2, 3, 5, 10),
  exponential = round(V0*exp(k*c(1,2,3,5,10)*365.25), 1),
  logistic    = round(K/(1 + (K/V0 - 1)*exp(-k*c(1,2,3,5,10)*365.25)), 1),
  Gompertz    = round(K*exp(log(V0/K)*exp(-kG*c(1,2,3,5,10)*365.25)), 1)),
  row.names = FALSE,
  caption = "Predicted tumour volume (cubic centimetres). At two years the three models agree to within a few percent; at ten years they differ by nearly an order of magnitude."))

## 
## 
## Table: (\#tab:growth-models)Predicted tumour volume (cubic centimetres). At two years the three models agree to within a few percent; at ten years they differ by nearly an order of magnitude.
## 
## | years| exponential| logistic| Gompertz|
## |-----:|-----------:|--------:|--------:|
## |     1|        13.3|     13.1|     12.9|
## |     2|        22.0|     21.4|     19.6|
## |     3|        36.5|     34.6|     28.5|
## |     5|       100.6|     84.9|     53.1|
## |    10|      1265.1|    360.0|    148.3|

h1 <- ggplot(gm, aes(t, V, colour = m)) +
  geom_line(linewidth = 0.9) +
  geom_hline(yintercept = K, linetype = "dotted", colour = "grey50") +
  scale_y_log10() + coord_cartesian(ylim = c(5, 1500)) +
  scale_colour_manual(values = bpad_pal[1:3]) +
  labs(title = "volume (log scale)", x = "years", y = expression("V (cm"^3*")"))
h2 <- ggplot(gm, aes(t, d_of(V), colour = m)) +
  geom_line(linewidth = 0.9) +
  scale_colour_manual(values = bpad_pal[1:3]) +
  labs(title = "equivalent spherical diameter", x = "years", y = "d (cm)")
bpad_grid(h1, h2, ncol = 2)

Three growth laws with matched initial specific growth rate, at parameters realistic for renal masses (initial volume 8 cubic centimetres, volume doubling time 500 days, carrying capacity 500 cubic centimetres). Left: on a logarithmic axis the three curves are indistinguishable for the first two years, which is the entire duration of most surveillance studies, and diverge eightfold by ten years. Right: the corresponding equivalent spherical diameters, showing why diameter-based surveillance is so insensitive: a doubling of diameter requires an eightfold increase in volume.

Figure 27: Three growth laws with matched initial specific growth rate, at parameters realistic for renal masses (initial volume 8 cubic centimetres, volume doubling time 500 days, carrying capacity 500 cubic centimetres). Left: on a logarithmic axis the three curves are indistinguishable for the first two years, which is the entire duration of most surveillance studies, and diverge eightfold by ten years. Right: the corresponding equivalent spherical diameters, showing why diameter-based surveillance is so insensitive: a doubling of diameter requires an eightfold increase in volume.

cat(sprintf("volume doubling time = %.0f days (k = %.5f /day)\n", VDT, k))

## volume doubling time = 500 days (k = 0.00139 /day)

cat(sprintf("volume 8 -> 64 cm3 is an 8-fold increase, but diameter only %.2f -> %.2f cm\n",
            d_of(8), d_of(64)))

## volume 8 -> 64 cm3 is an 8-fold increase, but diameter only 2.48 -> 4.96 cm

cat(sprintf("logistic inflection at day %.0f; Gompertz inflection at day %.0f\n",
            log(K/V0 - 1)/k, log(log(K/V0))/kG))

## logistic inflection at day 2971; Gompertz inflection at day 4234

From equation to surveillance protocol. Because $d \propto V^{1/3}$, a $26\%$ increase in diameter is required to signal a doubling of volume. Section 8.5.6 showed that a one-voxel boundary shift already moves diameter by roughly $\pm8\%$. A surveillance protocol that declares progression on a $20\%$ diameter increase is therefore operating close to its own measurement noise, and the growth-rate estimate inherits that noise multiplied by three (Eq. (19) applied to Eq. (47)). This is why volumetric criteria, despite being noisier per measurement, can detect real growth earlier: they do not throw away the cube.

## Propagate segmentation uncertainty into the estimated doubling time.
## Two measurements, each with a relative volume uncertainty from Section 8.5.6.
V1 <- 8; V2 <- 12; dt_days <- 365
rel_V <- 0.25                                    # +/-25%, the one-voxel band
k_hat <- log(V2/V1)/dt_days
cat(sprintf("point estimate: k = %.5f /day, VDT = %.0f days\n", k_hat, log(2)/k_hat))

## point estimate: k = 0.00111 /day, VDT = 624 days

set.seed(9); B <- 20000
V1s <- V1*exp(rnorm(B, 0, rel_V)); V2s <- V2*exp(rnorm(B, 0, rel_V))
k_mc <- log(V2s/V1s)/dt_days
vdt_mc <- log(2)/k_mc
vdt_mc <- vdt_mc[k_mc > 0]                       # negative k means apparent shrinkage
cat(sprintf("Monte Carlo: P(apparent shrinkage) = %.3f\n", mean(k_mc <= 0)))

## Monte Carlo: P(apparent shrinkage) = 0.129

cat(sprintf("VDT 95%% interval among growing draws: [%.0f, %.0f] days\n",
            quantile(vdt_mc, 0.025), quantile(vdt_mc, 0.975)))

## VDT 95% interval among growing draws: [226, 6962] days

cat("\nA 50% volume increase over one year is compatible with doubling times spanning\n")

## 
## A 50% volume increase over one year is compatible with doubling times spanning

cat("an order of magnitude once segmentation uncertainty is propagated. Reporting a\n")

## an order of magnitude once segmentation uncertainty is propagated. Reporting a

cat("single doubling time from two scans without this interval is overconfident.\n")

## single doubling time from two scans without this interval is overconfident.

8.8.8 Repeated measurements: mixed-effects and joint models

With several measurements per patient, observations are correlated (Section 8.1.3) and a linear model that ignores this understates standard errors. A linear mixed-effects model adds patient-specific random effects,

\[\begin{equation} y_{ij} = \beta_0 + \beta_1 t_{ij} + b_{0i} + b_{1i}t_{ij} + \varepsilon_{ij}, \qquad \begin{pmatrix}b_{0i}\\ b_{1i}\end{pmatrix} \sim \mathcal{N}(0, G), \quad \varepsilon_{ij}\sim\mathcal{N}(0,\sigma^2), \tag{49} \end{equation}\]

separating the population trajectory ($\beta$, the fixed effects) from individual deviations ($b_i$, the random effects). The random intercept absorbs each patient’s baseline level and the random slope each patient’s own growth rate, which is exactly the decomposition a surveillance clinician reasons with informally.

if (has_pkg("nlme")) {
  set.seed(14)
  np <- 12; nt <- 5
  b0 <- rnorm(np, log(10), 0.55); b1 <- rnorm(np, 0.20, 0.09)
  dat <- do.call(rbind, lapply(seq_len(np), function(i) {
    tt <- sort(runif(nt, 0, 4))
    data.frame(id = factor(i), t = tt,
               logV = b0[i] + b1[i]*tt + rnorm(nt, 0, 0.10))
  }))
  lmm  <- nlme::lme(logV ~ t, random = ~ t | id, data = dat,
                    control = nlme::lmeControl(opt = "optim"))
  pooled <- lm(logV ~ t, data = dat)
  print(knitr::kable(data.frame(
    model = c("mixed effects (random intercept and slope)", "pooled ordinary regression"),
    slope = round(c(nlme::fixef(lmm)["t"], coef(pooled)["t"]), 4),
    std_error = round(c(sqrt(diag(vcov(lmm)))["t"],
                        summary(pooled)$coefficients["t", 2]), 4),
    truth = 0.20),
    row.names = FALSE,
    caption = "Both models estimate the population slope well; only the mixed model reports an honest standard error, because the pooled model treats 60 correlated observations as 60 independent ones."))
  cat(sprintf("standard-error ratio (pooled / mixed) = %.2f\n",
              summary(pooled)$coefficients["t", 2]/sqrt(diag(vcov(lmm)))["t"]))

  pop <- data.frame(t = seq(0, 4, length.out = 100))
  pop$V <- exp(nlme::fixef(lmm)[1] + nlme::fixef(lmm)[2]*pop$t)
  m1 <- ggplot(dat, aes(t, exp(logV))) +
    geom_line(aes(group = id), alpha = 0.45, colour = bpad_pal[1]) +
    geom_point(alpha = 0.6, colour = bpad_pal[1], size = 1.3) +
    geom_line(data = pop, aes(t, V), colour = bpad_pal[2], linewidth = 1.1) +
    scale_y_log10() +
    labs(title = "individual trajectories and population mean",
         x = "years", y = expression("V (cm"^3*", log)"))
  m2 <- ggplot(dat, aes(t, logV)) +
    geom_point(aes(colour = id), alpha = 0.7, size = 1.6, show.legend = FALSE) +
    geom_abline(intercept = coef(pooled)[1], slope = coef(pooled)[2],
                linetype = "dashed", colour = "grey25", linewidth = 0.9) +
    geom_abline(intercept = nlme::fixef(lmm)[1], slope = nlme::fixef(lmm)[2],
                colour = bpad_pal[2], linewidth = 0.9) +
    labs(title = "pooled (dashed) versus mixed-effects (solid) fit",
         x = "years", y = "log V")
  bpad_grid(m1, m2, ncol = 2)
} else {
  cat("nlme unavailable; mixed-effects demonstration skipped\n")
}

## 
## 
## Table: (\#tab:mixed-effects)Both models estimate the population slope well; only the mixed model reports an honest standard error, because the pooled model treats 60 correlated observations as 60 independent ones.
## 
## |model                                      |  slope| std_error| truth|
## |:------------------------------------------|------:|---------:|-----:|
## |mixed effects (random intercept and slope) | 0.2533|    0.0274|   0.2|
## |pooled ordinary regression                 | 0.1626|    0.0654|   0.2|
## standard-error ratio (pooled / mixed) = 2.39

Why repeated measurements need a random-effects structure. Simulated longitudinal log-volume trajectories for twelve patients with a common population slope and substantial between-patient variation in both intercept and slope. Left: the individual trajectories with the population mean. Right: pooling all points into one ordinary regression (dashed) yields a similar slope but a standard error roughly half as large as the mixed model, so the pooled analysis is overconfident by a factor of about two.

Figure 28: Why repeated measurements need a random-effects structure. Simulated longitudinal log-volume trajectories for twelve patients with a common population slope and substantial between-patient variation in both intercept and slope. Left: the individual trajectories with the population mean. Right: pooling all points into one ordinary regression (dashed) yields a similar slope but a standard error roughly half as large as the mixed model, so the pooled analysis is overconfident by a factor of about two.

A joint model goes one step further, linking the longitudinal submodel of Eq. (49) to a survival submodel through the shared random effects,

\[\begin{equation} h_i(t) = h_0(t)\exp\{X_i^{\mathsf T}\gamma + \alpha\, m_i(t)\}, \tag{50} \end{equation}\]

where $m_i(t) = \beta_0 + \beta_1 t + b_{0i} + b_{1i}t$ is the patient’s underlying (noise-free) trajectory and $\alpha$ measures how strongly the current biomarker level drives the hazard. Joint models are the correct tool when the biomarker is measured with error and when its measurement schedule is itself informative, both of which hold for imaging surveillance: patients are scanned more often precisely when there is concern.

Section 8.8 summary.

Survival, hazard, and cumulative hazard are linked by $S = e^{-H}$, a first-order ODE solution from Chapter 1; a binary endpoint discards the follow-up that makes these estimable.
Kaplan–Meier handles censoring without discarding patients; late estimates are imprecise because the risk set shrinks.
Cox hazard ratios are relative and time-transformation invariant; the proportional-hazards assumption must be tested, not assumed.
Radical nephrectomy looks harmful unadjusted and neutral adjusted, because tumour size is both the indication and the prognostic factor: confounding by indication.
Immortal time creates a large apparent benefit from an exposure with no effect; a landmark analysis recovers the truth.
Under competing events $1-\mathrm{KM}$ overstates absolute risk substantially; use the cumulative incidence function.
Exponential, logistic, and Gompertz growth agree over typical surveillance intervals and diverge later, so growth-model choice is not identifiable from short follow-up.
Because $d\propto V^{1/3}$, diameter-based progression criteria operate near the segmentation-noise floor; repeated measurements require random effects for honest standard errors.

Checkpoint 8.8. A surveillance study reports that patients whose tumours were re-imaged within six months had better five-year survival, and concludes that closer surveillance improves outcome. Name the bias, describe the landmark analysis you would run instead, and state what result would persuade you that the effect is real.

Interim Summary: Sections 8.1–8.8

Eight sections in, the through-line is a single claim: the hard part of biomedical imaging AI is not the algorithm. Every quantitative result in this instalment supports it.

The pipeline is a composition of physical, computational, and statistical operators, Eq. (4), and prediction error decomposes across all of them.
Screening features on all the data before cross-validating turned pure noise into an AUC of 0.82 with a small standard deviation. Leakage does not look like noise; it looks like a good result.
Real clinical records contain truncated laboratory values, semantic string sentinels, deterministically structured missingness, and three near-identically-named surgical fields whose confusion produces a silently constant predictor.
Mask-derived tumour geometry agrees with the recorded clinical size to about 0.2 cm, which validates the physics pipeline; a one-voxel boundary shift nonetheless moves the volume biomarker by roughly $\pm25\%$, which is the reproducibility floor for everything downstream.
Texture features change by orders of magnitude with the discretization setting alone.
The development and external strata of a single public cohort differ in both slice thickness and tumour size, so technical and clinical shift are entangled.
The malignancy classifier has modest discrimination and essentially no Brier skill against the base rate. That is a finding about the cohort, not a failure of the code, and no model family in the comparison escaped it.
The eGFR-decline model does have real skill, because its outcome is mechanistically determined by a recorded intervention, and its coefficient magnitude matches a physiological prediction.
Radical nephrectomy is associated with worse survival unadjusted and with nothing once tumour size is controlled: confounding by indication, visible in two lines of output.
Immortal time manufactured a large protective hazard ratio from an exposure simulated to have no effect; $1-\mathrm{KM}$ overstated cause-specific risk by more than $80\%$ under competing events.

Decision guide for Sections 8.1–8.8.

Outcome is binary and complete $\rightarrow$ logistic regression; report AUC and calibration and net benefit (8.7.1).
Outcome is continuous $\rightarrow$ linear model with residual diagnostics; check whether the analysis population is the whole cohort (8.7.6).
Outcome is time-to-event $\rightarrow$ Kaplan–Meier and Cox; test proportional hazards (8.8.2, 8.8.3).
Competing events present $\rightarrow$ cumulative incidence function, never $1-\mathrm{KM}$ (8.8.6).
Exposure occurs after the time origin $\rightarrow$ landmark analysis (8.8.5).
Predictors are many and correlated $\rightarrow$ ridge for prediction, grouped representatives for interpretation; never invert $X^{\mathsf T}X$ (8.2.3, 8.6.3).
Features come from masks $\rightarrow$ quantify the one-voxel perturbation band before claiming a difference (8.5.6).
Sites, scanners, or protocols differ $\rightarrow$ measure the shift with standardized mean differences before modelling it away (8.6.6).
A hyperparameter is tuned $\rightarrow$ nested resampling (8.3.5, 8.7.7).
A variable was chosen because of the imaging $\rightarrow$ it cannot be a predictor in a model that informs the same choice (8.7.9).

Where this instalment connects to the rest of BPAD

Chapter 1 supplies every mathematical tool used here: the design-matrix algebra and the SVD behind Sections 8.2.3 and 8.6.4; the gradient and chain rule behind Eqs. (15) and (16); Bayes’ theorem behind Eq. (18); the delta method and Monte Carlo propagation behind Sections 8.2.9 and 8.5.6; the Gaussian, Poisson, and Rician noise models of Section 8.2.10; the sampling theorem behind Section 8.5.8; and the growth ODEs of Section 8.8.7.
Chapter 2 contributes the pulse-oximetry pigmentation bias as the clearest physically grounded example of a fairness problem with a wavelength-level explanation, which is the template for the device-equity discussion to come.
Chapters 3 through 6 determine the noise model, the contrast mechanism, and the irreducible error floor $\sigma^2$ for each modality, and therefore what any model can achieve.
Chapter 5 supplies the dose–noise trade-off that bounds low-dose CT model performance, and dual-energy material decomposition as the multi-channel unmixing problem of Section 8.2.2.
Chapter 6 supplies the Poisson likelihood that makes MLEM the statistically correct reconstruction, the same likelihood-matching argument made in Section 8.2.2.
Chapter 7 supplies every feature this chapter consumes: registration, segmentation, filtering, and quantification, together with the point spread function that Section 8.5.5 treats as a texture confound.
Sections 8.9 onward take up unsupervised learning and segmentation, deep learning and image-native AI, evaluation and clinical utility in depth, translation and deployment, and reproducible computational practice, all on this same cohort.

Problems for Sections 8.1–8.8

Problem 8.1. Effective sample size

A study uses 250 slices from each of 160 patients, split at the slice level, and reports a 95% confidence interval computed from $n = 40{,}000$. Assuming an intraclass correlation of 0.85, use Eq. (6) to compute $n_{\mathrm{eff}}$, and give the factor by which the reported interval is too narrow.

Problem 8.2. Conditioning

A radiomic design matrix has singular values spanning $\sigma_1 = 42.0$ to $\sigma_p = 0.021$. (a) Compute $\kappa(X)$ and $\kappa(X^{\mathsf T}X)$. (b) Given double-precision machine epsilon of $2.2\times10^{-16}$, roughly how many significant digits are lost by solving the normal equations rather than using QR? (c) Name the alternative you would use and the tuning decision it introduces.

Problem 8.3. Bias, variance, and the noise floor

For the demonstration in Section 8.2.7, the irreducible noise variance is $\sigma^2 = 0.09$. Using the tabulated bias$^2$ and variance at degrees 1, 3, and 12, verify that bias$^2$ + variance + $\sigma^2$ approximately reproduces the test MSE plus $\sigma^2$, and state which term dominates at each degree.

Problem 8.4. Predictive value across settings

A malignancy model has $Se = 0.85$ and $Sp = 0.75$. Compute PPV and NPV at prevalences 0.914 (this chapter’s surgical series), 0.20 (incidental small masses), and 0.05 (a hypothetical screening setting). State which of the four quantities is unchanged across the three settings and why.

Problem 8.5. Rician bias

A diffusion-weighted acquisition has $\sigma = 4$ arbitrary units. (a) Compute the expected magnitude signal in a voxel whose true signal is zero. (b) Explain why fitting a mono-exponential decay to magnitude data at high $b$-value systematically overestimates the apparent diffusion coefficient’s asymptote. (c) Name the section of Chapter 4 that treats the correction.

Problem 8.6. Diagnosing leakage

For each pipeline, state whether it leaks and, if so, exactly what information crosses the boundary: (a) impute missing values with the column median computed on all rows, then split; (b) standardize inside each training fold, then apply to the validation fold; (c) choose the classification threshold to maximize accuracy on the locked test set; (d) select the top 20 of 900 features by univariate $p$-value inside each training fold; (e) run five model families, report the best test AUC.

Problem 8.7. Events per predictor

The cohort has 21 deaths. Using the 10-to-20 events-per-predictor guidance, state the maximum number of covariates for a Cox model, and explain why the Section 8.8.3 model’s reported events-per-parameter figure is a warning rather than a reassurance.

Problem 8.8. Structured missingness

Ischaemia time is "not_applicable" for all 70 radical nephrectomies and numeric for 130 of 140 partial nephrectomies. (a) Classify this missingness. (b) Explain why median imputation is indefensible. (c) Explain why including ischaemia time as a predictor in a model that also contains the procedure indicator is problematic even after correct handling.

Problem 8.9. Voxel geometry

A label volume has dimensions $611\times512\times512$ with pixdim $= (0.5, 0.92, 0.92)$ mm, and the tumour label occupies 18,450 voxels. (a) Compute the tumour volume in cm$^3$. (b) Compute the equivalent spherical diameter. (c) Which axis is the slice axis, and how did you determine it from the header alone? (d) If the metadata spreadsheet instead lists spacing as $(0.92, 0.92, 0.5)$, what error results from using the spreadsheet order with the array dimensions?

Problem 8.10. Propagating the segmentation floor

A one-voxel dilation increases a tumour volume from 7.23 to 9.16 cm$^3$. (a) Compute the relative change in volume and in equivalent diameter. (b) Verify the relation $\sigma_d/d = \tfrac{1}{3}\sigma_V/V$. (c) A trial defines progression as a 20% diameter increase. Comment on the adequacy of that threshold.

Problem 8.11. Texture reproducibility

Using the table in Section 8.5.5, compute the ratio of GLCM contrast at 128 grey levels to that at 8 levels. State what must be reported for a texture feature to be reproducible, and explain why harmonizing voxel spacing alone is insufficient.

Problem 8.12. Collinearity remedies

A design contains tumour volume, surface area, equivalent diameter, and bounding-box diagonal, with VIFs of 41, 33, 58, and 29. (a) Explain why lasso would produce an unstable feature set. (b) Propose a composite index and state what it measures. (c) State which remedy you would choose if the goal were a publishable odds ratio for tumour size, and which if the goal were the best possible risk prediction.

Problem 8.13. PCA and supervision

PC1 of the development design explains a large share of variance and loads mainly on size and acquisition variables, yet the malignancy classes are not separated in PC space. (a) Explain why this is expected. (b) State a situation in which discarding PC12 would discard the entire signal. (c) Name two supervised alternatives.

Problem 8.14. Copula invariance

Show algebraically that Spearman’s rank correlation is invariant under any strictly increasing transformation of either variable, and state one imaging situation in which this invariance is a practical advantage.

Problem 8.15. Net benefit

At threshold $t = 0.90$, a model flags 100 of 161 patients, of whom 95 are truly malignant. (a) Compute net benefit using Eq. (39). (b) Compute net benefit for “treat everyone” at the same threshold. (c) Interpret the comparison for a surgeon.

Problem 8.16. Hazard and survival

A cohort has a constant hazard of 0.08 per year. (a) Write $S(t)$. (b) Compute five-year survival. (c) Compute the median survival time. (d) If the hazard instead doubles after year 3, sketch $S(t)$ and state whether a Cox model with a time-fixed covariate could represent this.

Problem 8.17. Confounding by indication

Using the Section 8.8.4 output, (a) state the unadjusted and fully adjusted hazard ratios for radical nephrectomy; (b) explain the direction of the change in terms of the confounding structure; (c) describe the study design that would answer the causal question, and one reason it may be infeasible.

Problem 8.18. Immortal time

In the Section 8.8.5 simulation the true hazard ratio is exactly 1. (a) Report the naive and landmark estimates. (b) Explain, in terms of the risk sets, why the naive estimate is biased downward. (c) State the cost of the landmark approach.

Problem 8.19. Competing risks

With cause-1 hazard 0.06 and competing hazard 0.20 per year, (a) compute the true cause-1 cumulative incidence at 5 years using $F_1(t) = \frac{h_1}{h_1+h_2}\{1-e^{-(h_1+h_2)t}\}$; (b) explain why $1-e^{-h_1 t}$ is larger; (c) state which estimator a clinician quoting five-year risk should use.

Problem 8.20. Growth models and doubling time

A renal mass measures 2.5 cm at baseline and 2.9 cm at 18 months. (a) Convert both to volumes assuming sphericity. (b) Compute the specific growth rate and volume doubling time. (c) Using a $\pm8\%$ diameter uncertainty, propagate to an interval for the doubling time and comment on whether growth has been demonstrated.

Solutions to Selected Problems

Solution 8.1

$n_{\mathrm{eff}} = nm/\{1+(m-1)\rho\} = 160\times250/\{1+249\times0.85\} = 40{,}000/212.65 = 188$. Because interval width scales as $n^{-1/2}$, the reported interval is too narrow by $\sqrt{40{,}000/188} = \sqrt{212.6} = 14.6$, a factor of nearly fifteen. A nominal 95% interval is therefore closer to a 20% interval, and every $p$-value from that analysis is meaningless.

Solution 8.2

$\kappa(X) = 42.0/0.021 = 2000$; $\kappa(X^{\mathsf T}X) = 4\times10^{6}$.
Digits lost are approximately $\log_{10}\kappa$: about 3.3 digits for QR on $X$, and about 6.6 for the normal equations. With roughly 16 available digits, the normal equations leave about 9 reliable digits, which is survivable here but fails outright once $\kappa(X)$ reaches $10^{8}$.
Truncated SVD, Eq. (12), or ridge regression. Both introduce a threshold ($\tau$ or $\lambda$) that is a hyperparameter and must be chosen inside resampling (Section 8.7.7).

Solution 8.4

Using Eq. (18) with $Se = 0.85$, $Sp = 0.75$:

Prevalence	PPV	NPV
0.914	0.973	0.398
0.20	0.459	0.949
0.05	0.152	0.990

The unchanged quantities are sensitivity and specificity, and therefore the AUC, because all three are computed conditionally on the true class and never involve $\pi$. PPV falls by a factor of six across the three settings while the AUC is constant, which is precisely the reporting gap identified in Worked Example 8.1.

Solution 8.5

$\mathbb{E}[M] = \sigma\sqrt{\pi/2} = 4\times1.2533 = 5.01$ units, although the true signal is zero.
At high $b$-value the true signal approaches the noise floor, so the measured magnitude flattens at approximately $\sigma\sqrt{\pi/2}$ rather than decaying to zero. A mono-exponential fit interprets that floor as residual signal, which flattens the fitted decay and biases the apparent diffusion coefficient downward while inflating the estimated asymptote.
Chapter 4’s treatment of Rician noise and noise-floor correction in magnitude diffusion data.

Solution 8.6

Leaks: the median is a function of the held-out rows, so information about the test distribution enters the training features.
Does not leak: this is the correct order, and it is what fit_preprocessor()/apply_preprocessor() implement.
Leaks: threshold selection is a fitting decision, and performing it on the locked test set converts that set into a tuning set.
Does not leak: screening is inside the fold. This is the “nested” arm of Section 8.3.6, which correctly returned AUC $\approx 0.5$ on pure noise.
Leaks: selecting the maximum over five families using test performance is multiplicity on the test set. The reported figure is a maximum of five noisy estimates and is biased upward; report all five, or reserve a further partition for the selection.

Solution 8.8

It is neither MCAR, MAR, nor MNAR in the usual sense: the value is undefined by design, a deterministic function of the procedure, because no vascular clamp is applied during a radical nephrectomy.
Median imputation would assign a clamping time to an operation in which no clamp existed, fabricating a physical quantity. The imputed value would then correlate perfectly with the procedure and act as a proxy for it.
Even handled correctly, ischaemia time is available only within partial nephrectomies, so a model containing both it and the procedure indicator is fitting an interaction that is estimable in only one stratum. The defensible choices are to model the partial-nephrectomy subgroup separately, or to drop the variable from any model spanning both procedures.

Solution 8.9

$V = 18{,}450\times0.5\times0.92\times0.92 = 7806$ mm$^3 = 7.81$ cm$^3$.
$d_{\mathrm{eq}} = 2(3\times7806/4\pi)^{1/3} = 2\times12.30 = 24.6$ mm $= 2.46$ cm.
Axis 1: it is the anisotropic axis, since axes 2 and 3 share a spacing of 0.92 mm while axis 1 differs at 0.5 mm. The code identifies it as which.max(abs(sp - median(sp))) rather than assuming axis 3.
Using the spreadsheet order pairs a spacing of 0.92 mm with the 611-element axis and 0.5 mm with a 512-element axis. The voxel volume is unchanged (the product is the same), so the volume would still be right, but every axis-specific quantity, including the slice-area profile, the slice thickness, the bounding-box extents, and the in-plane pixel area used for texture, would be wrong. This is why geometry and dimensions must come from the same header.

Solution 8.10

$\Delta V/V = (9.16-7.23)/7.23 = +26.7\%$. Diameters are $2(3V/4\pi)^{1/3}$, giving 2.399 cm and 2.596 cm, so $\Delta d/d = +8.2\%$.
$\tfrac{1}{3}\times26.7\% = 8.9\%$, close to the observed 8.2%; the small discrepancy is the second-order term omitted by the first-order delta method, and it shrinks as the perturbation shrinks.
A 20% diameter threshold corresponds to a 73% volume increase ($1.20^3 = 1.73$), which is well above the $\pm8\%$ diameter noise band, so the threshold is defensible for a single reader on a single protocol. It is not defensible across protocols, because Section 8.6.6 showed that slice thickness shifts the geometry systematically rather than randomly, and a systematic shift is not averaged away by a threshold.

Solution 8.15

TP $= 95$, FP $= 5$, $n = 161$, $t = 0.90$, so $t/(1-t) = 9$. Then $\mathrm{NB} = 95/161 - (5/161)\times9 = 0.590 - 0.280 = 0.310$.
Treating everyone flags all 161, of whom $0.914\times161 = 147$ are malignant and 14 are not: $\mathrm{NB} = 147/161 - (14/161)\times9 = 0.913 - 0.783 = 0.130$.
At a threshold implying that nine unnecessary operations are as bad as one missed cancer, the model’s net benefit (0.310) exceeds operating on everyone (0.130). At a lower threshold, where a missed cancer is weighted far more heavily, treat-everyone wins and the model adds nothing. The clinical question is therefore not “is the model accurate?” but “what is my exchange rate?”, and the decision curve makes that explicit.

Solution 8.19

$F_1(5) = \frac{0.06}{0.26}\{1-e^{-0.26\times5}\} = 0.2308\times(1-0.2725) = 0.168$.
$1-e^{-0.06\times5} = 1-0.741 = 0.259$, which is 54% larger. The naive expression assumes every patient remains at risk of cause 1 for the whole five years, whereas in reality many die of the competing cause first and can never experience cause 1.
The cumulative incidence function, Eq. (45), estimated by Aalen–Johansen, or equivalently a Fine–Gray model. Quoting $1-\mathrm{KM}$ to a patient overstates their absolute cause-specific risk, and in an elderly surgical cohort the overstatement can exceed a factor of two.

Solution 8.20

$V_1 = \frac{4}{3}\pi(1.25)^3 = 8.18$ cm$^3$; $V_2 = \frac{4}{3}\pi(1.45)^3 = 12.77$ cm$^3$.
$k = \ln(12.77/8.18)/547.5\ \text{d} = 0.4453/547.5 = 8.13\times10^{-4}$ per day, so $\mathrm{VDT} = \ln 2/k = 852$ days.
A $\pm8\%$ diameter uncertainty is a $\pm24\%$ volume uncertainty on each measurement. Propagating both into $k = \ln(V_2/V_1)/\Delta t$ gives $\sigma_k \approx \sqrt{2}\times0.24/547.5 = 6.2\times10^{-4}$ per day, comparable to $k$ itself. The 95% interval for $k$ therefore includes zero, and the doubling time is unbounded above. Growth has not been demonstrated: two measurements 18 months apart with this measurement uncertainty cannot distinguish an 852-day doubling time from no growth at all. A third timepoint, or a volumetric rather than diameter-based measurement, is required.

Continue to part 2 of BPAD Chapter 8.

References for Sections 8.1–8.8

Hastie, T., Tibshirani, R., and Friedman, J. (2009). The Elements of Statistical Learning, 2nd ed. Springer.
James, G., Witten, D., Hastie, T., and Tibshirani, R. (2021). An Introduction to Statistical Learning, 2nd ed. Springer.
Steyerberg, E. W. (2019). Clinical Prediction Models, 2nd ed. Springer.
Riley, R. D. et al. (2020). Calculating the sample size required for developing a clinical prediction model. BMJ, 368, m441.
Collins, G. S. et al. (2024). TRIPOD+AI statement: updated guidance for reporting clinical prediction models that use regression or machine learning methods. BMJ, 385, e078378.
Mongan, J., Moy, L., and Kahn, C. E. (2020). Checklist for Artificial Intelligence in Medical Imaging (CLAIM). Radiology: Artificial Intelligence, 2(2), e200029.
Zwanenburg, A. et al. (2020). The Image Biomarker Standardisation Initiative. Radiology, 295(2), 328–338.
Lambin, P. et al. (2017). Radiomics: the bridge between medical imaging and personalized medicine. Nature Reviews Clinical Oncology, 14(12), 749–762.
Heller, N. et al. (2021). The KiTS19 challenge data and results. Medical Image Analysis, 67, 101821.
Bland, J. M. and Altman, D. G. (1986). Statistical methods for assessing agreement between two methods of clinical measurement. Lancet, 327(8476), 307–310.
Therneau, T. M. and Grambsch, P. M. (2000). Modeling Survival Data: Extending the Cox Model. Springer.
Putter, H., Fiocco, M., and Geskus, R. B. (2007). Tutorial in biostatistics: competing risks and multi-state models. Statistics in Medicine, 26(11), 2389–2430.
Fine, J. P. and Gray, R. J. (1999). A proportional hazards model for the subdistribution of a competing risk. JASA, 94(446), 496–509.
Suissa, S. (2008). Immortal time bias in pharmacoepidemiology. American Journal of Epidemiology, 167(4), 492–499.
Vickers, A. J. and Elkin, E. B. (2006). Decision curve analysis: a novel method for evaluating prediction models. Medical Decision Making, 26(6), 565–574.
Rizopoulos, D. (2012). Joint Models for Longitudinal and Time-to-Event Data. CRC Press.
Nelsen, R. B. (2006). An Introduction to Copulas, 2nd ed. Springer.
Gudbjartsson, H. and Patz, S. (1995). The Rician distribution of noisy MRI data. Magnetic Resonance in Medicine, 34(6), 910–914.
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 and 7.

Symbol	Meaning
\(X \in \mathbb{R}^{n\times p}\)	design matrix, \(n\) observations by \(p\) predictors
\(y,\ \widehat y\)	outcome and prediction
\(\theta,\ \beta\)	model parameters; \(\beta\) for linear predictors
\(\ell(y,\widehat y),\ \Omega(\theta),\ \lambda\)	loss, penalty, penalty strength
\(\sigma(\eta) = (1+e^{-\eta})^{-1}\)	logistic (sigmoid) link
\(\kappa(A)\)	condition number of \(A\)
\(I \in \mathbb{R}^{n_x\times n_y\times n_z}\)	image volume (tensor)
\(M\)	segmentation mask; \(M_v \in \{0,1,2\}\) for background, kidney, tumour
\(z = g(I,M;\psi)\)	extracted features, with extraction choices \(\psi\)
\(\Delta_{\mathrm{acq}},\Delta_{\mathrm{recon}},\Delta_{\mathrm{seg}},\Delta_{\mathrm{model}},\Delta_{\mathrm{shift}}\)	error contributions by pipeline stage
\(D,\ T\)	development and deployment distributions
\(R_D(f),\ R_T(f)\)	risk under development and deployment
\(\widetilde T = \min(T,C),\ \delta = \mathbb{I}(T\le C)\)	observed time and event indicator
\(S(t),\ h(t),\ H(t)\)	survival, hazard, cumulative hazard
\(F_k(t)\)	cause-specific cumulative incidence function
\(\pi,\ Se,\ Sp\)	prevalence, sensitivity, specificity
\(V(t),\ K,\ k\)	tumour volume, carrying capacity, growth rate

Task	Loss \(\ell\)	Typical penalty \(\Omega\)
Linear regression	squared error \((y-\widehat y)^2\)	none, ridge, lasso
Logistic regression	binary cross-entropy	none, ridge, lasso
Support vector machine	hinge \(\max(0, 1-y\widehat y)\)	squared norm
Cox model	negative partial log-likelihood	ridge, lasso
Neural classifier	cross-entropy	weight decay, dropout, augmentation
Segmentation network	cross-entropy, Dice, or a sum	weight decay, augmentation

Model	Statement	Signal dependence	Modality
Gaussian	\(X = \mu + \epsilon\), \(\epsilon\sim\mathcal{N}(0,\sigma^2)\)	none	thermal and amplifier noise (CT, optical)
Poisson	\(N\sim\mathrm{Poisson}(\lambda)\), \(\mathbb{E}=\mathrm{Var}=\lambda\)	\(\sigma = \sqrt{\lambda}\)	photon counting (PET, SPECT, low-dose CT)
Rician	\(M = \sqrt{S_{\mathrm{re}}^2 + S_{\mathrm{im}}^2}\)	biased high at low SNR	magnitude MRI

Biomedical Physics with Applications to Disease (BPAD)

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

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

August 05, 2026