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