Centiles of Standardized Survival Functions

Background

In the previous tutorial I used stpm2_standsurv to obtain standardized survival functions. In this tutorial I show the first of a number of different measures of the standardized survival function where I obtain centiles of the standardized survival function.

As a reminder a centile of a survival function can be obtanined by solving $S(t) = \alpha$ for $t$. For example, for the median survival time we set $\alpha = 0.5$, i.e. the 50th (per)centile. For simple parametric distributions, such as the Weibull, we can solve for $t$ analytically, but for more complex models the centile is obtained through iterative root finding techniques. In stpm2 I have used Brent’s root finder when evaluating centiles.

The centile of a standardized survival function is obtained by solving the following equation for t.

$$ E\left(S(t | X=x,Z\right) = \alpha $$

This is done through root finding (using Brent’s root finder) by solving,

$$ \frac{1}{N}\sum_{i=1}^N {S(t | X=x,Z)} - \alpha = 0 $$

Variances can be obtained using M-estimation .

Example

I use a colon cancer example. I first load and stset the data

. //use https://www.pclambert.net/data/colon, clear
. use c:/cansurv/data/colon, clear
(Colon carcinoma, diagnosed 1975-94, follow-up to 1995)

. drop if stage==0
(2,356 observations deleted)

. stset surv_mm, f(status=1,2) scale(12) exit(time 120)

     failure event:  status == 1 2
obs. time interval:  (0, surv_mm]
 exit on or before:  time 120
    t for analysis:  time/12

------------------------------------------------------------------------------
     13,208  total observations
          0  exclusions
------------------------------------------------------------------------------
     13,208  observations remaining, representing
      8,866  failures in single-record/single-failure data
 43,950.667  total analysis time at risk and under observation
                                                at risk from t =         0
                                     earliest observed entry t =         0
                                          last observed exit t =        10

I drop those with missing stage information (stage == 0). I am investigating all cause survival (status=1,2).

I fit a model that includes stage, sex and age (using a rectricted cubic splines). I assume proportional hazards, but if I relax this assusmption the syntax for stpm2_standsurv would be identical. Stage is classified as localised, regional and distant and is modelled using two dummy covariates with localised as the reference category.

. tab stage, gen(stage)

   Clinical |
   stage at |
  diagnosis |      Freq.     Percent        Cum.
------------+-----------------------------------
  Localised |      6,274       47.50       47.50
   Regional |      1,787       13.53       61.03
    Distant |      5,147       38.97      100.00
------------+-----------------------------------
      Total |     13,208      100.00

. gen female = sex==2

. rcsgen age, df(3) gen(agercs) center(60)
Variables agercs1 to agercs3 were created

. stpm2 stage2 stage3 female agercs*, scale(hazard) df(4) nolog eform

Log likelihood = -19665.932                     Number of obs     =     13,208

------------------------------------------------------------------------------
             |     exp(b)   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
xb           |
      stage2 |   1.758926   .0617187    16.09   0.000     1.642026    1.884149
      stage3 |   5.656322   .1390555    70.48   0.000      5.39024    5.935539
      female |   .8634681     .01891    -6.70   0.000     .8271894     .901338
     agercs1 |   .9958187   .0053355    -0.78   0.434     .9854161    1.006331
     agercs2 |   1.000015   .0000127     1.22   0.223     .9999906     1.00004
     agercs3 |   .9999654   .0000155    -2.24   0.025     .9999351    .9999957
       _rcs1 |   3.484444   .0392205   110.90   0.000     3.408415    3.562169
       _rcs2 |     1.1908   .0101285    20.53   0.000     1.171113    1.210817
       _rcs3 |   .9620673   .0050498    -7.37   0.000     .9522207    .9720158
       _rcs4 |   1.015331    .003491     4.43   0.000     1.008512    1.022197
       _cons |   .1615451   .0046341   -63.55   0.000     .1527132    .1708879
------------------------------------------------------------------------------
Note: Estimates are transformed only in the first equation.

The is a clear effect of stage with a hazard ratio of 5.66 for distant stage versus localised stage. Remember that I am modelling all cause survival and one would expect a cause-specific hazard ratio to be higher. The all-cause mortality rate for females is 14% lower than males.

I will now predict the standardized survival functions, one where I force all subjects to be male and one where I force everyone to be female.

. range tt 0 10 100
(13,108 missing values generated)

. stpm2_standsurv, at1(female 0) at2(female 1) timevar(tt) atvar(ms_male ms_female) ci

. twoway  (line ms_male ms_female tt, sort) ///
>         , yline(0.5, lpattern(dash) lcolor(black)) ///
>         yline(0.5, lpattern(dash) lcolor(black)) ///
>         xtitle("Years since diagnosis") ///
>         ytitle("S(t)", angle(h)) ///
>         ylabel(0(0.2)1, format(%3.1f) angle(h)) ///
>         legend(order(1 "Male" 2 "Female") ring(0) pos(1) cols(1))

The graph of the two standardised survival functions can be seen below.

As expected (given the hazard ratio) females have better survival than males. I have added a horizontal reference line at $S(t)=0.5$. Where this line crosses the survival curves gives the median survival time. Reading from the graph, this is just under 2 years for the males and just under 2.5 years for females. Using the centile option of stpm2_standsurv will estimate these values more accurately with 95% confidenec intervals. We are also interested in contrasts of the centiles, so use of the contrast option will calculate either a difference or ratio of the median survival times with a 95% confidence interval.

. stpm2_standsurv, at1(female 0) at2(female 1) timevar(tt) centile(50) ///
>                 atvar(med_male med_female) contrast(difference) ci

. list med_male* in 1,  ab(15)    

     +-----------------------------------------+
     |  med_male   med_male_lci   med_male_uci |
     |-----------------------------------------|
  1. | 1.9801987      1.8728516      2.0936987 |
     +-----------------------------------------+

. list med_female* in 1, ab(15)   

     +----------------------------------------------+
     | med_female   med_female_lci   med_female_uci |
     |----------------------------------------------|
  1. |  2.4249751        2.3000843        2.5566472 |
     +----------------------------------------------+

. list _contrast* in 1, ab(18)

     +----------------------------------------------------+
     | _contrast2_1   _contrast2_1_lci   _contrast2_1_uci |
     |----------------------------------------------------|
  1. |    .44477636          .30074772            .588805 |
     +----------------------------------------------------+

The median survival time is 1.98 years for males with a 95% CI (1.87 to 2.09). The median for females is 2.42 years (95% CI, 2.30 to 2.56). As I used the contrast option I also get the difference in the median of the standardised survival curves with a 95% CI. Thus the time at which 50% of females have died is 0.44 years more than the time at which 50% of males have died, 95% CI (0.30 to 0.59).

It is possible to predict for multiple centiles by passing a numlist to the centiles option. For example, the code below calculates centiles between 10 and 60 at 10 unit intervals.

. stpm2_standsurv, at1(female 0) at2(female 1) timevar(tt) centile(10(10)60) ///
>                 atvar(cen_males cen_females) contrast(difference) ci ///
>                 centvar(centiles) contrastvar(cendiff)

. list centile cen_males cen_females cendiff in 1/6, sep(0) noobs

  +----------------------------------------------+
  | centiles   cen_males   cen_fem~s     cendiff |
  |----------------------------------------------|
  |       10   .14919254   .16940399   .02021145 |
  |       20   .32365079   .38459929    .0609485 |
  |       30   .63821705   .78393263   .14571558 |
  |       40   1.1648032   1.4192444   .25444124 |
  |       50   1.9801987   2.4249751   .44477636 |
  |       60   3.4239223    4.277573   .85365069 |
  +----------------------------------------------+

We can then plot the difference in these various centiles.

. twoway (rarea cendiff_lci cendiff_uci centile, sort color(red%30)) ///
>         (line cendiff centile, color(black)) ///
>         , xtitle(centile) xlabel(,format(%3.0f)) ///
>         ytitle("Difference in centile") ylabel(0(0.2)1.2,format(%3.1f) angle(h)) ///
>         legend(off)

There are probably more innovative ways of presenting such data.

Acknowledgement

I would like to acknowledge David Druker of StataCorp who I discussed these ideas with at two Nordic Stata User group meetings. David has written a command that estimates centiles of standardized distributions using a two parameter gamma distribution which is available here.

References

Stefanski, L. & Boos, D. The Calculus of M-Estimation. The American Statistician 2002;56:29-38