New package stcomlist available from SSC

[Phil Clayton]

Thanks to Kit Baum, a new package stcomlist is now available from SSC.

stcomlist is intended to be a "competing risks" version of the official sts list command. It uses the SSC program stcompet to calculate cumulative incidence functions (CIFs) for each competing outcome, optionally in groups, and reports the CIFs along with confidence intervals.


Example:

. webuse hypoxia
(Hypoxia study)

. stset dftime, failure(failtype==1)

     failure event:  failtype == 1
obs. time interval:  (0, dftime]
 exit on or before:  failure

------------------------------------------------------------------------------
        109  total observations
          0  exclusions
------------------------------------------------------------------------------
        109  observations remaining, representing
         33  failures in single-record/single-failure data
    353.129  total analysis time at risk and under observation
                                                at risk from t =         0
                                     earliest observed entry t =         0
                                          last observed exit t =     8.454

. stcomlist, compet1(2) at(1 5) by(pelvicln)

            failure:  failtype == 1
 competing failures:  failtype == 2

    Time       CIF         SE     [95% Conf. Int.]
--------------------------------------------------
pelvicln=E
       1    0.1304     0.0702     0.0327    0.2972
       5    0.2609     0.0916     0.1062    0.4469
pelvicln=N
       1    0.1256     0.0415     0.0587    0.2191
       5    0.2217     0.0523     0.1290    0.3302
pelvicln=Y
       1    0.4091     0.1048     0.2085    0.6007
       5    0.5568     0.1079     0.3261    0.7364


            failure:  failtype == 2
 competing failures:  failtype == 1

    Time       CIF         SE     [95% Conf. Int.]
--------------------------------------------------
pelvicln=E
       1    0.0435     0.0425     0.0031    0.1824
       5    0.2174     0.0860     0.0791    0.3993
pelvicln=N
       1    0.0312     0.0217     0.0059    0.0965
       5    0.1582     0.0496     0.0763    0.2667
pelvicln=Y
       1    0.1364     0.0732     0.0341    0.3087
       5    0.1364     0.0732     0.0341    0.3087

stcomlist can be installed from SSC using:

ssc install stcomlist

[Giorgio Roccaro]

This is a great development. I have two questions -

1) Is there a way to use this on a competing risk model following multiple imputation of several covariates (ie: using an mi estimate precursor command)?

2) This gives the estimated CIF at requested time points for a particular variable of interest. Can this be used to determine the adjusted CIF at selected time points for the selected exposure of interest (incorporating several selected covariates)?

Thank you

[Phil Clayton]

Unfortunately no, stcomlist is only designed to work with stompet which estimates empirical (unadjusted) CIFs.

You can calculate the estimated CIF from an adjusted model using stcrreg. The CIF for a particular value of covariates (beta) with a linear prediction of xbeta is:
CIF = 1 - exp(-exp(xbeta)*[base cumulative subhazard function])

Thus, you might try something like:

stcrreg ...
predict double basecsh, basecshazard
predict double xb, xb
gen ci=1 - exp(-exp(xb)*basecsh)

See help stcrreg postestimation and the accompanying PDF manual for more on this.

[bln tfte]

Table of cumulative incidence - the cif lists:
Hi Phil,
I developed cif curves after I run the stcrreg as a post estimation. But I needed to present the cumulative incidence at years 1 2 5 and 10 after start of follow-up. You said stcomlist cannot be used and I found some difference between cif curve and table values. Any more help to the above?

[Phil Clayton]

It might be helpful to show an example of what you're doing. I would expect a difference between the table (unadjusted CIF) and the results from a stcrreg postestimation (adjusted CIF). You should be able to calculate the adjusted CIF at particular timepoints using the code I provided in post #3. The CIF at, say, 5 years is just the highest CIF in the dataset with a time <5 years (you can do eg "sum cif if _t<5" and look at the maximum value).

[bln tfte]

More information on what I am doing -stcrreg and stcurve: Thank so far.
I have two groups: cardiac deaths and non-cardiac deaths. I did this: stset survtime_USE, failure(Failure=1) scale(365.25) and run 'stcrreg' (stcrreg i.age i.sex i.countrybirth i.alcohol......, compet(Failure==2)) followed by 'stcurve' for both groups using the command: stcurve, cif at1(countrybirth==1) at2(countrybirth==2. I have cif plot for 'cardiac deaths' by country of birth and the same again for non-cardiac deaths. I want to attach a table showing the cummulative incidence at 1, 2, 5 and 10 years after first admission. So, I want to fill in the table.

Country of origin	1-year cumulative incidence estimates (95% CI)	2-year cumulative incidence estimates (95% CI)	5-year cumulative incidence Estimates (95% CI)	10-year cumulative incidence estimates (95% CI)
Migrant=1	xxx	xxx	xxx	xxx
Native =2	xxx	xxx	xxx	xxx

[Phil Clayton]

Perhaps this example will help:

* example from Stata documentation
webuse hypoxia, clear
stset dftime, failure(failtype==1)
stcrreg ifp tumsize i.pelnode, compete(failtype==2)
stcurve, cif at1(pelnode=0) at2(pelnode=1) name(stcurve)

* now recreate the graph and show the CIF at times 2,4,6,8
predict double basecsh, basecshazard
margins pelnode, predict(xb) atmeans post
gen cif_0=1 - exp(-exp(_b[0.pelnode])*basecsh)
gen cif_1=1 - exp(-exp(_b[1.pelnode])*basecsh)
line cif_0 cif_1 _t, sort c(J J)
foreach time of numlist 2(2)8 {
    sum cif_0 if _t<`time', meanonly
    display "time=`time' pelnode=0 CIF=`=r(max)'"
    sum cif_1 if _t<`time', meanonly
    display "time=`time' pelnode=1 CIF=`=r(max)'"
}

This shows you how to get adjusted CIFs at different times in different groups. I don't know how to get CIs and I don't think you actually can, because (similar to a Cox model) the Fine-Gray model is only a semi-parametric model - even if you apply the CIs from the margins output you'd still not have accounted for the uncertainty in the baseline cumulative subhazard function.

[bln tfte]

It worked perfectly well, Phill.
Very much Appreciate it!
Bln

[Allen Buxton]

stcomlist is a great idea, thank you. I wish Stata would calculate crude cumulative incidence, its variance, and Gray's test with the same results as the R package cmprsk and SAS proc lifetest.

Please note: In the case of an event on exactly an at() time, stcomlist lists the incorrect cumulative incidence because the at() point is exactly on the jump. Here is an illustration of this from chapter 10 of Marubini, E. and Valsecchi, M. G., Analysing Survival Data from Clinical Trials and Observational Studies, John Wiley &amp; Sons, 1995.

As an option, I have the following that has allowed me the use of R package cmprsk to do Gray's test and get estimates and confidence intervals with a form to plot these estimates from Stata Re: (Gray R. J. (1988), A class of K-sample tests for comparing the cumulative incidence of a competing risk, Ann. Stat. 16:1141-1154). Stata does the Fine and Gray regression with stcrreg. stgrays.ado is described here: https://www.statalist.org/forums/for...Here%20is%20an

---

.
. scalar sc_at= "at(1 12(24)240)"

. scalar sc_at= "at(15 30 60 113 120 240)"

.
. use dtatbl10_1.dta,clear

.
. // Table 10.3 crude cumulative incidence stcompet
. stset t,failure(event==1)

Survival-time data settings

         Failure event: event==1
Observed time interval: (0, t]
     Exit on or before: failure

--------------------------------------------------------------------------
         70  total observations
          0  exclusions
--------------------------------------------------------------------------
         70  observations remaining, representing
         20  failures in single-record/single-failure data
      4,018  total analysis time at risk and under observation
                                                At risk from t =         0
                                     Earliest observed entry t =         0
                                          Last observed exit t =       240

. stcomlist , compet1(2) `=sc_at' by(group)

            failure:  event == 1
 competing failures:  event == 2

    Time       CIF         SE     [95% Conf. Int.]
--------------------------------------------------
group=1:A
      15    0.0571     0.0392     0.0103    0.1672
      30    0.0857     0.0473     0.0220    0.2057 The CIF is 0.1143 at 30 weeks, not 0.0857.
      60    0.1727     0.0642     0.0700    0.3135
     113    0.2402     0.0745     0.1126    0.3940
     120    0.2739     0.0783     0.1359    0.4317
     240    0.3144     0.0833     0.1632    0.4778
group=2:B
      15    0.0286     0.0282     0.0022    0.1268
      30    0.1143     0.0538     0.0362    0.2423
      60    0.2057     0.0693     0.0906    0.3530
     113    0.2914     0.0832     0.1436    0.4569 The CIF is 0.3390 at 113 weeks, not 0.2914.
     120    0.3390     0.0874     0.1780    0.5077
     240    0.3390     0.0874     0.1780    0.5077


            failure:  event == 2
 competing failures:  event == 1

    Time       CIF         SE     [95% Conf. Int.]
--------------------------------------------------
group=1:A
      15    0.1429     0.0591     0.0522    0.2774 The CIF is 0.1714 at 15 weeks, not 0.1429.
      30    0.1714     0.0637     0.0696    0.3113
      60    0.2896     0.0773     0.1510    0.4438
     113    0.4562     0.0878     0.2811    0.6154
     120    0.4562     0.0878     0.2811    0.6154
     240    0.6046     0.0886     0.4104    0.7525 The CIF is 0.6856 at 240 weeks, not 0.6046.
group=2:B
      15    0.2286     0.0710     0.1076    0.3764
      30    0.4286     0.0836     0.2643    0.5831
      60    0.4895     0.0851     0.3164    0.6419
     113    0.6133     0.0888     0.4175    0.7607
     120    0.6133     0.0888     0.4175    0.7607
     240    0.6133     0.0888     0.4175    0.7607

. stset , clear

.
. // Table 10.3 crude cumulative incidence and
. // R cmprsk cuminc Gray's test
. frame copy default dtaA, replace
(note: frame dtaA not found)

. stset t, fail(event=1)

... output omitted

. qui stgrays group event _t , `=sc_at'

.
. cwf default

. frame default : use "C:\data\stgraystest.dta" , clear

. format status %40.0g

. list , noobs clean

      status        stat          pv   df  
    1:type 1   .22254759    .6371056    1  
    2:type 2   1.7557283   .18515823    1  

. scalar strpv1  = string(pv[1] , "%5.4f")

. scalar strpv2  = string(pv[2] , "%5.4f")

.
. frame default : use "C:\data\stgrayslist.dta" , clear

. drop if missing(est)
(2 observations deleted)

. format group_status %40.0g

. format est se lci* uci* %6.4f

. bysort group_status: list time est se lci* uci*  , noobs clean

------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
-> group_status = 1:A 1:type 1

    time      est       se    lci95    uci95  
      15   0.0571   0.0398   0.0099   0.1693  
      30   0.1143   0.0546   0.0354   0.2445  **
      60   0.1727   0.0653   0.0687   0.3161  
     113   0.2402   0.0765   0.1098   0.3981  
     120   0.2739   0.0808   0.1322   0.4366  
     240   0.3144   0.0872   0.1571   0.4852  

------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
-> group_status = 2:B 1:type 1

    time      est       se    lci95    uci95  
      15   0.0286   0.0286   0.0021   0.1292  
      30   0.1143   0.0549   0.0351   0.2453  
      60   0.2057   0.0713   0.0880   0.3575  
     113   0.3390   0.0980   0.1611   0.5269  **
     120   0.3390   0.0980   0.1611   0.5269  

------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
-> group_status = 1:A 2:type 2

    time      est       se    lci95    uci95  
      15   0.1714   0.0647   0.0684   0.3136  **
      30   0.1714   0.0647   0.0684   0.3136  
      60   0.2896   0.0787   0.1488   0.4466  
     113   0.4562   0.0902   0.2763   0.6194  
     120   0.4562   0.0902   0.2763   0.6194  
     240   0.6856   0.1190   0.3947   0.8579  **

------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
-> group_status = 2:B 2:type 2

    time      est       se    lci95    uci95  
      15   0.2286   0.0722   0.1059   0.3789  
      30   0.4286   0.0855   0.2608   0.5862  
      60   0.4895   0.0872   0.3120   0.6453  
     113   0.6133   0.0958   0.4008   0.7700  
     120   0.6133   0.0958   0.4008   0.7700  


. frame default : use "C:\data\stgraysplot.dta" , clear

.
. #delimit;
delimiter now ;
. twoway
> (connected est time if inlist(group_status,1), sort connect(stairstep) msymbol(none) lpattern("solid") lcolor(black)  lwidth(thin))
> (connected est time if inlist(group_status,2), sort connect(stairstep) msymbol(none) lpattern("dash") lcolor(black) lwidth(thin))
> (connected est time if inlist(group_status,3), sort connect(stairstep) msymbol(none) lpattern("solid") lcolor(blue) lwidth(thin))
> (connected est time if inlist(group_status,4), sort connect(stairstep) msymbol(none) lpattern("dash") lcolor(blue) lwidth(thin))
> ,
> xlabel(0(24)240, labsize(*1.2)) xmtick(##6)
> ylabel(0(0.1)0.7, labsize(*1.2) angle(0)) ymtick(##5)
> ytitle("crude cumulative incidence")
> xtitle("Weeks")
> title(`""', color(black) size(medium))
> xline( 15, lcolor(red%40) lwidth(vthin) lpattern("dash"))
> xline( 30, lcolor(red%40) lwidth(vthin) lpattern("dash"))
> xline(113, lcolor(red%40) lwidth(vthin) lpattern("dash"))
> legend(label(1 "group:A 1:type 1") label(2 "group:B 1:type 1") label(3 "group:A 2:type 2") label(4 "group:B 2:type 2")
>        rows(1) bm(zero) region(lw(none)) position(6) ring(1))
> text(.36 160 "Gray's test event type 1, p=`=strpv1'" , place(e) nobox just(left) margin(l+4 t+1 b+1) width(20) size(*.8))
> text(.70 160 "Gray's test event type 2, p=`=strpv2'" , place(e) nobox just(left) margin(l+4 t+1 b+1) width(20) size(*.8))
> note(`"R cmprsk"')
> scheme(stcolor)
> ;

. #delimit cr
delimiter now cr
. cwf dtaA