meologit - problems with post-estimation prediction

[JasonP]

I'm having the some problems with the meologit command and wondered whether anyone had advice on a work-around?

Essentially, my analysis executes fine and converges nicely. However, any attempt to generate predictions of the marginal probabilities or random effects returns the error:

predict re*, remeans
(calculating posterior means of random effects)
(using 20 quadrature points)
could not compute empirical Bayes' means;
missing values were returned by the evaluator

I have tried all of the integration methods and (where appropriate) a variety of integration points, but the same error is returned each time.

Anyone else having this problem or suggestions for fixing it?

Thanks in advance!

[JasonP]

Actually, I may have answered my own question (in part), because my model is not using adaptive quadrature - instead it is reverting automatically to ghermite rather than mvaghermite as requested (although I cannot see anywhere why this should affect the post-estimation prediction?).

Is there a way to force adaptive quadrature with meologit - like the use of from() in the gsem procedure?

Alternatively, has anyone any recommendations on how to get post-estimation predictions for a model that does not use adaptive quadrature? Or am I barking up the wrong tree?

[Marcos Almeida]

Please see it the comments apply to your problem: there has been some recent discussion on that on the thread http://www.statalist.org/forums/foru...edict-meologit

[JasonP]

Thanks Marcos. I guess that post still leaves my question unanswered. Why will meologit not use adaptive quadrature? Is there no way to force it by feeding parameters from a non adaptive model? Why the exact same model converges using adaptive quadrature in GLAMM?

[Marcos Almeida]

Hi, Jason,

You didn't show your commands to perform meologit. But I wonder if there is something related to the options.

Below, an extract of the Stata manual ( here, the link: http://www.stata.com/manuals13/memeologit.pdf) on meologit and the integration methods:

intmethod(intmethod) specifies the integration method to be used for the random-effects model.
mvaghermite performs mean and variance adaptive Gauss–Hermite quadrature; mcaghermite
performs mode and curvature adaptive Gauss–Hermite quadrature; ghermite performs nonadaptive
Gauss–Hermite quadrature; and laplace performs the Laplacian approximation, equivalent to mode
curvature adaptive Gaussian quadrature with one integration point.
The default integration method is mvaghermite unless a crossed random-effects model is fit, in
which case the default integration method is laplace. The Laplacian approximation has been
known to produce biased parameter estimates; however, the bias tends to be more prominent in
the estimates of the variance components rather than in the estimates of the fixed effects.meologit — Multilevel mixed-effects ordered logistic regression 5
For crossed random-effects models, estimation with more than one quadrature point may be
prohibitively intensive even for a small number of levels. For this reason, the integration method
defaults to the Laplacian approximation. You may override this behavior by specifying a different
integration method.

Hopefully it helps.

Best,

Marcos

[JasonP]

Thanks Marcos, the model is really quite simple and my code is as follows. meologit risk || locat:, intpoints(20) intmethod (mvaghermite) 'RISK' is a four point ordinal variable. 'LOCAT' has 11 groups The dataset has 4406 observations clustered in 11 locations. I wonder whether there is simply insufficient variance at level 2, but as I said, the exact same model converges with adaptive quadrature in GLLAMM. This makes me suspicious about the GLLAMM results, but I don't know why. Thanks!!

[Marcos Almeida]

Hi, Jason,


From "what is new" in Stata 13 (here, the link:http://page2rss.com/701012f4d48dc18fbea1e3056ab48af7 ), we get the following information:

5. predict after gsem and me would incorrectly return zero-valued empirical Bayes' mean estimates when intmethod(mvaghermite) was not a viable integration method. In this case, gsem and me will now report a note below the coefficient table acknowledging that the results were produced using nonadaptive quadrature, and predict will issue the error message "could not compute the empirical Bayes' means"

Also, regarding the clustering form of your model, I wonder if you have tryed to add "vce(robust)".

Hopefully these suggestions will help you, because this is the furtherst I can go on this matter.


Best,

Marcos

[JasonP]

Thanks Marcos - I'm hoping someone from StataCorp can chime in with some further information and direction.

I found an earlier post on the Statalist archive about the me commands having difficulty with groups of large sample sizes -
http://hsphsun3.harvard.edu/cgi-bin/...rticle-73.html

According to this thread, the me commands do not use numerical derivatives (like GLLAMM does) and thus have difficulty with small groups of large sample size (n=+1000).

Still, it'd be great is StataCorp had any advice on working around this problem since the solution suggested in that thread (using evaltype(gf0)) does not resolve the issue in my case. Even if the advice is to "stick with GLLAMM" for this special case, is useful.

[Nina Rahman]

Hi Jason,

I was wondering if you managed to find a solution to your problem? I'm encountering a similar problem.

My code is as follows:
meologit y x1 x2 x3 x4 x5 x6... x13 || s024:, covariance(unstructured) intmethod(mcaghermite)

With y being ordinal, 8 level 1 independent variables, and 5 level 2 independent variables. I have 14,640 observations in 69 groups.

The output seems to be fine, but when I try to run post-estimation tests, the following error appears:

And I am really unsure of what to do next. Does anyone have any idea how I can go about this?

Thank you very much!

[Joseph Coveney]

And I am really unsure of what to do next. Does anyone have any idea how I can go about this?

My suggestion is to start out with the following code snippet and report back to the list.

generate byte missing = 0
forvalues i = 1/13 {
      quietly replace missing = 1 if missing(x`i')
}
quietly replace missing = 1 if missing(y)
tabulate missing
assert missing(outcome) if missing

Here's why:

.ÿversionÿ14.0

.ÿ
.ÿclearÿ*

.ÿsetÿmoreÿoff

.ÿsetÿseedÿ`=date("2015-06-19",ÿ"YMD")'

.ÿquietlyÿsetÿobsÿ69

.ÿ
.ÿgenerateÿbyteÿs024ÿ=ÿ_n

.ÿgenerateÿdoubleÿs024_uÿ=ÿrnormal()

.ÿforvaluesÿiÿ=ÿ1/8ÿ{
ÿÿ2.ÿÿÿÿÿÿÿÿÿgenerateÿdoubleÿx`i'ÿ=ÿruniform()
ÿÿ3.ÿ}

.ÿquietlyÿexpandÿ`=floor(14640ÿ/ÿ69)'

.ÿforvaluesÿiÿ=ÿ9/13ÿ{
ÿÿ2.ÿÿÿÿÿÿÿÿÿgenerateÿdoubleÿx`i'ÿ=ÿruniform()
ÿÿ3.ÿ}

.ÿgenerateÿdoubleÿlatÿ=ÿs024_uÿ+ÿrnormal()

.ÿgenerateÿbyteÿyÿ=ÿ0

.ÿquietlyÿcentileÿlat,ÿcentile(25ÿ50ÿ75)

.ÿforvaluesÿiÿ=ÿ1/3ÿ{
ÿÿ2.ÿÿÿÿÿÿÿÿÿquietlyÿreplaceÿyÿ=ÿyÿ+ÿ1ÿifÿlatÿ>=ÿr(c_`i')
ÿÿ3.ÿ}

.ÿ
.ÿmeologitÿyÿc.x*ÿ||ÿs024:,ÿcovariance(unstructured)ÿintmethod(mcaghermite)ÿ///
>ÿÿÿÿÿÿÿÿÿnolrtestÿnolog

Mixed-effectsÿologitÿregressionÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿNumberÿofÿobsÿÿÿÿÿ=ÿÿÿÿÿ14,628
Groupÿvariable:ÿÿÿÿÿÿÿÿÿÿÿÿs024ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿNumberÿofÿgroupsÿÿ=ÿÿÿÿÿÿÿÿÿ69

ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿObsÿperÿgroup:
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿminÿ=ÿÿÿÿÿÿÿÿ212
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿavgÿ=ÿÿÿÿÿÿ212.0
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿmaxÿ=ÿÿÿÿÿÿÿÿ212

Integrationÿmethod:ÿmcaghermiteÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿIntegrationÿpts.ÿÿ=ÿÿÿÿÿÿÿÿÿÿ7

ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿWaldÿchi2(13)ÿÿÿÿÿ=ÿÿÿÿÿÿ10.74
Logÿlikelihoodÿ=ÿ-16190.324ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿProbÿ>ÿchi2ÿÿÿÿÿÿÿ=ÿÿÿÿÿ0.6326
------------------------------------------------------------------------------
ÿÿÿÿÿÿÿÿÿÿÿyÿ|ÿÿÿÿÿÿCoef.ÿÿÿStd.ÿErr.ÿÿÿÿÿÿzÿÿÿÿP>|z|ÿÿÿÿÿ[95%ÿConf.ÿInterval]
-------------+----------------------------------------------------------------
ÿÿÿÿÿÿÿÿÿÿx1ÿ|ÿÿ-.1570915ÿÿÿ.7369334ÿÿÿÿ-0.21ÿÿÿ0.831ÿÿÿÿ-1.601455ÿÿÿÿ1.287271
ÿÿÿÿÿÿÿÿÿÿx2ÿ|ÿÿÿ.3056336ÿÿÿ.7277655ÿÿÿÿÿ0.42ÿÿÿ0.675ÿÿÿÿ-1.120761ÿÿÿÿ1.732028
ÿÿÿÿÿÿÿÿÿÿx3ÿ|ÿÿÿ.2129247ÿÿÿ.7506526ÿÿÿÿÿ0.28ÿÿÿ0.777ÿÿÿÿ-1.258327ÿÿÿÿ1.684177
ÿÿÿÿÿÿÿÿÿÿx4ÿ|ÿÿÿ1.308952ÿÿÿ.7212736ÿÿÿÿÿ1.81ÿÿÿ0.070ÿÿÿÿ-.1047179ÿÿÿÿ2.722623
ÿÿÿÿÿÿÿÿÿÿx5ÿ|ÿÿÿÿÿ.68202ÿÿÿ.6938127ÿÿÿÿÿ0.98ÿÿÿ0.326ÿÿÿÿ-.6778279ÿÿÿÿ2.041868
ÿÿÿÿÿÿÿÿÿÿx6ÿ|ÿÿÿ.0753138ÿÿÿ.6549596ÿÿÿÿÿ0.11ÿÿÿ0.908ÿÿÿÿ-1.208383ÿÿÿÿ1.359011
ÿÿÿÿÿÿÿÿÿÿx7ÿ|ÿÿÿÿ1.12136ÿÿÿ.6765177ÿÿÿÿÿ1.66ÿÿÿ0.097ÿÿÿÿ-.2045905ÿÿÿÿÿ2.44731
ÿÿÿÿÿÿÿÿÿÿx8ÿ|ÿÿ-.0831879ÿÿÿ.7576147ÿÿÿÿ-0.11ÿÿÿ0.913ÿÿÿÿ-1.568085ÿÿÿÿÿ1.40171
ÿÿÿÿÿÿÿÿÿÿx9ÿ|ÿÿ-.0342379ÿÿÿÿ.056055ÿÿÿÿ-0.61ÿÿÿ0.541ÿÿÿÿ-.1441037ÿÿÿÿ.0756278
ÿÿÿÿÿÿÿÿÿx10ÿ|ÿÿÿ.0740932ÿÿÿ.0560015ÿÿÿÿÿ1.32ÿÿÿ0.186ÿÿÿÿ-.0356678ÿÿÿÿ.1838542
ÿÿÿÿÿÿÿÿÿx11ÿ|ÿÿÿÿ.000982ÿÿÿ.0556496ÿÿÿÿÿ0.02ÿÿÿ0.986ÿÿÿÿ-.1080892ÿÿÿÿ.1100531
ÿÿÿÿÿÿÿÿÿx12ÿ|ÿÿÿ.0209931ÿÿÿ.0558528ÿÿÿÿÿ0.38ÿÿÿ0.707ÿÿÿÿ-.0884764ÿÿÿÿ.1304626
ÿÿÿÿÿÿÿÿÿx13ÿ|ÿÿ-.0012397ÿÿÿ.0557312ÿÿÿÿ-0.02ÿÿÿ0.982ÿÿÿÿ-.1104709ÿÿÿÿ.1079914
-------------+----------------------------------------------------------------
ÿÿÿÿÿÿÿ/cut1ÿ|ÿÿÿ.2333423ÿÿÿ1.016606ÿÿÿÿÿ0.23ÿÿÿ0.818ÿÿÿÿÿ-1.75917ÿÿÿÿ2.225854
ÿÿÿÿÿÿÿ/cut2ÿ|ÿÿÿ1.952038ÿÿÿ1.016726ÿÿÿÿÿ1.92ÿÿÿ0.055ÿÿÿÿ-.0407088ÿÿÿÿ3.944785
ÿÿÿÿÿÿÿ/cut3ÿ|ÿÿÿ3.534489ÿÿÿÿ1.01698ÿÿÿÿÿ3.48ÿÿÿ0.001ÿÿÿÿÿ1.541245ÿÿÿÿ5.527732
-------------+----------------------------------------------------------------
s024ÿÿÿÿÿÿÿÿÿ|
ÿÿÿvar(_cons)|ÿÿÿ2.653856ÿÿÿ.4578374ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿ1.892474ÿÿÿÿ3.721558
------------------------------------------------------------------------------

.ÿ
.ÿpredictÿoutcome,ÿmu
(predictionsÿbasedÿonÿfixedÿeffectsÿandÿposteriorÿmeansÿofÿrandomÿeffects)
(usingÿ7ÿquadratureÿpoints)

.ÿ
.ÿcountÿifÿmissing(outcome)
ÿÿ0

.ÿ
.ÿexit

endÿofÿdo-file


.

I also suggest that you omit covariance(unstructured) intmethod(mcaghermite) from your modeling command.

[Nina Rahman]

My suggestion is to start out with the following code snippet and report back to the list.

generate byte missing = 0
forvalues i = 1/13 {
quietly replace missing = 1 if missing(x`i')
}
quietly replace missing = 1 if missing(y)
tabulate missing
assert missing(outcome) if missing

Hi Joseph,

When I use the default integration method, my model doesn't converge.

Anyway, based on your suggestion I ran my model again:
meologit y x1 i.x2 x3 x4 i.x5 i.x6 x7 x8 x9 x10 x11 x12 x13, || s024:, intmethod(mcaghermite)

The independent variables with the i. prefix are binary variables.

I then ran your code and obtained the following output:

generate byte missing = 0

.
. forvalues i = 1/13 {
2.
. quietly replace missing = 1 if missing(x`i')
3.
. }

.
. quietly replace missing = 1 if missing(y)

.
. tabulate missing

missing | Freq. Percent Cum.
------------+-----------------------------------
0 | 14,640 61.77 61.77
1 | 9,061 38.23 100.00
------------+-----------------------------------
Total | 23,701 100.00

.
. assert missing(outcome) if missing
outcome not found
r(111);



Sorry, I didn't mention earlier that I actually have 23701 observations, but only 14640 valid observations after listwise deletion.
Would you be able to assist? Many thanks!

[Joseph Coveney]

When I use the default integration method, my model doesn't converge.

That's probably the first clue, there. From similar posts on the list with others' having trouble fitting mixed-effects ordered-categorical regression models, I'm guessing that you've got what you're going to get with meologit and that set of predictors. You might try what others have attempted, namely, the user-written gllamm, but I would be suspicious of the fit with a nonlinear model with so many predictors and so much difficulty in attaining even a semblance of convergence. You might want to look into whether your response variable is skewed overall, or in one category or another of the binary predictors, or at one end of the range of one or another of the continuous predictors, whether there is near linear dependence among the binary predictors, whether the random effects variance is essentially zero, and so on. Simplifying the model (omitting one or more predictors) might provide further clues as to troublemakers. Collapsing some response categories might help—examine the cut-points in your regression output.

[Nina Rahman]

That's probably the first clue, there. From similar posts on the list with others' having trouble fitting mixed-effects ordered-categorical regression models, I'm guessing that you've got what you're going to get with meologit and that set of predictors. You might try what others have attempted, namely, the user-written gllamm, but I would be suspicious of the fit with a nonlinear model with so many predictors and so much difficulty in attaining even a semblance of convergence. You might want to look into whether your response variable is skewed overall, or in one category or another of the binary predictors, or at one end of the range of one or another of the continuous predictors, whether there is near linear dependence among the binary predictors, whether the random effects variance is essentially zero, and so on. Simplifying the model (omitting one or more predictors) might provide further clues as to troublemakers. Collapsing some response categories might help—examine the cut-points in your regression output.

I will look into that. Thank you for your help, Joseph.