Cross validation of a survival model- what to make of “random effects” of parameter estimates?

Question

This is a question surrounding k-fold cross validation for time to event data. I am interested in what to do with the knowledge that certain variables fail to perform as well within some of the training/validation steps as they did in the full sample training/derivations sets. Essentially, if one is reestimating the model with each fold of the cross validation, how is one to interpret folds where the candidate variables fail to reach statistical significance?

Specifically, I have a dataset of heart transplant recipients and the goal is to create a model of post-transplant survival based on candidate heart donor and recipient/donor matching variables (such as difference in weights). The dataset contains ~19,200 recipient/donor pairs. I have estimated a Cox model using automated variable selection in STATA (done for various reasons). The final multivariable model contains 8 characteristics. However, upon cross validation (1000 random samples leaving out 10% of the sample as a validation sample in each loop) 1 of the parameters loses significance in 61% of training/derivation folds. While I could measure a fit estimate in the 10% leave-out sample, knowledge that the variable has lost significance seems to obviate the need to evaluate fit. While the data could still be overfit, a much larger problem seems to be occurring among the candidate variables.

For example from the pseudo data encoded below (results will vary depending on the randomly generated covariates), the estimates from the full dataset are:

• risk1 HR 0.91 (p = 0.044)
• risk2 HR 1.01 (p =0.152) -deleted from analysis
• risk3 HR 2.31 (p<0.0001)

However, when run through cross validation 1000-folds:

• risk1 mean p-value across folds = 0.078
• risk3 mean p-value across folds <0.0001

Should I pause after discovering this discrepancy and eject variable risk1? If the model has variation in the fit parameter, such as Harrell's c-index, how does one parse the cause of the poor fit, if not by looking at the parameters in each fold?

Thanks

*output data sets

clear
input _k_
0 
end
save k_pseudo
save somersd_pseudo

*pseudo data

clear
set obs 1000
gen t1 = rpoisson(5)
gen risk1 = rnormal(0,1)
gen risk2 = rpoisson(42)
gen risk3 = rbinomial(1,0.5)
gen death = rbinomial(1,0.75) if risk3 == 1
replace death = rbinomial(1,0.3) if risk3 == 0 
replace t1 = t1 - 2.5 if death == 1 & risk1 < 0
save k_fold_data, replace

stset t1 , failure(death ==1)
sts graph, survival by(risk3)
sts gen tx_surv = s, by(risk3)
stcox risk1 risk2 risk3

*k-fold cross validation

foreach i of numlist 1/100 {
    capture drop uniform
    capture drop sample_10
    capture drop censind
    capture drop xb_`i'

    gen uniform = runiform()
    xtile sample_10 = uniform, nq(10)

    qui stset t1 if sample_10 != 1 , failure(death ==1)
    qui stcox risk1 risk2 risk3
    qui regsave risk1 risk2 risk3 using k_pseudo, addlabel(k_fold, `i') pval tstat   
    append
    qui predict xb_`i' if _st !=1, xb
    }
    .

    stset t1, failure(death==1)
    capture drop censind
    gen censind = _d

    foreach x of numlist 1/100{ 
        qui somersd _t xb_`x', cenind(censind) tdist transf(c)
        qui regsave using somersd_pseudo, addlabel(k_fold, `x') pval tstat append
        }
        .

*results- regression parameters and p-values

use k_pseudo, clear
keep if _k_ != 0

*parameter values
    encode var, gen(vars)
    gen stat_sig = 0 
        replace stat_sig = 1 if pval <=0.05

    tab vars stat_sig
    mean coef pval, over(vars)

*results- fit in k_fold cross validation

use somersd_pseudo, clear
keep if _k_ !=0

*statistically significant c-index
gen stat_sig = 0 
    replace stat_sig = 1 if pval <=0.05

tab stat_sig
univar coef pval tstat

Answer 1

Think of cross-validation as evaluating model-fitting strategies, not particular models. In this case, variable selection is part of your model-fitting strategy. So the cross-validated error estimate for this strategy is valid only for a fit of the exact same strategy on the entire dataset. If a variable is not kept in many CV folds, but is kept when running the same procedure on the full dataset, that may just mean that having 10% more data was important for pushing it over the significance threshold, for instance.

However, the variation in which features get kept over CV folds can serve as a useful indicator of the robustness of the fitting procedure. Depending on the automated feature selection method used (I don't know Stata), your p-values may be distorted and the selected features may not be very robust. If this is a worry for you, I recommend Frank Harrell's Regression Modeling Strategies as a reference for sound variable selection methodology; section 4.3 has good info. If you want a sparse model, you could consider doing $L^1$-regularized Cox regression instead.