Direct modelling of marginal relative survival models
Background
When using relative survival we may be intested in both estimating relative survival conditional on specific covariate patterns, for example,
for a male aged 70 diagnosed in 2018 with localized cancer, or marginal relative survival where we may interested in an average effect in a population or making comparisons
where we average over the same covariate (confounder) patterns. My standsurv command can be used to estimate marginal relative survival after fitting a model conditional on covariates.
The estimand of interest, is marginal relative survival. Consider a set of covariates, $\mathbf{X}_i$, for the $i^{th}$ individual that may affect the rate of death from the cancer under study and the rate of death from other causes. The all cause rate of death, $h(t|\mathbf{X}_i)$, can be partitioned into two components,
$$ h(t|\mathbf{X}_i) = h^*(t|\mathbf{X}_i) + \lambda(t|\mathbf{X}_i) $$
where $h^*(t|\mathbf{X}_i)$ is the expected mortaliity rate and $\lambda(t|\mathbf{X}_i)$ is the excess mortality rate for the $i^{th}$ individual. The relative survival for covariate pattern, $\mathbf{X}_i$ is,
$$ R(t|\mathbf{X}_i) = \int_0^t {\lambda(u|\mathbf{X}_i) du} $$
The marginal relative survival involves taking the expectation of $R(t|\mathbf{X}_i)$ over covariate pattern, $\mathbf{X}_i$,
$$ R^m(t|\mathbf{X}) = E_{\mathbf{X}}\left[R(t|\mathbf{X})\right] \tag{Equation 1} $$
Note that in the above for simplicity, I assume the same covariates act on the expected and excess mortality rates, but this is not a requirement
Example
I use the Melanoma data, restricting to those diagnosed in the later calendar perdiod, 1985-1994.
I restrict follow-up to 10 years after diagnosis using the exit() option.
. use https://pclambert.net/data/melanoma.dta if year8594 == 1
(Skin melanoma, diagnosed 1975-94, follow-up to 1995)
. stset surv_mm, failure(status=1,2) id(id) exit(time 120.5) scale(12)
Survival-time data settings
ID variable: id
Failure event: status==1 2
Observed time interval: (surv_mm[_n-1], surv_mm]
Exit on or before: time 120.5
Time for analysis: time/12
--------------------------------------------------------------------------
4,744 total observations
0 exclusions
--------------------------------------------------------------------------
4,744 observations remaining, representing
4,744 subjects
1,401 failures in single-failure-per-subject data
22,003.417 total analysis time at risk and under observation
At risk from t = 0
Earliest observed entry t = 0
Last observed exit t = 10.04167
I will first estimate the non-parametric estimate of marginal relative survival using stpp, so we have something to compare our model based estimates to.
. stpp R_pp using https://pclambert.net/data/popmort.dta, ///
> agediag(age) datediag(dx) pmother(sex)
. preserve
. rename _t PP_t
. keep R_pp* PP_t
. save PP, replace
file PP.dta saved
. restore
I have saved the Pohar Perme estimates in a dataset, so I can merge them in after using mrsprep.
I will now fit some relative survival models, but first I need to merge in the expected mortality rates at the event/censoring times.
. // conditional model (no covariates)
. gen _age = floor(min(age + _t,99))
. gen _year = floor(year(dx + _t*365.24))
. merge m:1 _age _year sex using https://pclambert.net/data/popmort.dta, ///
> keep(match master)
Result Number of obs
-----------------------------------------
Not matched 0
Matched 4,744 (_merge==3)
-----------------------------------------
Now I will fit a flexible parametric relative survival model with no covariate using stpm2. I will then predict the estimated relative survival.
. stpm2, scale(hazard) df (5) bhazard(rate)
Iteration 0: log likelihood = -4364.5376
Iteration 1: log likelihood = -4279.6772
Iteration 2: log likelihood = -4279.5867
Iteration 3: log likelihood = -4279.5864
Log likelihood = -4279.5864 Number of obs = 4,744
------------------------------------------------------------------------------
| Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
xb |
_rcs1 | .7654828 .0324456 23.59 0.000 .7018906 .8290749
_rcs2 | .1834916 .029522 6.22 0.000 .1256295 .2413537
_rcs3 | .0572084 .017442 3.28 0.001 .0230227 .091394
_rcs4 | .0242512 .0122991 1.97 0.049 .0001453 .048357
_rcs5 | .0120989 .008138 1.49 0.137 -.0038512 .028049
_cons | -2.037147 .0457164 -44.56 0.000 -2.12675 -1.947545
------------------------------------------------------------------------------
. range tt_cond 0 10 101
(4,643 missing values generated)
. predict s_cond, surv timevar(tt) ci
I can now compare the model based and the non-parametric estimates.
. twoway (rarea R_pp_lci R_pp_uci _t, sort connect(stairstep) color(${C1ci})) ///
> (line R_pp _t, sort connect(stairstep) color(${C1})) ///
> (line s_cond* tt_cond, color(${C2} ${C2} ${C2}) lpattern(solid dash dash)) ///
> , legend(order(2 "Pohar Perme" 3 "stpm2 model without covariates") ///
> ring(0) cols(1) pos(7)) ///
> ylabel(0.6(0.1)1, format(%3.1f)) ///
> ytitle("Marginal relative survival") ///
> xtitle("Years from diagnosis") ///
> name(int_stand, replace)
It can be seen that there is disagrement between the non-parametric estimate and the model based estimate. This is not good and differs in what we would expect to see in a standard survival model.
For example, if I fit an all-cause stpm2 survival model without covariates I get the following graph comparing the model based estmates with a Kaplan Meier estimate.
There is now near perfect agreement. I explain in the next section why there is disagreement between the model based and non-parametric estimate.
Why is there disagreement when a model with no covariates is fitted.
Consider the relative survival model fitted when not including any covariates.
$$ h(t|\mathbf{X}_i) = h^*(t|\mathbf{X}_i) + \lambda(t) $$
In this model the excess mortality is assumed to be exactly the same for each individual. In this model the all cause mortality rate varies between individuals only through variation in expected (other cause) mortality rates and the excess (cancer) mortality rate is assumed to be the same for all individuals. This is different from the definition in Equation 1 where relative survival is allowed to vary between individuals. Assuming that the excess mortality is the same over age, sex etc is a very strong assumption and almost certainly not true.
Regression standardization
Regression standardization can be used in the relative survival framework. This means that we should include all covariates that affect expected mortality rates in the model. In the case of the Melanoma data this is age, sex and calendar year
I will fit a model that uses restricted cubic splines to model the effect of age at diagnosis and also relax the proportional hazards assumption for the effect of age by allowing an interaction with time. The model will include sex and calendar years as these both impact the expected mortality rates. I will allow the effect of sex to be time-dependent (non-proportional), and model the effect of year of diagnosis using restricted cubic splines. A key point here is that various modelling choices need to be made, for example, I have chosen not to include interactions between any of the covariates. Different modelling choices will result in different estimates.
. gen female = sex==2
. rcsgen age , df(3) gen(agercs) orthog
Variables agercs1 to agercs3 were created
. rcsgen yydx , df(3) gen(yearrcs) orthog
Variables yearrcs1 to yearrcs3 were created
.
. stpm2 agercs* female yearrcs*, scale(hazard) df(5) bhazard(rate) ///
> tvc(agercs* female) dftvc(3)
Iteration 0: log likelihood = -4291.7959
Iteration 1: log likelihood = -4221.2221
Iteration 2: log likelihood = -4216.8452
Iteration 3: log likelihood = -4216.714
Iteration 4: log likelihood = -4216.7071
Iteration 5: log likelihood = -4216.7071
Log likelihood = -4216.7071 Number of obs = 4,744
---------------------------------------------------------------------------------
| Coefficient Std. err. z P>|z| [95% conf. interval]
----------------+----------------------------------------------------------------
xb |
agercs1 | .4014643 .0509245 7.88 0.000 .3016541 .5012744
agercs2 | -.0241269 .0519829 -0.46 0.643 -.1260115 .0777577
agercs3 | -.077221 .0497178 -1.55 0.120 -.1746662 .0202241
female | -.4660131 .0871942 -5.34 0.000 -.6369105 -.2951157
yearrcs1 | .0378156 .0466302 0.81 0.417 -.0535778 .129209
yearrcs2 | -.0840698 .0434742 -1.93 0.053 -.1692778 .0011382
yearrcs3 | .0226007 .0427815 0.53 0.597 -.0612495 .1064508
_rcs1 | .8485124 .0520516 16.30 0.000 .7464932 .9505316
_rcs2 | .2333051 .0517937 4.50 0.000 .1317914 .3348189
_rcs3 | .0467025 .0240512 1.94 0.052 -.0004371 .093842
_rcs4 | .0323118 .0137245 2.35 0.019 .0054122 .0592113
_rcs5 | .0194916 .0086331 2.26 0.024 .0025711 .036412
_rcs_agercs11 | .0222324 .0489672 0.45 0.650 -.0737416 .1182063
_rcs_agercs12 | .0379927 .0464296 0.82 0.413 -.0530077 .128993
_rcs_agercs13 | .0350472 .0215183 1.63 0.103 -.0071278 .0772223
_rcs_agercs21 | -.0939309 .0459412 -2.04 0.041 -.183974 -.0038878
_rcs_agercs22 | -.0120687 .0433049 -0.28 0.780 -.0969448 .0728074
_rcs_agercs23 | .0225369 .020735 1.09 0.277 -.0181029 .0631767
_rcs_agercs31 | -.0494584 .0434172 -1.14 0.255 -.1345546 .0356378
_rcs_agercs32 | .0334398 .0398716 0.84 0.402 -.0447071 .1115867
_rcs_agercs33 | -.0040932 .0201338 -0.20 0.839 -.0435547 .0353684
_rcs_female1 | -.0637989 .067838 -0.94 0.347 -.196759 .0691611
_rcs_female2 | -.0909384 .0607976 -1.50 0.135 -.2100995 .0282228
_rcs_female3 | .0296021 .0326255 0.91 0.364 -.0343426 .0935468
_cons | -1.720422 .0584158 -29.45 0.000 -1.834914 -1.605929
---------------------------------------------------------------------------------
. standsurv, surv timevar(tt) ci atvar(mrs_cond)
After fitting the model I have used standsurv to obtain the estimate of marginal relative survival using regression standardization. This predicts a relative survival
function for each individual conditional on their observed covariate pattern and takes the average of these curves. In this case there are 4,744
individuals in the study and so the estimated marginal relative survival is an average of 4,744 different survival curves.
I can now compare the model based estimate, based on regression standardization, and the non-parametric Pohar estimate.
. twoway (rarea R_pp_lci R_pp_uci _t, sort connect(stairstep) color(${C1ci})) ///
> (line R_pp _t, sort connect(stairstep) color(${C1})) ///
> (line mrs_cond* tt_cond, color(${C2} ${C2} ${C2}) lpattern(solid dash dash)) ///
> , legend(order(2 "Pohar Perme" 3 "Regression standardization after stpm2 model") ///
> ring(0) cols(1) pos(7)) ///
> ylabel(0.6(0.1)1, format(%3.1f)) ///
> ytitle("Marginal relative survival") ///
> xtitle("Years from diagnosis") ///
> name(int_stand_standsurv, replace)
There is now good agreement between the estimate based on regression standardizion and the non-parametric estimate. I will now move on to describing how to model marginal relative survival directly, so we can make fewer modelling decisions if we are only interested in the estimation of marginal relative survival.
Using mrsprep to prepare data for fitting a marginal modelling
In order to directly model marginal relative survival I will run mrsprep. This does two things,
- It calculates time-dependent weights which are the inverse of the expected survival. These are needed as the estimand of interest is in the net world, where it is not possible to die from causes other than the cancer under study. However, we have data in the real world and as follow-up time increases we have fewer at risk and fewer deaths than we would see in the net world. The weights are based on the same idea as the weights used in the Pohar Perme non-parametric estimate.
- At each event time it calculates the weighted mean mortality (hazard) rate for those still at risk. The weights are based on the inverse of expected survival among those at risk. The weighted mean is needed as the marginal relative survival is of interest. See the paper for more details.
The code for mrsprep is shown below.
. mrsprep using https://pclambert.net/data/popmort.dta ///
> , pmother(sex) agediag(age) datediag(dx) ///
> breaks(0(0.2)10)
mrsprep needs the filename of where the expected mortality rates are stored.
It requires the name of the variable for age at diagnosis and the name
of the variable for date of diagnosis. It also needs the name of variables
other than age and calendar year that the expected mortality rates are stratified by, in this
case this is just sex. The final option is breaks(0(0.2)10). This splits the
time scale into intervals, each of width 0.2 years. The weights are calculated at the mid-point of
each interval. This is an approximation, greater precision can be obtained without
narrower intervals, but the expanded dataset becomes larger. See the paper for
a sensitivity analysis for different interval widths.
Below is a listing for the first two individuals in the dataset.
. list id age sex tstart tstop wt meanhazard event if inlist(id,51,574), ///
> noobs sepby(id) abbrev(13)
+--------------------------------------------------------------------------+
| id age sex tstart tstop wt meanhazard_wt event |
|--------------------------------------------------------------------------|
| 51 86 1 0 .2 1.017644 999 0 |
| 51 86 1 .2 .4 1.0556412 999 0 |
| 51 86 1 .4 .6 1.0968963 999 0 |
| 51 86 1 .6 .8 1.1378526 999 0 |
| 51 86 1 .8 1 1.1783593 999 0 |
| 51 86 1 1 1.2 1.2222716 999 0 |
| 51 86 1 1.2 1.375 1.2674326 .02960256 1 |
|--------------------------------------------------------------------------|
| 574 69 2 0 .2 1.0039818 999 0 |
| 574 69 2 .2 .4 1.0119931 999 0 |
| 574 69 2 .4 .6 1.0206144 999 0 |
| 574 69 2 .6 .8 1.0298603 999 0 |
| 574 69 2 .8 1 1.0391899 999 0 |
| 574 69 2 1 1.2 1.0486292 999 0 |
| 574 69 2 1.2 1.4 1.0581798 999 0 |
| 574 69 2 1.4 1.6 1.0673769 999 0 |
| 574 69 2 1.6 1.8 1.0762101 999 0 |
| 574 69 2 1.8 2 1.085564 999 0 |
| 574 69 2 2 2.2 1.0958287 999 0 |
| 574 69 2 2.2 2.4 1.1065721 999 0 |
| 574 69 2 2.4 2.6 1.1174208 999 0 |
| 574 69 2 2.6 2.8 1.1283759 999 0 |
| 574 69 2 2.8 3 1.1394383 999 0 |
| 574 69 2 3 3.2 1.1510196 999 0 |
| 574 69 2 3.2 3.4 1.1631333 999 0 |
| 574 69 2 3.4 3.6 1.1753744 999 0 |
| 574 69 2 3.6 3.8 1.1877444 999 0 |
| 574 69 2 3.8 4 1.2002446 999 0 |
| 574 69 2 4 4.2 1.2135969 999 0 |
| 574 69 2 4.2 4.4 1.2278268 999 0 |
| 574 69 2 4.4 4.6 1.2422236 999 0 |
| 574 69 2 4.6 4.8 1.2567893 999 0 |
| 574 69 2 4.8 5 1.2715257 999 0 |
| 574 69 2 5 5.2 1.2863285 999 0 |
| 574 69 2 5.2 5.4 1.3011962 999 0 |
| 574 69 2 5.4 5.6 1.3162357 999 0 |
| 574 69 2 5.6 5.8 1.331449 999 0 |
| 574 69 2 5.8 6 1.3468382 999 0 |
| 574 69 2 6 6.2 1.3633453 999 0 |
| 574 69 2 6.2 6.4 1.381007 999 0 |
| 574 69 2 6.4 6.6 1.3988974 999 0 |
| 574 69 2 6.6 6.8 1.4170196 999 0 |
| 574 69 2 6.8 7 1.4353766 999 0 |
| 574 69 2 7 7.2 1.4547646 999 0 |
| 574 69 2 7.2 7.2916667 1.4696508 .03708329 1 |
+--------------------------------------------------------------------------+
Individual 51 is 86 years old at diagnosis and male. They have 7 rows of data with each row corresponding
to a different time interval. The start of the interval is given by tstart and the end of the interval by tstop.
Each time interval is 0.2 years, execpt the last interval in which they die (event==1) at 1.375 years.
For each interval there is an associated weight (wt), which is the inverse of the expected survival at the midpoint of the interval.
As the expected survival decreases over time, the weights increase over time. The meanhazard_wt gives the weighted mean
expected mortality rate at each individuals event time. Note that for any censored time it is set to 999. When fitting
relative survival models using stpm2 or strcs the expected mortality rate at the event time is needed, but is not required
for any censored times. However, having a missing value would exclude these rows from the analysis and so we feed it
a value that is actually not used when we fit the model. Individual 574 is younger than Individual 1 and so the weights are lower
at the same time points, e.g. 1.039 vs 1.222 at 1 year.
Having restructured the data we can now use stset where we need to give the end of each interval (tstop),
the start of the interval (tstart). The weights are passed using [iweights=wt].
. stset tstop [iweight=wt], enter(tstart) failure(event==1)
Survival-time data settings
Failure event: event==1
Observed time interval: (0, tstop]
Enter on or after: time tstart
Exit on or before: failure
Weight: [iweight=wt]
--------------------------------------------------------------------------
112,229 total observations
0 exclusions
--------------------------------------------------------------------------
112,229 observations remaining, representing
1,401 failures in single-record/single-failure data
21,994.417 total analysis time at risk and under observation
At risk from t = 0
Earliest observed entry t = 0
Last observed exit t = 10
The marginal model van now be fitted using stpm2.
As there are time-dependent weights cluster robust standard errors are used using vce(cluster id).
. stpm2, scale(hazard) df(5) bhazard(meanhazard_wt) vce(cluster id)
note: delayed entry models are being fitted
Iteration 0: log pseudolikelihood = -5670.0926
Iteration 1: log pseudolikelihood = -5558.0634
Iteration 2: log pseudolikelihood = -5557.3568
Iteration 3: log pseudolikelihood = -5557.3542
Iteration 4: log pseudolikelihood = -5557.3542
Log pseudolikelihood = -5557.3542 Number of obs = 112,229
(Std. err. adjusted for 4,744 clusters in id)
------------------------------------------------------------------------------
| Robust
| Coefficient std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
xb |
_rcs1 | .983882 .0430693 22.84 0.000 .8994676 1.068296
_rcs2 | .1594691 .0243798 6.54 0.000 .1116855 .2072526
_rcs3 | .0581832 .0161444 3.60 0.000 .0265407 .0898257
_rcs4 | .0197788 .0125402 1.58 0.115 -.0047995 .044357
_rcs5 | .0184667 .0097193 1.90 0.057 -.0005827 .0375162
_cons | -2.208731 .0511667 -43.17 0.000 -2.309016 -2.108446
------------------------------------------------------------------------------
. range tt 0 10 101
(112,128 missing values generated)
. predict s_mrsprep, surv timevar(tt) ci
After fitting the model, the marginal relative survival has been predicted. This can now be compared to the Pohar Perme non-parametric estimate.
. merge 1:1 _n using PP, nogenerate
Result Number of obs
-----------------------------------------
Not matched 107,485
from master 107,485
from using 0
Matched 4,744
-----------------------------------------
. twoway (rarea R_pp_lci R_pp_uci PP_t, sort connect(stairstep) color(${C1ci})) ///
> (line R_pp PP_t, sort connect(stairstep) color(${C1})) ///
> (line s_mrsprep* tt, color(${C2} ${C2} ${C2}) lpattern(solid dash dash)) ///
> , legend(order(2 "Pohar Perme" 3 "Marginal stpm2 model") ///
> ring(0) cols(1) pos(7)) ///
> ylabel(0.6(0.1)1, format(%3.1f)) ///
> ytitle("Marginal relative survival") ///
> xtitle("Years from diagnosis") ///
> name(int_stand_standsurv, replace)
There is now good agreement between the model based and the non parametric Pohar Perme estimtor. Here the marginal estimate obtained from regression standardization and them arginal model are very similar.
The estimate here is an internally standardized estimate, over the observed covariate distribution and thus would not be comparable to another study with a different age/sex distribution or if separate analysis were performed for males and females. See the example of external age standardization and modelling covarites for further extensions.
Note that mrsprep makes uses of frames.
. frame
(current frame is mrs_data)
. frames dir
* default 4744 x 74; Skin melanoma, diagnosed 1975-94, follow-up to 1995
* mrs_data 112229 x 37
Note: Frames marked with * contain unsaved data.
It is possible to switch to the orginal data using frame change default.
References
Lambert PC, Syriopoulou E, Rutherford MJ. Direct modelling of age standardized marginal relative survival through incorporation of time-dependent weights. BMC Medical Research Methodology 2021;21:84
