# How can I assess internal validation (discrimination, calibration) of a Fine-and-Gray competing risk model, fitted using a MFP-algoritm in Stata?

I'm a post-doc at Karolinska Institute, and I'm working on developing a Fine-and-Gray competing risk model to predict endometrial cancer recurrence/progression with death as a competing event.

I have used a multiple fractional polynominal algorithm that simultaneously performs backwards elimination and selects fractional polynominals for continuous variables https://mfp.imbi.uni-freiburg.de/. Since this algoritm is very computationally intensive I've expanded my dataset and assigned weights by using the stcrprep command by P C Lambert (https://pubmed.ncbi.nlm.nih.gov/30542252/). By doing so the native Stata stcox command becomes a Fine-and-Gray model. However, I’ve been unable to get standard commands for assessing calibration, such as stcoxcal by P Royston (https://journals.sagepub.com/doi/pdf/10.1177/1536867X1401400403) to work.

One way to assess calibration is by using pseudo-values (described by Gerds https://onlinelibrary.wiley.com/doi/abs/10.1002/sim.6152), and recommended in Validation of prediction models in the presence of competing risks: a guide through modern methods by van Gerloven et al (https://www.bmj.com/content/bmj/377/bmj-2021-069249.full.pdf). I’ve got the code to work, almost. I can draw one calibration plot, but I have been unable to resample the dataset and repeat the calibration procedure. Since the calibration plot is working, for the first run, I’m thinking that the problem lies in my programming of the bootstrap resampling. Thus, a coding issue, not necessarily a statistical issue. What I want the code to do is resample the data, fit the model, spit out the results, store those in a separate file, repeat X-hundreds of times, and then analyse that file to assess calibration.

A few weeks ago I posted my problem on Statalist, but still with zero replies statalistpost. To that post I attached a sample dataset with 10% of my data.

        clear all
frames reset

************************************************
*CHOOSE WHICH DATASET TO USE********************
*10% OF DATA (FOR TRYING OUT THE CODE) OR FULL**
************************************************
*use "the database file pathway" <---ADD YOURS HERE
*drop if excl_fu == 1

************************************************
*DRAW CALIBRATION PLOT**************************
************************************************

stset date_of_primary_event, id(idnum) failure(primary_end_point == 1) origin(date_of_surgery) scale(365.24)

// generate pseduo obs for calibration plot
stpci primary_end_point, at(3) competingvalues(2)  //Calibration assessment at 3 years.
frame put id pseudo, into(pseudo)

// run stcrprep, same as above.
stset date_of_primary_event, id(idnum) failure(primary_end_point == 1, 2) origin(date_of_surgery) scale(365.24) //Both failcodes has to be included in the failure() option

stcrprep, events(primary_end_point) keep(esmo_pre_op_sub extrauterine age colorscore vascularpattern_dicho emborder echogenicity bmi tumor_size) trans(1 2) //Here I'm omitting the byg() option, since I want the same output as for stcrreg. See Lambert for details.

gen event = primary_end_point == failcode

// fit model to event of interest
stset tstop [pw=weight_c], failure(event) enter(tstart) noshow
mfpa, cycles(3) dfdefault(3) select(0.157, i.esmo_pre_op_sub:1) alpha(0.25) : ///
stcox i.esmo_pre_op_sub c.age i.colorscore i.extrauterine ///
i.vascularpattern_dicho i.emborder i.echogenicity c.bmi c.tumor_size if failcode==1 ///
, vce(cluster idnum) //After using stcrprep, stcox will do the same as stcrreg did before transforming the data.

// predict baseline survival
predict S0 if failcode==1, basesurv

// baseline at 3 years
summ  S0 if _t<=3
global S0_3 r(min)'

// only use first row for each individual as we have expanded data
bysort id (_t): gen first = _n==1
predict xb if first, xb

// Gen predicted CIF at 3 years
gen CIF3 = 1-${S0_3}^exp(xb) hist CIF3, freq // merge with pseudo obs frame put id CIF3 xb if first, into(predcif) frame change pseudo frlink 1:1 id, frame(predcif) gen(l) frget CIF3 xb, from(l) lowess pseudo CIF3, gen(lowess) nograph //Smoothing 'pseudo' and 'CIF3' using LOWESS egen xbgrp = cut(xb), group(10) //Cutting xb into 10 equally spaced groups, generating 'xbgrp' frame put xbgrp CIF3 pseudo, into(summary) frame summary: collapse (mean) CIF3 pseudo, by(xbgrp) //Put variable 'xbgrp', 'CIF3' and 'pseudo' into a frame "summary", collapse the data within each group defined by 'xbgrp', calculate the mean of 'CIF3' and 'pseudo' within each group. frame summary: tw (scatter pseudo CIF3) (function y=x, xscale(range(0 1)) yscale(range(0 1))), aspectratio(1) addplot: line lowess CIF3, sort *************************************************** **BOOTSTRAP THE CALIBRATION PLOT PROCEDURE********* *************************************************** clear // Load your dataset use "the database file pathway", clear //<---change file directory frame create originaldata // Define the number of bootstrap iterations local B = 10 <---just a low number when setting up the code. // Initialize an empty frame to store bootstrap results frame create bootstrap_results // Define a loop to iterate through each bootstrap sample forval i = 1/B' { // Draw a bootstrap sample frame change originaldata bsample gen unique_id = _n // After the bsample one observation might have been drawn multiple times, stcrprep only allows single records, so we will trick stcrprep that each row is unique by adding a new identifier variable. // ************************************************ // * DRAW CALIBRATION PLOT ************************ // ************************************************ // Set survival time stset date_of_primary_event, id(unique_id) failure(primary_end_point == 1) origin(date_of_surgery) scale(365.24) // Generate pseudo observations for calibration plot stpci primary_end_point, at(3) competingvalues(2) // Calibration assessment at 3 years frame put id unique_id pseudo, into(pseudo_i') // Run stcrprep stset date_of_primary_event, id(unique_id) failure(primary_end_point == 1, 2) origin(date_of_surgery) scale(365.24) stcrprep, events(primary_end_point) keep(esmo_pre_op_sub extrauterine age colorscore vascularpattern_dicho emborder echogenicity bmi tumor_size) trans(1 2) gen event = primary_end_point == failcode // Fit Fine-and-Gray model stset tstop [pw=weight_c], failure(event) enter(tstart) noshow mfpa, cycles(3) dfdefault(3) select(0.157, i.esmo_pre_op_sub:1) alpha(0.25) : /// stcox i.esmo_pre_op_sub c.age i.colorscore i.extrauterine /// i.vascularpattern_dicho i.emborder i.echogenicity c.bmi c.tumor_size if failcode==1 /// , vce(cluster unique_id) // Predict baseline survival predict S0 if failcode==1, basesurv summ S0 if _t<=3 global S0_3 r(min)' bysort unique_id (_t): gen first = _n==1 predict xb if first, xb gen CIF3 = 1 -${S0_3}^exp(xb)

// Merge with pseudo observations
frame put unique_id CIF3 xb if first, into(predcif_i')
frame change pseudo_i'
frlink 1:1 unique_id, frame(predcif_i') gen(l)
frget CIF3 xb, from(l)

// Smooth 'pseudo' and 'CIF3' using LOWESS
lowess pseudo CIF3, gen(lowess) nograph

// Cut 'xb' into equally spaced groups
egen xbgrp = cut(xb), group(10)

// Store results in separate frames
frame put xbgrp CIF3 pseudo, into(summary_i')
frame summary_i': collapse (mean) CIF3 pseudo, by(xbgrp)
frame change summary_i'
frlink 1:1 xbgrp, frame(summary_i') gen(bootstrap_results)

}

}
`

Is it obvious to someone what I'm doing wrong? Can you help me fix it?

Do you have any other suggestions (and code!) on how to assess discrimination and calibration with resampling techniques in my setting?

Any help is greatly appreciated, and I will of course post any helpful replies also to Statalist.

BR Rasmus Green

• As this seems primarily to be a matter of coding rather than statistical matters, this question is off-topic here. I'd suggest posting this instead to Stack Overflow.
– EdM
Commented Mar 29 at 15:59
• Well, since I'm not even 100% sure that this is the right way of assessing calibration, and that I here haven't even discussed how to assess discrimination (which I also need help with), I guess it can fall into primarily a coding issue if my approach indeed is valid. I tried posting to Stack Overflow, but my question somehow triggered the spam filter, so it got blocked. Commented Mar 29 at 16:29
• Do you really want/need a Fine-Gray model instead of a simple two-state competing risks Cox model? See Sections 3 and 4 of the R competing risks vignette. Also, is there some advantage of the fractional polynomial approach over something like a restricted cubic regression spline for the continuous predictors?
– EdM
Commented Mar 29 at 17:22