Semi-parametric multi-state modelling
Multi-state models with a Cox model for every transition: fitted with merlin, and turned into exact transition probabilities, length of stay and contrasts with pendragon.
A multi-state model is usually fitted one transition at a time, and each transition model can be any survival model we like. In this tutorial every one of them is a Cox model, which makes the multi-state model semi-parametric: no transition’s baseline hazard is given a shape. The work is split between two commands:
merlinfits each transition’s Cox model, withfamily(cox), on that transition’s rows of the data. The baseline cumulative hazard is Breslow’s, a step function that jumps at the transition’s event times;pendragoncombines the fitted transitions into transition probabilities, the Aalen–Johansen estimator over the event times of all the transitions. It is exact: there is no simulation and no approximation on a time grid;- from the same estimate,
pendragongives length of stay, and differences and ratios between covariate patterns.
Let’s take a look, going from data in wide format to fitted models, transition probabilities, a stacked plot, length of stay and contrasts.
Setting up the data
We’ll simulate a dataset shaped like the Rotterdam breast cancer data, a standard example for multi-state models, so that we know the answers the models should find. We start with 1,500 patients and their baseline covariates: age at surgery, tumour size in three classes, the number of positive lymph nodes, nodes, the log of the progesterone receptor measurement plus one, pr_1, and whether they received hormonal therapy, hormon.
. clear
. set seed 4430
. set obs 1500
Number of observations (_N) was 0, now 1,500.
. gen pid = _n
. gen age = min(max(round(rnormal(55,10)), 25), 90) // age at surgery
. gen size = irecode(runiform(), 0.45, 0.9) + 1 // tumour size class
. label define size 1 "<=20 mm" 2 ">20-50 mm" 3 ">50 mm"
. label values size size
. label variable size "Tumour size"
. gen nodes = rpoisson(3) // positive nodes
. gen pr_1 = runiform(0,8) // log(PR + 1)
. gen hormon = runiform() < 0.2 // hormonal therapyEach patient then follows an illness-death process. From surgery (state 1) they can relapse (state 2) or die without relapse (state 3), and after a relapse they can go on to die. Each of the three transitions has a Weibull hazard with proportional covariate effects,
with t the time in years since surgery for every transition, a clock-forward (Markov) model; with , the hazard of death after relapse is in fact constant. The true values are:
- Transition 1, relapse: and , with log hazard ratios of −0.01 per year of age, 0.4 and 0.7 for the two larger size classes, 0.08 per positive node, −0.05 for
pr_1and −0.2 for hormonal therapy. - Transition 2, death without relapse: and , with 0.12 per year of age, 0.2 and 0.4 for size and 0.04 per node, and no effect of
pr_1or hormonal therapy. ( is tiny because age enters the linear predictor uncentred.) - Transition 3, death after relapse: and , with 0.01 per year of age, 0.15 and 0.3 for size, 0.03 per node, −0.1 for
pr_1and 0.1 for hormonal therapy.
Everyone is followed for 20 years. avalon msm simulates the whole process from the transition-specific hazards and the transition matrix,
. gen byte size2 = size==2 // indicators, for the simulation
. gen byte size3 = size==3
. matrix tmat = (.,1,2\.,.,3\.,.,.) // illness-death
. avalon msm time state event, transmatrix(tmat) maxtime(20) ///
> hazard1(dist(weibull) lambda(0.1) gamma(0.9) /// relapse
> covariates(age -0.01 size2 0.4 size3 0.7 ///
> nodes 0.08 pr_1 -0.05 hormon -0.2)) ///
> hazard2(dist(weibull) lambda(0.00001) gamma(1.2) /// death
> covariates(age 0.12 size2 0.2 size3 0.4 ///
> nodes 0.04)) ///
> hazard3(dist(weibull) lambda(0.2) gamma(1) /// death
> covariates(age 0.01 size2 0.15 size3 0.3 /// after
> nodes 0.03 pr_1 -0.1 hormon 0.1)) // relapseand returns the time at which each state was entered, and which state it was. The Rotterdam data store this information in wide format, with times in months, in four variables, so we build four to match:
. gen rfi = state1==2 // relapsed
. gen rf = 12*time1 // months to relapse, death or censori
> ng
. gen osi = state1==3 | state2==3 // died
. gen os = 12*cond(rfi, time2, time1) // months to death or censoring
. drop time* state* event* size2 size3which gives us our dataset in wide format,
. list pid rf rfi os osi if inlist(pid,1,2,4), sepby(pid) noobs
+---------------------------------------+
| pid rf rfi os osi |
|---------------------------------------|
| 1 21.32643 1 127.8385 1 |
|---------------------------------------|
| 2 49.45266 0 49.45266 1 |
|---------------------------------------|
| 4 240 0 240 0 |
+---------------------------------------+where pid is our patient identifier, rf contains the time of relapse (or of censoring or death, if there was no relapse), rfi the relapse event indicator, os overall survival time and osi our survival indicator. Patient 1 relapsed at 21.3 months and died at 127.8 months, patient 2 died without a relapse at 49.5 months, and patient 4 was still alive and relapse-free when follow-up ended at 240 months. We can use pendragon set to reshape our wide dataset into the stacked format, with a row for each transition a patient is at risk of.
. pendragon set, id(pid) states(rfi osi) times(rf os)
transition-long data: 3889 rows, 3 transitions, 2007 events
Transition frequencies (_status==1):
to: to: to:
start rfi osi
from:start 0 889 315
from:rfi 0 0 803
from:osi 0 0 0pendragon set replaces the data in memory with the stacked data, counts the transitions that happened (889 relapses, 315 deaths without relapse and 803 deaths after relapse), and creates variables for use in subsequent analyses, similar to stset,
. list pid _start _stop _from _to _status _trans if inlist(pid,1,2,4), ///
> sepby(pid) noobs
+--------------------------------------------------------------+
| pid _start _stop _from _to _status _trans |
|--------------------------------------------------------------|
| 1 0 21.326426 1 2 1 1 -> 2 |
| 1 0 21.326426 1 3 0 1 -> 3 |
| 1 21.326426 127.83852 2 3 1 2 -> 3 |
|--------------------------------------------------------------|
| 2 0 49.452663 1 2 0 1 -> 2 |
| 2 0 49.452663 1 3 1 1 -> 3 |
|--------------------------------------------------------------|
| 4 0 240 1 2 0 1 -> 2 |
| 4 0 240 1 3 0 1 -> 3 |
+--------------------------------------------------------------+where _trans is labelled with the move it stands for. pendragon set also returns a default transition matrix, if one was not provided. We need to store this for later use. In this case it’s the illness-death transition matrix, as it will assume an upper triangular transition matrix, with a common initial state.
. mat tmat = r(transmatrix)
. mat list tmat
tmat[3,3]
to: to: to:
start rfi osi
from:start . 1 2
from:rfi . . 3
from:osi . . .We can now stset our dataset, using the variables pendragon set created,
. stset _stop, enter(_start) failure(_status==1) scale(12)
Survival-time data settings
Failure event: _status==1
Observed time interval: (0, _stop]
Enter on or after: time _start
Exit on or before: failure
Time for analysis: time/12
--------------------------------------------------------------------------
3,889 total observations
0 exclusions
--------------------------------------------------------------------------
3,889 observations remaining, representing
2,007 failures in single-record/single-failure data
31,893.215 total analysis time at risk and under observation
At risk from t = 0
Earliest observed entry t = 0
Last observed exit t = 20We also changed the timescale from months into years by using scale(12). Tumour size at diagnosis, size, is a three-level factor variable, for which we now create dummy indicator variables,
. tab size, gen(sz)
Tumour size | Freq. Percent Cum.
------------+-----------------------------------
<=20 mm | 1,696 43.61 43.61
>20-50 mm | 1,772 45.56 89.17
>50 mm | 421 10.83 100.00
------------+-----------------------------------
Total | 3,889 100.00Fitting the transition models
Neither merlin nor pendragon takes factor variables, so you must create your own dummies. We can now fit our multi-state model, in this case transition-specific Cox models, using merlin’s family(cox). For transition 1,
. merlin (_t age sz2 sz3 nodes pr_1 hormon if _trans==1, ///
> family(cox, failure(_d) ltruncated(_t0)))
Fitting full model:
Iteration 0: Log likelihood = -5975.4624
Iteration 1: Log likelihood = -5925.8487
Iteration 2: Log likelihood = -5925.6435
Iteration 3: Log likelihood = -5925.6435
Fixed effects regression model Number of obs = 1,500
Log likelihood = -5925.6435
------------------------------------------------------------------------------
| Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
_t: |
age | -.0067375 .0036984 -1.82 0.068 -.0139863 .0005113
sz2 | .3912494 .0722216 5.42 0.000 .2496977 .5328011
sz3 | .6268919 .1104953 5.67 0.000 .4103251 .8434587
nodes | .0754635 .01858 4.06 0.000 .0390473 .1118797
pr_1 | -.0792186 .0147572 -5.37 0.000 -.1081422 -.0502949
hormon | -.215271 .0861233 -2.50 0.012 -.3840696 -.0464723
------------------------------------------------------------------------------
. estimates store m1merlin does not read the stset settings itself, so we hand it the variables stset made: _t, the time in years, as the response, _d, the event indicator, in failure(), and _t0, the time the patient became at risk of the transition, in ltruncated(). Everyone is at risk of the two transitions out of state 1 from surgery, but only at risk of death after relapse from the time they relapse, so transition 3 needs that delayed entry; the if inside the model picks out the transition’s rows. We can reassure ourselves that merlin’s Cox model agrees with official Stata’s stcox,
. stcox age sz2 sz3 nodes pr_1 hormon if _trans==1, breslow
Failure _d: _status==1
Analysis time _t: _stop/12
Enter on or after: time _start
Iteration 0: Log likelihood = -5975.4624
Iteration 1: Log likelihood = -5925.8487
Iteration 2: Log likelihood = -5925.6435
Iteration 3: Log likelihood = -5925.6435
Refining estimates:
Iteration 0: Log likelihood = -5925.6435
Cox regression with no ties
No. of subjects = 1,500 Number of obs = 1,500
No. of failures = 889
Time at risk = 14,493.0121
LR chi2(6) = 99.64
Log likelihood = -5925.6435 Prob > chi2 = 0.0000
------------------------------------------------------------------------------
_t | Haz. ratio Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
age | .9932851 .0036736 -1.82 0.068 .9861111 1.000511
sz2 | 1.478827 .1068033 5.42 0.000 1.283637 1.703698
sz3 | 1.871783 .2068232 5.67 0.000 1.507307 2.324392
nodes | 1.078384 .0200364 4.06 0.000 1.03982 1.118378
pr_1 | .923838 .0136333 -5.37 0.000 .8974999 .950949
hormon | .8063229 .0694432 -2.50 0.012 .681084 .954591
------------------------------------------------------------------------------which indeed it does: the same log likelihood, −5925.6435, and hazard ratios that are the exponentiated merlin coefficients. Phew. For transition 2,
. merlin (_t age sz2 sz3 nodes pr_1 hormon if _trans==2, ///
> family(cox, failure(_d) ltruncated(_t0)))
Fitting full model:
Iteration 0: Log likelihood = -2106.833
Iteration 1: Log likelihood = -1932.6525
Iteration 2: Log likelihood = -1931.8796
Iteration 3: Log likelihood = -1931.8795
Fixed effects regression model Number of obs = 1,500
Log likelihood = -1931.8795
------------------------------------------------------------------------------
| Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
_t: |
age | .118193 .0064388 18.36 0.000 .1055731 .1308128
sz2 | .2543679 .1200156 2.12 0.034 .0191417 .4895941
sz3 | .3989257 .2005832 1.99 0.047 .0057899 .7920616
nodes | .0338848 .0314819 1.08 0.282 -.0278186 .0955882
pr_1 | -.0062765 .0249233 -0.25 0.801 -.0551254 .0425723
hormon | .1328961 .1342732 0.99 0.322 -.1302746 .3960668
------------------------------------------------------------------------------
. estimates store m2and for transition 3,
. merlin (_t age sz2 sz3 nodes pr_1 hormon if _trans==3, ///
> family(cox, failure(_d) ltruncated(_t0)))
Fitting full model:
Iteration 0: Log likelihood = -4035.8988
Iteration 1: Log likelihood = -4005.1951
Iteration 2: Log likelihood = -4005.1765
Iteration 3: Log likelihood = -4005.1765
Fixed effects regression model Number of obs = 889
Log likelihood = -4005.1765
------------------------------------------------------------------------------
| Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
_t: |
age | .0075543 .0039804 1.90 0.058 -.0002472 .0153558
sz2 | .1466675 .0765161 1.92 0.055 -.0033014 .2966363
sz3 | .3255602 .1164242 2.80 0.005 .097373 .5537474
nodes | .0355739 .0201757 1.76 0.078 -.0039698 .0751175
pr_1 | -.1041699 .0153823 -6.77 0.000 -.1343186 -.0740212
hormon | .1800177 .0907928 1.98 0.047 .002067 .3579684
------------------------------------------------------------------------------
. estimates store m3Each time we store the model results using estimates store. We can then pass the model objects to pendragon to obtain a whole range of predictions. It’s that simple.
Because we simulated the data, we can check the estimates against the truth. For relapse, the log hazard ratios are −0.0067 for age (we simulated −0.01), 0.391 and 0.627 for the two larger size classes (0.4 and 0.7), 0.075 for nodes (0.08), −0.079 for pr_1 (−0.05) and −0.215 for hormonal therapy (−0.2). For death without relapse, the age effect is 0.118 (0.12), and the estimates for pr_1, −0.006, and hormonal therapy, 0.133, are compatible with their true value of zero. For death after relapse, the pr_1 effect is −0.104 (−0.1) and the age effect 0.0076 (0.01). All eighteen are within two standard errors of the truth; the furthest is the pr_1 effect on relapse, −0.079 against −0.05, very nearly two standard errors away.
Transition probabilities
First we create a time variable, at which to calculate predictions,
. range tvar 0 15 100
(3,789 missing values generated)pendragon always writes transition probabilities, and from(1) says we start in state 1 at time 0 (the default, made explicit). I’m predicting for a patient aged 50 through the at() option, and zeros sets every covariate not named in at() to 0. Without zeros, an unlisted covariate keeps each row’s own value, which would leave each transition with many linear predictors rather than one, so on Cox transitions pendragon stops with an error rather than guess. So, with zeros,
. pendragon , transmatrix(tmat) models(m1 m2 m3) timevar(tvar) ///
> at(age 50) zeros from(1)
pendragon: 3 Cox transitions, 3 states (Aalen-Johansen over the Breslow event ti
> mes, exact)
wrote prob_1_1 prob_1_2 prob_1_3which gives us some new variables, prob_1_1 to prob_1_3, the probability of being in each state given that we started in state 1,
. list prob_1_* tvar if inlist(_n,34,67,100), noobs
+------------------------------------------+
| prob_1_1 prob_1_2 prob_1_3 tvar |
|------------------------------------------|
| .73147985 .10645278 .16206738 5 |
| .53608422 .10823676 .35567902 10 |
| .38302865 .09293704 .52403432 15 |
+------------------------------------------+With Cox transition models these are exact. Each transition’s cumulative hazard is Breslow’s, a step function that jumps at that transition’s event times, and the transition probabilities are the product-limit (Aalen–Johansen) estimator over the event times of all three transitions, the estimate R’s mstate forms with msfit and probtrans. There is no simulation, so no Monte Carlo error, and no approximation on the time grid: the 100 values of tvar only say where the estimate is written. After a transition’s last event its cumulative hazard stays flat, and beyond the latest follow-up time in the data there is no risk set to estimate from, so pendragon writes missing values there rather than extrapolate; it returns that time in r(horizon). Our follow-up runs to 20 years, so predicting to 15 is safe. Our patient aged 50 has a probability of 0.731 of being alive and relapse-free at 5 years, 0.536 at 10 years and 0.383 at 15 years. As we know the true model, we can also calculate what that probability really is: it’s the chance of avoiding both transitions out of state 1, , where is the cumulative hazard of transition k,
. // true P(alive, relapse-free) at 5, 10, 15 years: age 50, others 0
. forvalues t = 5(5)15 {
2. display %6.4f exp(-0.1*`t'^0.9*exp(-0.01*50) ///
> - 0.00001*`t'^1.2*exp(0.12*50))
3. }
0.7513
0.5794
0.45020.751, 0.579 and 0.450. The estimates sit below the truth, most at 15 years, 0.383 against 0.450. That is this particular dataset rather than the method: our patient has no positive nodes and a pr_1 of 0, at the edge of the data, so the prediction leans on the estimated covariate effects, and the pr_1 effect on relapse came out steeper than the truth, which raises the predicted relapse hazard for a patient with pr_1 of 0.
We can use pendragon graph occupancy, with stack, to obtain a stacked plot of the transition probabilities. Its twopts() option passes twoway options on to the graph, so we can put the legend under the plot and name the states,
. pendragon graph occupancy, timevar(tvar) from(1) stack ///
> twopts(title("Aged 50 at surgery") ///
> xtitle("Years since surgery") ///
> ytitle("Probability") ///
> legend(pos(6) order(3 "Dead" 2 "Relapsed" ///
> 1 "Alive, relapse-free")))
Length of stay
As well as transition probabilities, we can calculate the mean length of time spent in each state, as a function of follow-up time, by adding the los option,
. pendragon , transmatrix(tmat) models(m1 m2 m3) timevar(tvar) ///
> at(age 50) zeros from(1) los
pendragon: 3 Cox transitions, 3 states (Aalen-Johansen over the Breslow event ti
> mes, exact)
wrote prob_1_1 los_1_1 prob_1_2 los_1_2 prob_1_3 los_1_3. list los_1_* tvar if inlist(_n,34,67,100), noobs
+------------------------------------------+
| los_1_1 los_1_2 los_1_3 tvar |
|------------------------------------------|
| 4.2714425 .41269799 .31585954 5 |
| 7.3913748 .98273687 1.6258883 10 |
| 9.6855851 1.4589964 3.8554185 15 |
+------------------------------------------+So over the first 15 years after surgery, a patient aged 50 spends on average 9.69 years alive and relapse-free and 1.46 years alive after a relapse, and the remaining 3.86 years are spent dead (the truth is 10.19, 1.25 and 3.56). Length of stay is the area under each transition probability curve, and here it is exact too: the curves are step functions, so each area is a sum of rectangles.
Contrasts between covariate patterns
We can get useful contrasts between covariate patterns through use of the at1() and at2() options (you can use more at#()s if you wish). Here, we’ll calculate both the difference and the ratio, for the transition probabilities and for length of stay, for a patient aged 60 compared to a patient aged 50. at1() is the default reference group (you can change that with atreference()), so each contrast is at2() minus, or over, at1(),
. pendragon , transmatrix(tmat) models(m1 m2 m3) timevar(tvar) ///
> at1(age 50) at2(age 60) zeros from(1) los ///
> difference ratio
pendragon: 3 Cox transitions, 3 states (Aalen-Johansen over the Breslow event ti
> mes, exact)
contrast: at2 vs the reference at1
wrote diff_prob_at2_1_1 diff_prob_at2_1_2 diff_prob_at2_1_3 ratio_prob_at2_1
> _1 ratio_prob_at2_1_2 ratio_prob_at2_1_3 diff_los_at2_1_1 diff_los_at2_1_2 dif
> f_los_at2_1_3 ratio_los_at2_1_1 ratio_los_at2_1_2 ratio_los_at2_1_3The variable names say what was contrasted: diff_prob_at2_1_3 is the difference in the probability of being in state 3, having started in state 1, between at2() and the reference. The differences first,
. list diff_prob_at2_1_* tvar if inlist(_n,34,67,100), noobs ab(17)
+------------------------------------------------------------------+
| diff_prob_at2_1_1 diff_prob_at2_1_2 diff_prob_at2_1_3 tvar |
|------------------------------------------------------------------|
| -.03201961 -.01436989 .04638951 5 |
| -.05490025 -.02056624 .07546649 10 |
| -.05867058 -.02225497 .08092555 15 |
+------------------------------------------------------------------+Being ten years older at surgery lowers the probability of being alive and relapse-free, by 0.032 at 5 years and 0.059 at 15, and raises the probability of being dead, by 0.046 at 5 years and 0.081 at 15. The age effect on death without relapse (0.118 per year, on the log hazard scale) far outweighs the slightly lower relapse hazard at older ages (−0.0067 per year). Because we know the true model, we can check the first of these,
. // true difference in P(alive, relapse-free), age 60 minus age 50, others 0
. forvalues t = 5(5)15 {
2. display %7.4f exp(-0.1*`t'^0.9*exp(-0.01*60) ///
> - 0.00001*`t'^1.2*exp(0.12*60)) ///
> - exp(-0.1*`t'^0.9*exp(-0.01*50) ///
> - 0.00001*`t'^1.2*exp(0.12*50))
3. }
-0.0295
-0.0564
-0.0724−0.0295, −0.0564 and −0.0724. The estimates, −0.0320, −0.0549 and −0.0587, match the first two closely and fall short at 15 years, for the same reason as the probabilities themselves.
The same call wrote the ratios,
. list ratio_prob_at2_1_* tvar if inlist(_n,34,67,100), noobs ab(18)
+---------------------------------------------------------------------+
| ratio_prob_at2_1_1 ratio_prob_at2_1_2 ratio_prob_at2_1_3 tvar |
|---------------------------------------------------------------------|
| .95622626 .86501155 1.2862359 5 |
| .89759026 .80998839 1.2121758 10 |
| .84682456 .76053708 1.154428 15 |
+---------------------------------------------------------------------+so at 15 years, a patient aged 60 is 1.15 times as likely to have died as a patient aged 50 (truly 1.19), and 0.85 times as likely to be alive and relapse-free (truly 0.84). And it contrasted length of stay too,
. list diff_los_at2_1_* tvar if inlist(_n,34,67,100), noobs ab(16)
+---------------------------------------------------------------+
| diff_los_at2_1_1 diff_los_at2_1_2 diff_los_at2_1_3 tvar |
|---------------------------------------------------------------|
| -.06789399 -.0430417 .1109357 5 |
| -.29082811 -.1358013 .42662941 10 |
| -.57814356 -.23908156 .81722513 15 |
+---------------------------------------------------------------+Over 15 years, the patient aged 60 spends 0.58 years less alive and relapse-free and 0.24 years less alive after relapse than the patient aged 50 (truly 0.61 and 0.26 years less), so 0.82 years more of the 15 are spent dead (truly 0.87).
Confidence intervals are not yet available for Cox transition models: pendragon refuses the ci option on that path rather than ignore it. We’re working on this. Drop us an email at contact@reddooranalytics.se for any feedback, feature requests and bug reports.
MethodMulti-state models