A simulation to study the efficiency gain of random effects vs. semiparametric Poisson.

cluster_effects_sim.R

library(tidyverse)
library(sandwich)
library(lmtest)
library(lme4)

patients <- 3e3

seq_len(patients) %>%
    map_df(function(id) {
        rnorm(1, 0, 1) -> patient_effect
        N <- 2
        rbinom(N, 1, 0.5) -> trt
        tibble(id,
               trt,
               y = rpois(N, exp(1 + 0.5 * trt +
                                    patient_effect)))
    }) -> d

glm(y ~ trt,
    data = d,
    family = "poisson") -> fit
coeftest(fit, vcovHAC(fit))

glmer(y ~ trt + (1 | id),
      data = d,
      family = "poisson") %>%
    summary()