Mixed models in STAN vs. glmer

[Larry Hunsicker]

Greetings, STAN-list users.  Another question from a STAN (and Bayesian) neophyte concerning mixed model analysis.

I have been trying to understand better the differences in assumptions and conclusions from fully Bayesian analyses (STAN) and the frequentist analyses (with which I am more familiar) of mixed models.  To be explicit, I have been looking at analyses of the BLOCKER data as in the BUGS, and now the STAN models, and using glmer from the R package nlme4.  The estimates of "d" in the two sets of analyses are virtually identical -- -0.25 using STAN, and -0.26 using glmer.  But the estimates of the SE of d are quite different, about 0.12 for STAN (taking the sqrt of the estimated sigmasqu_delta), but only 0.05 using glmer.  In this case, the conclusions from the two approaches are the same.  But it appears that STAN is estimating a relatively much larger SEM than glmer, and it is easy to see that with a d a bit closer to 0, the conclusions might be quite different.

I am using RSTAN 2.6.0, R version 3.2.0, RStudio 0.99.441, and lme4 1.1-7, on a Windows 8.1 machine with all updates.  I used blocker.stan blocker.dat, and blocker.init, unchanged from the files in the STAN examples.   For the glmer analysis, I restructured blocker.dat so that there was one case per record, with three columns:  death (1 vs 0), study number, and treatment group.  My glmer function was:

m1 <- glmer(death ~ group + (1 | study), data = bldat2, family = binomial)
//group is treated vs. control (1 vs 0);  study is the number of the specific cited study.

I also checked a variant of blocker.stan, using normal distributions instead of t distributions (4 df) for d.  The results of this modified blocker.stan were essentially identical to the results from the original blocker.stan.

My set of questions includes:
a)  What are the differences in assumptions that lead to the different estimates of the SE of d?  Is glmer missing a source of variance that STAN is finding, and if so, where?  Alternatively, is the difference due to error introduced by approximations introduced in the calculations?
b)  Given that the model results are somewhat different, which is more likely to represent "reality"?
c)  My understanding is that the estimates of SE(d) in STAN come directly from sampling the posterior estimates of d, but that with glmer, the SE is estimated as the partial slope of the inverted information matrix, which, as I understand it, is only asymptotically equivalent to a normal distribution.  Does this imply that the STAN estimate is more accurate than the glmer estimate?  If I am correct in my understanding of this, the implication is that one should do most of one's mixed model analysis using STAN rather than glmer.  Am I way off?
d)  If someone could point me to a useful reference that will answer my questions above, I would be deeply grateful to you.  (I am slowly making my way through Gelman's book, "Bayesian Data Analysis.") 

Many thanks in advance
.
Larry Hunsicker
L. G. Hunsicker, M.D.
Professor Emeritus, University of Iowa College of Medicine

[Bob Carpenter]

> On Jun 9, 2015, at 10:09 AM, Larry Hunsicker <lhuns...@mchsi.com> wrote:
>
> Greetings, STAN-list users. Another question from a STAN (and Bayesian) neophyte concerning mixed model analysis.
>
> I have been trying to understand better the differences in assumptions and conclusions from fully Bayesian analyses (STAN) and the frequentist analyses (with which I am more familiar) of mixed models. To be explicit, I have been looking at analyses of the BLOCKER data as in the BUGS, and now the STAN models, and using glmer from the R package nlme4. The estimates of "d" in the two sets of analyses are virtually identical -- -0.25 using STAN, and -0.26 using glmer. But the estimates of the SE of d are quite different, about 0.12 for STAN (taking the sqrt of the estimated sigmasqu_delta), but only 0.05 using glmer. In this case, the conclusions from the two approaches are the same. But it appears that STAN is estimating a relatively much larger SEM than glmer, and it is easy to see that with a d a bit closer to 0, the conclusions might be quite different.
>
> I am using RSTAN 2.6.0, R version 3.2.0, RStudio 0.99.441, and lme4 1.1-7, on a Windows 8.1 machine with all updates. I used blocker.stan blocker.dat, and blocker.init, unchanged from the files in the STAN examples. For the glmer analysis, I restructured blocker.dat so that there was one case per record, with three columns: death (1 vs 0), study number, and treatment group. My glmer function was:
>
> m1 <- glmer(death ~ group + (1 | study), data = bldat2, family = binomial)
> //group is treated vs. control (1 vs 0); study is the number of the specific cited study.
>
> I also checked a variant of blocker.stan, using normal distributions instead of t distributions (4 df) for d. The results of this modified blocker.stan were essentially identical to the results from the original blocker.stan.
>
> My set of questions includes:
> a) What are the differences in assumptions that lead to the different estimates of the SE of d? Is glmer missing a source of variance that STAN is finding, and if so, where?

nlme does max marginal likelihood --- marginalizing out low-level
parameters, then optimizing the hierarchical parameters, then plugging
in point estimates for the hierarchical parameters in estimating the
low-level params. It's not a well-founded process probabilistically.

Andrew or Ben or Rob will know the answer for sure, but my guess is
that the uncertainty in the hierarchical parameters is being missed in
the nlme analysis.

> Alternatively, is the difference due to error introduced by approximations introduced in the calculations?

No. Neither calculation is "approximate" in the meaningful sense.
nlme's whole approach is an approximation, but the calculation is
not (although it uses an iterative algorithm). MCMC is only approximate
in that its error goes down with sqrt of number of independent draws.

> b) Given that the model results are somewhat different, which is more likely to represent "reality"?

That's a loaded question asking whether to do classical point
estimates or full Bayes.

Full Bayes is the way to go. It captures all the uncertainty
in your estimates and lets you use them downstream for inference.

We think of max marginal likelihood as an approximation procedure.

> c) My understanding is that the estimates of SE(d) in STAN come directly from sampling the posterior estimates of d, but that with glmer, the SE is estimated as the partial slope of the inverted information matrix, which, as I understand it, is only asymptotically equivalent to a normal distribution.

There is that, too.

> Does this imply that the STAN estimate is more accurate than the glmer estimate?

They're computing different things; see above.

> If I am correct in my understanding of this, the implication is that one should do most of one's mixed model analysis using STAN rather than glmer. Am I way off?

We'd agree, of course.

> d) If someone could point me to a useful reference that will answer my questions above, I would be deeply grateful to you. (I am slowly making my way through Gelman's book, "Bayesian Data Analysis.")

I'd recommend the lme4 doc:

http://cran.r-project.org/web/packages/lme4/vignettes/lmer.pdf

to actually see what's going on.

- Bob

[Larry Hunsicker]

Thanks, Bob, for a very helpful answer.  It is both encouraging and discouraging:  the former in that I have a better sense of where I have to go, and discouraging in that I now have to question the accuracy of my past analyses using glmer.
I would take mild issue with only one of your answers.  "They're computing different things."  I would offer that they are both intended to provide an answer to the same question, in this case "What is the size of the effect of using beta blockers, and how accurately do we know this size?"  (I am a physician, after all.) They are computing the answer to the same question using different methods.  Both are approximations, but one (STAN) appears to be a giving a somewhat more accurate approximation, and the accuracy of the approximation is under my control (by how many chains and iterations I request).  
The history of statistics is a record of approximations acceptable because they were tractable.  But with today's computational resources, there is no good excuse for using a less good approximation simply because it runs a bit faster.
I'll spend my free time (after my day job) reading more of Gelman's book rather than trying to digest the innards of glmer.
Again, thanks.
Larry Hunsicker

[Andrew Gelman]

> On Jun 9, 2015, at 11:51 AM, Bob Carpenter <ca...@alias-i.com> wrote:
>
>
>> On Jun 9, 2015, at 10:09 AM, Larry Hunsicker <lhuns...@mchsi.com> wrote:
>>
>
>> Alternatively, is the difference due to error introduced by approximations introduced in the calculations?
>
> No. Neither calculation is "approximate" in the meaningful sense.
> nlme's whole approach is an approximation, but the calculation is
> not (although it uses an iterative algorithm). MCMC is only approximate
> in that its error goes down with sqrt of number of independent draws.
>

nmle is approximate in that, for nonlinear models, it is not actually maximizing the intergrated likelihood, it is maximizing a normal approximation to the integrated likelihood. According to Sophia, this can create problems sometimes.

[Bob Carpenter]

> On Jun 9, 2015, at 1:48 PM, Larry Hunsicker <lhuns...@mchsi.com> wrote:
>
> Thanks, Bob, for a very helpful answer. It is both encouraging and discouraging: the former in that I have a better sense of where I have to go, and discouraging in that I now have to question the accuracy of my past analyses using glmer.
> I would take mild issue with only one of your answers. "They're computing different things." I would offer that they are both intended to provide an answer to the same question, in this case "What is the size of the effect of using beta blockers, and how accurately do we know this size?" (I am a physician, after all.)

The fact that you said "and how accurately" means you're at least
part statistician!

> They are computing the answer to the same question using different methods.

I think I understand what you mean --- it's like saying Newtonian physics
and Einsteinian physics are trying to compute the answer to the same
question using different methods. They are different mathematical
systems that provide different answers if computed correctly, but they're trying
to get at the same physical phenomenon.

Once you take that view, the world opens up into a huge array of possible
models if you conflate every possible approach to answering the effect of
using beta blockers. It's almost certainly not linear, for instance, so
all of these approximations to reality using (generalized) linear models are
highly questionable.

What you want to do is get into testing the predictions of your model on
held-out data (or using cross validation). I think that'll give you the means
you need to properly compare models.

> Both are approximations, but one (STAN) appears to be a giving a somewhat more accurate approximation, and the accuracy of the approximation is under my control (by how many chains and iterations I request).

Stan computes posterior expectations up to MCMC precision. All computation
on a computer is approximate in that we have double-precision arithmetic,
not real numbers. The computer approximation is better for some operations (like
multiplication) than it is for others (like subtraction).

> The history of statistics is a record of approximations acceptable because they were tractable. But with today's computational resources, there is no good excuse for using a less good approximation simply because it runs a bit faster.

Our point exactly! Of course, we keep taking on bigger and harder
problems where we can't get the solutions we want using MCMC in
less than a week, which becomes problematic for doing science
again, and forces us to think about approximations like variational
inference and max marginal likelihood.

> I'll spend my free time (after my day job) reading more of Gelman's book rather than trying to digest the innards of glmer.

I'd strongly suggest the Gelman and Hill book for this kind of thing.
It's what convinced me to quit my day job and join the circus :-)

- Bob