Prediction and forest plots

metaGLMM authors

A fitted model

dat <- data.frame(
  study = paste0("Study ", 1:8),
  estimate = c(-0.80, -0.30, 0.10, 0.70, 1.10, -0.50, 0.40, 1.30),
  vi = c(0.04, 0.05, 0.03, 0.06, 0.04, 0.05, 0.03, 0.05),
  ni = c(100, 90, 120, 80, 110, 95, 130, 85)
)

fit <- metaGLMM(
  estimate ~ 1,
  data = dat,
  vi = dat$vi,
  ni = dat$ni,
  tau2 = NA,
  tau2_var = TRUE,
  family = gaussian(link = "identity"),
  fast = TRUE
)
stopifnot(is.finite(fit$tau), fit$tau > 0)

Prediction for a future study

predict.metaGLMM() accepts either newdata or a numeric contrast matrix, but not both. Neither is required for this intercept-only example: the confidence interval targets the pooled mean, while the prediction interval targets the underlying true effect in one future study.

pred_conf <- predict(fit, type = "link", interval = "confidence")
pred_pred <- predict(fit, type = "link", interval = "prediction",
                     random_scale = 1)
pred_response <- predict(fit, type = "response", interval = "confidence")

pred_conf
#>              fit      lower     upper
#> overall 0.248429 -0.2419502 0.7388083
pred_pred
#>              fit     lower    upper
#> overall 0.248429 -1.164695 1.661553
pred_response
#>              fit      lower     upper
#> overall 0.248429 -0.2419502 0.7388083
stopifnot(diff(as.numeric(pred_pred[1, c("lower", "upper")])) >
            diff(as.numeric(pred_conf[1, c("lower", "upper")])))

A confidence interval uses X V_beta X'. A prediction interval adds random_scale^2 * tau2 for a future study-specific effect; the observed study’s sampling variance vi is not added. Set random_scale = 0 for the fixed mean, or use another loading when the random effect has a known scale.

L <- matrix(1, nrow = 1L, ncol = 1L,
            dimnames = list("future study" = "(Intercept)"))
predict(fit, contrast = L, type = "link",
        interval = "prediction", random_scale = 1)
#>                  fit     lower    upper
#> (Intercept) 0.248429 -1.164695 1.661553

For a log or logit link, type = "exp" returns the exponentiated prediction and interval endpoints. It is restricted to those link scales so that transformations remain explicit.

Forest plots

The conservative forest interface accepts study-level estimates, variances, and labels explicitly. This is the safest choice for arm-level or otherwise complex input.

plot_data <- forest(
  fit,
  estimate = dat$estimate,
  vi = dat$vi,
  labels = dat$study,
  parm = "(Intercept)",
  ci_methods = c("wald", "profile", "SBC"),
  ci_args = list(renge.c = 10, silent = TRUE),
  xlab = "Effect estimate"
)

plot_data[plot_data$kind != "study", ]
#>          kind      label estimate      lower     upper     method
#> 9      pooled    Wald CI 0.248429 -0.2419502 0.7388083       wald
#> 10     pooled Profile CI 0.248429 -0.3070376 0.8050757    profile
#> 11     pooled     SBC CI 0.248429 -0.3719478 0.8701520        sbc
#> 12 prediction     95% PI 0.248429 -1.1646953 1.6615534 prediction

The lower panel shows the selected fixed-effect intervals as filled diamonds, followed by the plug-in prediction interval. ci_methods is case-insensitive and preserves the requested order. Profile and SBC rows require a single coefficient selected by parm; arbitrary linear contrasts are available for the Wald row only. The prediction row remains the Wald-type plug-in interval described above, rather than a profile- or Bartlett-corrected prediction interval. The right column reports each estimate and its limits numerically.

The same display can be requested through the generic plot method:

plot(fit, type = "forest",
     estimate = dat$estimate,
     vi = dat$vi,
     labels = dat$study,
     parm = "(Intercept)",
     ci_methods = c("wald", "profile", "SBC"),
     ci_args = list(renge.c = 10, silent = TRUE),
     xlab = "Effect estimate")

Both calls use base graphics and request Wald, profile, and SBC confidence intervals. A prediction interval is included when the fitted heterogeneity is finite. Automatic extraction is available only when a one-row-per-study structure (or a safely matched two-arm structure) is evident; otherwise pass the vectors as above.

Plain data exchange

as_metafor_data() returns an ordinary data.frame with the conventional yi, vi, and slab columns. It is useful for exchanging study-level summaries without changing the fitted object’s class.

exchange_dat <- as_metafor_data(fit)
stopifnot(is.data.frame(exchange_dat))
head(exchange_dat[c("yi", "vi", "slab")])
#>     yi   vi slab
#> 1 -0.8 0.04    1
#> 2 -0.3 0.05    2
#> 3  0.1 0.03    3
#> 4  0.7 0.06    4
#> 5  1.1 0.04    5
#> 6 -0.5 0.05    6