RMST of Standardized Survival Functions

Here I will show another useful measure from standardized survival functions. There have been several papers promoting the use of restricted mean survival time (RMST) in clinical trials. The arguments are (i) ease of interpretation (though I am not convinced a restricted mean is that easy to explain) and (ii) providing a simple summary in the presence of non-proportional hazards. See Royston and Parmar (2013) for a description of the use of the measure in RCTs.

The restricted mean survival time at time $t^{*}$ is defined as, $E [m i n (t, t^{*})]$ i.e. it is the mean up to some point $t^{*}$ . The treatment effect in a RCT can be defined as the difference in RMST between the randomized arms at time $t^{*}$ . The RMST can be estimated by calculating the area under the survival curve between 0 and $t^{*}$ . In an observational study where we need to take account of potential confounders, we can define the RMST of the standardized survival function as

$R M S T (t^{*} | X = x, Z) = \int_{0}^{t^{*}} E [S (t | X = x, Z)]$

and is estimated by

$\hat{R M S T} (t^{*} | X = x, Z) = \int_{0}^{t^{*}} \frac{1}{N} \sum_{i = 1}^{N} S (t | X = x, Z = z_{i})]$

Contrasts between exposure groups can either be differences or ratios,

$\hat{R M S T} (t^{*} | X = 1, Z) - \hat{R M S T} (t^{*} | X = 0, Z)$

$\frac{\hat{R M S T} (t^{*} | X = 1, Z)}{\hat{R M S T} (t^{*} | X = 0, Z)}$

Standardized RMST and contrasts is implemented in standsurv using the rmst option.

Example

I will use the Rotterdam Breast cancer data. The code below loads and stset’s the data and then fits a model using stpm3.

. use https://www.pclambert.net/data/rott3, clear
(Rotterdam breast cancer data (augmented with cause of death))

. stset os, f(osi==1) scale(12) exit(time 120)

Survival-time data settings

         Failure event: osi==1
Observed time interval: (0, os]
     Exit on or before: time 120
     Time for analysis: time/12

--------------------------------------------------------------------------
      2,982  total observations
          0  exclusions
--------------------------------------------------------------------------
      2,982  observations remaining, representing
      1,171  failures in single-record/single-failure data
 20,002.424  total analysis time at risk and under observation
                                                At risk from t =         0
                                     Earliest observed entry t =         0
                                          Last observed exit t =        10

. stpm3 i.hormon age enodes pr_1, scale(lncumhazard) df(4) eform nolog tvc(i.hormon) dftvc(3)

                                                        Number of obs =  2,982
                                                        Wald chi2(4)  = 615.95
Log likelihood = -2666.5968                             Prob > chi2   = 0.0000

-----------------------------------------------------------------------------------
                  |     exp(b)   Std. err.      z    P>|z|     [95% conf. interval]
------------------+----------------------------------------------------------------
xb                |
           hormon |
             yes  |   .8499802   .0963501    -1.43   0.152     .6806444    1.061445
              age |   1.013249   .0024115     5.53   0.000     1.008534    1.017987
           enodes |   .1132408    .011008   -22.41   0.000     .0935963    .1370084
             pr_1 |   .9061179   .0119267    -7.49   0.000     .8830411    .9297979
------------------+----------------------------------------------------------------
time              |
             _ns1 |  -27.09524   2.109681   -12.84   0.000    -31.23014   -22.96034
             _ns2 |   8.647725   1.122097     7.71   0.000     6.448455    10.84699
             _ns3 |  -1.072205   .0477674   -22.45   0.000    -1.165827   -.9785823
             _ns4 |  -.6930019   .0518048   -13.38   0.000    -.7945373   -.5914664
                  |
hormon#c._ns_tvc1 |
             yes  |   5.425507    3.92237     1.38   0.167    -2.262197    13.11321
                  |
hormon#c._ns_tvc2 |
             yes  |  -3.309698   2.096769    -1.58   0.114    -7.419291    .7998943
                  |
hormon#c._ns_tvc3 |
             yes  |  -.1256484    .195217    -0.64   0.520    -.5082667      .25697
                  |
            _cons |   .7984459   .1615956     4.94   0.000     .4817244    1.115168
-----------------------------------------------------------------------------------
Note: Estimates are transformed only in the first equation.

I have made the effect of our exposure, hormon, time-dependent using the tvc option to illustrate that we can have interactions etc with our exposure in our model. This is an interaction with time, i.e. non proportional hazards.

I first calculate the standardized survival curves where everyone is forced to be exposed and then unexposed.

. range timevar 0 10 100
(2,882 missing values generated)

. standsurv, surv timevar(timevar) ci frame(surv, replace)          ///
>            at1(hormon 0) at2(hormon 1) atvar(S_hormon0 S_hormon1)

. 
. frame surv {
.   twoway (area S_hormon0 timevar, sort fcolor(%30)) ///
>          , legend(off)                              ///
>          ylabel(0(0.1)1, format(%3.1f))             ///
>          ytitle("S(t)")                             ///
>          xtitle("Years from surgery")               ///
>          title("No treatment")                      ///
>          name(hormon0, replace)
.   
.   twoway (area S_hormon1 timevar, sort pstyle(p2line) fcolor(%30)) ///
>          , legend(off)                                             ///
>          ylabel(0(0.1)1, format(%3.1f))                            ///
>          ytitle("S(t)")                                            ///
>          xtitle("Years from surgery")                              ///
>          title("Treatment")                                        ///
>          name(hormon1, replace)
.                 
.   graph combine hormon0 hormon1, nocopies               
. }

The RMST at 10 years for each of the standardized survival functions is the area under the standardized survival curve, shown by the shaded areas in the graphs above.

I will now run standsurv again with the rmst option to estimate these togther with the difference in RMST. I only want the RMST at 10 years so create a variable t_rmst10 with only one observation, equal to 10.

. gen t10 = 10 in 1
(2,981 missing values generated)

. standsurv, rmst timevar(t10) ci frame(rmst, replace)          ///
>            at1(hormon 0) at2(hormon 1) atvar(rmst_h0 rmst_h1) ///
>            contrast(difference) contrastvar(rmstdiff)

I will first list the standardized RMST in both treatment groups.

. frame rmst: list t10 rmst_h0* rmst_h1* in 1, noobs abb(12) 

  +-------------------------------------------------------------------------------------+
  | t10     rmst_h0   rmst_h0_lci   rmst_h0_uci     rmst_h1   rmst_h1_lci   rmst_h1_uci |
  |-------------------------------------------------------------------------------------|
  |  10   7.5444253     7.4318217      7.658735   7.9386172      7.687098      8.198366 |
  +-------------------------------------------------------------------------------------+

The RMST at 10 years is 7.54 years in those not taking treatment and 7.94 years in those taking treatment. The 95% confidence intervals are also shown. As I used the contrast(difference) option I can look at the difference in RMST at 10 years.

. frame rmst: list t10    rmstdiff* in 1, noobs abb(12)

  +-----------------------------------------------+
  | t10    rmstdiff   rmstdiff_lci   rmstdiff_uci |
  |-----------------------------------------------|
  |  10   .39419189        .116234      .67214979 |
  +-----------------------------------------------+

The difference is 0.39 years (95% CI 0.12 to 0.67).

The RMST will vary by the choice of $t^{*}$ . A range of values of $t^{*}$ can be given and then plotted.

. range t_rmst 0 10 50
(2,932 missing values generated)

. standsurv, rmst timevar(t_rmst) ci frame(rmst2, replace)      ///
>            at1(hormon 0) at2(hormon 1) atvar(rmst_h0 rmst_h1) ///
>            contrast(difference) contrastvar(rmstdiff)

We can plot how the RMST changes and the difference in RMST changes as a function of $t^{\*}$ .

. frame rmst2 {
.   twoway (line rmst_h0 rmst_h1 t_rmst) ///
>          , legend(order(1 "No treatment" 2 "Treatment") cols(1) pos(11)) ///
>          ytitle("RMST (years)") ///
>          xtitle("Years from surgery") ///
>          name(RMST,replace)
.   
.   twoway (rarea rmstdiff_lci rmstdiff_uci t_rmst, color(blue%20)) ///
>          (line rmstdiff t_rmst, lcolor(blue))                     ///
>          , legend(off)                                            ///
>          ylabel(, format(%3.1f))                                  ///
>          ytitle("Difference in RMST (years)")                     ///
>          xtitle("Years from surgery")                             ///
>                    name(RMSTdiff, replace)
.                 
.   graph combine RMST RMSTdiff, nocopies
. }

References

Royston, P. Parmar, M. K. B. Restricted mean survival time: an alternative to the hazard ratio for the design and analysis of randomized trials with a time-to-event outcome. BMC medical research methodology 2013;13:152