stcrprep - computational benefits

When using stcrprep there are some computational benefits when compared to using Stata’s inbuilt stcrreg. One reason for this is that everytime you fit a model using stcrreg you the probability of censoring weights are calculated and the data must be expanded (in the background) when maximising the likelihood. When using stcrprep you only need to do this once.

I have run some timings. If I fit a simple model to the embt1 data with risk score as the only covariate (2 dummy variables) then these are the timings no my current work laptop (Intel i5 - running Stata 15 MP2).

First I load and stset the data.

. use https://www.pclambert.net/data/ebmt1_stata.dta, clear

. stset time, failure(status==1) scale(365.25) id(patid) noshow

Now, stcrreg can be used

. timer clear

. timer on 1

. stcrreg i.score, compete(status==2) nolog noshow

Competing-risks regression                       No. of obs       =      1,977
                                                 No. of subjects  =      1,977
Failure event  : status == 1                     No. failed       =        456
Competing event: status == 2                     No. competing    =        685
                                                 No. censored     =        836

                                                 Wald chi2(2)     =       9.87
Log pseudolikelihood = -3333.3217                Prob > chi2      =     0.0072

                              (Std. Err. adjusted for 1,977 clusters in patid)
------------------------------------------------------------------------------
             |               Robust
          _t |        SHR   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       score |
Medium risk  |   1.271221   .1554323     1.96   0.050     1.000333    1.615465
  High risk  |   1.769853   .3238535     3.12   0.002     1.236465    2.533337
------------------------------------------------------------------------------

. timer off 1

. timer list
   1:     23.32 /        1 =      23.3150

This takes 23.3 seconds to fit.

I now reload and stset the data, but this time declaring both status=1 and status=2 as events.

. use https://www.pclambert.net/data/ebmt1_stata.dta, clear

. stset time, failure(status==1,2) scale(365.25) id(patid)

We can now run stcrprep.

. timer on 2

. stcrprep, events(status) keep(score) trans(1)   

. timer off 2

. timer list 2
   2:      6.22 /        1 =       6.2240

This takes 6.2 seconds to run. However, this only restructures the data and calculates the weights. To fit the model, we first generate the event indicator and use stset.

. gen      event = status == failcode

. stset tstop [iw=weight_c], failure(event) enter(tstart) 

We use stcox to fit the model.

. timer on 3

. stcox i.score

         failure _d:  event
   analysis time _t:  tstop
  enter on or after:  time tstart
             weight:  [iweight=weight_c]

Iteration 0:   log likelihood = -3338.1244
Iteration 1:   log likelihood = -3333.4173
Iteration 2:   log likelihood = -3333.3113
Iteration 3:   log likelihood = -3333.3112
Refining estimates:
Iteration 0:   log likelihood = -3333.3112

Cox regression -- Breslow method for ties

No. of subjects =       72,880                  Number of obs    =      72,880
No. of failures =          456
Time at risk    =   6026.27434
                                                LR chi2(2)       =        9.63
Log likelihood  =   -3333.3112                  Prob > chi2      =      0.0081

------------------------------------------------------------------------------
          _t | Haz. Ratio   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       score |
Medium risk  |   1.271235   .1593392     1.91   0.056     .9943389    1.625238
  High risk  |   1.769899   .3219273     3.14   0.002     1.239148     2.52798
------------------------------------------------------------------------------

. timer off 3

. timer list
   1:     23.32 /        1 =      23.3150
   2:      6.22 /        1 =       6.2240
   3:      1.41 /        1 =       1.4060

This takes 1.4 seconds to run giving a combined total of 7.6 seconds. What is important is that if we want to fit other models (including other covariates etc), then we do not need to run stcrprep again.

To assess the time on larger data I have expanded the data by 20 times and added a small random number to each time, so that there are no ties. I used the following code.

expand 20
replace time = time + runiform()*0.0001
replace patid = _n

This leads to 19,770 indviduals in the analysis. The fact that there are no ties is perhaps a little unrealistic in a dataset this size, but this is still a usefull assessment of computational speed. The same analysis as above on this larger dataset gave the following times.

I think this really highlights the benfits of restructuring the data and using stcox in terms of computational time. Unless there is need to recalculate the probability of censoring weights, there is no need to do this every time you fit a model. Thus, in this case an stcrreg model takes almost 35 minutes, whilst the same model using stcox after using stcrprep takes only 46 seconds.

It is worthwhile noting that Stata’s implementation of Fine and Grays proportional subhazards model using stcrreg seems particularly slow. If I fit the model in R using crr the model fitted to the expanded data it only takes 370 seconds compared to 2066 in Stata.

There are other benefits with using stcox to fit the subhazards model, mainly because we can now use many of the other commands and extensions associated with stcox. I will discuss these in other tutorials.