Calculating optimism-corrected estimates of performance after bootstrapping

[Annemarie Pantke]

Hello everyone,

I have seen a few posts from the past that had covered similar questions before, but as I unfortunately wasn’t able to solve my problem with the help of those threads I thought I’d start a new topic. I have developed a prediction model (Cox) and would now like to obtain the optimism-corrected measures of performance (p), which in my analysis would be Harrell’s C (C-index), the Brier score as well as the calibration slope. The formula which I’m using for that follows the "regular bootstrap" as referred to by Steyerberg et al. (2001): p_corrected = p_apparent - p_optimism

Here is an outline of the steps:

1. Train a model on the original dataset and record the value of a performance metric of interest.
2. Generate a bootstrap sample.
3. Develop a model using the bootstrap sample (applying the same predictors) and record the corresponding performance metric for the bootstrap-sample-derived model.
4. Apply the bootstrap model to the original dataset and obtain the performance metric.
5. Estimate optimism by taking the mean of the differences between the values calculated in step 3 (the apparent performance of the bootstrap-sample-derived model) and step 4 (the bootstrap-sample-derived model's performance when tested on the original sample).
6. Calculate the optimism-corrected value of the performance metric as the difference between the values calculated in step 1 (the naive value) and step 5 (the estimated optimism).

The code that I have so far (for the C-index) goes as follows:

capture program drop optimism
program define optimism, rclass
     preserve
     bsample
     stcox i.risk ib4.cd4baseline_group i.vlbaseline_group age i.SEX, nohr
     estat concordance
     return scalar c = r(C)
end

stcox i.risk ib4.cd4baseline_group i.vlbaseline_group age i.SEX, nohr
estat concordance
local base_harrell = r(C)
tempfile sim_results
simulate C = r(c), reps(200) seed(12345) saving(`sim_results'): optimism

use `sim_results', clear
gen diff = C - `base_harrell'
summ diff

I believe that what this code did so far is to calculate the difference between the bootstrap C-index and the original one. Here, I’m not quite sure how to run the bootstrap model in the original dataset and obtain the C-index. I assume the coefficients from each bootstrap sample would need to be saved (matrix b = e(b) ? ) and then applied to the original sample, but unfortunately I can’t figure out how to change the code accordingly so that all the steps as indicated above are carried out properly.

If anyone has some advice or help to share, that would be very much appreciated - many thanks in advance!



Reference:
Steyerberg EW, Harrell FE Jr, Borsboom GJ, Eijkemans MJ, Vergouwe Y, Habbema JD. Internal validation of predictive models: Efficiency of some procedures for logistic regression analysis. J Clin Epidemiol. 2001 Aug;54(8):774-81. doi: 10.1016/s0895-4356(01)00341-9.

[Juan Vicente-Valor]

Hi, did Annemarie or anyone else manage to find a solution? I have the same need (to run code to calculate performance estimates after bootstrapping), but I haven’t been able to find anything.

[Annemarie Pantke]

Hi Juan,

I was actually able to solve this issue with the help from a kind professor from Mahidol University in Thailand who had sent me the entire code for it. It's been a few years now since I did the analyses, but the following code worked for calculating the optimism-corrected performance metrics. I'm copying you here the one i used for the C-index and calibration slope (with my dataset and variables, so those would need to be replaced accordingly of course ) Hope this helps!

use $Datadir/AIDSFile_AP_followup1999-2008.dta, clear    
    
**stset
stset censored, failure(aids) origin(beodate_my) enter(beodate_my) id(IdPatmain) scale(12)


*Fit model
stcox i.risk2 ib4.cd4baseline_group ib5.vlbaseline_group age, nohr base

predict lp, xb
summarize lp

*C-index
estat concordance
global cstat_orig = r(C)    

*Calibration slope
stcox lp, nohr
global cslope_orig = _b[lp]    
    
    
    
*We first create an empty matrix to hold the bs results. Setting a seed means we can reproduce the same results if we need to run the bootstrapping again at a later date
matrix results = J(500,4,.)
set seed 132166
    

    
qui forvalues i=1/500 {
*then load the original sample data
use $Datadir/AIDSFile_AP_followup2009-2018_BS.dta, clear

bsample

nois _dots `i' 0

stset censored, failure(aids) origin(beodate_my) enter(beodate_my) id(IdPatmain) scale(12)

*Fit model to the bootstrap sample
stcox i.risk2 i.ib4.cd4baseline_group ib5.vlbaseline_group age, nohr base
    
*predict probabilities & lp from the bootstrap model in the bs sample
predict lp, xb


**calculate the apparent performance of the bootstrap model in the bs sample

*first calculate the C-statistic
estat concordance
matrix results[`i',1] = r(C)

*calculate the bootstrap apparent C-slope
stcox lp
matrix results[`i',2] = _b[lp]

/* test performance of the bootstrap model in original sample
    NB: to do this we must first fit the model to the bootstrap data again, so
    that Stata has the coefficients stored in memory */
stcox i.risk2 ib4.cd4baseline_group ib5.vlbaseline_group age, nohr base 


*then load the original sample data
use $Datadir/AIDSFile_AP_followup2009-2018_BS.dta, clear

*predict probabilities & lp from the bs model in the original dataset
predict lp, xb

*stset the data as always before fitting survival models in Stata
stset censored, failure(aids) origin(beodate_my) enter(beodate_my) id(IdPatmain) scale(12)

*calculate the test performance in original data C-statistic
stcox lp
estat concordance
matrix results[`i',3] = r(C)

*calculate the test performance in original data C-slope
stcox lp
matrix results[`i',4] = _b[lp]

}    
    
    
    

*rename columns of the matrix used to store bs results
mat colnames results = cstat_boot cslope_boot cstat_test cslope_test 

*clear and load the matrix results into stata dataset
clear
svmat results, n(col)



*We now have a dataset of 500 estimates of apparent and test performance

*We can calculate the optimism
gen optimism_c = cstat_boot - cstat_test
gen optimism_slope = cslope_boot - cslope_test


*summarise the apparent performance of the bootstrap model in the bs samples
summarize cstat_boot cslope_boot 

*summarise the performance of the bs model in the original sample
su cstat_test cslope_test 


*summarise the optimism in the performance

*NB: this is simply the apparent bs performance minus the test performance
su optimism_c 
global cstat_opt = r(mean)
su optimism_slope
global cslope_opt = r(mean)


/* now calculate the optimism adjusted performance of the original model by
    subtracting the optimism from the original models apparent performance in
    the original dataset 
*/
di "Optimism adj C-statistic = " $cstat_orig - $cstat_opt
di "Optimism adj C-slope = " $cslope_orig - $cslope_opt