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

## ----family-------------------------------------------------------------------
nb_size <- 8

negative_binomial_rate <- metaGLMM_family(
  name = "negative-binomial-rate",
  link = "log",
  loglik = function(y, eta, vi, ni) {
    count <- round(y * ni)
    dnbinom(
      x = count,
      size = nb_size,
      mu = exp(eta) * ni,
      log = TRUE
    )
  },
  validate = function(y, vi, ni) {
    count <- y * ni
    if (any(y < 0) || any(abs(count - round(count)) > 1e-7)) {
      stop("y * ni must define non-negative integer counts.")
    }
    if (any(!is.finite(vi)) || any(vi <= 0)) {
      stop("vi must contain finite positive display variances.")
    }
    TRUE
  },
  initialize = function(y, X, vi, ni) {
    value <- rep(0, ncol(X))
    intercept <- match("(Intercept)", colnames(X))
    if (!is.na(intercept)) {
      value[intercept] <- log(sum(y * ni) / sum(ni))
    }
    value
  }
)
negative_binomial_rate

## ----data-fit-----------------------------------------------------------------
nb_dat <- data.frame(
  study = paste0("Study ", 1:8),
  events = c(8, 30, 5, 55, 18, 90, 9, 70),
  exposure = c(20, 35, 18, 42, 30, 50, 25, 45)
)
nb_dat$rate <- nb_dat$events / nb_dat$exposure
nb_dat$vi <- 1 / pmax(0.5, nb_dat$events) + 1 / nb_size
nb_dat

qmc_points <- qnorm((seq_len(512) - 0.5) / 512)

nb_fit <- metaGLMM(
  rate ~ 1,
  data = nb_dat,
  vi = nb_dat$vi,
  ni = nb_dat$exposure,
  tau2 = NA,
  tau2_var = TRUE,
  start_tau2 = 0.3,
  family = negative_binomial_rate,
  rstdnorm = qmc_points,
  fast = FALSE
)

summary(nb_fit)
exp(coef(nb_fit))

nb_ci <- list(
  Wald = confint(nb_fit, method = "wald"),
  Profile = confint(nb_fit, method = "profile"),
  SBC = confint(nb_fit, method = "SBC")
)
do.call(rbind, nb_ci)
stopifnot(is.finite(nb_fit$tau), nb_fit$tau > 0)

## ----forest-------------------------------------------------------------------
forest(
  nb_fit,
  labels = nb_dat$study,
  type = "exp",
  xlab = "Negative-binomial rate",
  ci_methods = c("Wald", "profile", "SBC")
)

## ----expected-error, eval=FALSE-----------------------------------------------
# try(
#   metaGLMM(
#     rate ~ 1,
#     data = nb_dat,
#     vi = nb_dat$vi,
#     ni = nb_dat$exposure,
#     tau2 = NA,
#     tau2_var = TRUE,
#     family = negative_binomial_rate,
#     fast = TRUE
#   )
# )