Red Door Analytics
Resources · Tutorial

Fitting your first model in merlin

The same Weibull model fitted twice, in streg and in merlin, so you can see which parts of the output are identical, which are only presented differently, and why the two log-likelihoods are not the same number.

Tutorial3 min readStata

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

cancer.dta ships with Stata: 48 patients in a drug trial, of whom 31 died. drug is coded 1 for placebo and 2 or 3 for the two active arms, so it needs two indicators.

Stata
sysuse cancer, clear
gen byte drug2 = drug==2
gen byte drug3 = 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.

Stata
. stset studytime, failure(died)

Survival-time data settings

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

--------------------------------------------------------------------------
         48  total observations
          0  exclusions
--------------------------------------------------------------------------
         48  observations remaining, representing
         31  failures in single-record/single-failure data
        744  total analysis time at risk and under observation
                                                At risk from t =         0
                                     Earliest observed entry t =         0
                                          Last observed exit t =        39

. streg age drug2 drug3, dist(weibull) nohr

        Failure _d: died
  Analysis time _t: studytime

Fitting constant-only model:
Iteration 0:  Log likelihood = -61.342985
Iteration 1:  Log likelihood =   -60.6269
Iteration 2:  Log likelihood = -60.624022
Iteration 3:  Log likelihood = -60.624022

Fitting full model:
Iteration 0:  Log likelihood = -60.624022  
Iteration 1:  Log likelihood = -53.956776  
Iteration 2:  Log likelihood = -42.173112  
Iteration 3:  Log likelihood = -42.090795  
Iteration 4:  Log likelihood = -42.090672  
Iteration 5:  Log likelihood = -42.090672  

Weibull PH regression

No. of subjects =  48                                   Number of obs =     48
No. of failures =  31
Time at risk    = 744
                                                        LR chi2(3)    =  37.07
Log likelihood = -42.090672                             Prob > chi2   = 0.0000

------------------------------------------------------------------------------
          _t | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
         age |   .1172843   .0364704     3.22   0.001     .0458037     .188765
       drug2 |  -1.768649   .4874724    -3.63   0.000    -2.724077   -.8132205
       drug3 |  -2.547726   .5144271    -4.95   0.000    -3.555985   -1.539468
       _cons |  -10.58208   2.299523    -4.60   0.000    -15.08907   -6.075102
-------------+----------------------------------------------------------------
       /ln_p |   .5573333   .1402154     3.97   0.000     .2825163    .8321504
-------------+----------------------------------------------------------------
           p |    1.74601   .2448175                      1.326463    2.298256
         1/p |   .5727343   .0803062                      .4351126    .7538844
------------------------------------------------------------------------------

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 =       -744  
Iteration 1:  Log likelihood = -129.09057  
Iteration 2:  Log likelihood =  -114.3412  
Iteration 3:  Log likelihood = -110.49148  
Iteration 4:  Log likelihood = -110.28014  
Iteration 5:  Log likelihood = -110.26737  
Iteration 6:  Log likelihood = -110.26736  
Iteration 7:  Log likelihood = -110.26736  

Fixed effects regression model                              Number of obs = 48
Log likelihood = -110.26736
------------------------------------------------------------------------------
             | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
studytime:   |            
         age |   .1172843   .0364704     3.22   0.001     .0458037     .188765
       drug2 |  -1.768649   .4874724    -3.63   0.000    -2.724077   -.8132205
       drug3 |  -2.547726   .5144271    -4.95   0.000    -3.555985   -1.539468
       _cons |  -10.58208   2.299523    -4.60   0.000    -15.08907   -6.075102
  log(gamma) |   .5573333   .1402154     3.97   0.000     .2825163    .8321504
------------------------------------------------------------------------------

Reading the output

Every coefficient, standard error, z statistic and confidence limit matches streg's to all six printed digits. That is the reassurance worth having: for a model both commands can fit, merlin is not an approximation to streg — it is the same likelihood maximised the same way.

Two things are only presented differently. streg reports the shape parameter as /ln_p and again as p; merlin reports that same quantity as log(gamma), with the identical estimate 0.5573333 and standard error 0.1402154. And merlin labels the equation with the response variable, where streg labels it _t.

The log-likelihoods, though, genuinely differ: -42.090672 from streg and -110.26736 from merlin. The gap is 68.176683, which is exactly the sum of log(t) over the 31 failure times. 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.

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.

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.

Need this applied to your own data?

Tell us what you're modelling and we'll point you to the closest worked example — or build one with you.

Get in touch Start a project