Let’s say we fit a logistic regression model for the purposes of predicting the probability of low infant birth weight, which is an infant weighing less than 2.5 kg. Below we fit such a model using the “birthwt” data set that comes with the MASS package in R. (This is an example model and not to be used as medical advice.)

We first subset the data to select four variables:

- low: indicator of birth weight less than 2.5 kg (0 = no, 1 = yes)
- ht: history of hypertension (0 = no, 1 = yes)
- ptl: previous premature labor (0 = no, 1 = yes)
- lwt: mother’s weight in pounds at last menstrual period

Then we model “low” as a simple additive function of the other three variables. The models says if we know the mother’s history of hypertension, whether or not she had previous premature labor, and her weight at last menstrual period, then we can use that information to estimate the probability of a low birth weight infant.

library(MASS) data("birthwt") d <- birthwt[,c("low", "ht", "ptl", "lwt")] # convert ptl to indicator d$ptl <- ifelse(d$ptl < 1, 0, 1) m <- glm(low ~ ht + ptl + lwt, family = binomial, data = d) summary(m) ## ## Call: ## glm(formula = low ~ ht + ptl + lwt, family = binomial, data = d) ## ## Deviance Residuals: ## Min 1Q Median 3Q Max ## -1.7420 -0.8018 -0.6895 0.9647 2.2460 ## ## Coefficients: ## Estimate Std. Error z value Pr(>|z|) ## (Intercept) 1.017367 0.853337 1.192 0.23317 ## ht 1.893971 0.721090 2.627 0.00863 ** ## ptl 1.406770 0.428501 3.283 0.00103 ** ## lwt -0.017280 0.006787 -2.546 0.01090 * ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## (Dispersion parameter for binomial family taken to be 1) ## ## Null deviance: 234.67 on 188 degrees of freedom ## Residual deviance: 210.12 on 185 degrees of freedom ## AIC: 218.12 ## ## Number of Fisher Scoring iterations: 4

The model certainly seems like it’s doing something. The coefficients are “significant” which means they’re reliably positive or negative. But is the model any good? Does it do a good job of estimating an appropriate probability for the event of interest?

One statistic to help us assess this is Somers’ D, which is a rank correlation between predicted probabilities and observed responses. It ranges between -1 and 1. When D = 1, the predictions are perfectly discriminating. When D = 0, the model is making random predictions. The basic idea is that a good model will return higher probabilities more often than not when the binary response is 1 versus 0.

We can calculate Somers’ D using the `somers2()`

function in the Hmisc package. In the output it is listed as “Dxy”. This function also lists another statistic, the c-index, which is also known as AUC or area-under-the-curve. (See this StatLab article to learn more about AUC.) A c-index value close to 1 is best. A value close to 0.5 indicates the model is guessing at random. Somers’ D and the c-index are related as follows:

\[D_{xy} = 2 \times (c – 0.5)\] To use the `somers2()`

function, set the `x`

argument to the predicted probabilities and the `y`

argument to the observed responses.

library(Hmisc) D_orig <- somers2(x = predict(m, type = "response"), y = d$low) D_orig ## C Dxy n Missing ## 0.7189048 0.4378096 189.0000000 0.0000000

The results are not bad. Somers’ D is well above 0 and the c-index is about 0.72. Our model may be good enough for clinical use.

But both indices are *optimistic*. We’re using the same data we used to fit the model to evaluate the model. How will our model work with new data *not used in fitting the model*? To answer this we use *model validation*.

A common approach to model validation is data-splitting. This is where we randomly split our data into training and test sets, build a model using the training data, and then see how the model performs on the test data. For our “birthwt” data set, this would mean holding out some portion of our data (say 30%), building a model with the remaining data, and then evaluating how well the model works on the data held out. Hopefully the model would discriminate between low and normal birth weights with suitable probabilities.

The problem with the data-splitting approach is that it greatly reduces the sample size for model building. In designed experiments and field studies, data can be expensive and/or hard to obtain. We would prefer to not randomly “hold out” some portion of data when we run our analysis. Plus any model validation statistic that we calculate is for the model built on a subset of the data, not all the data. In addition a lucky split could produce an unusually good model. Cross-validation alleviates some of these problems, but again it does not validate the full model built using all the data. (Harrell 2015, p. 112)

A more efficient approach is *bootstrap model validation*. The basic steps are as follows:

- Resample the data with replacement.
- Refit the same model to the resampled data.
- Evaluate the performance of the refit model on the resampled data set (using a index such as Somers’ D).
- Evaluate the performance of the refit model on the original data set
- Take the difference between steps 3 and 4.
- Repeat steps 1 – 5 about 200 times.
- Take the average of the 200 differences and subtract from the original index estimate to get a
*bias-corrected*estimate of how the final model will perform on future data.

Let’s demonstrate by iterating through steps 1-5 above one time. First we resample the number of rows with replacement and save to `i`

. (We use `set.seed(11)`

so the interested reader can follow along and reproduce these results.) Then we use `i`

to subset the data when we refit the model. Next we calculate two versions of Somers’ D: one using the resampled data (which we call “D_train”) and another using the original data (which we call “D_test”). Notice we use subsetting brackets to extract just the “Dxy” value. Finally we take the difference between the two statistics and subtract from the original estimate. Notice the “Corrected” Somers’ D is slightly lower.

set.seed(11) i <- sample(nrow(d), size = nrow(d), replace = TRUE) m2 <- glm(low ~ ht + ptl + lwt, family = binomial, data = d[i,]) D_train <- somers2(x = predict(m2, type = "response"), y = m2$y)["Dxy"] D_test <- somers2(x = predict(m2, newdata = d, type = "response"), y = d$low)["Dxy"] D_diff <- D_train - D_test cbind("Original" = D_orig["Dxy"], "Corrected" = D_orig["Dxy"] - D_diff) ## Original Corrected ## Dxy 0.4378096 0.4366723

Now let’s do this 200 times. We could create a loop to do this but we’ll use the `boot()`

function from boot package that comes with R. We first need to create a function that takes at least two arguments, one for the data frame and one for resampling the data frame. We name the function `somersd`

and call the first argument `data`

. This argument will take the data frame we use for resampling. We call the second argument `index`

. This argument will take the resampled row numbers of the data frame that `boot()`

will use to create the bootstrap sample. You can see the body of the function is basically our code from above, with `d`

and `i`

replaced with `data`

and `index`

. The function returns the difference between the two Somers’ D estimates. Recall the first estimate is for the performance of the bootstrapped model on the bootstrapped data, and the second estimate is for the performance of the bootstrapped model on the original data.

somersd <- function(data, index){ m <- glm(low ~ ht + ptl + lwt, family = binomial, data = data[index,]) # train D_train <- somers2(predict(m, type = "response"), m$y)["Dxy"] # test D_test <- somers2(predict(m, newdata = data, type = "response"), d$low)["Dxy"] D_train - D_test }

With our function created, we’re ready to proceed with bootstrapping. We first load the boot package so we can access the `boot()`

function. The first argument is the data frame we want to bootstrap. In this case this is our original data set, `d`

. The second argument is the statistic we want to calculate on each bootstrap sample. This is where our custom `somersd`

function comes in. The third argument, `R`

, is how many bootstrap replicates we wish to perform. Finally, we save the result to an object called `sd.out`

. We also use `set.seed(222)`

so you can replicate the results.

library(boot) set.seed(222) sd.out <- boot(data = d, statistic = somersd, R = 200)

The bootstrap replicates are in a vector called `t`

in the `sd.out`

object. We can take the mean of that vector and subtract it from our original Somers’ D estimate to get a bias-corrected estimate. The estimate is a bit lower reflecting the fact our model won’t perform as well on future data.

D_orig["Dxy"] - mean(sd.out$t) ## Dxy ## 0.4247608

The rms package provides the `validate()`

function that automatically performs bootstrap model validation for Somers’ D and a number of other indices. The only catch is we need to use the modeling functions in the rms package in order to use `validate()`

. In this case we need to use the `lrm()`

function, which stands for “logistic regression modeling”. It uses the same formula notation as the `glm()`

function but does not have a `family`

argument. To use `validate()`

with the model object, we need to also include the arguments `x = TRUE, y = TRUE`

when fitting the model. This simply saves the original data with the model object. Below we tell `validate()`

to run 200 bootstrap replicates for our fitted rms model.

library(rms) m_rms <- lrm(low ~ ht + ptl + lwt, data = d, x = TRUE, y = TRUE) set.seed(333) v <- validate(m_rms, B = 200) v ## index.orig training test optimism index.corrected n ## Dxy 0.4378 0.4496 0.4294 0.0201 0.4177 200 ## R2 0.1713 0.1933 0.1557 0.0377 0.1336 200 ## Intercept 0.0000 0.0000 -0.0554 0.0554 -0.0554 200 ## Slope 1.0000 1.0000 0.9007 0.0993 0.9007 200 ## Emax 0.0000 0.0000 0.0327 0.0327 0.0327 200 ## D 0.1246 0.1434 0.1121 0.0313 0.0933 200 ## U -0.0106 -0.0106 0.0052 -0.0158 0.0052 200 ## Q 0.1352 0.1540 0.1069 0.0471 0.0880 200 ## B 0.1852 0.1800 0.1897 -0.0097 0.1950 200 ## g 0.8943 0.9722 0.8331 0.1390 0.7553 200 ## gp 0.1751 0.1807 0.1638 0.0168 0.1583 200

The first row shows Somers’ D as “Dxy”. The bootstrapped validated value is listed in the “index.corrected” column. Notice it’s similar to the value we calculated using `boot()`

. Of course it’s slightly different because bootstrapping is based on randomly selected data sets. The “index.orig” column shows our original estimate of Somers’ D, 0.4378. The “training” and “test” columns shows the average of the two Somers’ D estimates made during bootstrapping. The “training” estimate is the average bootstrap model performance on the bootstrapped data. The “test” estimate is the average bootstrap model performance on the original unsampled data. The “optimism” column is the difference between training and test. The “index.corrected” value is calculated by subtracting “optimism” from “index.orig”.

The `validate()`

output lists 10 more indices. See this extremely helpful blog post for a summary of what they mean and how to interpret them. The “R2” index, as you may have guessed, is a version of the traditional R-squared associated with linear models. It ranges from 0 to 1 and is interpreted as the proportion of variance explained. Clearly values near 1 would be preferred. The “B” index is the Brier score. In its most common formulation, it ranges between 0 and 1, with a score near 0 being desirable. A score of 0.25 is associated with random guessing. See this StatLab article to learn more about Brier Scores.

Based on these bootstrapped indices, we’re a little more modest about our model’s potential performance in a clinical setting. Whereas before we may have been excited about “significant” coefficients, we now see that our model’s predictive ability is decent but not great. Hopefully you now have a better understanding of how to use the bootstrap to validate a statistical model.

**References**

- Angelo Canty and Brian Ripley (2021). boot: Bootstrap R (S-Plus) Functions. R package version 1.3-28.
- Davison, A. C. & Hinkley, D. V. (1997) Bootstrap Methods and Their Applications. Cambridge University Press, Cambridge. ISBN 0-521-57391-2
- Harrell Jr FE (2022).
*rms: Regression Modeling Strategies*. R package version 6.3-0, https://CRAN.R-project.org/package=rms. - Hosmer, D.W. and Lemeshow, S. (1989) Applied Logistic Regression. New York: Wiley
- Harrell Jr F (2015).
*Regression Modeling Strategies, 2nd edition*. Springer. - Harrell Jr F (2022).
*Hmisc: Harrell Miscellaneous*. R package version 4.7-1, https://CRAN.R-project.org/package=Hmisc. - R Core Team (2022). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/. Version 4.2.1.
- Venables, W. N. & Ripley, B. D. (2002) Modern Applied Statistics with S. Fourth Edition. Springer, New York. ISBN 0-387-95457-0

For questions or clarifications regarding this article, contact the UVA Library StatLab: statlab@virginia.edu

View the entire collection of UVA Library StatLab articles.

*Clay Ford
Statistical Research Consultant
University of Virginia Library
September 09, 2022*