an example where glmer is much faster than Stan

143 views
Skip to first unread message

Andrew Gelman

unread,
Jul 14, 2015, 11:14:22 PM7/14/15
to stan...@googlegroups.com, Jalaj Bhandari
Hi all. I came across this example in my applied work. I’ve scrambled the data and replaced with fake outcomes so as not to violate any privacy with the dataset. Attached on this message are a zipfile with the data, the Stan code, and the R code.

Here’s the story: I have a relatively simple hierarchical model (nonnested, with 3 batches of varying intercepts of dimension 51, 15, and 51*15) and 150,000 data points. Running Stan for 4 chains in parallel, 100 iterations each, takes about 3 times as long as glmer. But really most people would run Stan for at least 1000 iterations, so Stan is much slower than glmer here.

This is a big motivation for us to get MML working in Stan, also a great example for trying out the new ADVI on a real live problem!

In the attached code I’m fitting the model to a random subset of 20,000 data points, just so it will run faster and you can see what’s going on.

For what it’s worth, glmer and Stan give similar estimates for the regression coefficients but they differ on the estimated variance parameters, where Stan seems to be doing much better. That’s no surprise because glmer has no regularization on these variance parameters. Our MML and ADVI will allow these, so they should outperform glmer.

So our goal is to have a quick approx algorithm that runs fast. The example is a motivator.

P.S. Stan runs faster on these simulated data than on the real data. Folk theorem in action, I guess.

nbme_fake.zip

Ben Goodrich

unread,
Dec 6, 2015, 1:26:56 AM12/6/15
to stan development mailing list, gel...@stat.columbia.edu

stan_glmer with meanfield ADVI ran on the whole fake dataset in about 15 minutes with a QR reparameterization, which is only about twice as long as glmer with its checks turned off. The results look okay

cbind(glmer = fixef(fit1), stan_glmer = fixef(fit2))
                  glmer  stan_glmer
(Intercept) -4.13102216 -3.22480000
S2          
-0.01073689 -0.01429655
year        
0.15051490  0.15833861
female      
-0.65045473 -0.77372205

> cbind(glmer = VarCorr(fit1), stan_glmer = VarCorr(fit2))
                glmer      stan_glmer
specialty
:state 0.1509191  0.07504937
state          
0.05464908 0.04413998
specialty      
0.0636179  0.1038228

Ben

Reply all
Reply to author
Forward
0 new messages