Custom families

metaGLMM authors

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().

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
#> 
#> Family: negative-binomial-rate 
#> Link function: log

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

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
#>     study events exposure      rate        vi
#> 1 Study 1      8       20 0.4000000 0.2500000
#> 2 Study 2     30       35 0.8571429 0.1583333
#> 3 Study 3      5       18 0.2777778 0.3250000
#> 4 Study 4     55       42 1.3095238 0.1431818
#> 5 Study 5     18       30 0.6000000 0.1805556
#> 6 Study 6     90       50 1.8000000 0.1361111
#> 7 Study 7      9       25 0.3600000 0.2361111
#> 8 Study 8     70       45 1.5555556 0.1392857

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)
#> Aggregate-data generalized linear mixed-effects meta-analysis
#> 
#> Call:
#> metaGLMM(formula = rate ~ 1, data = nb_dat, vi = nb_dat$vi, ni = nb_dat$exposure, 
#>     tau2 = NA, family = negative_binomial_rate, tau2_var = TRUE, 
#>     rstdnorm = qmc_points, fast = FALSE, start_tau2 = 0.3)
#> 
#> Family:negative-binomial-rate(log)
#> Random-effect structure: row_intercept 
#> Integration: qmc 
#> 
#> Fixed effects:
#>             Estimate Std. Error z value Pr(>|z|)
#> (Intercept)  -0.2089     0.2377  -0.879     0.38
#> 
#> Heterogeneity:
#>   tau^2: 0.2522 
#>   tau:   0.5022 
#> 
#> Plug-in prediction interval for a future underlying true effect:
#>      fit  lower  upper
#>  -0.2089 -1.298 0.8802
#> 
#> Fit diagnostics:
#>   Optimizer: L-BFGS-B 
#>   Convergence code: 0 
#>   Boundary tau^2: no 
#>   Positive-definite Hessian: yes
exp(coef(nb_fit))
#> (Intercept) 
#>   0.8114576

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)
#>                  lower     upper
#> (Intercept) -0.6748997 0.2570535
#> (Intercept) -0.7702812 0.2864727
#> (Intercept) -0.8377501 0.3449623
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.

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.

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.