Bug in models with >2 levels of nesting and no fixed effects

[ellessenne]

Hi,

stmixed seems to have some issues when trying to fit a multilevel model with three levels of nesting and no fixed effects:

. webuse jobhistory
(Job history data, Event History Analysis with Stata, Blossfeld et al. 2007)

. gen t = tend - tstart
. stset t, failure(failure == 1)

Survival-time data settings

         Failure event: failure==1
Observed time interval: (0, t]
     Exit on or before: failure

--------------------------------------------------------------------------
        600  total observations
          0  exclusions
--------------------------------------------------------------------------
        600  observations remaining, representing
        458  failures in single-record/single-failure data
     40,782  total analysis time at risk and under observation
                                                At risk from t =         0
                                     Earliest observed entry t =         0
                                          Last observed exit t =       428

. stmixed || birthyear: || id: , dist(weibull)
: invalid name
r(198);

Compared to mestreg:

. mestreg || birthyear: || id: , dist(weibull)

        Failure _d: failure==1
  Analysis time _t: t

Fitting fixed-effects model:

Iteration 0:  Log likelihood = -5549636.6  
Iteration 1:  Log likelihood = -2590.9107  
Iteration 2:  Log likelihood = -2507.8507  
Iteration 3:  Log likelihood = -2504.9718  
Iteration 4:  Log likelihood = -2504.8994  
Iteration 5:  Log likelihood = -2504.8994  

Refining starting values:

Grid node 0:  Log likelihood = -2531.2799

Fitting full model:

Iteration 0:  Log likelihood = -2531.2799  (not concave)
Iteration 1:  Log likelihood = -2503.6304  (not concave)
Iteration 2:  Log likelihood = -2499.6646  (not concave)
Iteration 3:  Log likelihood = -2497.9739  
Iteration 4:  Log likelihood = -2495.0353  
Iteration 5:  Log likelihood = -2491.8732  
Iteration 6:  Log likelihood = -2490.6464  
Iteration 7:  Log likelihood = -2490.4876  
Iteration 8:  Log likelihood = -2490.4823  
Iteration 9:  Log likelihood = -2490.4822  

Mixed-effects Weibull PH regression             Number of obs     =        600

        Grouping information
        -------------------------------------------------------------
                        |     No. of       Observations per group
         Group variable |     groups    Minimum    Average    Maximum
        ----------------+--------------------------------------------
              birthyear |         12          3       50.0         99
                     id |        201          1        3.0          9
        -------------------------------------------------------------

Integration method: mvaghermite                 Integration pts.  =          7

                                                Wald chi2(0)      =          .
Log likelihood = -2490.4822                     Prob > chi2       =          .
------------------------------------------------------------------------------
          _t |     Hazard   Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
       _cons |   .0135981   .0030285   -19.30   0.000     .0087882    .0210405
-------------+----------------------------------------------------------------
       /ln_p |  -.0299377   .0439552                     -.1160883    .0562129
-------------+----------------------------------------------------------------
birthyear    |
   var(_cons)|   .0695696   .0454997                      .0193072    .2506797
-------------+----------------------------------------------------------------
birthyear>id |
   var(_cons)|   .1868868   .0836493                      .0777297    .4493354
------------------------------------------------------------------------------
Note: Estimates are transformed only in the first equation to hazard.
LR test vs. Weibull model: chi2(2) = 28.83                Prob > chi2 = 0.0000

Note: LR test is conservative and provided only for reference.

Compared to merlin:

. merlin (t M1[birthyear]@1 M2[birthyear>id]@1, family(weibull, failure(failure)))

Fitting fixed effects model:

Fitting full model:

Iteration 0:  Log likelihood = -2494.4923  
Iteration 1:  Log likelihood =  -2493.098  (not concave)
Iteration 2:  Log likelihood = -2491.5986  
Iteration 3:  Log likelihood = -2490.5647  
Iteration 4:  Log likelihood = -2490.4804  
Iteration 5:  Log likelihood = -2490.4801  
Iteration 6:  Log likelihood = -2490.4801  

Mixed effects regression model                             Number of obs = 600
Log likelihood = -2490.4801
------------------------------------------------------------------------------
             | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
t:           |            
M1[birthye~] |          1          .        .       .            .           .
M2[birthye~] |          1          .        .       .            .           .
       _cons |  -4.298231   .2228186   -19.29   0.000    -4.734947   -3.861515
  log(gamma) |  -.0298658   .0439616    -0.68   0.497     -.116029    .0562974
-------------+----------------------------------------------------------------
birthyear:   |            
      sd(M1) |   .2637157   .0862376                      .1389264    .5005957
-------------+----------------------------------------------------------------
id:          |            
      sd(M2) |   .4325904   .0967619                      .2790486    .6706159
------------------------------------------------------------------------------

Given that merlin just works, I suspect it is a matter of parsing the model formula with multiple || separators – unless I am doing something silly here?