Why is the mixed command slow? Can it be sped up?

[Paul Burkander]

I'm attempting to do a simulation comparing resutls from a random effects model with a random intercept to some other regression models, but the mixed command slows down my simulation too much to be useful. In my experience running mixed effects commands in SAS is relatively quick, and the LMER command in R is relatively quick. Why is mixed so slow, and can it be sped up?

For context, I'm using stata 14. I have a very simple model that I run with "mixed y T || cluster_id:"

Thanks.

[Joseph Coveney]

If your sample size is much larger than, say, 200 cluster IDs, then you could probably get away with fewer. With about 200 and 10 observations each (see code below), mixed takes about 350 milliseconds to converge on my aging laptop. Even with a thousand iterations for the simulation, you're looking at less than six minutes. In the code, I have called for mixed not to compute standard errors of the random effect variances, which also saves time, if you don't need them.

In at least resampling operations, there is often a time advantage to using starting values from the first converged model fit for all successive iterations. Depending upon the nature of the simulations, a similar time advantage might obtain. mixed doesn't allow for starting values as it first uses expectation-maximization to get initial values, but you can fit the same maximum likelihood model with meglm, which does accept starting values. Under similar circumstances as above, it takes about a quarter of a second to converge on average. Again, any advantage is only potential, and will depend upon how good the starting values are for each successive dataset that you simulate.

version 15.1

clear *

set seed `=strreverse("1495995")'
quietly set obs 200

generate int pid = _n
generate double pid_u = rnormal()

quietly expand 10
generate double prd = runiform() - 0.5

generate double rsp = pid_u + rnormal()

timer clear 1
timer on 1
forvalues i = 1/10 {
    quietly mixed rsp c.prd || pid: , nostderr nolrtest
}
timer off 1

timer clear 2
timer on 2
quietly meglm rsp c.prd || pid: , family(gaussian) link(identity) nolrtest
tempname B
matrix define `B' = e(b)
forvalues i = 2/10 {
    quietly meglm rsp c.prd || pid: , family(gaussian) link(identity) nolrtest from(`B')
}
timer off 2

timer list

exit

[Paul Burkander]

Great - thanks for the useful tips!