When I produce cumulative incidence curves (CIC) for the Fine and Gray model using stcrreg and stcurve, I noticed that the estimated plateau level of the CIC (for large durations) highly depends on the degree of censoring – see the graph below. E.g. if 60% of the spells are censored, the CIC attains a maximum of 0.18, if instead 40% are censored, the maximum CIC is 0.29. Censoring is random. In theory, however, the estimated plateau level of the CIC should NOT depend on the degree of censoring. Indeed, the estimated plateau level of the CIC is invariant, when I do the same exercise in R using "cmprsk". So, either, I do a basic coding error or there is an issue with the Stata implementation.

Array

Here is the Stata code generating the data plus the graph:

Code:
set seed 1000

foreach level in 40 60 80 100 { // level of censoring

clear
set obs 1000
gen i=_n
gen ra=round(runiform())
gen rb=1-ra
gen udur=-ln(uniform())

gen female=rnormal(0,1)
gen deutsch=rnormal(0,1)
gen VT_HE=rnormal(0,1)

*add censoring to data:
qui gen ran=runiform()
qui gen ran2=runiform()
qui replace ra=0 if ran>0.`level'
qui replace rb=0 if ran>0.`level'
qui replace udur=udur*ran2 if ran>0.`level'

stset udur, failure(ra==1)
stcrreg VT_HE deutsch female, compete(ra==0) iter(100)

stcurve, cif legend(off) xla(, grid) outfile("FG_`level'.dta", replace)
}

foreach level in 40 60 80 100 { // level of censoring
use "FG_`level'.dta", clear
rename ci1 ci`level'
save "FG_`level'.dta", replace
}

use "FG_40.dta", clear
foreach level in 60 80 100 {
joinby _t using FG_`level'.dta, unmatched(both)
drop _merge
}

sort _t
twoway (line ci40 _t)(line ci60 _t)(line ci80 _t)(line ci100 _t), legend(order(1 "60" 2 "40" 3 "20" 4 "0"))
graph export statalist.png, replace


foreach level in 40 60 80 100 {
erase "FG_`level'.dta"
}
Has anybody come across this problem before?
Thank you in advance!