---
title: "Poisson rate outcomes"
author: "metaGLMM authors"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Poisson rate outcomes}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>",
                      fig.width = 6, fig.height = 4)
set.seed(20260821)
library(metaGLMM)
```

## Rates from aggregate counts

For a Poisson analysis, use a rate as the response and the corresponding
exposure as `ni`.  A small positive variance approximation is useful for
the aggregate-data representation, especially when an observed count is
zero.

```{r data}
poisson_dat <- data.frame(
  study = paste0("Study ", 1:12),
  sex = factor(rep(c("Female", "Male"), 6)),
  year = c(0, 0, 2, 2, 6, 6, 6, 6, 7, 7, 8, 8),
  events = c(13, 1, 16, 3, 38, 3, 13, 4, 5, 10, 11, 23),
  exposure = c(9.04, 8.46, 23.75, 23.75, 18.47, 18.47,
               9.38, 9.19, 2.27, 3.34, 2.72, 26.00)
)
poisson_dat$rate <- poisson_dat$events / poisson_dat$exposure
poisson_dat$vi <- 1 / pmax(0.5, poisson_dat$events)
poisson_dat
```

## Basic single-group meta-analysis and forest plot

An intercept-only model estimates the overall log rate while allowing each
study to have its own random effect.  This is the direct single-group analogue
of a conventional random-effects meta-analysis.

```{r basic-fit}
poisson_mean_fit <- metaGLMM(
  rate ~ 1,
  data = poisson_dat,
  vi = poisson_dat$vi,
  ni = poisson_dat$exposure,
  tau2 = NA,
  family = poisson(link = "log"),
  tau2_var = TRUE,
  fast = TRUE,
  ghq_Q = 40L
)

summary(poisson_mean_fit)
coef(poisson_mean_fit)
confint(poisson_mean_fit, method = "wald")
stopifnot(is.finite(poisson_mean_fit$tau), poisson_mean_fit$tau > 0)
```

Each forest row is one observed study rate.  `type = "exp"` converts the
study estimates, pooled intervals, and prediction interval from the log-rate
scale to the rate scale.

```{r basic-forest}
poisson_studies <- as_metafor_data(
  poisson_mean_fit, labels = poisson_dat$study
)
head(poisson_studies)

forest(
  poisson_mean_fit,
  labels = poisson_dat$study,
  type = "exp",
  xlab = "Rate",
  ci_methods = c("Wald", "profile", "SBC")
)
```

## Meta-regression with moderators

The formula can include factors, continuous moderators, and interactions.  The
coefficient names are the names generated by `model.matrix()`.

```{r fit}
poisson_fit <- metaGLMM(
  rate ~ sex + year + sex:year,
  data = poisson_dat,
  vi = poisson_dat$vi,
  ni = poisson_dat$exposure,
  tau2 = NA,
  family = poisson(link = "log"),
  tau2_var = TRUE,
  fast = TRUE,
  ghq_Q = 40L
)

summary(poisson_fit)
coef(poisson_fit)
exp(coef(poisson_fit))
```

The exponentiated coefficients are rate ratios relative to the factor
reference level and at `year = 0`.  Center `year` before fitting when another
reference year is more meaningful.  Because this model has several fixed
effects, a pooled forest summary would require an explicitly chosen coefficient
or contrast; the basic intercept-only fit above is the clearer study-rate
display.

## Offsets

An offset is part of the ordinary formula interface.  For example, when the
response is defined on a count scale, the model matrix contains the exposure
offset as follows:

```{r offset, eval=FALSE}
count_formula <- events ~ sex + year + offset(log(exposure))
model.matrix(count_formula, data = poisson_dat)
```

Choose either a rate response with exposure supplied through `ni` or a count
response with a formula offset according to the sampling convention of the
analysis; keep the response and the exposure definition consistent.