Centiles of Standardized Survival Functions

Background

In a previous tutorial I used 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 obtained 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 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.

There 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 two 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)

. 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 standsurv will estimate these values more accurately with 95% confidence 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.

. standsurv, at1(female 0) at2(female 1) 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.8875488      2.0773963 |
     +-----------------------------------------+

. list med_female* in 1, ab(15)   

     +----------------------------------------------+
     | med_female   med_female_lci   med_female_uci |
     |----------------------------------------------|
  1. |  2.4249751        2.3197847        2.5349353 |
     +----------------------------------------------+

. list _contrast* in 1, ab(18)

     +----------------------------------------------------+
     | _contrast2_1   _contrast2_1_lci   _contrast2_1_uci |
     |----------------------------------------------------|
  1. |    .44477636          .31389716          .57565556 |
     +----------------------------------------------------+

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.

. standsurv, at1(female 0) at2(female 1) 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    .3236508   .38459929   .06094849 |
  |       30   .63821705   .78393263   .14571558 |
  |       40   1.1648032   1.4192444   .25444123 |
  |       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

Avatar
Paul C Lambert
Professor of Biostatistics

Related