How to save a Kaplan-Meier survival point estimate to a local macro?

[Jonathan Tward]

Hello STATA Community,

I am currently using version 15.1 of STATA.

I have some straightforward survival time data, (time in years, censor (1=failed, 0=censored) , group of interest). example as follows:

* Example generated by -dataex-. To install: ssc install dataex
clear
input double pfsyrs byte pfscensor str17 first_deftx_type
   4.802739726027397 1 "RT"     
   3.147945205479452 0 "Surgery"
                   0 0 "RT"     
   1.704109589041096 0 "RT"     
  4.5698630136986305 0 "RT"     
   3.734246575342466 0 "Surgery"
  3.2246575342465755 0 "Surgery"
   7.043835616438356 0 "Surgery"
   4.816438356164384 0 "Surgery"
   1.189041095890411 0 "Surgery"
   .6958904109589041 0 "Surgery"
   -.273972602739726 0 "RT"     
  .15616438356164383 0 "RT"     
  2.1890410958904107 0 "RT"     
  .11232876712328767 1 "Surgery"
   .4986301369863014 1 "RT"     
-.052054794520547946 0 "RT"     
  11.734246575342466 0 "RT"     
                   0 1 "Surgery"
  11.542465753424658 0 "RT"     
  12.126027397260273 0 "RT"     
   4.594520547945206 0 "RT"     
  1.6383561643835616 0 "RT"     
   .3863013698630137 1 "Surgery"
                   0 1 "RT"     
   .8493150684931506 1 "RT"     
  1.6821917808219178 1 "RT"     
   7.682191780821918 0 "RT"     
   .8712328767123287 0 "RT"     
  1.4712328767123288 0 "Surgery"
   1.515068493150685 0 "Surgery"
  1.1863013698630136 0 "Surgery"
   4.183561643835616 1 "RT"     
   10.35068493150685 0 "RT"     
  2.6082191780821917 1 "RT"     
   .6520547945205479 0 "Surgery"
  3.4547945205479453 0 "RT"     
   .1178082191780822 1 "Surgery"
   .8246575342465754 1 "Surgery"
  3.0027397260273974 0 "Surgery"
   .7068493150684931 0 "Surgery"
   2.802739726027397 0 "Surgery"
   4.561643835616438 1 "RT"     
    8.24931506849315 0 "RT"     
                 -.4 0 "RT"     
  .23835616438356164 0 "RT"     
 -.03287671232876712 0 "RT"     
  2.8054794520547945 0 "RT"     
  1.5068493150684932 0 "RT"     
   5.175342465753425 0 "Surgery"
  1.9205479452054794 0 "Surgery"
end

I can easily display the 5 and 10 year survival point estimates for the groups using the "sts list" command:

. stset pfsyrs, failure(pfscensor) scale(1)

     failure event:  pfscensor != 0 & pfscensor < .
obs. time interval:  (0, pfsyrs]
 exit on or before:  failure

------------------------------------------------------------------------------
      4,929  total observations
        568  observations end on or before enter()
------------------------------------------------------------------------------
      4,361  observations remaining, representing
      1,012  failures in single-record/single-failure data
 15,887.025  total analysis time at risk and under observation
                                                at risk from t =         0
                                     earliest observed entry t =         0
                                          last observed exit t =  89.83836



. sts list, by(first_deftx_type) failure at(5 10)

         failure _d:  pfscensor
   analysis time _t:  pfsyrs

              Beg.                      Failure       Std.
    Time     Total     Fail             Function     Error     [95% Conf. Int.]
-------------------------------------------------------------------------------
RT 
       5       442       51              0.0686    0.0097     0.0520    0.0903
      10       147       28              0.1511    0.0180     0.1194    0.1902
Surgery 
       5       800       95              0.0577    0.0061     0.0468    0.0711
      10       235       45              0.1550    0.0155     0.1272    0.1883
-------------------------------------------------------------------------------
Note: Failure function is calculated over full data and evaluated at indicated
      times; it is not calculated from aggregates shown at left.

However, what I would like to do is actually just save the value of the 5-year failure function for any particular group (in this case, for the "RT" group it is equal to 0.0686) as a local scalar so that I can use this discrete value in a program. Does anyone know how to save a failure function value from a discrete specified time (like 5 years) for a KM estimate in a local scalar?

Many Thanks,

Jonathan Tward

[William Lisowski]

The output of help sta list tells us of the saving() option.

. webuse stan3
(Heart transplant data)

. stset, noshow

. tempfile stsout

. sts list, at(10 40 to 160) by(posttran) saving(`"`stsout'"')

              Beg.                      Survivor      Std.
    Time     Total     Fail             Function     Error     [95% Conf. Int.]
-------------------------------------------------------------------------------
posttran=0 
      10        74       12              0.8724    0.0346     0.7858    0.9256
      40        31       11              0.6781    0.0601     0.5446    0.7801
      70        17        2              0.6126    0.0704     0.4603    0.7339
     100        11        1              0.5616    0.0810     0.3900    0.7022
     130        10        1              0.5054    0.0903     0.3199    0.6646
     160         7        1              0.4422    0.0986     0.2480    0.6204
posttran=1 
      10        16        1              0.9286    0.0688     0.5908    0.9896
      40        43        6              0.7391    0.0900     0.5140    0.8716
      70        45        9              0.6002    0.0841     0.4172    0.7423
     100        40        9              0.4814    0.0762     0.3271    0.6198
     130        38        1              0.4687    0.0752     0.3174    0.6063
     160        36        1              0.4561    0.0742     0.3076    0.5928
-------------------------------------------------------------------------------
Note: Survivor function is calculated over full data and evaluated at indicated
      times; it is not calculated from aggregates shown at left.

. use `"`stsout'"', clear

. list, sepby(posttran)

     +----------------------------------------------------------------------------+
     | posttran   time   begin   fail   survivor    std_err         lb         ub |
     |----------------------------------------------------------------------------|
  1. |        0     10      74     12   .8724239   .0345523   .7858152   .9256319 |
  2. |        0     40      31     11   .6780517   .0601395   .5445634   .7800649 |
  3. |        0     70      17      2   .6126396   .0704132   .4602586   .7338921 |
  4. |        0    100      11      1   .5615863   .0809651    .390011   .7021749 |
  5. |        0    130      10      1   .5054277   .0902676   .3199082   .6646314 |
  6. |        0    160       7      1   .4422492   .0986462   .2480202   .6203719 |
     |----------------------------------------------------------------------------|
  7. |        1     10      16      1   .9285714   .0688303   .5907668   .9896199 |
  8. |        1     40      43      6   .7390694   .0899507   .5140094    .871645 |
  9. |        1     70      45      9     .60024   .0840877   .4172456   .7422516 |
 10. |        1    100      40      9   .4813914   .0761764   .3271319   .6198243 |
 11. |        1    130      38      1   .4687232   .0752178   .3173997   .6063312 |
 12. |        1    160      36      1    .456055    .074244   .3076466   .5927731 |
     +----------------------------------------------------------------------------+

.

From here you can then create a scalar in the usual ways.

[Jonathan Tward]

William,

Thank you for your reply, but what you are referring to is saving the output to a file, and I still don't know how to extract a single piece of data from the list you suggested. the STATA help file for 'list' does not refer to any postestimation variables.

I am trying to store the specific value into a local scalar, so I can then pass it as an argument in STATA programming.

Perhaps I can state my question a different way for clarity,

"Does anyone know of a function that will compute the failure rate for stset data given the inputs of a single timepoint of interest?"

so from a programming point of view this hypothetical function that I can't identify using my example data might look like this to compute the 5-year failure risk for the "RT" group....

stset pfsyrs, failure(pfscensor) scale(1)


local five_year_failure = SOMEFUNCTION (timepoint_of_interest=5 if first_deftx_type=="RT")

display "the computed 5 year failure probability for RT is `five_year_failure' "

would echo:

the computed 5 year failure probability for RT is 0.0686

The basic need for this function is that sts list does not have any post estimation matrices that stores the results.

Anyone have a way to do this?

[William Lisowski]

The basic need for this function is that sts list does not have any post estimation matrices that stores the results.

Instead it creates a dataset with the results, which you can use to create scalars and macros as in the code below.

. webuse stan3
(Heart transplant data)

. stset, noshow

. tempfile stsout

. sts list, failure at(10 40 to 160) by(posttran) saving(`"`stsout'"')

              Beg.                      Failure       Std.
    Time     Total     Fail             Function     Error     [95% Conf. Int.]
-------------------------------------------------------------------------------
posttran=0 
      10        74       12              0.1276    0.0346     0.0744    0.2142
      40        31       11              0.3219    0.0601     0.2199    0.4554
      70        17        2              0.3874    0.0704     0.2661    0.5397
     100        11        1              0.4384    0.0810     0.2978    0.6100
     130        10        1              0.4946    0.0903     0.3354    0.6801
     160         7        1              0.5578    0.0986     0.3796    0.7520
posttran=1 
      10        16        1              0.0714    0.0688     0.0104    0.4092
      40        43        6              0.2609    0.0900     0.1284    0.4860
      70        45        9              0.3998    0.0841     0.2577    0.5828
     100        40        9              0.5186    0.0762     0.3802    0.6729
     130        38        1              0.5313    0.0752     0.3937    0.6826
     160        36        1              0.5439    0.0742     0.4072    0.6924
-------------------------------------------------------------------------------
Note: Failure function is calculated over full data and evaluated at indicated
      times; it is not calculated from aggregates shown at left.

. use `"`stsout'"' if posttran==1 & time==10, clear

. list

     +----------------------------------------------------------------------------+
     | posttran   time   begin   fail    failure    std_err         lb         ub |
     |----------------------------------------------------------------------------|
  1. |        1     10      16      1   .0714286   .0688303   .0103802   .4092331 |
     +----------------------------------------------------------------------------+

. local ten_year_failure = failure[1]

. display "the computed 10 year failure probability for posttran=1 is `ten_year_failure'"
the computed 10 year failure probability for posttran=1 is .0714285746216774

. scalar tyf = failure[1]

. display tyf
.07142857

.

[Jonathan Tward]

William,

PERFECT! thank you.