MELOGIT initial values not feasible error

[Tim Anderson]

Hello,

I am attempting to fit a fairly straightforward melogit model to examine the relationship between the dichotomous outcome of receiving a medical test and the binary predictor of year (2009 vs 2013), with clustering by hospital and a sample size of ~500,000. (This is a simplified version of my final model which has multiple covariates, but I receive the same error with this simplified model)

I am looking at two different dichotomous tests as outcomes, test1 and test2. When I run the same melogit model for test1 I receive a result but for the test2 I receive an "initial values not feasible error." There is no missing data for either test, receiving test2 is much more common, as shown below using frequency tables.

I have tried running a simple logit model to generate start values, as was suggested in a prior post (https://www.stata.com/statalist/arch.../msg00906.html) but this did not change the error.

Are there other strategies anyone would suggest for getting around the initial values not feasible error in this scenario?

Appreciate your help!
Tim

. tab test1 year
| Calendar year
test1 | 2009 2013 | Total
-----------+----------------------+----------
0 | 233,925 247,522 | 481,447
1 | 29,669 45,642 | 75,311
-----------+----------------------+----------
Total | 263,594 293,164 | 556,758

. tab test2 year
| Calendar year
test2 | 2009 2013 | Total
-----------+----------------------+----------
0 | 48,884 46,859 | 95,743
1 | 214,710 246,305 | 461,015
-----------+----------------------+----------
Total | 263,594 293,164 | 556,758

melogit test1 i.year || hospital_id:, or

Fitting fixed-effects model:

Iteration 0: log likelihood = -219971.48
Iteration 1: log likelihood = -219518.87
Iteration 2: log likelihood = -219517.63
Iteration 3: log likelihood = -219517.63

Refining starting values:

Grid node 0: log likelihood = -195609.94

Fitting full model:

Iteration 0: log likelihood = -195609.94
Iteration 1: log likelihood = -195490.36
Iteration 2: log likelihood = -195418.4
Iteration 3: log likelihood = -195331.55
Iteration 4: log likelihood = -195262.55
Iteration 5: log likelihood = -195243.63
Iteration 6: log likelihood = -195241.22
Iteration 7: log likelihood = -195241.19
Iteration 8: log likelihood = -195241.19

Mixed-effects logistic regression Number of obs = 556,758
Group variable: hospital_id Number of groups = 676

Obs per group:
min = 3
avg = 823.6
max = 6,786

Integration method: mvaghermite Integration pts. = 7

Wald chi2(1) = 2242.83
Log likelihood = -195241.19 Prob > chi2 = 0.0000
------------------------------------------------------------------------------
test1 | Odds Ratio Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
|
year |
2013 | 1.496314 .0127331 47.36 0.000 1.471565 1.52148
|
_cons | .040664 .0025562 -50.94 0.000 .0359504 .0459957
-------------+----------------------------------------------------------------
hospital_id |
var(_cons)| 2.329856 .1600997 2.036279 2.665758
------------------------------------------------------------------------------
LR test vs. logistic model: chibar2(01) = 48552.88 Prob >= chibar2 = 0.0000


melogit test2 i.year || hospital_id:, or

Fitting fixed-effects model:

Iteration 0: log likelihood = -255276.2
Iteration 1: log likelihood = -255226.3
Iteration 2: log likelihood = -255226.29

Refining starting values:

Grid node 0: log likelihood = -224213.99

Fitting full model:

initial values not feasible
r(1400);

[Weiwen Ng]

Try the same command with the -noestimate- option, and see what the allegedly infeasible starting values were. In this syntax, you'd have them reported on the log-odds scale.

melogit test2 i.year || hospital_id:, noestimate

Reading the section of the -meglm- manual which discusses how to search for better starting values, it may be that Stata already tried its default starting values and found them infeasible. The problem, the default for -meglm- and -melogit- appears to be fitting a fixed effect-only model, and then using those coefficients and a random effects variance of 1 as starting values. If so, then fitting the logistic model and saving the coefficients didn't help. In any case, try this and see what the allegedly infeasible starting values were. If the variance somehow started near 0 or was something very large like 10, then perhaps forcing that starting value back to something reasonable would be a good idea. You could use the option -startvalues(iv)- as a first step, which is recommended by the manual. Or you could use the variance of 2.23 from your other model as a prior, perhaps. Or maybe, because most people appear to get test 2, you could use something more like 1 / 2.23 as the starting value.

As a worked example of the latter, here's some code using a stock Stata dataset. Note: it's easier to read your code when you use the code delimiters; you can hit the hex (#) button in the post control panel, or you can manually type the code delimiters - they're CODE and \CODE, each enclosed in square brackets ([ and ]), with the code written between the delimiters.

webuse bangladesh
logit c_use urban age child*
melogit c_use urban age child* || district:, noestimate

You'll see that Stata is already using a model with only fixed effects for start values in this example. It decided to start with the variance for the random intercept for district at 1.

melogit c_use urban age child* || district:
melogit c_use urban age child* || district:, startgrid(2)

The second model makes Stata start its estimation with the variance of the random effects set at 2, instead of the default 1. In this case, we wound up with the same maximum likelihood estimates. The final model log likelihoods are identical, which is good. In fact, I got the same results when I set the starting value for the variance to 20000, which would be ridiculous in real life.

[Steve Samuels]

Weiwen gives good advice. In addition, meqrlogit and the difficult option might help.

Another possible cause of your problem is that the model is too simple. Estimating the single random effect could be difficult because part of variation is due to differing hospital/patient characteristics.

[Tim Anderson]

Thank you Weiwen and Steve, I appreciate your thoughtful advice. I took your advice first to try the noestimate approach, and obtained the following:

[CODE]
melogit test2 i.year || hosted:, noestimate

Fitting fixed-effects model:

Iteration 0: log likelihood = -255276.2
Iteration 1: log likelihood = -255226.3
Iteration 2: log likelihood = -255226.29

Refining starting values:

Posting starting values:

Mixed-effects logistic regression Number of obs = 556,758
Group variable: hospid Number of groups = 676

Obs per group:
min = 3
avg = 823.6
max = 6,786

Integration method: mvaghermite Integration pts. = 7

( 1) [test2]2009b.year = 0
------------------------------------------------------------------------------
comb_ekg | Coef. Legend
-------------+----------------------------------------------------------------
|
year |
2013 | .1795896 _b[test2:2013.year]
|
_cons | 1.479838 _b[test2:_cons]
-------------+----------------------------------------------------------------
hosted |
var(_cons)| 1 _b[var(_cons[hospid]):_cons]
------------------------------------------------------------------------------
Note: The above coefficient values are starting values and not the result of
a fully fitted model.[

/CODE]

Though I am not sure how to interpret this result. I then tried each of the 4 alternative startvalues options which you suggested as well as the startgrid option set at both startgrid(2) and startgrid(1/2.23). However, each of these approaches gave me the same "initial values not feasible" error. Appreciate any additional advice you might be able to provide.

Thanks,
Tim

[Tim Anderson]

Sorry the prior post did not format correctly.

Thank you Weiwen and Steve, I appreciate your thoughtful advice. I took your advice first to try the noestimate approach, and obtained the following:

melogit test2 i.year || hospid:, noestimate

Fitting fixed-effects model:

Iteration 0: log likelihood = -255276.2  
Iteration 1: log likelihood = -255226.3  
Iteration 2: log likelihood = -255226.29  

Refining starting values:

Posting starting values:

Mixed-effects logistic regression Number of obs = 556,758
Group variable: hospid Number of groups = 676

Obs per group:
min = 3
avg = 823.6
max = 6,786

Integration method: mvaghermite Integration pts. = 7

( 1) [test2]2009b.year = 0
------------------------------------------------------------------------------
comb_ekg | Coef. Legend
-------------+----------------------------------------------------------------
|
year |
2013 | .1795896 _b[test2:2013.year]
|
_cons | 1.479838 _b[test2:_cons]
-------------+----------------------------------------------------------------
hosted |
var(_cons)| 1 _b[var(_cons[hospid]):_cons]
------------------------------------------------------------------------------
Note: The above coefficient values are starting values and not the result of
a fully fitted model.

Though I am not sure how to interpret this result. I then tried each of the 4 alternative startvalues options which you suggested as well as the startgrid option set at both startgrid(2) and startgrid(1/2.23). However, each of these approaches gave me the same "initial values not feasible" error. Appreciate any additional advice you might be able to provide.

Thanks,
Tim

[Joseph Coveney]

Stata's having trouble estimating the variance component, You could try starting values with a very small number for it, getting over the initial values not feasible hurdle (it's essentially fitting a fixed-effects-only model at that point, which you know will be feasible), and then allowing the maximization to seek larger values for it.

Also, in order to speed things up while you see whether it's going to work, you could do the first run with the intmethod(laplace).

Something like the following.

melogit test2 i.year || hospid:, noestimate
matrix define B = e(b)

// Choose something arbitrarily small, but larger than -epsdouble()-, for the variance component
matrix define B[1, 4] = 1e-8

// See whether it gets you over the hurdle so that maximization can proceed
melogit test2 i.year || hospid: , from(B)  startvalues(iterate(0)) intmethod(laplace)

If that works, then you might try going back to the default (mvaghermite) integration method.

[Tim Anderson]

Thanks Joseph!

I tried your approach which did run with the laplace integration method option, though I had to remove the start values(iterate(0)) term, as it caused repeated errors.

However, then when I tried to switch back to the default (mvaghermite) approach I receive a discontinuity error.

#omitting noestimate model output for space
melogit comb_ekg i.year  || ahaid:, noestimate

matrix define B = e(b)

matrix define B[1, 4] = 1e-8

melogit test2 i.year || hospid: , from(B) intmethod(laplace)

Fitting fixed-effects model:

Iteration 0:   log likelihood =  -255276.2  
Iteration 1:   log likelihood =  -255226.3  
Iteration 2:   log likelihood = -255226.29  

Refining starting values:

Grid node 0:   log likelihood =  -255226.2

Fitting full model:

Iteration 0:   log likelihood =  -255226.2  
Iteration 1:   log likelihood = -246558.26  
Iteration 2:   log likelihood = -240642.05  
Iteration 3:   log likelihood =  -235970.9  
Iteration 4:   log likelihood =  -232358.6  
Iteration 5:   log likelihood = -229666.23  
Iteration 6:   log likelihood = -227734.88  
Iteration 7:   log likelihood =  -226397.8  
Iteration 8:   log likelihood = -225503.01  
Iteration 9:   log likelihood = -224925.69  
Iteration 10:  log likelihood = -224569.91  
Iteration 11:  log likelihood = -224364.76  
Iteration 12:  log likelihood = -224258.48  
Iteration 13:  log likelihood = -224212.89  
Iteration 14:  log likelihood = -224199.25  
Iteration 15:  log likelihood = -224197.29  
Iteration 16:  log likelihood = -224197.22  
Iteration 17:  log likelihood = -224197.22  

Mixed-effects logistic regression               Number of obs     =    556,758
Group variable:           hospid                 Number of groups  =        676

                                                Obs per group:
                                                              min =          3
                                                              avg =      823.6
                                                              max =      6,786

Integration method:     laplace

                                                Wald chi2(1)      =     779.26
Log likelihood = -224197.22                     Prob > chi2       =     0.0000
------------------------------------------------------------------------------
    test2    |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
             |
        year |
       2013  |   .2177583   .0078007    27.92   0.000     .2024692    .2330474
             |
       _cons |   1.385414    .036566    37.89   0.000     1.313746    1.457082
-------------+----------------------------------------------------------------
hospid       |
   var(_cons)|   .8442908   .0479567                      .7553405     .943716
------------------------------------------------------------------------------
LR test vs. logistic model: chibar2(01) = 62058.13    Prob >= chibar2 = 0.0000

Then, without the laplace option

melogit test2 i.year || hospid: , from(B)

Fitting fixed-effects model:

Iteration 0:   log likelihood =  -255276.2  
Iteration 1:   log likelihood =  -255226.3  
Iteration 2:   log likelihood = -255226.29  

Refining starting values:

Grid node 0:   log likelihood =  -255226.2

Fitting full model:

Iteration 0:   log likelihood =  -255226.2  
Iteration 1:   log likelihood = -229241.91  
Iteration 2:   log likelihood = -226377.45  
Iteration 3:   log likelihood = -225270.05  
Iteration 4:   log likelihood = -224778.74  
Iteration 5:   log likelihood = -224482.54  
cannot compute an improvement -- discontinuous region encountered

[Weiwen Ng]

...
Then, without the laplace option

melogit test2 i.year || hospid: , from(B)

Fitting fixed-effects model:
...

Iteration 0: log likelihood = -255226.2
Iteration 1: log likelihood = -229241.91
Iteration 2: log likelihood = -226377.45
Iteration 3: log likelihood = -225270.05
Iteration 4: log likelihood = -224778.74
Iteration 5: log likelihood = -224482.54 
cannot compute an improvement -- discontinuous region encountered

At this point, I am running short of good ideas. But, just to check, did you try any of the other two available integration methods (options -intmethod(mcaghermite)- and -intemethod(ghermite)-)?

[Joseph Coveney]

At this point, I'd try something that Steve Samuels mentioned upthread, either meqrlogit or the difficult option (or both).

melogit comb_ekg i.year  || ahaid:, noestimate // Huh?
matrix define B = e(b)
matrix define B[1, 4] = 1e-8
meqrlogit test2 i.year || hospid: , from(B)

Another alternative: because the iterations converged with the first test, you might be able to exploit that in order to gain leverage in attaining convergence (i.e., estimating the variance) with the second. If your dataset is arranged like I think it is, then you could fit the two simultaneously without any dataset manipulation.

gsem (test1 <- i.year M[hospid], logit)  (test2 <- i.year M[hospid], logit), nocnsreport nodvheader

My guess, though, is that in order to get it to work you'll have to impose the not unreasonable constraint of a single, common variance. Something like this.

constraint define 1 _b[test1:M[hospid]] = 1
constraint define 2 _b[test2:M[hospid]] = 1
gsem (test1 <- i.year M[hospid], logit)  (test2 <- i.year M[hospid], logit), constraints(1/2) nocnsreport nodvheader

[Tim Anderson]

Hello,

Thank you again. I did end up using the meqrlogit approach which worked without the difficult option. I have one follow up question related to post-estimation with meqrlogit as it appears to be different than when using melogit.

When I have used melogit in the past I have often used the margins post-estimation command to obtain marginal predicted mean probabilities (in the case for each year), which generally takes the following form:

melogit test2 i.year || hospid:, or
(output omitted for space)

margins year

Adjusted predictions                            Number of obs     =  1,006,900
Model VCE    : OIM

Expression   : Pr(test2), predict()

------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        year |
       2009  |   .4118406   .0007168   574.52   0.000     .4104356    .4132456
       2013  |   .2903434   .0006203   468.08   0.000     .2891277    .2915592
------------------------------------------------------------------------------

However, now that I am using meqrlogit when I run the post estimation margins command I get a very different type of result:

meqrlogit test1 i.year || hospid:, or
(output omitted for space)

margins year

Predictive margins                              Number of obs     =  1,001,215

Expression   : Linear prediction, fixed portion, predict(xb)

------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        year |
       2009  |  -4.669686   .0864848   -53.99   0.000    -4.839194   -4.500179
       2013  |  -3.252255   .0862166   -37.72   0.000    -3.421237   -3.083274
------------------------------------------------------------------------------

It appears that the output is similar to if I had specified the xb option and thus am only obtaining linear predictors for the fixed effects. When I tried to specify the mu option for predicted means I received only an error:

margins year, mu
option mu not allowed

I tried reading through the Stata meqrlogit post estimation manual which seems to indicate that it is possible to use margins following meqrlogit regression to obtain predicted mean probabilities so I am not sure what to try next. Appreciate any suggestions you might have.

Sincerely,
Tim

[Joseph Coveney]

Try something like margins year, expression(invlogit(predict(xb))).

As in

sysuse auto
logit foreign i.rep78
margins rep78
margins rep78, expression(invlogit(predict(xb)))

[Tim Anderson]

Thanks Joseph,

I tried your approach and it ran, but the resulting values were very far from the results that I received when running the same model as a regular logit model omitting the random effects. The ORs for the logit vs meqrlogit approaches are similar, thus I am concerned the post estimation approach might not be right? Output below:

Here is the logit model using the regular margin command and the margins with expression specification - output are identical

logit test1 i.year, or

Iteration 0:   log likelihood = -287971.89  
Iteration 1:   log likelihood = -276508.41  
Iteration 2:   log likelihood = -275818.47  
Iteration 3:   log likelihood = -275816.71  
Iteration 4:   log likelihood = -275816.71  

Logistic regression                             Number of obs     =  1,006,900
                                                LR chi2(1)        =   24310.36
                                                Prob > chi2       =     0.0000
Log likelihood = -275816.71                     Pseudo R2         =     0.0422

------------------------------------------------------------------------------
         obs | Odds Ratio   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        year |
       2013  |   3.444853   .0297141   143.40   0.000     3.387104    3.503586
             |
       _cons |   .0403076   .0003042  -425.49   0.000     .0397158    .0409083
------------------------------------------------------------------------------

. margins year, expression(invlogit(predict(xb)))

Adjusted predictions                            Number of obs     =  1,006,900
Model VCE    : OIM

Expression   : invlogit(predict(xb))

------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        year |
       2009  |   .0387459   .0002811   137.84   0.000      .038195    .0392968
       2013  |   .1219243   .0004471   272.69   0.000     .1210479    .1228006
------------------------------------------------------------------------------

. margins year

Adjusted predictions                            Number of obs     =  1,006,900
Model VCE    : OIM

Expression   : Pr(obs), predict()

------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        year |
       2009  |   .0387459   .0002811   137.84   0.000      .038195    .0392968
       2013  |   .1219243   .0004471   272.69   0.000     .1210479    .1228006
------------------------------------------------------------------------------

Here is the meqrlogit model using the regular margins and the margins with expression specification

meqrlogit test1 i.year ||hospid:, or

Refining starting values: 

Iteration 0:   log likelihood = -222474.12  
Iteration 1:   log likelihood = -221660.93  
Iteration 2:   log likelihood = -221391.76  

Performing gradient-based optimization: 

Iteration 0:   log likelihood = -221391.76  
Iteration 1:   log likelihood = -221211.81  
Iteration 2:   log likelihood = -221203.88  
Iteration 3:   log likelihood = -221203.83  
Iteration 4:   log likelihood = -221203.83  

Mixed-effects logistic regression               Number of obs     =  1,006,900
Group variable: hospid                          Number of groups  =       1048

                                                Obs per group:
                                                              min =          1
                                                              avg =      960.8
                                                              max =     14,436

Integration points =   7                        Wald chi2(1)      =   21890.93
Log likelihood = -221203.83                     Prob > chi2       =     0.0000

------------------------------------------------------------------------------
         obs | Odds Ratio   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        year |
       2013  |   3.896306   .0358153   147.96   0.000     3.826738    3.967139
             |
       _cons |    .011363   .0009801   -51.91   0.000     .0095957    .0134558
------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
hospid: Identity             |
                  var(_cons) |    6.88235   .3896002      6.159585    7.689925
------------------------------------------------------------------------------
LR test vs. logistic model: chibar2(01) = 1.1e+05     Prob >= chibar2 = 0.0000

. margins year

Adjusted predictions                            Number of obs     =  1,006,900

Expression   : Linear prediction, fixed portion, predict(xb)

------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        year |
       2009  |  -4.477394   .0862507   -51.91   0.000    -4.646442   -4.308345
       2013  |  -3.117365   .0860058   -36.25   0.000    -3.285933   -2.948797
------------------------------------------------------------------------------

. margins year, expression(invlogit(predict(xb)))

Adjusted predictions                            Number of obs     =  1,006,900

Expression   : invlogit(predict(xb))

------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        year |
       2009  |   .0112353   .0009582    11.73   0.000     .0093573    .0131133
       2013  |   .0423966   .0034918    12.14   0.000     .0355529    .0492404
------------------------------------------------------------------------------

So using logit approaches the marginal predicted probabilities are 3.87% and 11.2% in 2009 and 2013
While for the meqrlogit approach the marginal predicted probabilities are 1.12% and 4.24% in 2009 and 2013

Could this still be an issue specific with meqrlogit?

Thanks,
Tim

[Joseph Coveney]

They're remarkably similar, actually, given that you have a relatively rare success and that the hospital variance is 7. The latter strikes me as large for a logistic regression model.

[Tim Anderson]

Sorry, I was using a different test example than the original post (I have two related tests outcomes which both have the same problem with starting values when using melogit)

For the test that I posted the meqrlogit results above...the unadjusted 2x2 of this test is:

| Calendar year
test1 | 2009 2013 | Total
-----------+----------------------+----------
0 | 453,115 470,228 | 923,343
| 96.13 87.81 | 91.70
-----------+----------------------+----------
1 | 18,264 65,293 | 83,557
| 3.87 12.19 | 8.30
-----------+----------------------+----------
Total | 471,379 535,521 | 1,006,900
| 100.00 100.00 | 100.00

This sensibly fits closely with the basic logistic regression predicted probability results (3.87% and 12.19%) as the model has no other covariates, but is quite different from the meqrlogit results (1.12% and 4.24%).