## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>",
                      fig.width = 6, fig.height = 4)
set.seed(20260821)
library(metaGLMM)

## ----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-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)

## ----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")
)

## ----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))

## ----offset, eval=FALSE-------------------------------------------------------
# count_formula <- events ~ sex + year + offset(log(exposure))
# model.matrix(count_formula, data = poisson_dat)

