Skip to content
Prev 3951 / 5636 Next

[R-meta] Question regarding metarate calling the Poisson model for meta-analysis

Dear Roel,

I can't quite follow exactly what you are tying to do or have done.

"We wanted to use the random effects, to account for different cumulative time periods of follow-up." This sounds like you want some kind of model that includes multiple random effects, at multiple levels (e.g., a random effect at the estimate level (for overdispersion) and a random effect at some higher level to account for clustering / a multilevel structure). rma.glmm() is not set up for that. In this case, you can just use glmer() from lme4 directly.

"The pooled results however, showed confidence limits that differed significantly when we compared these to exact Poisson 95% CIs that we calculated separately." Are you talking about CIs for the individual rates or for the pooled result? If you already have an exact CI for the pooled result, then why bother with metarate(), rma.glmm(), or glmer()? But if you are talking about the individual rates: The CIs are constructed are typically Wald-type CIs (possibly on a transformed scale and then back-transformed). So these can differ from exact CIs computed directly based on a Poisson distribution. For example:

summary(escalc(measure="IR", xi=10, ti=500))[c("ci.lb", "ci.ub")]
summary(escalc(measure="IRLN", xi=10, ti=500), transf=exp)[c("ci.lb", "ci.ub")]
summary(escalc(measure="IRS", xi=10, ti=500), transf=\(x) x^2)[c("ci.lb", "ci.ub")]

And comparing those against the exact CI:

round(c(poisson.test(10, 500)$conf.int), 4)

Note that these CIs are typically just used in a forest plot for visualization of the individual studies/estimates. They are not used in the actual meta-analysis. So even though the CIs may not be 'exact', they do convey approximately the uncertainty in the different estimates and how these uncertainties differ from each other.

If you prefer to show exact CIs in a forest plot, you can do that as follows:

dat <- dat.nielweise2008
dat <- escalc(measure="IRLN", xi=x1i, ti=t1i, data=dat, slab=paste0(authors, ", ", year))

res <- rma(yi, vi, data=dat)

par(mfrow=c(1,2))

forest(res, refline=NA, transf=exp, digits=3, psize=1, header=TRUE)

# compute the 'exact' CIs
cis <- tapply(dat, dat$study, FUN = \(x) poisson.test(x$x1i, x$t1i)$conf.int)
cis <- do.call(rbind, cis)

# set up the forest plot using the exact CIs for the individual rates
dat <- escalc(measure="IR", xi=x1i, ti=t1i, data=dat, slab=paste0(authors, ", ", year))
forest(dat$yi, ci.lb=cis[,1], ci.ub=cis[,2], refline=NA, digits=3, psize=1, header=TRUE, ylim=c(-1.5,res$k+3))
abline(h=0)

# add the estimate from the model
pred <- predict(res, transf=exp)
addpoly(pred$pred, ci.lb=pred$ci.lb, ci.ub=pred$ci.ub, row=-1)

I put the two forest plots side-by-side to illustrate how they differ. Note that one cannot construct the forest plot on some transformed scale such that the CIs will be symmetric when using the exact CIs, so you will end up with these skewed-looking CIs, as opposed to:

forest(res, refline=NA, atransf=exp, digits=3, psize=1, header=TRUE)

Best,
Wolfgang