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