---
title: "Custom families"
author: "metaGLMM authors"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Custom families}
  %\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)
```

## The family contract

`metaGLMM_family()` defines an R-level conditional likelihood.  Supply a
name, a link accepted by `stats::make.link()`, and a vectorized callback
with arguments `(y, eta, vi, ni)`.  The callback returns one finite
contribution for every value of `eta`.  Optional `validate` and `initialize`
hooks can add family-specific checks or deterministic starting values.

The example below defines a negative-binomial likelihood for overdispersed
study counts.  Each response is an observed rate `y`, `ni` is its exposure,
and `y * ni` recovers the count.  The conditional count mean is
`exp(eta) * ni`.  A fixed negative-binomial size is captured by the callback;
the between-study log-rate variance remains the `tau2` estimated by
`metaGLMM()`.

```{r 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
```

The `loglik` callback is vectorized in `eta`: during QMC integration a single
study count is evaluated at all random-effect draws.  The validation function
checks the aggregate count contract before optimization, and the initialization
function supplies a deterministic pooled log-rate start.

## Fit an overdispersed rate model with QMC

```{r 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)
```

The exponentiated intercept is the pooled conditional rate.  The fitted
`tau2` describes additional between-study variation on the log-rate scale,
while `nb_size` controls negative-binomial variation within the conditional
count model.  In this interface `nb_size` is prespecified rather than estimated;
record its choice and consider sensitivity analyses when it is not externally
known.

Custom families support Wald, profile, and SBC intervals.  The simplified
Bartlett correction uses `vi` and the heterogeneity estimate under each profile
constraint; for a custom family, supply `vi` as an appropriate study-level
sampling variance on the link scale.

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

## Validation and limitations

The callback must return a finite numeric vector of the same length as `eta`.
Check response domains and variance conventions before fitting; family
validation errors are reported before numerical optimization.

Custom families are supported through `fast = FALSE` in this release.  A
request for `fast = TRUE` with a custom family is rejected because compiled
quadrature for arbitrary R callbacks is not part of the interface.

```{r 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
  )
)
```

For production analyses, keep the callback definition, the QMC sequence, and
the fitted object's integration settings with the analysis record.