Competing risks in merlin

Fitting your first model in merlin fitted one model and checked it against streg, which left the syntax looking like a lot of brackets for no return. This is where the return starts: a second set of brackets is a second model, fitted at the same time as the first.

The data

Rather than borrow a dataset, we simulate one, so that we know what the models should find: 400 women treated for cervical cancer, followed for up to five years to their first recurrence of disease. Each has a tumour diameter between 1 and 7 cm, a 40% chance that the tumour is hypoxic, and a 50% chance of chemoradiotherapy rather than radiotherapy alone.

Stata

. clear

. set seed 1906

. set obs 400
Number of observations (_N) was 0, now 400.

. gen tumsize = runiform(1, 7)            // tumour diameter, cm

. gen byte hypoxic = runiform() < 0.4     // hypoxic tumour

. gen byte chemo = runiform() < 0.5       // chemoradiotherapy

A recurrence is either pelvic or distant, and each has its own cause-specific hazard, a Weibull,

h_{k} (t) = λ_{k} γ_{k} t^{γ_{k} - 1} \exp (β_{k 1} tumsize + β_{k 2} hypoxic + β_{k 3} chemo), k = 1, 2

For pelvic recurrence ( $k = 1$ ), $λ_{1} = 0.02$ and $γ_{1} = 1.2$ , with log hazard ratios of 0.3 per cm of tumour, 0.5 for a hypoxic tumour and −0.5 for chemoradiotherapy. For distant recurrence ( $k = 2$ ), $λ_{2} = 0.01$ and $γ_{2} = 1.5$ , with 0.1, 0.8 and −0.2. merlin estimates each constant as $\log λ_{k}$ and each shape as $\log γ_{k}$ , so the true values there are −3.91 and 0.18 for pelvic recurrence and −4.61 and 0.41 for distant. avalon, the merlin family's simulator, does competing risks with its msm subcommand, and simulates from the two hazards together, one in hazard1() and one in hazard2():

Stata

. avalon msm dftime state event,                           /// new variables
>         hazard1(dist(weibull) lambda(0.02) gamma(1.2)    /// cause 1: pelvic
>             covariates(tumsize 0.3 hypoxic 0.5 chemo -0.5)) ///
>         hazard2(dist(weibull) lambda(0.01) gamma(1.5)    /// cause 2: distant
>             covariates(tumsize 0.1 hypoxic 0.8 chemo -0.2)) ///
>         maxtime(5)                                       //  5-year follow-up

avalon msm follows each woman to her first recurrence. It draws the time from the sum of the two hazards, then which kind of recurrence it was in proportion to their sizes at that moment. So each recurrence is classed as pelvic or distant, never both, and a woman whose first recurrence is distant is not followed on for a pelvic one. avalon msm numbers each new variable by the transition it records, so the first recurrence is in dftime1 and state1. state1 records where she ended up: 2 for pelvic, 3 for distant, and 1 for still free of recurrence at five years. We rename dftime1 to dftime, recode state1 as failtype and make an indicator for each cause.

Stata

. rename dftime1 dftime

. gen byte failtype = state1 - 1

. label define failtype 0 "none" 1 "pelvic" 2 "distant"

. label values failtype failtype

. gen byte pelvic  = failtype==1

. gen byte distant = failtype==2

Stata

. tab failtype

   failtype |      Freq.     Percent        Cum.
------------+-----------------------------------
       none |        196       49.00       49.00
     pelvic |        143       35.75       84.75
    distant |         61       15.25      100.00
------------+-----------------------------------
      Total |        400      100.00

So 143 pelvic recurrences, 61 distant, and 196 women censored without either. dftime is the time to whichever came first, or to censoring at five years.

One cause at a time

A cause-specific hazard model asks about one cause and treats the competing cause as censoring: a woman whose first recurrence was distant contributes follow-up time to the pelvic model, and then leaves it without a pelvic event. That is what failure(pelvic) does here, because pelvic is 0 for those women.

Stata

. merlin (dftime tumsize hypoxic chemo, family(weibull, failure(pelvic)))

Fitting full model:

Iteration 0:  Log likelihood = -1488.1383
Iteration 1:  Log likelihood = -466.80394
Iteration 2:  Log likelihood = -457.36779
Iteration 3:  Log likelihood = -444.18011
Iteration 4:  Log likelihood = -443.90034
Iteration 5:  Log likelihood = -443.89944
Iteration 6:  Log likelihood = -443.89944

Fixed effects regression model                             Number of obs = 400
Log likelihood = -443.89944
------------------------------------------------------------------------------
             | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
dftime:      |
     tumsize |   .3265336   .0521662     6.26   0.000     .2242897    .4287776
     hypoxic |   .7426681    .169327     4.39   0.000     .4107933    1.074543
       chemo |     -.4689   .1690237    -2.77   0.006    -.8001805   -.1376196
       _cons |  -4.326558   .3306674   -13.08   0.000    -4.974654   -3.678462
  log(gamma) |   .2924924   .0752393     3.89   0.000      .145026    .4399588
------------------------------------------------------------------------------

. scalar ll_pelvic = e(ll)

. matrix b_pelvic = e(b)

Distant recurrence is the same command with the other indicator, so now it is a pelvic recurrence that counts as censoring. The two lines after each fit keep its log-likelihood and its coefficients for a comparison further down.

Stata

. merlin (dftime tumsize hypoxic chemo, family(weibull, failure(distant)))

Fitting full model:

Iteration 0:  Log likelihood = -1488.1383
Iteration 1:  Log likelihood = -260.75084
Iteration 2:  Log likelihood = -246.47051
Iteration 3:  Log likelihood = -246.18135
Iteration 4:  Log likelihood = -246.18102
Iteration 5:  Log likelihood = -246.18102

Fixed effects regression model                             Number of obs = 400
Log likelihood = -246.18102
------------------------------------------------------------------------------
             | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
dftime:      |
     tumsize |   .0736668   .0747558     0.99   0.324    -.0728518    .2201854
     hypoxic |   .9062612   .2589698     3.50   0.000     .3986897    1.413833
       chemo |  -.4591817   .2584978    -1.78   0.076    -.9658282    .0474647
       _cons |  -4.168497   .4579717    -9.10   0.000    -5.066105   -3.270889
  log(gamma) |   .3106665     .11616     2.67   0.007      .082997     .538336
------------------------------------------------------------------------------

. scalar ll_distant = e(ll)

. matrix b_distant = e(b)

Both models find what we put in. For pelvic recurrence the log hazard ratios are 0.327 per cm, 0.743 for hypoxia and −0.469 for chemoradiotherapy (we simulated 0.3, 0.5 and −0.5); for distant recurrence they are 0.074, 0.906 and −0.459 (0.1, 0.8 and −0.2). The constants and log shapes are −4.33 and 0.292 for pelvic (−3.91 and 0.18) and −4.17 and 0.311 for distant (−4.61 and 0.41). The furthest from its true value is the pelvic log shape, 1.5 standard errors away, with the effect of hypoxia on pelvic recurrence next at 1.4; neither is remarkable among ten estimates.

Both causes in one command

Now the second set of brackets. Each pair carries its own response, its own covariates and its own family, and merlin fits them together.

Stata

. merlin (dftime tumsize hypoxic chemo, family(weibull, failure(pelvic)))  ///
>        (dftime tumsize hypoxic chemo, family(weibull, failure(distant)))

Fitting full model:

Iteration 0:  Log likelihood = -2976.2765
Iteration 1:  Log likelihood = -727.55478
Iteration 2:  Log likelihood = -705.26932
Iteration 3:  Log likelihood = -690.39384
Iteration 4:  Log likelihood = -690.08138
Iteration 5:  Log likelihood = -690.08047
Iteration 6:  Log likelihood = -690.08047

Fixed effects regression model                             Number of obs = 400
Log likelihood = -690.08047
------------------------------------------------------------------------------
             | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
dftime:      |
     tumsize |   .3265336   .0521662     6.26   0.000     .2242897    .4287776
     hypoxic |   .7426681    .169327     4.39   0.000     .4107933    1.074543
       chemo |     -.4689   .1690237    -2.77   0.006    -.8001805   -.1376196
       _cons |  -4.326558   .3306674   -13.08   0.000    -4.974654   -3.678462
  log(gamma) |   .2924924   .0752393     3.89   0.000      .145026    .4399588
-------------+----------------------------------------------------------------
dftime:      |
     tumsize |   .0736668   .0747558     0.99   0.324    -.0728518    .2201854
     hypoxic |   .9062612   .2589698     3.50   0.000     .3986897    1.413833
       chemo |  -.4591817   .2584978    -1.78   0.076    -.9658282    .0474647
       _cons |  -4.168497   .4579717    -9.10   0.000    -5.066105   -3.270889
  log(gamma) |   .3106665     .11616     2.67   0.007      .082997     .538336
------------------------------------------------------------------------------

. scalar ll_joint = e(ll)

. matrix b_joint = e(b)

The two blocks appear in the order you wrote them. Both are labelled dftime, because merlin labels an equation by its response variable and both causes are measured on the same clock — the failure indicator is what separates them.

The joint fit's iteration 0, −2976.2765, is twice the −1488.1383 that each single-cause fit started from; the last digit differs only because each is rounded separately. merlin starts every parameter at zero, where every hazard is 1, so each cause contributes minus the total follow-up time, whoever had which event, as in the first tutorial; the joint fit counts it once for each cause.

Is it really the same model?

It should be. With nothing shared between the two sets of brackets, the joint likelihood is the product of the two separate ones, so the log-likelihoods should add up.

Stata

. display "pelvic alone   = " %11.6f ll_pelvic
pelvic alone   = -443.899444

. display "distant alone  = " %11.6f ll_distant
distant alone  = -246.181024

. display "their sum      = " %11.6f ll_pelvic + ll_distant
their sum      = -690.080468

. display "both together  = " %11.6f ll_joint
both together  = -690.080468

They do, to every digit displayed. The coefficients should agree too. The joint fit stores the pelvic parameters first and the distant ones second, in the order of the brackets, so the two separate vectors side by side line up with it:

Stata

. matrix b_separate = b_pelvic, b_distant

. mata: d = st_matrix("b_joint") :- st_matrix("b_separate")

. mata: d'
                   1
     +----------------+
   1 |             0  |
   2 |             0  |
   3 |             0  |
   4 |             0  |
   5 |             0  |
   6 |   7.37455e-12  |
   7 |   3.32442e-11  |
   8 |  -3.15460e-11  |
   9 |  -1.28090e-10  |
  10 |   2.86743e-11  |
     +----------------+

. mata: max(abs(d))
  1.28090e-10

Across all ten parameters the largest disagreement between the joint fit and the two separate fits is 1.28e-10, which is within the optimiser's convergence tolerance rather than a difference between the models. The five pelvic parameters agree exactly. The distant ones differ because the fits stopped at different iterations: the joint fit judges convergence on all ten parameters at once and took six iterations, like the pelvic fit, while the distant fit on its own stopped after five. Fitting the two causes together has not changed either of them.

Then why write it this way?

Because nothing is shared yet. The moment the two causes share a parameter — the same coefficient constrained across both, or a random effect that carries a patient's frailty into each — the joint likelihood stops factorising and two separate commands can no longer produce it. The second set of brackets is what makes that possible to write down; this tutorial is the check that it costs nothing when you do not need it.

MethodCompeting risks

The data

One cause at a time

Both causes in one command

Is it really the same model?

Then why write it this way?

Need this applied to your own data?