dear statalist members,

I am trying to estimate the posterior means of two crossed random effects from a logistic random intercept model (estimated with melogit). Although the model is quite huge (some 24000 observations, around 25 individual-level covariates and two crossed random effects), estimation nevertheless converges without problems after some time (around 50 hours). However, when I try to calculate posterior means of both random effects(with predict re_*, reffects) this takes an appearingly endless amount of time (currently 15 days).

To check how fast this type of postestimation works with a much smaller problem, I have created an artificial data set with just 200 observations, two individual-level covariates and two crossed random effects for two grouping variables with 12 and 18 categories. Model estimation runs smoothly and the results conform to what I would expect, given the structure of the artificial data (all syntax and output below). However, the calculation of the posterior means again takes a lot of time and has not come to an end by now (after approximately 3 days). Although the user is warned that "computing empirical Bayes means for a crossed-effects model is very time consuming", calculation seems extremely slow, given the apparently limited size of the estimation problem (?).

I should probably mention that I have no problems to estimate posterior modes (with predict …, reffects rmodes).

It would appreciate very much any comments and ideas to the following points:
1. Is there some fundamental issue with the calculation of posterior means for crossed random effects, I should be aware of?
2. Has anybody successfully used "predict …, reffects" after a logistic mixed model with crossed random effects? What's the experience regarding the amount of time needed?
3. Are there any alternative (and preferably faster) estimation methods for posterior means in this type of situation?

Thanks a lot for any suggestions!
Stefan


My do-file:

Code:
clear
set obs 200

gen x1 = rnormal(0,1)
gen x2 = rnormal(0,1)

gen u1 = rnormal(0,1)
gen u2 = rnormal(0,1)
correlate x* u*

xtile re1 = u1, nq(12)
xtile re2 = u2, nq(18)

by re1, sort: egen postmean1 = mean(u1)
by re2, sort: egen postmean2 = mean(u2)

gen e = rnormal(0,1.5)

gen z = (0.7*x1)+(-0.5*x2)+postmean1+postmean2+e
gen prob = 1/(1+exp(-1*z))
sum prob
recode prob (min/0.6=1)(*=0), gen(d)

melogit d x1 x2 || _all:R.re1 || re2:, diff

* started ca 15:30 pm 21. 12. 2019
predict pre*, reffects

Model-Output:

Array