How to use add_ci in presence of observation-level random effect

[tomwenseleers]

I have a binomial GLMM of the form fit = glmer(cbind(strain, TOTAL-strain) ~ 0+(1|REGION/obs) +
COUNTRY*SAMPLE_DATE, family=binomial(logit), data=data) Now I would like to produce confidence intervals, taking into account my fixed factors COUNTRY & SAMPLE_DATE and random factor intercept REGION, but ignoring the observation-level random effect (put in to take into account overdispersion). I would also like to be able to extrapolate my model beyond the original sample dates in my dataset. Is there any way to achieve something like this? Or is the only option for me to set includeRanef = FALSE and then ignore the random intercepts per region? (which would not be ideal for me)

[jthaman]

1. About ignoring random effects, the two options are
   • Switch from a GLMM model to a GEE model, where the model coefficients have a more direct population-level interpretation. You might look at geepack R package.
   • Use a package that is more clever than ciTools about this. GLMMadaptive is a reputable package that can successfully marginalize out the random effects from a GLMM using the function marginal_coefs(). I don't know about confidence intervals, and I'm not sure if the marginalization works for the nested structure you have.

I can't endorse one strategy over the other, but the GEE approach is maybe better established in the literature.

2. On extrapolation, I'm not sure. I might try to use a spline term for the time variable and try to only make very near-term extrapolations. Extrapolation is a pretty uncommon goal for GLMMs.

[tomwenseleers]

Thanks for the advice! I also tried the emmeans package by the way. I might give geepack a go, thanks for that!

[jthaman]

Did the emmeans approach bear any fruit? I recall that the package has an allergic reaction to random effects (in the reference grid).

[tomwenseleers]

If I used extra options with bias.adjust = TRUE, sigma = total.SD to correct for bias in the marginal means due to the presence of random effects with total.SD = sqrt(sum(sapply(as.data.frame(VarCorr(fit))$sdcor, function (x) x^2))) the emmeans plots in the end looked OK. That's an option that's been recently introduced in emmeans, https://github.com/rvlenth/emmeans/issues/231. And switching to a model glmer(cbind(strain, TOTAL-strain) ~ 0+(1|obs) + REGION+COUNTRY*SAMPLE_DATE, family=binomial(logit), data=data) was also an option... In my case it's mainly the slopes in function of SAMPLE_DATE by COUNTRY that I am interested in, and that I can also look at using emtrends(). So I think in the end this will do...