# 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 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 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.

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