Announcement

Collapse
No announcement yet.
X
  • Filter
  • Time
  • Show
Clear All
new posts

  • New package stcomlist available from SSC

    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:
    Code:
    . 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:
    Code:
    ssc install stcomlist

  • #2
    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

    Comment


    • #3
      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:

      Code:
      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.

      Comment


      • #4
        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?
        Last edited by bln tfte; 22 Aug 2019, 23:57.

        Comment


        • #5
          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).

          Comment


          • #6
            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

            Comment


            • #7
              Perhaps this example will help:

              Code:
              * 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.

              Comment


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

                Comment


                • #9
                  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


                  Code:
                  .
                  . 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
                  
                  ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
                  -&gt; 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  
                  
                  ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
                  -&gt; 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  
                  
                  ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
                  -&gt; 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  **
                  
                  ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
                  -&gt; 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
                  &gt; (connected est time if inlist(group_status,1), sort connect(stairstep) msymbol(none) lpattern("solid") lcolor(black)  lwidth(thin))
                  &gt; (connected est time if inlist(group_status,2), sort connect(stairstep) msymbol(none) lpattern("dash") lcolor(black) lwidth(thin))
                  &gt; (connected est time if inlist(group_status,3), sort connect(stairstep) msymbol(none) lpattern("solid") lcolor(blue) lwidth(thin))
                  &gt; (connected est time if inlist(group_status,4), sort connect(stairstep) msymbol(none) lpattern("dash") lcolor(blue) lwidth(thin))
                  &gt; ,
                  &gt; xlabel(0(24)240, labsize(*1.2)) xmtick(##6)
                  &gt; ylabel(0(0.1)0.7, labsize(*1.2) angle(0)) ymtick(##5)
                  &gt; ytitle("crude cumulative incidence")
                  &gt; xtitle("Weeks")
                  &gt; title(`""', color(black) size(medium))
                  &gt; xline( 15, lcolor(red%40) lwidth(vthin) lpattern("dash"))
                  &gt; xline( 30, lcolor(red%40) lwidth(vthin) lpattern("dash"))
                  &gt; xline(113, lcolor(red%40) lwidth(vthin) lpattern("dash"))
                  &gt; 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")
                  &gt;        rows(1) bm(zero) region(lw(none)) position(6) ring(1))
                  &gt; 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))
                  &gt; 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))
                  &gt; note(`"R cmprsk"')
                  &gt; scheme(stcolor)
                  &gt; ;
                  
                  . #delimit cr
                  delimiter now cr
                  . cwf dtaA
                  Click image for larger version

Name:	note02d_stataforumCIstcompet.png
Views:	1
Size:	105.9 KB
ID:	1736178
                  Attached Files

                  Comment

                  Working...
                  X