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