5
$\begingroup$

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

Improve this question
$\endgroup$
4
  • $\begingroup$ Why you need to treat the effect of group as random, instead of fixed? $\endgroup$ – user158565 Nov 19 '18 at 1:56
  • $\begingroup$ 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) $\endgroup$ – Bakaburg Nov 19 '18 at 12:11
  • $\begingroup$ 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. $\endgroup$ – user158565 Nov 19 '18 at 23:03
  • 1
    $\begingroup$ 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 $\endgroup$ – Bakaburg Nov 20 '18 at 14:19
1
$\begingroup$

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.

Improve this answer
$\endgroup$
3
  • $\begingroup$ 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. $\endgroup$ – Bakaburg Nov 20 '18 at 21:03
  • 1
    $\begingroup$ 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. $\endgroup$ – Dimitris Rizopoulos Nov 21 '18 at 5:29
  • $\begingroup$ 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 $\endgroup$ – Bakaburg Nov 21 '18 at 22:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.