Hazard of Standardized Survival Functions

This will be a short tutorial as the ideas are very simple. I have previously discussed standardized survival functions. In survival analysis we know that there is a simple mathematical transformation from hazard to survival function and vice versa. The idea here is to transform to a hazard function from the standardized survival function. Recall that a standardized survival funnction; $S_s(t|X=x,Z)$ is estimated by

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

If we apply the usual transformation from survival to hazard to function ($h(t) = \frac{-d}{dt}\log[S(t)]$) we get

$$ h_s(t|X=x,Z) = \frac{1}{N} \frac{\sum_{i=1}^{N}S(t|X=x,Z=z_i)h(t|X=x,Z=z_i)}{\sum_{i=1}^{N}S(t|X=x,Z=z_i)} $$

This is a weighted average of the $N$ individual hazard functions with weights equal to $S(t|X=x,Z=z_i)$, i.e. the predicted survival function for individual $i$ when forced to take a specific value of the exposure variable, $X$, but their observed values of confounding variables, $Z$.

This is implemented in stpm2_standsurv using the hazard 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 stpm2.

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

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

     failure event:  osi == 1
obs. time interval:  (0, os]
 exit on or before:  time 120
    t 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

. stpm2 hormon age enodes pr_1, scale(hazard) df(4) eform nolog tvc(hormon) dftvc(3)

Log likelihood = -2666.5999                     Number of obs     =      2,982

--------------------------------------------------------------------------------
               |     exp(b)   Std. Err.      z    P>|z|     [95% Conf. Interval]
---------------+----------------------------------------------------------------
xb             |
        hormon |   .8019893   .0741703    -2.39   0.017     .6690322     .961369
           age |   1.013249   .0024115     5.53   0.000     1.008534    1.017987
        enodes |   .1132406    .011008   -22.41   0.000     .0935961    .1370082
          pr_1 |   .9061179   .0119267    -7.49   0.000      .883041    .9297979
         _rcs1 |   2.644573   .0814503    31.58   0.000     2.489656    2.809129
         _rcs2 |   1.209479   .0379393     6.06   0.000      1.13736    1.286172
         _rcs3 |      1.014   .0162037     0.87   0.384     .9827339    1.046262
         _rcs4 |   .9961807   .0072731    -0.52   0.600     .9820273    1.010538
  _rcs_hormon1 |   1.003465   .0756175     0.05   0.963     .8656822    1.163176
  _rcs_hormon2 |    .891054    .056664    -1.81   0.070      .786637    1.009331
  _rcs_hormon3 |   1.025052   .0390804     0.65   0.516     .9512477    1.104583
         _cons |   1.103353   .1771893     0.61   0.540     .8054133    1.511508
--------------------------------------------------------------------------------
Note: Estimates are transformed only in the first equation.

I have made the effect of our exposure, hormon, time-dependent using the tvc option.

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)

. stpm2_standsurv, at1(hormon 0) at2(hormon 1) timevar(timevar) ci contrast(difference)

. 
. twoway  (rarea _at1_lci _at1_uci timevar, color(red%25)) ///
>         (rarea _at2_lci _at2_uci timevar, color(blue%25)) ///
>         (line _at1 timevar, sort lcolor(red)) ///
>         (line _at2  timevar, sort lcolor(blue)) ///
>         , legend(order(1 "No hormonal treatment" 2 "Hormonal treatment") ring(0) cols(1) pos(1)) ///
>         ylabel(0.5(0.1)1,angle(h) format(%3.1f)) ///
>         ytitle("S(t)") ///
>         xtitle("Years from surgery")

If I run stpm2_standsurv again with the hazard option I get the corresponding hazard functions of the standardized curves. This is the marginal hazard ratio (as a function of time).

. capture drop _at* _contrast*

. stpm2_standsurv, at1(hormon 0) at2(hormon 1) timevar(timevar) hazard ci contrast(ratio) per(1000)

Plot the standardized hazard functions.

. twoway (rarea _at1_lci _at1_uci timevar, color(red%30)) ///
>         (rarea _at2_lci _at2_uci timevar, color(blue%30)) ///
>         (line _at1 timevar, color(red)) ///
>         (line _at2 timevar, color(blue)) ///
>     , legend(off) ///
>     ylabel(,angle(h) format(%3.1f)) ///
>     xtitle("Years from surgery")         

I can’t explain the lower and then higher hazard for those on hormon therapy. Perhaps better adjustment for confounders would change this.

I can also plot the ratio of these two hazard functions with a 95% confidence interval.

. twoway (rarea _contrast2_1_lci _contrast2_1_uci timevar, color(red%30)) ///
>     (line _contrast2_1 timevar, color(red)) ///
>     , yscale(log) ///
>     ylabel(0.5 1 2 4 8 20 40, angle(h) format(%3.1f)) ///
>     xtitle("Years from surgery") ///
>     legend(off) ///
>     yscale(log) 

If I had used the difference argument of the contrast() option I would have obtained the absolute difference in the standardized hazard functions.

I am still thinking about the usefulness of this - in general I prefer the idea of standardized survival functions rather than the corresponding hazard function. However, it is harder to see how the risk of events changes over follow-up time with a cumulative measure (i.e. standardized survival).