5  Describing, Resampling, Validating, and Simplifying the Model

5.1 Describing the Fitted Model

5.1.1 Interpreting Effects

  • Regression coefficients if 1 d.f. per factor, no interaction A
  • Not standardized regression coefficients
  • Many programs print meaningless estimates such as effect of increasing age\(^2\) by one unit, holding age constant
  • Need to account for nonlinearity, interaction, and use meaningful ranges
  • For monotonic relationships, estimate \(X\hat{\beta}\) at quartiles of continuous variables, separately for various levels of interacting factors
  • Subtract estimates, anti-log, e.g., to get inter-quartile-range B odds or hazards ratios. Base C.L. on s.e. of difference. See Figure 21.4.
  • Partial effect plot: Plot effect of each predictor on \(X\beta\) or some transformation. See Figure 21.2. See also Karvanen & Harrell (2009). Fox and Weisberg have other excellent approaches to displaying predictor effects
  • Nomogram. See Figure 21.5
  • Use regression tree to approximate the full model


Forming contrasts through differences in predicted values is the most general and easy way to estimate effects of interest in a model. This is because

  • One needs only to specify two or more sets of predictor settings to contrast
  • Predictor coding, nonlinear terms, and interactions are automatically taken into account
  • Unspecified predictors can default to reference values such as medians for continuous predictors and modes for categorical ones. The values used for unspecified predictors only matter if they interact with predictors that are changing in the contrast.
  • Approximate methods are available for getting simultaneous confidence intervals over a whole sequence of contrasts (e.g., simultaneous confidence bands for an entire range of ages) when contrasts produce more than one estimate
  • The approach extends to any order of differences. E.g., interaction contrasts are formed with double differences.
  • Contrasts allow assessment of linearity over a given predictor range in a spline fit
  • Contrasts allow assessment of “partial interactions”, i.e., whether there is an interaction between one predictor and a subset of the levels of another predictor
  • The R predict function provides the \(X\) (design matrix) that goes along with each requested prediction, and simple matrix algebra provides needed covariance matrices once differences in design matrices are formed

On the latter point, a single difference contrast such as one comparing 30 year old males with 40 year old males would generate two design matrices \(X_{1}, X{2}\), and the contrast \(X_{1}\beta - X_{2}\beta = (X_{1} - X_{2})\beta\).

Consider a series of examples in which predicted values from a fitted model are symbolized by a function \(f\) of the predictor settings, with continuous predictors \(A\) and \(B\) and a categorical predictor \(C\) with levels \(i, j, k\). Let the reference (detault) values of the predictors be \(a\) (median of \(A\)), \(b\) (median of \(B\)), and \(c\) (modal (most frequent) category of \(C\)). The term “average slope” used below will be the slope if the predictor has a linear effect, otherwise it is the average slope over the indicated interval. See here for more information.

Single Difference Contrasts
  • Average slope of \(A\) over \(A \in [1, 2]\): \(f(2, b, c) - f(1, b, c)\)
  • Difference between categories \(j\) and \(k\) in \(C\): \(f(a, b, j) - f(a, b, k)\)
  • Effect of an entire range of predictor values against a single reference value: \(f(1:10, b, c) - f(5, b, c)\) (a simultaneous confidence band would be useful here)
Double Difference Contrasts
  • Linearity in \(A\) over \(A \in [1, 3]\): \(f(3, b, c) - f(2, b, c) - [f(2, b, c) - f(1, b, c)] = f(3, b, c) - 2f(2, b, c) + f(1, b, c)\)
  • Interaction between \(A\) and \(C\) with respect to only two levels of \(C\): \(f(2, b, j) - f(1, b, j) - [f(2, b, k) - f(1, b, k)]\)
  • Interaction between \(A\) and \(B\) when \(A\) varies from 1 to 3 and \(B\) varies from 5 to 6: \(f(3, 5, c) - f(1, 5, c) - [f(3, 6, c) - f(1, 6, c)]\)
  • Effect of \(C=i\) vs. \(C=j\) over a regular grid of values of \(B\): \(f(a, 5:15, i) - f(a, 5:15, j)\) (with possible simultaneous band; this is often a better way of getting hazard ratios in a Cox proportional hazards model, as the user has control of the reference value)
Triple Difference Contrasts
  • Third-order interactions
  • Assessment of nonlinearity in the \(A\times B\) interaction (differences in double differences)

The R rms package contrast function (full name contrast.rms) provides predicted values, differences, and double difference contrasts. The latter is specified by providing four lists of predictor settings.

5.1.2 Indexes of Model Performance

Error Measures

  • Central tendency of prediction errors C
    • Mean absolute prediction error: mean \(|Y - \hat{Y}|\)
    • Mean squared prediction error
      • Binary \(Y\): Brier score (quadratic proper scoring rule)
    • Logarithmic proper scoring rule (avg. log-likelihood)
  • Discrimination measures D
    • Pure discrimination: rank correlation of \((\hat{Y}, Y)\)

      • Spearman \(\rho\), Kendall \(\tau\), Somers’ \(D_{xy}\)
      • \(Y\) binary \(\rightarrow\) \(D_{xy} = 2\times (C - \frac{1}{2})\)
        \(C\) = concordance probability = area under receiver operating characteristic curve \(\propto\) Wilcoxon-Mann-Whitney statistic
    • Mostly discrimination: \(R^{2}\)

      • \(R^{2}_{\mathrm{adj}}\)—overfitting corrected if model pre-specified
    • Brier score can be decomposed into discrimination and calibration components

    • Discrimination measures based on variation in \(\hat{Y}\)

      • regression sum of squares
      • \(g\)–index
      • see here for related information
  • Calibration measures E
    • calibration–in–the–large: average \(\hat{Y}\) vs. average \(Y\)
    • high-resolution calibration curve (calibration–in–the–small). See ?fig-titanic-calibrate.
    • calibration slope and intercept
    • maximum absolute calibration error
    • mean absolute calibration error
    • 0.9 quantile of calibration error

See Van Calster et al. (2016) for a nice discussion of different levels of calibration stringency and their relationship to likelihood of errors in decision making.


  • Based on Gini’s mean difference F
    • mean over all possible \(i \neq j\) of \(|Z_{i} - Z_{j}|\)
    • interpretable, robust, highly efficient measure of variation
  • \(g =\) Gini’s mean difference of \(X_{i}\hat{\beta} = \hat{Y}\)
  • Example: \(Y=\) systolic blood pressure; \(g = 11\)mmHg is typical difference in \(\hat{Y}\)
  • Independent of censoring etc.
  • For models in which anti-log of difference in \(\hat{Y}\) G represent meaningful ratios (odds ratios, hazard ratios, ratio of medians):
    \(g_{r} = \exp(g)\)
  • For models in which \(\hat{Y}\) can be turned into a probability estimate (e.g., logistic regression):
    \(g_{p} =\) Gini’s mean difference of \(\hat{P}\)
  • These \(g\)–indexes represent e.g. “typical” odds ratios, “typical” risk differences
  • Can define partial \(g\)

5.1.3 Relative Explained Variation

  • In OLS linear models \(R^2\) is an excellent measure of explained variation in \(Y\)
  • Fraction of variance of \(Y\) that is explained by \(X\)
  • SSR = \(\sum_{i=1}^{n} (\hat{Y}_{i} - \bar{\hat{Y}})^2\) = regression sum of squares
    • proportional to variance of \(\hat{Y}\)
  • SSE = \(\sum_{i=1}^{n} (Y_{i} - \hat{Y}_{i})^2\) = error sum of squares = residual sum of squares
    • proportional to error variance
  • SST = \(\sum_{i=1}^{n} (Y_{i} - \bar{Y})^2\) = total sum of squares
    • proportion to variance of \(Y\)
  • SST = SSR + SSE
  • \(R^{2} =\) SSR/SST = 1 - SSE/SST
  • For each part of the model can compute its partial SSR = loss in SSR after refitting the model upon removing the part of the model being assessed
  • Partial \(R^{2}\) = (partial SSR)/SST
  • Relative partial \(R^{2} = \frac{\text{partial } R^{2}}{R^{2}} =\) partial SSR divided by full model SSR
  • Relative partial \(R^{2}\) = portion of explainable (by \(X\)) variation in \(Y\) that is explained by one or more terms in the model
  • Relative explained variation: REV
  • REVs do not need to sum to 1.0
    • collinearities
    • computing REVs for overlapping model terms
  • When a predictor interacts with another predictor, the predictor’s REV is taken from the combined predictive information in
    • the predictor main effect (including any nonlinear terms)
    • any interaction terms involving that effect

Extension to Other Models

  • For models fitted using maximum likelihood estimation, can compute the partial likelihood ratio \(\chi^2\) statistic for each set and divide that by the total model LR \(\chi^2\)
  • See also Chapter 9 and the example below
  • This is an excellent measure but requires more computation and doesn’t extend to models not having a likelihood or to Bayesian modeling
  • Need a general purpose REV measure that works on virtually all models using all estimation techniques
  • Letting \(\hat{Y} = X\hat{\beta}\) be the estimated linear predictor from the model
  • Use the same design matrix \(X\) used to get \(\hat{Y}\)
  • Predict \(\hat{Y}\) from \(X\) using OLS
  • The \(R^{2}\) from this model is 1.0
  • SSE = 0; ignore this and take variance-covariance matrix as \((X'X)^{-1}\) since SSR is used only on a relative basis
  • REV = \(\frac{\text{partial SSR }}{\text{SSR}}\)
  • Identical to relative \(R^2\) if original model was OLS
  • Easily lends itself to uncertainty intervals
    • frequentist: bootstrap the original model fit and save the repeated coefficient estimates (done by R rms package bootcov function); compute bootstrap confidence intervals for REVs; see Section 5.2 and example below
    • Bayesian: compute REV for each posterior draw of \(\beta\) and from the repeated REVs compute highest posterior density uncertainty intervals
  • Implemented in R rms package rexVar function

See Luchman (2014) for a mice review of log-likelihood-based relative importance measures.


  • Ordinal regression model (proportional odds model for continuous Y)
  • Allow for quadratic effects for two predictors which are allowed to interact
  • Three other predictors that are irrelevant are modeled linearly
  • Bootstrap
  • Through uncertainty intervals for REVs show that when the sample size is inadequate for fitting an accurate reproducible model it’s also inadequate for measuring variable importance
n <- 60
x1 <- rnorm(n)
x2 <- rnorm(n)
x3 <- rnorm(n)
x4 <- rnorm(n)
x5 <- rnorm(n)
y  <- x1 + x2 + rnorm(n)
d  <- data.frame(x1, x2, x3, x4, y)
f  <- orm(y ~ pol(x1, 2) * pol(x2, 2) + x3 + x4 + x5,
          x=TRUE, y=TRUE)
a <- anova(f, test='LR')
Makes print methods output html for prettier printing
x=TRUE, y=TRUE are needed for later bootstrapping
Computes likelihood ratio \(\chi^2\) test statistics for all partial effects in the model plus for the whole model
Likelihood Ratio Statistics for y
χ2 d.f. P
x1 (Factor+Higher Order Factors) 40.74 6 <0.0001
All Interactions 4.71 4 0.3183
Nonlinear (Factor+Higher Order Factors) 3.58 3 0.3102
x2 (Factor+Higher Order Factors) 29.88 6 <0.0001
All Interactions 4.71 4 0.3183
Nonlinear (Factor+Higher Order Factors) 3.84 3 0.2788
x3 0.08 1 0.7786
x4 1.07 1 0.3018
x5 0.15 1 0.6957
x1 × x2 (Factor+Higher Order Factors) 4.71 4 0.3183
Nonlinear 4.24 3 0.2363
Nonlinear Interaction : f(A,B) vs. AB 4.24 3 0.2363
f(A,B) vs. Af(B) + Bg(A) 1.74 1 0.1876
Nonlinear Interaction in x1 vs. Af(B) 3.53 2 0.1715
Nonlinear Interaction in x2 vs. Bg(A) 3.42 2 0.1807
TOTAL NONLINEAR 4.56 5 0.4721
TOTAL 47.65 11 <0.0001
plot(a, what='proportion chisq', sort='ascending')
Plots ratios of partial \(\chi^2\) to total model LR \(\chi^2\). ascending is really a misnomer.
Figure 5.1: Relative LR \(\chi^2\) explained. Interaction effects are added to main effects.
g <- bootcov(f, B=300)
Did not converge in 15 iterations
plot(rexVar(g, data=d))
Figure 5.2: Relative explained variation due to each predictor. Interaction effects are added to main effects. Intervals are 0.95 bootstrap percentile confidence intervals.

5.2 The Bootstrap

  • If know population model, use simulation or analytic derivations H to study behavior of statistical estimator
  • Suppose \(Y\) has a cumulative dist. fctn. \(F(y) = \Pr\{Y \leq y\}\)
  • We have sample of size \(n\) from \(F(y)\),
    \(Y_{1}, Y_{2}, \ldots, Y_{n}\)
  • Steps:
    1. Repeatedly simulate sample of size \(n\) from \(F\)
    2. Compute statistic of interest
    3. Study behavior over \(B\) repetitions
  • Example: 1000 samples, 1000 sample medians, compute their sample variance I
  • \(F\) unknown \(\rightarrow\) estimate by empirical dist. fctn. \[ F_{n}(y) = \frac{1}{n}\sum_{i=1}^{n} [Y_{i} \leq y]. \]
  • Example: sample of size \(n=30\) from a normal distribution with mean 100 and SD 10
x <- rnorm(30, 100, 20)
xs <- seq(50, 150, length=150)
cdf <- pnorm(xs, 100, 20)
plot(xs, cdf, type='l', ylim=c(0,1),
     ylab=expression(paste("Pr(", X <= x, ")")))
lines(ecdf(x), cex=.5)
Figure 5.3: Empirical and population cumulative distribution function
  • \(F_{n}\) corresponds to density function placing probability J \(\frac{1}{n}\) at each observed data point (\(\frac{k}{n}\) if point duplicated \(k\) times)
  • Pretend that \(F \equiv F_{n}\)
  • Sampling from \(F_{n} \equiv\) sampling with replacement from observed data \(Y_{1},\ldots,Y_{n}\)
  • Large \(n\) \(\rightarrow\) selects \(1-e^{-1} \approx 0.632\) of original data points in each bootstrap sample at least once
  • Some observations not selected, others selected more than once
  • Efron’s bootstrap \(\rightarrow\) K general-purpose technique for estimating properties of estimators without assuming or knowing distribution of data \(F\)
  • Take \(B\) samples of size \(n\) with replacement, choose \(B\) L so that summary measure of individual statistics \(\approx\) summary if \(B=\infty\)
  • Bootstrap based on distribution of observed differences between a resampled parameter estimate and the original estimate telling us about the distribution of unobservable differences between the original estimate and the unknown parameter

Example: Data \((1,5,6,7,8,9)\), obtain 0.80 confidence interval for M population median, and estimate of population expected value of sample median (only to estimate the bias in the original estimate of the median).

spar(ps=9, mfrow=c(1,2))
y <- c(2,5,6,7,8,9,10,11,12,13,14,19,20,21)
y <- c(1,5,6,7,8,9)
n   <- length(y)
n2  <- n/2
n21 <- n2+1
B   <- 400
M <- double(B)
plot(0, 0, xlim=c(0,B), ylim=c(3,9),
     xlab="Bootstrap Samples Used",
     ylab="Mean and 0.1, 0.9 Quantiles", type="n")
for(i in 1:B) {
  s <- sample(1:n, n, replace=T)
  x <- sort(y[s])
  m <- .5*(x[n2]+x[n21])
  M[i] <- m
  points(i, mean(M[1:i]), pch=46)
  if(i>=10) {
    q <- quantile(M[1:i], c(.1,.9))
    points(i, q[1], pch=46, col='blue')
    points(i, q[2], pch=46, col='blue')
  1   3 3.5   4 4.5   5 5.5   6 6.5   7 7.5   8 8.5   9 
  2   7   6   2   1  30  45  59  72  70  45  48   8   5 
hist(M, nclass=length(unique(M)), xlab="", main="")
Figure 5.4: Estimating properties of sample median using the bootstrap

First 20 samples:

Bootstrap Sample Sample Median
1 5 5 7 8 9 6.0
1 1 5 7 9 9 6.0
6 7 7 8 9 9 7.5
1 1 5 6 8 9 5.5
1 6 7 7 8 8 7.0
1 5 6 8 8 9 7.0
1 6 8 8 9 9 8.0
5 5 6 7 8 9 6.5
1 5 6 7 7 8 6.5
1 5 6 8 9 9 7.0
1 5 7 7 8 9 7.0
1 5 6 6 7 8 6.0
1 6 6 7 8 9 6.5
5 6 7 7 8 9 7.0
1 5 6 8 8 8 7.0
1 1 6 6 7 8 6.0
5 5 5 8 8 9 6.5
5 6 6 6 7 7 6.0
1 5 7 9 9 9 8.0
1 1 5 5 5 7 5.0
  • Histogram tells us whether we can assume normality for the bootstrap N medians or need to use quantiles of medians to construct C.L.
  • Need high \(B\) for quantiles, low for variance (but see Booth & Sarkar (1998))
  • See Bradley Efron & Narasimhan (2020) for useful information about bootstrap confidence intervals and the latest R functions

5.3 Model Validation

5.3.1 Introduction

  • External validation (best: another country at another time); O also validates sampling, measurements
  • Internal
    • apparent (evaluate fit on same data used to create fit)
    • data splitting
    • cross-validation
    • bootstrap: get overfitting-corrected accuracy index
  • Best way to make model fit data well is to discard much of the data P
  • Predictions on another dataset will be inaccurate
  • Need unbiased assessment of predictive accuracy

Working definition of external validation: Validation of a prediction tool on a sample that was not available at publication time.

Alternate: Validation of a prediction tool by an independent research team.

One suggested hierarchy of the quality of various validation methods is as follows, ordered from worst to best.

  1. Attempting several validations (internal or external) and Q reporting only the one that “worked”
  2. Reporting apparent performance on the training dataset (no validation)
  3. Reporting predictive accuracy on an undersized independent test sample
  4. Internal validation using data splitting where at least one of the training and test samples is not huge and the investigator is not aware of the arbitrariness of variable selection done on a single sample
  5. Strong internal validation using 100 repeats of 10-fold R cross-validation or several hundred bootstrap resamples, repeating all analysis steps involving \(Y\) afresh at each re-sample and the arbitrariness of selected “important variables” is reported (if variable selection is used)
  6. External validation on a large test sample, done by the original research team
  7. Re-analysis by an independent research team using strong internal validation of the original dataset
  8. External validation using new test data, done by an independent research team
  9. External validation using new test data generated using different instruments/technology, done by an independent research team

Some points to consider:

  • Unless both sample sizes are huge, external validation can be S low precision
  • External validation can be costly and slow and may result in disappointment that would have been revealed earlier with rigorous internal validation
  • External validation is sometimes gamed; researchers disappointed in the validation sometimes ask for a “do over”; resampling validation is harder to game as long as all analytical steps using \(Y\) are repeated each time.
  • Instead of external validation to determine model applicability T at a different time or place, and being disappointed if the model does not work in that setting, consider building a unified model containing time and place as predictors
  • When the model was fully pre-specified, external validation tests the model
  • But when the model was fitted using machine learning, feature screening, variable selection, or model selection, the model developed using training data is usually only an example of a model, and the test sample validation could be called an example validation
  • When resampling is used to repeat all modeling steps for each resample, rigorous internal validation tests the process used to develop the model and happens to also provide a high-precision estimate of the likely future performance of the “final” model developed using that process, properly penalizing for model uncertainty.
  • Resampling also reveals the volatility of the model selection process

\(\rightarrow\) See BBR

Collins et al. (2016) estimate that a typical sample size needed for externally validating a time-to-event model is 200 events.

5.3.2 Which Quantities Should Be Used in Validation?

  • OLS: \(R^2\) is one good measure for quantifying drop-off in U predictive ability
  • Example: \(n=10, p=9\), apparent \(R^{2}=1\) but \(R^2\) will be close to zero on new subjects
  • Example: \(n=20, p=10\), apparent \(R^{2}=.9\), \(R^2\) on new data 0.7, \(R^{2}_{adj} = 0.79\)
  • Adjusted \(R^2\) solves much of the bias problem assuming \(p\) in its formula is the largest number of parameters ever examined against \(Y\)
  • Few other adjusted indexes exist
  • Also need to validate models with phantom d.f.
  • Cross-validation or bootstrap can provide unbiased estimate of V any index; bootstrap has higher precision
  • Two main types of quantities to validate
    1. Calibration or reliability: ability to make unbiased estimates of response (\(\hat{Y}\) vs. \(Y\))
    2. Discrimination: ability to separate responses
      OLS: \(R^2\); \(g\)-index; binary logistic model: ROC area, equivalent to rank correlation between predicted probability of event and 0/1 event
  • Unbiased validation nearly always necessary, to detect overfitting

5.3.3 Data-Splitting

  • Split data into training and test sets W
  • Interesting to compare index of accuracy in training and test
  • Freeze parameters from training
  • Make sure you allow \(R^{2} = 1-SSE/SST\) for test sample to be \(<0\)
  • Don’t compute ordinary \(R^2\) on \(X\hat{\beta}\) vs. \(Y\); this allows for linear recalibration \(aX\hat{\beta} + b\) vs. \(Y\)
  • Test sample must be large enough to obtain very accurate assessment of accuracy X
  • Training sample is what’s left
  • Example: overall sample \(n=300\), training sample \(n=200\), develop model, freeze \(\hat{\beta}\), predict on test sample (\(n=100\)), \(R^{2} = 1 - \frac{\sum(Y_{i}-X_{i}\hat{\beta})^{2}}{\sum(Y_{i}-\bar{Y})^{2}}\).
  • Disadvantages of data splitting: Y
    1. Costly in \(\downarrow n\) (Breiman, 1992; Roecker, 1991)
    2. Requires decision to split at beginning of analysis
    3. Requires larger sample held out than cross-validation
    4. Results vary if split again
    5. Does not validate the final model (from recombined data)
    6. Not helpful in getting CL corrected for var. selection
    7. Nice summary of disadvantages: Steyerberg (2018)

5.3.4 Improvements on Data-Splitting: Resampling

  • No sacrifice in sample size Z
  • Work when modeling process automated
  • Bootstrap excellent for studying arbitrariness of variable selection (Sauerbrei & Schumacher, 1992).
  • Cross-validation solves many problems of data splitting (B. Efron, 1983; Shao, 1993; van Houwelingen & le Cessie, 1990; Wu, 1986)
  • Example of \(\times\)-validation: A
    1. Split data at random into 10 tenths
    2. Leave out \(\frac{1}{10}\) of data at a time
    3. Develop model on \(\frac{9}{10}\), including any variable selection, pre-testing, etc.
    4. Freeze coefficients, evaluate on \(\frac{1}{10}\)
    5. Average \(R^2\) over 10 reps
  • Drawbacks:
    1. Choice of number of groups and repetitions
    2. Doesn’t show full variability of var. selection
    3. Does not validate full model
    4. Lower precision than bootstrap
    5. Need to do 50 repeats of 10-fold cross-validation to ensure adequate precision
  • Randomization method B
    1. Randomly permute \(Y\)
    2. Optimism = performance of fitted model compared to what expect by chance

5.3.5 Validation Using the Bootstrap

  • Estimate optimism of final whole sample fit without C holding out data
  • From original \(X\) and \(Y\) select sample of size \(n\) with replacement
  • Derive model from bootstrap sample
  • Apply to original sample
  • Simple bootstrap uses average of indexes computed on original sample
  • Estimated optimism = difference in indexes
  • Repeat about \(B=100\) times, get average expected optimism
  • Subtract average optimism from apparent index in final model
  • Example: \(n=1000\), have developed a final model that is D hopefully ready to publish. Call estimates from this final model \(\hat{\beta}\).
    • final model has apparent \(R^2\) (\(R^{2}_{app}\)) =0.4
    • how inflated is \(R^{2}_{app}\)?
    • get resamples of size 1000 with replacement from original 1000
    • for each resample compute \(R^{2}_{boot}\) = apparent \(R^2\) in bootstrap sample
    • freeze these coefficients (call them \(\hat{\beta}_{boot}\)), apply to original (whole) sample \((X_{orig}, Y_{orig})\) to get \(R^{2}_{orig} = R^{2}(X_{orig}\hat{\beta}_{boot}, Y_{orig})\)
    • optimism = \(R^{2}_{boot} - R^{2}_{orig}\)
    • average over \(B=100\) optimisms to get \(\overline{optimism}\)
    • \(R^{2}_{overfitting~corrected} = R^{2}_{app} - \overline{optimism}\)
  • Example: Chapter 8
  • Is estimating unconditional (not conditional on \(X\)) distribution E of \(R^2\), etc. (Faraway, 1992, p. 217)
  • Conditional estimates would require assuming the model one is trying to validate
  • Efron’s “\(.632\)” method may perform better (reduce bias further) for small \(n\) (B. Efron, 1983), (Bradley Efron & Tibshirani, 1993, p. 253), Bradley Efron & Tibshirani (1997)

Bootstrap useful for assessing calibration in addition to discrimination:

  • Fit \(C(Y|X) = X\beta\) on bootstrap sample F
  • Re-fit \(C(Y|X) = \gamma_{0} + \gamma_{1}X\hat{\beta}\) on same data
  • \(\hat{\gamma}_{0}=0, \hat{\gamma}_{1}=1\)
  • Test data (original dataset): re-estimate \(\gamma_{0}, \gamma_{1}\)
  • \(\hat{\gamma}_{1}<1\) if overfit, \(\hat{\gamma}_{0} > 0\) to compensate
  • \(\hat{\gamma}_{1}\) quantifies overfitting and useful for improving calibration (Spiegelhalter, 1986)
  • Use Efron’s method to estimate optimism in \((0,1)\), estimate \((\gamma_{0}, \gamma_{1})\) by subtracting optimism from \((0,1)\)
  • See also Copas (1987) and van Houwelingen & le Cessie (1990), p. 1318

See Freedman et al. (1988) for warnings about the bootstrap, and B. Efron (1983) for variations on the bootstrap to reduce bias.

Use bootstrap to choose between full and reduced models:

  • Bootstrap estimate of accuracy for full model G
  • Repeat, using chosen stopping rule for each re-sample
  • Full fit usually outperforms reduced model(Spiegelhalter, 1986)
  • Stepwise modeling often reduces optimism but this is not offset by loss of information from deleting marginal var.
Method Apparent Rank Correlation of Predicted vs. Observed Over-optimism Bias-Corrected Correlation
Full Model 0.50 0.06 0.44
Stepwise Model 0.47 0.05 0.42

In this example, stepwise modeling lost a possible \(0.50 - 0.47 = 0.03\) predictive discrimination. The full model fit will especially be an improvement when

  1. The stepwise selection deleted several variables which were H almost significant.
  2. These marginal variables have some real predictive value, even if it’s slight.
  3. There is no small set of extremely dominant variables that would be easily found by stepwise selection.

Other issues:

  • See van Houwelingen & le Cessie (1990) for many interesting ideas I
  • Faraway (1992) shows how bootstrap is used to penalize for choosing transformations for \(Y\), outlier and influence checking, variable selection, etc. simultaneously
  • Brownstone (1988), p. 74 feels that “theoretical statisticians have been unable to analyze the sampling properties of (usual multi-step modeling strategies) under realistic conditions” and concludes that the modeling strategy must be completely specified and then bootstrapped to get consistent estimates of variances and other sampling properties
  • See Blettner & Sauerbrei (1993) and Chatfield (1995) for more interesting examples of problems resulting from data-driven analyses.

5.3.6 Multiple Imputation and Resampling-Based Model Validation

When there are missing data, especially missing \(X\), there are simple approaches that are sometimes used for strong internal validation, including

  • When the fraction of incomplete observations is small, fit the model on complete cases and validate it with the bootstrap or cross-validation.
  • When the fraction of incomplete observations is moderately small, do single imputation and validate the model the ordinary way as in the last bullet.
  • For either solution, still use multiple imputation for coefficient estimation, confidence intervals, and statistical tests.

A more general approach is needed.

  • Spirit of multiple imputation: average the regression coefficients over \(k\) model fits for \(k\) completed datasets
  • For validation: do a bootstrap validation separately for each model fit on each completed dataset
  • Average the \(k\) estimates of model performance
  • For each fit use \(B\) bootstrap resamples for a total of \(kB\) iterations
  • May be OK to choose lower value of \(B\) than usual because of the averaging over \(k\) validations
  • For \(k=20\) consider using \(B=100\), with a lower \(B\) for very large dataset
  • Procedure is implemented using the Hmisc fit.mult.impute function (for Hmisc 5.0-1 or later) and the rms processMI function (for rms 6.6-0 or later)
  • Example: Simulate data making two predictors missing completely at random for \(\frac{1}{5}\) of the observations
  • Ordinary linear model
  • 5 noise variables added to create more overfitting
n <- 200
d <- data.frame(age = rnorm(n, 50, 10),
                sex = sample(c('f', 'm'), n, TRUE),
                x1  = runif(n),
                x2  = runif(n),
                x3  = runif(n),
                x4  = runif(n),
                x5  = runif(n)   )
mis <- function(x) ifelse(runif(n) < 0.2, NA, x)
d <- transform(d,
  y = 0.07 * (age - 50) + (sex == 'f') + rnorm(n),
  age = mis(age),
  sex = mis(sex) )
fracmiss <- with(d, mean(is.na(age) | is.na(sex)))
# Do inefficient complete case analysis
f <- ols(y ~ age + sex + x1 + x2 + x3 + x4 + x5, data=d)

Linear Regression Model

ols(formula = y ~ age + sex + x1 + x2 + x3 + x4 + x5, data = d)
Frequencies of Missing Values Due to Each Variable
  y age sex  x1  x2  x3  x4  x5 
  0  51  41   0   0   0   0   0 
Model Likelihood
Ratio Test
Obs 118 LR χ2 57.26 R2 0.384
σ 1.0919 d.f. 7 R2adj 0.345
d.f. 110 Pr(>χ2) 0.0000 g 0.935


     Min       1Q   Median       3Q      Max 
-2.95513 -0.63703  0.06559  0.64790  2.55598 
β S.E. t Pr(>|t|)
Intercept  -2.5095  0.7361 -3.41 0.0009
age   0.0787  0.0109 7.25 <0.0001
sex=m  -0.5208  0.2062 -2.53 0.0130
x1  -0.2570  0.3390 -0.76 0.4501
x2   0.4045  0.3322 1.22 0.2260
x3  -0.1084  0.3652 -0.30 0.7671
x4  -0.4420  0.3429 -1.29 0.2001
x5  -0.8222  0.3499 -2.35 0.0206
  • Fraction of missing data is 0.41 so let’s run 40 imputations
a <- aregImpute(~ y + age + sex + x1 + x2 + x3 + x4 + x5,
                data=d, n.impute=40, pr=FALSE)   # 2s
  • For each imputation, complete the dataset using predictive mean matching
  • Fit a linear model on that completed dataset and run bootstrap validation of both model performance metrics and a smooth calibration curve
  • Just do \(B=50\) resamples per imputaton because we are doing so many imputations
  • Fitting linear models is so fast that higher \(B\) could have been used
  • Specify a function to fit.mult.impute to tell it which extra analyses are needed per completed dataset
g <- function(fit) 
  list(validate  = validate(fit, B=50),
       calibrate = calibrate(fit, B=50) )
f <- fit.mult.impute(y ~ age + sex + x1 + x2 + x3 + x4 + x5,
                     ols, a, data=d, pr=FALSE,
                     fun=g, fitargs=list(x=TRUE, y=TRUE))   # 3s

Linear Regression Model

fit.mult.impute(formula = y ~ age + sex + x1 + x2 + x3 + x4 + 
    x5, fitter = ols, xtrans = a, data = d, fun = g, pr = FALSE, 
    fitargs = list(x = TRUE, y = TRUE))
Model Likelihood
Ratio Test
Obs 200 LR χ2 78.74 R2 0.325
σ 1.0623 d.f. 7 R2adj 0.300
d.f. 192 Pr(>χ2) 0.0000 g 0.814


     Min       1Q   Median       3Q      Max 
-2.89841 -0.76497 -0.05365  0.68108  2.68169 
β S.E. t Pr(>|t|)
Intercept  -2.4529  0.5472 -4.48 <0.0001
age   0.0710  0.0089 8.01 <0.0001
sex=m  -0.5044  0.1824 -2.77 0.0062
x1  -0.0454  0.2779 -0.16 0.8704
x2   0.4237  0.2560 1.66 0.0995
x3  -0.1978  0.2815 -0.70 0.4831
x4  -0.4577  0.2720 -1.68 0.0940
x5  -0.4373  0.2736 -1.60 0.1116
  • The numbers in the top box of the output are averaged over the imputations. So are the \(\beta\)s.
  • Average the bootstrapped performance measures and print
processMI(f, 'validate')
Index Original
Optimism Corrected
R2 0.3249 0.3444 0.2961 0.0483 0.2766 2000
MSE 1.0837 1.039 1.1299 -0.0909 1.1746 2000
g 0.8137 0.8313 0.7882 0.0431 0.7706 2000
Intercept 0 0 0.0235 -0.0235 0.0235 2000
Slope 1 1 0.9506 0.0494 0.9506 2000
  • For the estimated calibration curves, processMI by default shows all the individual curves from all multiple imputations (set plotall=FALSE to suppress this)
  • Use an option to also show the first three individual validations followed by the average
cal <- processMI(f, 'calibrate', nind=3)

n=200   Mean absolute error=0.091   Mean squared error=0.0139
0.9 Quantile of absolute error=0.186

n=200   Mean absolute error=0.108   Mean squared error=0.0171
0.9 Quantile of absolute error=0.181

n=200   Mean absolute error=0.095   Mean squared error=0.0123
0.9 Quantile of absolute error=0.19

n=200   Mean absolute error=0.12   Mean squared error=0.021
0.9 Quantile of absolute error=0.235
  • Run plot(cal) if you want to see a full size final average calibration

5.3.7 Validation of Data Reduction

A frequently asked question is whether or not unsupervised learning (data reduction) needs to be repeated inside an internal validation loop. See this for some excellent discussion. Here are are few thoughts.

  • To be safe, include fresh data reduction inside each loop
  • If you want to take a slight risk you can do the data reduction once on the whole dataset, because of the following
  • Principal components (especially minor components), variable and observation clustering, and redundancy analysis can be very unstable
  • They way in which they are unstable is uninformed by \(Y\), so this instability can work either for or against you
  • Typically we are allowed to use procedures that don’t always work in our favor
  • Contrast this with stepwise variable selection, which always attempts to work in our favor by maximizing \(R^2\)

5.4 Bootstrapping Ranks of Predictors

  • Order of importance of predictors not pre-specified J
  • Researcher interested in determining “winners” and “losers”
  • Bootstrap useful in documenting the difficulty of this task
  • Get confidence limits of the rank of each predictor in the scale of partial \(\chi^2\) - d.f.
  • Example using OLS
# Use the plot method for anova, with pl=FALSE to suppress actual
# plotting of chi-square - d.f. for each bootstrap repetition.
# Rank the negative of the adjusted chi-squares so that a rank of
# 1 is assigned to the highest.  It is important to tell
# plot.anova.rms not to sort the results, or every bootstrap
# replication would have ranks of 1,2,3,... for the stats.
n <- 300
d <- data.frame(x1=runif(n), x2=runif(n), x3=runif(n), x4=runif(n),
                x5=runif(n), x6=runif(n), x7=runif(n), x8=runif(n),
             x9=runif(n), x10=runif(n), x11=runif(n), x12=runif(n))
d$y <- with(d, 1*x1 + 2*x2 +  3*x3  +  4*x4  + 5*x5 + 6*x6 + 7*x7 +
               8*x8 + 9*x9 + 10*x10 + 11*x11 + 12*x12 + 9*rnorm(n))

f <- ols(y ~ x1+x2+x3+x4+x5+x6+x7+x8+x9+x10+x11+x12, data=d)
B <- 1000
ranks <- matrix(NA, nrow=B, ncol=12)
rankvars <- function(fit)
  rank(plot(anova(fit), sort='none', pl=FALSE))
Rank <- rankvars(f)
for(i in 1:B) {    # 6s
  j <- sample(1:n, n, TRUE)
  bootfit <- update(f, data=d, subset=j)
  ranks[i,] <- rankvars(bootfit)
lim <- t(apply(ranks, 2, quantile, probs=c(.025,.975)))
predictor <- factor(names(Rank), names(Rank))
w <- data.frame(predictor, Rank, lower=lim[,1], upper=lim[,2])
ggplot(w, aes(x=predictor, y=Rank)) + geom_point() + coord_flip() +
  scale_y_continuous(breaks=1:12) +
  geom_errorbar(aes(ymin=lim[,1], ymax=lim[,2]), width=0)
Figure 5.5: Bootstrap percentile 0.95 confidence limits for ranks of predictors in an OLS model. Ranking is on the basis of partial \(\chi^2\) minus d.f. Point estimates are original ranks

See this for a higher-dimension example.

5.5 Simplifying the Final Model by Approximating It

5.5.1 Difficulties Using Full Models

  • Predictions are conditional on all variables, standard errors K \(\uparrow\) when predict for a low-frequency category
  • Collinearity
  • Can average predictions over categories to marginalize, \(\downarrow\) s.e.

5.5.2 Approximating the Full Model

  • Full model is gold standard L
  • Approximate it to any desired degree of accuracy
  • If approx. with a tree, best c-v tree will have 1 obs./node
  • Can use least squares to approx. model by predicting \(\hat{Y} = X\hat{\beta}\)
  • When original model also fit using least squares, coef. of approx. model against \(\hat{Y} \equiv\) coef. of subset of variables fitted against \(Y\) (as in stepwise)
  • Model approximation still has some advantages M
    1. Uses unbiased estimate of \(\sigma\) from full fit
    2. Stopping rule less arbitrary
    3. Inheritance of shrinkage
  • If estimates from full model are \(\hat{\beta}\) and approx.
    model is based on a subset \(T\) of predictors \(X\), coef. of approx.
    model are \(W \hat{\beta}\), where
    \(W = (T'T)^{-1}T'X\)
  • Variance matrix of reduced coef.: \(W V W'\)

5.6 How Do We Break Bad Habits?

  • Insist on validation of predictive models and discoveries N
  • Show collaborators that split-sample validation is not appropriate unless the number of subjects is huge
    • Split more than once and see volatile results
    • Calculate a confidence interval for the predictive accuracy in the test dataset and show that it is very wide
  • Run simulation study with no real associations and show that associations are easy to find
  • Analyze the collaborator’s data after randomly permuting the \(Y\) vector and show some positive findings
  • Show that alternative explanations are easy to posit
    • Importance of a risk factor may disappear if 5 “unimportant” risk factors are added back to the model
    • Omitted main effects can explain apparent interactions
    • Uniqueness analysis: attempt to predict the predicted values from a model derived by data torture from all of the features not used in the model

5.7 Study Questions

Section 5.1

  1. What is a general way to solve the problem of how to understand a complex model?
  2. Why is it useful to estimate the effect of changing a continuous predictor on the outcome by varying the predictor from its \(25^{th}\) to its \(75^{th}\) percentile?

Section 5.2

  1. Describe briefly in very general terms what the bootstrap involves.
  2. The bootstrap may be used to estimate measures of many types of statistical quantities. Which type is being estimated when it comes to model validation?

Section 5.3

  1. What are some disadvantages of external validation?
  2. What would cause an independently validated \(R^2\) to be negative?
  3. A model to predict the probability that a patient has a disease was developed using patients cared for at hospital A. Sort the following in order from the worst way to validate the predictive accuracy of a model to the best and most stringent way. You can just specify the letters corresponding to the following, in the right order.
    1. compute the predictive accuracy of the model when applied to a large number of patients in another country
    2. compute the accuracy on a large number of patients in hospital B that is in the same state as hospital A
    3. compute the accuracy on 5 patients from hospital \(A\) that were collected after the original sample used in model fitting was completed
    4. use the bootstrap
    5. randomly split the data into two halves assuming that the dataset is not extremely large; develop the model on one half, freeze it, and evaluate the predictive accuracy on the other half
    6. report how well the model predicted the disease status of the patients used to fit the model
    7. compute the accuracy on patients in hospital \(A\) that were admitted after the last patient’s data used to fit the model were collected, assuming that both samples were large

Section 5.6

  1. What is a way to do a misleading interaction analysis?
  2. What would cause the most embarrassment to a researcher who analyzed a high-dimensional dataset to find “winning” associations?