# Build confidence intervals for random effects in intercept only models

I'm building a Poisson model for a rate between an outcome and an offset (to be precise between the observed and predicted values of a phenomenon). I want to get this rate + confidence intervals (CI) for each observation group. This is the model:

mod <- glmer(observed ~ 1 + offset(log(preds)) + (1 | group),
family = poisson(), DF)


What should I do to get the rate (not exponentiated) for each group? I know about the dotplot() function to plot the conditional means and confidence intervals of the groups, but I can't understand how it works and I would like to have the raw estimates to do my own plots.

I tried extracting the mean and standard error for each group using ranef() and than approximate a CI:

est = ranef(mod, condVar = TRUE)$$group se = attr(ranef(mod, condVar = TRUE)$$group, 'postVar')[,,]

upr = est + 1.96 * se
lwr = est - 1.96 * se


A second approach I used was through empirical Bayes simulation:

s <- sim(mod, n.sims = 1000)

lapply(group_names, function(g_name) {
q <- quantile(s@ranef$$h_code[,g_name,], c(0.025, .0975)) data.frame(group = code, lwr = q[1], upr = q[2], est = ranef(mod)$$group[g_name,])
}) %>% bind_rows


The CI of the first method are way shorter than in the second case and both look different than those produced by the dotplot() function. Furthermore, I have not clear whether I should use ranef() estimates as they are or join them somehow with fixef() estimates.

NOTE: I know how to solve this problem using full Bayesian estimation eg. using rstanarm. I was wondering if there was a faster approximate solution using the lme4 framework

• Why you need to treat the effect of group as random, instead of fixed? Nov 19, 2018 at 1:56
• Because there are 50+ groups and I'm interested in the ratio at the group level, not in the effect of the specific group on the global ratio. The alternative would be to run 50+ regressions but then I would lose the regularization brought by the exchange of information between groups (one of the characteristics of hierarchical models) Nov 19, 2018 at 12:11
• Unclear what you want. 1. rate or ratio? 2. "I'm interested in the ratio at the group level" = you want group level thing, "not in the effect of the specific group on the global ratio." = you do not want group level thing. Nov 19, 2018 at 23:03
• Uhm, I don't understand your doubts... Imagine that I want to estimate the healthcare-associated infection rate (actually is my case) for a group of hospitals. I could just do the ratio $n.\ infections/n.\ patients$ for each hospital (adjusting for the predicted risk related to patients' clinical gravity with an offset). Because smaller hospitals may have artificially high or low estimates due to sampling noise, one uses the shrinkage effect brought by mixed models to have more stable and precise estimates at the cost of a certain amount of bias. Here's a ref. bit.ly/2FyCZlS Nov 20, 2018 at 14:19

When you are interested in predictions conditional on the random effects, to my view it is easier to work with the hiearachical formulation of the mixed model that has a intrinsically Bayesian flavor. In particular, in your specific case, your are interested in the mean of the Poisson model conditional on the random effects, i.e., $$\mu_i = \exp(\beta + b_i),$$ with $$i$$ denoting the group, $$\beta$$ the fixed effect intercept, and $$b_i$$ the random intercept. You can derive a confidence interval for $$\mu_i$$ by using the following simulation scheme:

• Step I: Simulate a value $$\theta^*$$ from the approximate posterior distribution $$\mathcal N(\hat\theta, \hat\Sigma)$$, where $$\hat\theta$$ denotes the maximum likelihood estimatates for $$\beta$$ and $$\sigma_b$$ with $$\sigma_b$$ denoting the standard deviation of the random intercepts term $$b_i$$, and $$\hat\Sigma$$ the variance of $$\hat\theta$$.
• Step II: Simulate a value $$b_i^*$$ from the posterior distribution of the random effects $$[b_i \mid y_i, \theta^*]$$, where $$y_i$$ denotes the outcome data for group $$i$$ (note we condition on $$\theta^*$$ from the previous step).
• Step III: Calculate $$\mu_i^* = \exp(\beta^* + b_i^*)$$.

Step I accounts for the sampling variability of the maximum likelihood estimates, and Step II for the variability in the random effects.

Repeating Steps I-III $$L$$ times, you obtain a Monte Carlo sample for $$\mu_i$$ based on which you could obtain a 95% CI using the 2.5% and 97.5% percentile.

This procedure is implemented in the predict() method for mixed models fitted using the GLMMadaptive package. For an example, check the vignette Methods for MixMod Objects.

• Thanks Dimitris. I do use rstanarm when I want to go full Bayesian. But in this case I wanted to move in the lme4 framework, if just for speed reasons. Nov 20, 2018 at 21:03
• Note that in the procedure described above the model is estimated under maximum likelihood (i.e., as in lme4). Based on this fitted model, the confidence interval for the group-specific predictions is calculated on Bayesian grounds. When you want this type of predictions and even if you fitted under maximum likelihood, tackling the problem under the Bayesian approach is standard. E.g., function ranef() gives you the empirical Bayes estimates of the random effects, derived as the modes of the posterior distribution of the random effects given the observed data. Nov 21, 2018 at 5:29
• yep, I know it's an approximation of Bayes. I just would like to know which approximate method is better without resolving to full bayesian Nov 21, 2018 at 22:22