Fitting your first model in merlin

merlin feels unfamiliar before it feels useful. The quickest way past that is to fit something you have already fitted a hundred times and watch it come back with the numbers you expect — so this walks through a Weibull proportional hazards model in streg first, then the same model in merlin, and compares them line by line.

The data

Rather than borrow a dataset, we simulate one, so that we know the answers both commands should find. It is a small three-arm trial: 90 patients aged 45 to 70, 30 on placebo (drug 1) and 30 on each of two active drugs (2 and 3). Because drug has three levels, the models need two indicators for it, drug2 and drug3.

Stata

. clear

. set seed 4817

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

. gen byte drug = mod(_n, 3) + 1    // 1 placebo, 2 and 3 active, 30 per arm

. gen age = runiformint(45, 70)     // age at entry, in whole years

. gen byte drug2 = drug==2

. gen byte drug3 = drug==3

Survival times, in months, come from a Weibull proportional hazards model,

h (t) = λ γ t^{γ - 1} \exp (β_{1} age + β_{2} drug2 + β_{3} drug3)

with scale $λ = 0.001$ , shape $γ = 1.5$ , and log hazard ratios $β_{1} = 0.05$ per year of age, $β_{2} = - 1$ for drug 2 and $β_{3} = - 1.5$ for drug 3. avalon, the merlin family's simulator, draws the times with its standard subcommand, and anyone still alive at 36 months is censored there.

Stata

. avalon standard studytime died,                  /// new variables
>         distribution(weibull)                    /// Weibull baseline hazard
>         lambda(0.001) gamma(1.5)                 /// scale and shape
>         covariates(age 0.05 drug2 -1 drug3 -1.5) /// log hazard ratios
>         maxtime(36)                              //  censored at 36 months

Both commands estimate the constant as $\log λ$ and the shape as $\log γ$ , so the true values they are aiming at are $\log 0.001 = - 6.91$ and $\log 1.5 = 0.405$ , alongside 0.05, −1 and −1.5.

Stata

. tab drug died

           |         died
      drug |         0          1 |     Total
-----------+----------------------+----------
         1 |         1         29 |        30
         2 |        10         20 |        30
         3 |        15         15 |        30
-----------+----------------------+----------
     Total |        26         64 |        90

64 of the 90 died within the 36 months: 29 of the 30 on placebo, 20 on drug 2 and 15 on drug 3.

The model you already know

streg works from the data being stset, so declare the survival setup first. Asking for coefficients rather than hazard ratios keeps the two outputs directly comparable. The last line keeps streg's log-likelihood, which we come back to below.

Stata

. stset studytime, failure(died)

Survival-time data settings

         Failure event: died!=0 & died<.
Observed time interval: (0, studytime]
     Exit on or before: failure

--------------------------------------------------------------------------
         90  total observations
          0  exclusions
--------------------------------------------------------------------------
         90  observations remaining, representing
         64  failures in single-record/single-failure data
  1,835.212  total analysis time at risk and under observation
                                                At risk from t =         0
                                     Earliest observed entry t =         0
                                          Last observed exit t =        36

. streg age drug2 drug3, dist(weibull) nohr

        Failure _d: died
  Analysis time _t: studytime

Fitting constant-only model:
Iteration 0:  Log likelihood =  -129.4186
Iteration 1:  Log likelihood = -129.13548
Iteration 2:  Log likelihood = -129.13517
Iteration 3:  Log likelihood = -129.13517

Fitting full model:
Iteration 0:  Log likelihood = -129.13517
Iteration 1:  Log likelihood = -119.16215
Iteration 2:  Log likelihood = -111.92416
Iteration 3:  Log likelihood = -111.89654
Iteration 4:  Log likelihood = -111.89652
Iteration 5:  Log likelihood = -111.89652

Weibull PH regression

No. of subjects =         90                            Number of obs =     90
No. of failures =         64
Time at risk    = 1,835.2116
                                                        LR chi2(3)    =  34.48
Log likelihood = -111.89652                             Prob > chi2   = 0.0000

------------------------------------------------------------------------------
          _t | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
         age |   .0320407   .0175389     1.83   0.068    -.0023348    .0664162
       drug2 |  -1.112716   .3011575    -3.69   0.000    -1.702974    -.522458
       drug3 |   -1.74678   .3336523    -5.24   0.000    -2.400726   -1.092833
       _cons |  -5.258999   1.186898    -4.43   0.000    -7.585276   -2.932722
-------------+----------------------------------------------------------------
       /ln_p |   .2821586   .1071567     2.63   0.008     .0721353    .4921819
-------------+----------------------------------------------------------------
           p |   1.325989   .1420886                      1.074801    1.635882
         1/p |   .7541541   .0808127                      .6112912     .930405
------------------------------------------------------------------------------

. scalar ll_streg = e(ll)

The same model in merlin

merlin does not use the stset declaration. The response variable and its failure indicator are named inside the model statement instead, which is why this runs on studytime and died rather than on the _t and _d that stset created — the command would work just as well in a dataset that had never been stset.

Everything about one model sits inside one set of brackets: the response, then its covariates, then a family. For a single model that looks like ceremony. It is the reason you can write a second set of brackets later.

Stata

. merlin (studytime age drug2 drug3, family(weibull, failure(died)))

Fitting full model:

Iteration 0:  Log likelihood = -1835.2116
Iteration 1:  Log likelihood = -288.78602
Iteration 2:  Log likelihood = -262.89682
Iteration 3:  Log likelihood =  -261.3348
Iteration 4:  Log likelihood = -261.26401
Iteration 5:  Log likelihood = -261.26397
Iteration 6:  Log likelihood = -261.26397

Fixed effects regression model                              Number of obs = 90
Log likelihood = -261.26397
------------------------------------------------------------------------------
             | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
studytime:   |
         age |   .0320407   .0175389     1.83   0.068    -.0023348    .0664162
       drug2 |  -1.112716   .3011575    -3.69   0.000    -1.702974    -.522458
       drug3 |   -1.74678   .3336523    -5.24   0.000    -2.400726   -1.092833
       _cons |  -5.258999   1.186898    -4.43   0.000    -7.585276   -2.932722
  log(gamma) |   .2821586   .1071567     2.63   0.008     .0721353    .4921819
------------------------------------------------------------------------------

. scalar ll_merlin = e(ll)

Reading the output

Every coefficient, standard error, z statistic and confidence limit matches streg's to every printed digit. That is the reassurance worth having: for a model both commands can fit, merlin is not an approximation to streg — it maximises the same likelihood, up to a constant we come to below, and arrives at the same estimates.

Both are also close to the truth. The log hazard ratios are 0.032 per year of age (we simulated 0.05), −1.11 for drug 2 (−1) and −1.75 for drug 3 (−1.5); the constant is −5.26 (−6.91) and the log shape 0.282 (0.405). With 90 patients, none of the five estimates is more than one and a half standard errors from its true value.

Two things are only presented differently. streg reports the log of the shape parameter as /ln_p, and the shape itself as p (and 1/p); merlin reports that log shape as log(gamma), with the identical estimate 0.2821586 and standard error 0.1071567. And merlin labels the equation with the response variable, where streg labels it _t.

The log-likelihoods, though, genuinely differ: −111.89652 from streg and −261.26397 from merlin. The last block of the script measures the gap, and the sum of log(t) over the 64 failure times:

Stata

. gen double logt = log(studytime)

. quietly summarize logt if died

. display "sum of log(t) over the failures   = " %11.6f r(sum)
sum of log(t) over the failures   =  149.367442

. display "streg minus merlin log likelihood = " %11.6f ll_streg - ll_merlin
streg minus merlin log likelihood =  149.367442

. quietly summarize studytime

. display "total follow-up time, sum of t    = " %11.4f r(sum)
total follow-up time, sum of t    =   1835.2116

They are the same number, 149.367442. The two commands evaluate the same likelihood on different scales — merlin on the time scale, streg on the log-time scale — and changing the variable from t to log(t) contributes precisely that sum. The density of log(t) is the density of t multiplied by t, so each failure adds its own log(t); a censored patient contributes a survival probability rather than a density, and adds nothing.

ℓ_{streg} = ℓ_{merlin} + \sum_{i : d_{i} = 1} \log t_{i}

It depends on the data and not on any parameter, so it shifts the log-likelihood without moving a single estimate, and it cancels in any likelihood-ratio test between nested models. Compare log-likelihoods across the two commands and you will think something is wrong; compare differences of log-likelihoods and they agree.

The last line of that block explains the first line of merlin's iteration log. merlin starts every parameter at zero. The constant and the coefficients at zero give $λ \exp (X β) = 1$ , and log(gamma) at zero gives $γ = 1$ , so the hazard is 1 at every time and the cumulative hazard is just t. Each patient contributes $d_{i} \log 1 - t_{i}$ , and the log-likelihood at iteration 0 is minus the total follow-up time,

ℓ (0) = \sum_{i} {d_{i} \log h (t_{i}) - H (t_{i})} = - \sum_{i} t_{i}

which here is −1835.2116: the 1,835.2116 months that stset and streg report as time at risk. It does not depend on who died, only on how long everyone was followed.

What this buys you

On its own, nothing: this is a model streg already fits, in more keystrokes. What it buys is trust in the syntax — a response, its covariates, and a family, inside brackets — and the knowledge that the numbers are the ones you would have got anyway. The next step is the second set of brackets, because that is where merlin starts doing what streg cannot.

MethodSurvival analysis

Fitting your first model in merlin

The data

The model you already know

The same model in merlin

Reading the output

What this buys you

Need this applied to your own data?