xtmelogit - How to solve the problem "error obtaining starting values; try fitting a marginal model in order to diagnose the problem"

[Tiago Pereira]

Dear all,

I am been working on a study to examine the effect of two drugs (drug A = 1, drug B = 0) over time (time 0, 1 and 2 semesters) on depression (1 = yes, 0 = no). My Stata version is 14.2

Here is an example data set:

use http://atshaoc.com.br/depression.dta, clear
xtmelogit depression drug##time || patient_id: , intpoints(5)

I have tried a number of different approaches, such as:

logit depression drug##time
matrix a = e(b)
xtmelogit depression drug##time || patient_id: , intpoints(3) refineopts(dif) from(a)

unsuccessfully, however. Would you have any thoughts on how to estimate decent coefficients using any alternative method?

All the best,

Tiago

P.S. Edited for typos

[Jeph Herrin]

Do you really intend to treat -time- as categorical? If you treat it as continuous, this model runs fine:

. use http://atshaoc.com.br/depression.dta, clear

. 
. xtmelogit depression i.drug##c.time || patient_id:

Refining starting values: 

Iteration 0:   log likelihood = -126.95844  
Iteration 1:   log likelihood =   -125.679  
Iteration 2:   log likelihood = -125.66428  

Performing gradient-based optimization: 

Iteration 0:   log likelihood = -125.66428  
Iteration 1:   log likelihood = -125.66427  

Mixed-effects logistic regression               Number of obs     =        274
Group variable: patient_id                      Number of groups  =        100

                                                Obs per group:
                                                              min =          1
                                                              avg =        2.7
                                                              max =          3

Integration points =   7                        Wald chi2(3)      =      40.93
Log likelihood = -125.66427                     Prob > chi2       =     0.0000

------------------------------------------------------------------------------
  depression |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
      1.drug |  -.3969482    .690502    -0.57   0.565    -1.750307    .9564109
        time |   1.821564   .3789941     4.81   0.000      1.07875    2.564379
             |
 drug#c.time |
          1  |   1.408457   .5644762     2.50   0.013     .3021036     2.51481
             |
       _cons |  -2.527696   .5590665    -4.52   0.000    -3.623447   -1.431946
------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
patient_id: Identity         |
                   sd(_cons) |    1.32776   .4202793      .7139806    2.469179
------------------------------------------------------------------------------
LR test vs. logistic model: chibar2(01) = 6.28        Prob >= chibar2 = 0.0061

[Tiago Pereira]

Many thanks, Jeph.

Yes, time should be categorical, because I want to see if drug A is better (or not) than drug B at time points 12 and 24 monhs respectively. By considering time as a continuous variable and not as a factor, I believe the estimates seem ackward under a clinical point of view. Suggestions?

Thank you so much again.

All the best,

Tiago

[Andrew Musau]

You could give the old gllamm command a shot. Install from SSC. Like is the case for melogit, Newton-Raphson fails but adaptive quadrature combined with Newton-Raphson does generate results.

. use http://atshaoc.com.br/depression.dta, clear

. xi: gllamm depression i.drug*i.time, i(patient_id) nocons l(mlogit) f(binom) nip(5) adapt
i.drug            _Idrug_0-1          (naturally coded; _Idrug_0 omitted)
i.time            _Itime_0-2          (naturally coded; _Itime_0 omitted)
i.drug*i.time     _IdruXtim_#_#       (coded as above)

Running adaptive quadrature
Iteration 0:    log likelihood = -132.55331
Iteration 1:    log likelihood = -135.51525
Iteration 2:    log likelihood =  -122.5878
Iteration 3:    log likelihood = -121.74142
Iteration 4:    log likelihood = -121.63818
Iteration 5:    log likelihood = -121.31046
Iteration 6:    log likelihood = -121.33407
Iteration 7:    log likelihood = -121.35016
Iteration 8:    log likelihood =  -121.3331
Iteration 9:    log likelihood = -121.33398
Iteration 10:    log likelihood = -121.32524
Iteration 11:    log likelihood = -121.32447
Iteration 12:    log likelihood = -121.32447


Adaptive quadrature has converged, running Newton-Raphson
Iteration 0:   log likelihood = -121.32447  
Iteration 1:   log likelihood = -121.32447  (backed up)
Iteration 2:   log likelihood = -121.32442  
Iteration 3:   log likelihood = -121.32442  
 
number of level 1 units = 274
number of level 2 units = 100
 
Condition Number = 2904.4699
 
gllamm model 
 
log likelihood = -121.32442
 
-------------------------------------------------------------------------------
   depression |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
--------------+----------------------------------------------------------------
     _Idrug_1 |  -56.60878    568.473    -0.10   0.921    -1170.795    1057.578
     _Itime_1 |   4.283289   1.432229     2.99   0.003     1.476171    7.090407
     _Itime_2 |   3.449187   1.283677     2.69   0.007     .9332252    5.965148
_IdruXtim_1_1 |   56.68102    568.473     0.10   0.921    -1057.506    1170.868
_IdruXtim_1_2 |   58.32228   568.4739     0.10   0.918    -1055.866    1172.511
-------------------------------------------------------------------------------
 
 
Variances and covariances of random effects
------------------------------------------------------------------------------

 
***level 2 (patient_id)
 
    var(1): 18.273264 (11.386074)
------------------------------------------------------------------------------

[Tiago Pereira]

Many thanks, Andrew! Extremely smart alternative. I don't know how to get the marginal effects, such as we typically obtain:

contrast drug@time, nowald effects
margins drug##time

Any suggestions? I somehow feel those parameters are poorly estimated.

Thanks again for all your help so far.


All the best,

Tiago

[Andrew Musau]

Here is one suggestion.

http://statalist.blogspot.com/2006/0...l-effects.html

[Tiago Pereira]

Thank you, Andrew! Very useful!

All the best,

Tiago

[Tiago Pereira]

Dear all,

In theory, should the following models be similar in terms of z- and p-values?

use http://atshaoc.com.br/depression.dta, clear
* Fails
xtmelogit depression drug##time || patient_id:
* Runs OK
meprobit event i.time#i.trt || id:

All the best,

Tiago

[Tiago Pereira]

I mean:

meprobit event i.drug##i.time || patient_id: 
 xtmelogit depression i.drug##i.time || patient_id:

[Andrew Musau]

Why are you using different outcome variables? For a two-level random intercept model, xtmelogit and melogit will give you identical results. The marginal effects for melogit and meprobit will also be very similar. See the example below.

webuse bangladesh, clear
xtmelogit c_use urban age child* || district:
melogit c_use urban age child* || district:
margins, dydx(*)
meprobit c_use urban age child* || district:
margins, dydx(*)

So yes, I would say that in the case of estimating a two-level random intercept model, empirically, there is not much to choose between melogit and meprobit.

[Tiago Pereira]

Thank you so much, Andrew! Extremely useful. Many, many thanks!

All the best,

Tiago