[R-sig-ME] Getting the right uncertainties when fitting multilevel models

[Poe, John]

I'm just wondering if anyone who knows more about lme4 than i do has any
comments on the Gelman post today. He's pretty vague and brief in his
reply.



Cesare Aloisi writes:

I am writing you regarding something I recently stumbled upon in your book
Data Analysis Using Regression and Multilevel/Hierarchical Models which
confused me, in hopes you could help me understand it. This book has been
my reference guide for many years now, and I am extremely grateful for
everything I learnt from you.

On page 261, a 95% confidence interval for the intercept in a certain group
(County 26) is calculated using only the standard error of the “random
effect” (the county-level error). The string is as follows:

coef(M1)$county[26,1] + c(-2,2)*se.ranef(M1)$county[26]

My understanding is that, since the group-level prediction (call it y.hat_j
= coef(M1)$county[26,1]) is a linear combination of a global average and a
group-level deviation from the average (y.hat_j = beta_0 + eta_j), then the
variance of y.hat_j should be the sum of the covariances of beta_0 and
eta_j, not just the variance of eta_j, as the code on page 261 seems to
imply. In other words:

Var(y.hat_j) = Var(beta_0) + Var(eta_j) + 2Cov(beta_0, eta_j)

Admittedly, lme4 does not provide an estimate for the last term, the
covariance between “fixed” and “random” effects. Was the code used in the
book to simplify the calculations, or was there some deeper reason to it
that I failed to grasp?

My (Gelman's) reply: The short answer is that it’s difficult to get this
correct in lmer but very easy when using stan_lmer() in the rstanarm
package. That’s what I recommend, and that’s what we’ll be doing in the 2nd
edition of our book


http://andrewgelman.com/2017/09/23/getting-right-uncertainties-fitting-multilevel-models/

	[[alternative HTML version deleted]]