metaGLMM

metaGLMM fits likelihood-based random-effects meta-analyses as generalized linear mixed-effects models using aggregate study data. The same interface supports covariate adjustment, non-Gaussian responses, and a small set of explicit random-effect structures without requiring individual participant data.

The package is intended for reproducible analyses in base R. A fitted model is an S3 object with familiar methods such as summary(), coef(), vcov(), confint(), predict(), and plot().

Installation

Install the development version from GitHub with:

install.packages("remotes")
remotes::install_github("keisuke-hanada/metaGLMM", build_vignettes = TRUE)

The package imports only the numerical tools needed for fitting. The examples and vignettes use base R data manipulation and graphics.

A minimal analysis

The response is an aggregate study estimate, vi is its sampling variance, and ni is the relevant study size or exposure. The formula defines the fixed-effects model matrix.

library(metaGLMM)

dat <- data.frame(
  study = paste0("Study ", 1:6),
  estimate = c(-0.30, -0.08, 0.05, 0.18, 0.27, 0.41),
  moderator = c(0, 0, 1, 1, 0, 1),
  vi = c(0.08, 0.11, 0.09, 0.13, 0.10, 0.12),
  ni = c(90, 80, 110, 75, 95, 85)
)

fit <- metaGLMM(
  estimate ~ moderator,
  data = dat,
  vi = dat$vi,
  ni = dat$ni,
  tau2 = NA,
  family = gaussian(link = "identity"),
  tau2_var = TRUE,
  fast = TRUE
)

summary(fit)
coef(fit)             # fixed effects only
vcov(fit)             # fixed-effect covariance matrix
fit$tau2              # between-study variance
confint(fit, method = "wald")
predict(fit, newdata = data.frame(moderator = c(0, 1)),
        interval = "prediction")

The vignettes provide small, self-contained examples:

• Getting started introduces the data contract and the main methods.
• Binary outcomes fits grouped binomial summaries and plots within-study treatment contrasts.
• Poisson outcomes models single-group rates, draws a study-rate forest plot, and then introduces covariates.
• Gaussian outcomes covers weighted continuous study estimates, random effects, and a forest plot.
• Likelihood inference compares Wald, profile, and SBC intervals.
• Prediction and forest plots demonstrates prediction scales, intervals, and safe study-level plotting.
• Custom families shows the R-level family extension using QMC integration.

Supported families

The built-in stats families are gaussian(), binomial(), poisson(), and Gamma(). Their link functions are passed in the usual R style, for example binomial(link = "logit") or Gamma(link = "log"). The fixed-effect coefficient names always follow model.matrix(), including interactions, factor contrasts, and formulas without an intercept.

For binary data, y is a proportion and ni is the number of trials. For Poisson data, y is a rate and ni is the exposure. For Gaussian and Gamma analyses, provide the aggregate response, sampling variance, and the study size used by the likelihood. Values at a binomial or Poisson boundary can be represented with the legacy vi = Inf convention; finite variances are needed when an automatic study-level display is requested.

Custom families

Use metaGLMM_family() when a response requires an R-level conditional log-likelihood. The callback is vectorized in eta and returns one log-likelihood contribution per row.

nb_size <- 8
nb_rate_family <- metaGLMM_family(
  name = "negative-binomial-rate",
  link = "log",
  loglik = function(y, eta, vi, ni) {
    dnbinom(round(y * ni), size = nb_size,
            mu = exp(eta) * ni, log = TRUE)
  }
)

nb_dat <- data.frame(events = c(8, 30, 5, 55),
                     exposure = c(20, 35, 18, 42))
nb_dat$rate <- nb_dat$events / nb_dat$exposure
nb_dat$vi <- 1 / nb_dat$events + 1 / nb_size

qmc <- qnorm((seq_len(256) - 0.5) / 256)
custom_fit <- metaGLMM(rate ~ 1, nb_dat, vi = nb_dat$vi,
                       ni = nb_dat$exposure, tau2 = NA,
                       family = nb_rate_family, tau2_var = TRUE,
                       rstdnorm = qmc, fast = FALSE)

Custom families use the QMC path (fast = FALSE) in the initial release. Compiled quadrature for custom callbacks is intentionally not part of this interface. The custom-family vignette provides validation, initialization, and a larger negative-binomial example.

Inference and prediction

The official confint() methods are "wald", "profile", and "SBC" (case-insensitive). They target fixed effects; tau2 and tau remain available as fitted-object components and in summary().

predict() can return the link or response scale. A confidence interval uses the fixed-effect covariance, while a prediction interval additionally uses random_scale^2 * tau2. Observation-level sampling variance is not added to the prediction interval. For log and logit models, type = "exp" provides a convenient multiplicative response scale.

Forest plots and data exchange

forest(fit) and plot(fit, type = "forest") use base graphics and show study estimates, their sampling intervals, selected fixed-effect confidence intervals as diamonds, and a plug-in prediction interval when available. The default uses SBC when supported and Wald otherwise. Request several rows with ci_methods = c("wald", "profile", "SBC"). Automatic study extraction uses as_metafor_data() and remains deliberately conservative. For complex arm-level data, pass estimate, vi, and labels explicitly.

as_metafor_data(fit) returns an ordinary data.frame with yi, vi, and slab columns. It does not create a class from another package or add a dependency on that package.

Citation

Please cite metaGLMM and the methodological work that motivates the likelihood construction:

Hanada, K. and Sugimoto, T. (2026). Random-effects meta-analysis via generalized linear mixed models: A Bartlett-corrected approach for few studies. doi:10.1093/biomtc/ujag148.

@article{HanadaSugimoto2026,
  author  = {Hanada, Keisuke and Sugimoto, Tomoyuki},
  title   = {Random-effects meta-analysis via generalized linear mixed models:
             A Bartlett-corrected approach for few studies},
  year    = {2026},
  doi     = {10.1093/biomtc/ujag148}
}