---
title: "Prediction and forest plots"
author: "metaGLMM authors"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Prediction and forest plots}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>",
                      fig.width = 7, fig.height = 5)
set.seed(20260821)
library(metaGLMM)
```

## A fitted model

```{r fit}
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.

```{r prediction}
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
pred_pred
pred_response
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.

```{r contrast}
L <- matrix(1, nrow = 1L, ncol = 1L,
            dimnames = list("future study" = "(Intercept)"))
predict(fit, contrast = L, type = "link",
        interval = "prediction", random_scale = 1)
```

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.

```{r forest}
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", ]
```

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:

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

```{r exchange}
exchange_dat <- as_metafor_data(fit)
stopifnot(is.data.frame(exchange_dat))
head(exchange_dat[c("yi", "vi", "slab")])
```