Standardized survival functions using standsurv
Background
When we are performing data exploration on survival data we usually start with plotting Kaplan-Meier curves. In clinical trials with a survival outcome, one would nearly always expect to see a Kaplan-Meier curve plotted. They are simple to interpret (though there can be confusion when there are competing risks).
In observational studies, we expect that there will be confounding and we would usually adjust for these confounders in a Cox model. If you have read my other tutorials then you will know that I prefer fitting parametric models, but the choice is that not that important if all you want is an adjusted hazard ratio and that you are (i) happy with the proportional hazards assumption, (ii) believe you have included all relevant counfounders and (iii) made sensible modelling assumptions (non-linear effect, interactions etc).
Given we are happy with the model, an adjusted hazard ratio is reported. This is fine, but hazard ratios are more difficult to interpret and there are further problems when using hazard ratios as causal effects (Hernan 2010, Aaalen et al. 2015). Risks are much easier to interpret than rates and so quantifying the difference on the survival scale can be desirable.
Some statistical software implements something called “adjusted” survival curves, but it is not always clear what this means. For example, in Stata stcurve gives survival curves where certain covariates can be given specific values, but those not specified are given as mean values. Thus it gives a prediction for an individual who happens to have the mean values of each covariate. This is a prediction for an individual and may not reflect the average in the population. A more appropriate way is to average over the survival curves. For example, if we have 1000 individuals in our study we can predict a survival curve for each individual and then take the average of these 1000 curves. This is essentially what standsurv does.
Example
I use the Rotterdam breast cancer data and use all cause survival as the outcome. I restrict follow-up to 10 years after diagnosis using the option exit(time 10).
. 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 = 10The scale(12) option converts the times recorded in months to years.
I will explore differences between women who received hormonal treatment and those who did not. This is our exposure, but as this is observational study we know that any association we see may be due to confounding. But.. it is always good to start with a Kaplan-Meier plot
. sts graph, by(hormon) risktable ///
> legend(order(1 "No hormonal treatment" 2 "Hormonal treatment") ring(0) cols(1) pos(1)) ///
> ylabel(0(0.2)1,angle(h) format(%3.1f)) ///
> ytitle("S(t)") xtitle("Time since surgery")
Failure _d: osi==1
Analysis time _t: os/12
Exit on or before: time 120
This plot shows that those receiving hormal treatment had worse survival. If we fit a proportional hazards model with hormon as the only covariate we get the following,
. stpm3 i.hormon, df(4) scale(lncumhazard) eform nolog
Number of obs = 2,982
Wald chi2(1) = 25.19
Log likelihood = -2930.4853 Prob > chi2 = 0.0000
------------------------------------------------------------------------------
| exp(b) Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
xb |
hormon |
yes | 1.540779 .1326967 5.02 0.000 1.301464 1.824099
-------------+----------------------------------------------------------------
time |
_ns1 | -25.45166 1.863924 -13.65 0.000 -29.10488 -21.79844
_ns2 | 8.063828 .9957355 8.10 0.000 6.112222 10.01543
_ns3 | -.997684 .0432372 -23.07 0.000 -1.082427 -.9129407
_ns4 | -.6392856 .0466345 -13.71 0.000 -.7306874 -.5478837
_cons | -.5692268 .033194 -17.15 0.000 -.6342858 -.5041679
------------------------------------------------------------------------------
Note: Estimates are transformed only in the first equation.The hazard ratio indicates that there is a 54% higher mortality rate in those receiving hormonal therapy.
As this is an observational study we should not stop there and conclude that hormonal therapy is bad for you. We do not know if there are differences between women taking or not taking the hormonal therapy. If those who took teh treatment tended to be older and have more severe disease then we do not have a fair comparison.
In fact a simple tabulation shows this to be the case,
. tabstat age nodes, by(hormon)
Summary statistics: Mean
Group variable: hormon (Hormonal therapy)
hormon | age nodes
---------+--------------------
no | 54.09762 2.326523
yes | 62.54867 5.719764
---------+--------------------
Total | 55.05835 2.712274
------------------------------Those who received the hormonal therapy tended to be older and have more lymph node involvment. Thus, even if hormonal treatment did not have any effect on survival, we would expect to see a difference in such a simplistic analysis due to the type of people who receieved the treatment.
So, we now adjust for some covariates. To simplify things, I will assume proportional hazards and include the covariates age, nodes and progesterone receptor. Previous analyses of this data have found that transformation of the nodes variable (exp(-0.12*nodes)) and the progesterone variable (log(pr + 1)) model the non-linear effects of these variables fairly well and so I will use these transformed variables. The model is fitted below.
. stpm3 i.hormon age enodes pr_1, scale(lncumhazard) df(4) eform nolog
Number of obs = 2,982
Wald chi2(4) = 619.62
Log likelihood = -2668.4925 Prob > chi2 = 0.0000
------------------------------------------------------------------------------
| exp(b) Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
xb |
hormon |
yes | .7906432 .0715077 -2.60 0.009 .66221 .9439854
age | 1.013244 .0024119 5.53 0.000 1.008528 1.017983
enodes | .1132534 .0110135 -22.40 0.000 .0935998 .1370337
pr_1 | .9064855 .0119282 -7.46 0.000 .8834055 .9301685
-------------+----------------------------------------------------------------
time |
_ns1 | -25.90082 1.871965 -13.84 0.000 -29.5698 -22.23184
_ns2 | 7.980587 1.003724 7.95 0.000 6.013324 9.947851
_ns3 | -1.091126 .0461407 -23.65 0.000 -1.18156 -1.000691
_ns4 | -.70103 .0504635 -13.89 0.000 -.7999366 -.6021234
_cons | .801967 .161537 4.96 0.000 .4853603 1.118574
------------------------------------------------------------------------------
Note: Estimates are transformed only in the first equation.Things have now changed, the adjusted hazard ratio for hormon is 0.79 (95% CI 0.66 to 0.94) indicating a benficial effect. There is strong confounding here as we have gone from a significant harmful effect to a significant beneficial effect when adjusting for age, number of positive lymph nodes and progesterone receptor.
We could stop here, but would the fun be in that. Instead we will try to understand what this hazard ratio means in terms of survival. First I will replicate what stcurve does and then obtain the standardized survival functions.
The code below obtains the mean of the covariates age, enodes and pr_1 and puts these into a macro which can then be passed to the at options of stpm3’s predict command when predicting survival.
. foreach var in age enodes pr_1 {
2. summ `var', meanonly
3. local atopt `atopt' `var' `r(mean)'
4. }
. predict s0_covave s1_covave, surv ci timevar(0 10, step(0.1)) ///
> frame(f1, replace) ///
> at1(hormon 0 `atopt') ///
> at2(hormon 1 `atopt')
Predictions are stored in frame - f1I will then plot these curves
. frame f1 {
. twoway (line s0_covave tt, sort) ///
> (line s1_covave tt, sort) ///
> , legend(order(1 "No hormonal treatment" 2 "Hormonal treatment") ///
> ring(0) cols(1) pos(1)) ///
> ylabel(0(0.2)1,angle(h) format(%3.1f)) ///
> ytitle("S(t)") ///
> xtitle("Years since surgery")
. }
These are the survival functions for two “average” women (i.e. who had the average values of each covariate) with one receiving hormonal treatment and the other not receiving it. We can see that for such women the difference is in the way we expect, given the hazard ratio, with those on hormonal treatment having better survival. If I had categorial covariates then interpretation is awkward. For example, if we modelled sex (not appropriate for this data) then the prediction would be for someone who was part male and part female as we would be taking the average of a binary variable.
Standardized survival curves
We want to calculate the expected survival under two counterfactuals. One where everyone has hormonal treatment and one where everybody does not take it. See (Sjölander 2016) for a nice description of these issues and implementation in an R package, stdreg. We also want contrasts of these standardized curves, for example a difference in standardized survival curves.
\[ E\left(S(t | X=1,Z\right) - E\left(S(t | X=0,Z\right) \]
In the above, \(X\) is the exposure of interest and \(Z\) are the confounders. We are interested in the expectation over the distribution of \(Z\), with the key point being this distribution is forced to be the same for \(X=0\) and \(X=1\). If our model is sufficent for confounding control then the above formula gives an average causal effect.
To estimate the difference in the standardized curves we need to generate the two standardized survival curves. In each of these we predict as many survival curves as there are observations in the data set and then take the average of these curves. The only difference is that in one we make everyone be exposed (hormon=1) and in the other we make everybody be unexposed (hormon=0).
\[\frac{1}{N}\sum_{i=1}^{N}S(t|\mbox{hormon=1},\mbox{age}_i,\mbox{enodes}_i,\mbox{pr\_1}_i) - \frac{1}{N}\sum_{i=1}^{N}S(t|\mbox{hormon=0},\mbox{age}_i,\mbox{enodes}_i,\mbox{pr\_1}_i)\]
There are 2982 observation in the dataset (and in the model). Thus for each of the two standardized curves we need to predict 2982 survival curves and then take the average of these curves. We can make some computational efficiency savings by only estimating the survival curves a small number of time points. This could be a single time point, e.g. the survival at 5 years, or over a range. We often want to plot survival curves and about 50-100 time points is usually sufficient for plotting purposes.
To do all this (and more) we use the standsurv command. I will run the command and the explain the syntax.
. range tt 0 10 101
(2,881 missing values generated)
. standsurv S0 S1, surv timevar(tt) ci ///
> frame(f2, replace) ///
> at1(hormon 0) at2(hormon 1) ///
> contrast(difference) contrastvar(Sdiff)Each of the at() options creates a standardized survival curve. Here a covariate (or covariates) can be set to take specific values. Any covariates not specified keep their observed values. Thus we are just implementing the equation above. The timevar() option gives the name of the variable that gives the survival times in which to evaluate the survival function. I have defined a variable tt above to be 101 rows ranging from 0 to 10, i.e. in steps of 0.1 years. The ci option requests that confidence intervals be calculated. Standard errors are either obtained using the delta-method or M-estimation. The default is the delta-method (for standardized survival). The frame() options saves the predictions to a new frame. The standardized estmates will be named S0 and S1. The contrast() option asks for a comparison of the two survival curves with the difference argument asking to take differences in the standardized survival curves. By default at1() is the reference, i.e. the contrast will be at2()-at1(), but this can be changed using the atref() option. The contrast will be stored in the new variable Sdiff.
You do not have to specify the names of the new variables and if you do not, standsurv will create new variables, _at1, _at2, _contrast2_1. However, in general it is more useful to create sensible names for these varables. As the ci option was specified there will be lower and upper for the confidence intervals (95% by deafult) for each estimate.
Below I list the standardized curves at 10 years, followed by their difference.
. frame f2 {
. list S0* S1* if tt==10, noobs
+---------------------------------------------------------------------+
| S0 S0_lci S0_uci S1 S1_lci S1_uci |
|---------------------------------------------------------------------|
| .54380594 .52512603 .56315034 .60749077 .56414068 .654172 |
+---------------------------------------------------------------------+
. list Sdiff* if tt==10, noobs
+-----------------------------------+
| Sdiff Sdiff_lci Sdiff_uci |
|-----------------------------------|
| .06368483 .01748572 .10988394 |
+-----------------------------------+
. } Thus the average survival at 10 years when everyone is forced to be unexposed (not on hormonal treatment) is 0.54 and when everyone is exposed it is 0.61. The difference is 0.064. We have actually evaluated each function at 101 time points and so we can plot the estimates together with 95% confidence intervals.
. frame f2 {
. twoway (rarea S0_lci S0_uci tt, color(red%25)) ///
> (rarea S1_lci S1_uci tt, color(blue%25)) ///
> (line S0 tt, sort lcolor(red)) ///
> (line S1 tt, 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")
. }
And now we can plot the difference in standardized curves together with a 95% confidence interval.
. frame f2 {
. twoway (rarea Sdiff_lci Sdiff_uci tt, color(red%25)) ///
> (line Sdiff tt, sort lcolor(red)) ///
> , legend(off) ///
> ylabel(,angle(h) format(%3.2f)) ///
> ytitle("Difference in S(t)") ///
> xtitle("Years from surgery")
. }
The covariate distribution we are averaging over is a combination of those on and not on hormon treatment. It may also be of interest to restrict to the covariate distribution of the exposed or the unexposed. Assuming our model has controlled for confounding this will give the average causal effect (difference) in the exposed. All we need to do is to add an if statement.
. standsurv S0 S1 if hormon==1, surv timevar(tt) ci ///
> frame(f3, replace) ///
> at1(hormon 0) at2(hormon 1) ///
> contrast(difference) contrastvar(Sdiff)The resulting standardized curves can then be plotted.
. frame f3 {
. twoway (rarea S0_lci S0_uci tt, color(red%25)) ///
> (rarea S1_lci S1_uci tt, color(blue%25)) ///
> (line S0 tt, sort lcolor(red)) ///
> (line S1 tt, 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")
. } 
Note that these curves give higher survival. This is because on average those who received hormonal treatment were younger and had less severe disease.
We can also plot the difference in these standardized curves together with a 95% confidence interval.
. frame f3 {
. twoway (rarea Sdiff_lci Sdiff_uci tt, color(red%25)) ///
> (line Sdiff tt, sort lcolor(red)) ///
> , legend(off) ///
> ylabel(,angle(h) format(%3.2f)) ///
> ytitle("Difference in S(t)") ///
> xtitle("Years from surgery")
. }
References
Aalen O.O., Cook R.J. Røysland K. Does Cox analysis of a randomized survival study yield a causal treatment effect? Lifetime data analysis 2015;21:579-593.
Hernán M.A. The hazards of hazard ratios. Epidemiology 2010;21:13-15
Sjölander A. Regression standardization with the R package stdReg. European Journal of Epidemiology 2016;31:563–574