---
title: "Random-effects meta-analysis with known sampling variances"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Random-effects meta-analysis with known sampling variances}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
if (!"package:drmTMB" %in% search()) {
  library(drmTMB)
}
```

<style>
@media (max-width: 575.98px) {
  main h1 {
    hyphens: none;
    overflow-wrap: normal;
    word-break: normal;
  }

  main table {
    display: block;
    max-width: 100%;
    overflow-x: auto;
  }
}
</style>

A random-effects meta-analysis pools effect sizes from several studies while
acknowledging that the studies differ. Each study reports an effect size
together with its sampling variance, which is treated as **known**. The model
estimates two things on top of that known uncertainty: a pooled mean effect and
the between-study heterogeneity. This is a specialist route for data that
already are effect sizes with known sampling variances; for raw observations,
start with [Distributional regression with drmTMB](drmTMB.html) instead.

In `drmTMB` this is ordinary Gaussian distributional regression with a known
sampling covariance. There is no separate meta-analysis family. You fit
`family = gaussian()`, you supply the known per-study variances with
`meta_V()` inside the location formula, and the between-study heterogeneity is
the residual scale `sigma`.

Throughout, `Normal(a, b)` uses the variance (not the standard deviation) as
its second argument.

## The model

Write \(y_i\) for the observed effect size from study \(i\) and \(v_i\) for its
**known** sampling variance. The random-effects model is

\[
y_i \mid \mu, \tau, v_i
  \;\sim\;
  \operatorname{Normal}\!\left(\mu,\; v_i + \tau^2\right),
\qquad i = 1, \ldots, K.
\]

Each study sees the same pooled mean \(\mu\) but its own total variance
\(v_i + \tau^2\): the known sampling variance \(v_i\) that the primary study
already quantified, plus a shared between-study variance \(\tau^2\) that the
meta-analysis estimates. When \(\tau = 0\) every study is just a noisy
measurement of one common effect (a fixed-effect, or common-effect,
meta-analysis); when \(\tau > 0\) the true study effects themselves scatter
around \(\mu\).

The two unknowns map onto the two `drmTMB` distributional parameters:

| Meta-analysis quantity | Symbol | `drmTMB` parameter | How to read it |
|---|---|---|---|
| pooled effect | \(\mu\) | `mu` intercept | the average effect across studies |
| known sampling variance | \(v_i\) | `meta_V(V = vi)` | supplied, not estimated |
| between-study SD | \(\tau\) | `sigma` | how much true effects differ across studies |
| between-study variance | \(\tau^2\) | `sigma^2` | the heterogeneity variance \(\tau^2\) |

The matching R syntax is

```r
drmTMB(
  bf(yi ~ 1 + meta_V(V = vi), sigma ~ 1),
  family = gaussian(),
  data = dat
)
```

Two points are worth stating plainly.

First, `vi` must be a **variance**. If your dataset stores standard errors,
square them (`vi <- se^2`) before fitting.

Second, the residual scale `sigma` *is* the between-study heterogeneity \(\tau\).
Meta-analysts usually call this quantity `tau` and report `tau^2`. `drmTMB`
keeps the public parameter name `sigma` so that meta-analysis uses exactly the
same distributional grammar as every other Gaussian model in the package; you
recover the familiar `tau^2` simply by squaring `sigma`. The marker
`meta_V()` puts the *known* variance into the likelihood; `sigma` estimates the
*unknown* variance that is left over.

This is an implemented, source-tested interface and the example below checks
its ML fit against `metafor`. It does **not** currently have a registered
`meta_V()` capability-ledger cell, so this tutorial does not assign it an
evidence tier or make an interval-coverage claim. Treat fitted-model intervals
as the methods returned by the current fit, and state their source when
reporting them.

## A simulated dataset

We simulate `K = 30` studies. Each study has a true effect drawn around a
common mean, and we observe that true effect with known sampling error.

```{r}
set.seed(101)

K <- 30
mu_true <- 0.40 # pooled effect
tau_true <- 0.30 # between-study SD

# Known sampling variances: larger studies (smaller vi) and smaller studies.
vi <- runif(K, 0.02, 0.10)

# True study effects scatter around mu_true with SD tau_true.
theta_i <- rnorm(K, mean = mu_true, sd = tau_true)

# Observed effect sizes: each true effect seen with its known sampling error.
yi <- rnorm(K, mean = theta_i, sd = sqrt(vi))

dat <- data.frame(study = factor(seq_len(K)), yi = yi, vi = vi)
head(dat)
```

The data frame has one row per study: the effect size `yi` and its known
sampling variance `vi`.

## Fitting the model

```{r}
fit <- drmTMB(
  bf(yi ~ 1 + meta_V(V = vi), sigma ~ 1),
  family = gaussian(),
  data = dat
)

summary(fit)
```

Before reading the coefficients, confirm the fit is sound. A clean optimisation
and a positive-definite Hessian are necessary diagnostics for routine Wald
output, but they do not by themselves establish that a `sigma` interval is
finite, usable, or coverage-valid.

```{r}
is_converged(fit)                 # optimizer convergence
is_converged(fit, include_hessian = TRUE) # also requires a positive-definite Hessian
```

```{r}
diagnostics <- check_drm(fit)
diagnostics[, c("check", "status", "value", "message")]
```

If `is_converged(fit, include_hessian = TRUE)` is `FALSE`, use the status and
message in `diagnostics` to identify the failed check before interpreting Wald
intervals. A Hessian warning is an inference warning, not by itself proof that
the fitted pooled effect is unusable.

## The pooled effect

The pooled effect is the `mu` intercept. Its Wald confidence interval comes from
`confint()`.

```{r}
mu_hat <- coef(fit, "mu")[["(Intercept)"]]
mu_hat

confint(fit, parm = "mu:(Intercept)")[, c("parm", "lower", "upper")]
```

In this run the pooled estimate is about `r round(mu_hat, 3)`, which sits near
the simulated `mu_true = 0.40`. The interval reflects uncertainty in the mean
*after* both the known sampling variances and the estimated between-study
heterogeneity have been accounted for.

## Between-study heterogeneity

The between-study SD \(\tau\) is the residual scale `sigma`. Because the `sigma`
formula here is intercept-only, every study shares the same value, so we take
the first element. Squaring it gives the heterogeneity variance \(\tau^2\) that
meta-analysis reports.

```{r}
tau_hat <- sigma(fit)[1]
c(tau = unname(tau_hat), tau_squared = unname(tau_hat^2))
```

Heterogeneity is easier to communicate as a proportion. \(I^2\) is the share of
the total variation that is between-study rather than sampling noise. With the
usual "typical" within-study variance \(\tilde v\) of Higgins and Thompson
(2002),

\[
\tilde v =
  \frac{(K-1)\sum_i w_i}{\left(\sum_i w_i\right)^2 - \sum_i w_i^2},
\qquad w_i = 1 / v_i,
\qquad
I^2 = \frac{\tau^2}{\tau^2 + \tilde v}.
\]

```{r}
w <- 1 / dat$vi
v_typical <- ((K - 1) * sum(w)) / (sum(w)^2 - sum(w^2))
I2 <- tau_hat^2 / (tau_hat^2 + v_typical)

c(
  tau_squared = unname(tau_hat^2),
  typical_v = v_typical,
  I2_percent = unname(100 * I2)
)
```

An \(I^2\) of this size means a substantial fraction of the variation among the
observed effect sizes reflects genuine differences between studies, not just
within-study sampling error. The pooled mean is still meaningful, but it is a
mean of effects that really do differ.

## Cross-check against metafor

The same model can be fitted with `metafor::rma()` using maximum likelihood. It
should agree with `drmTMB`, because both fit the identical random-effects
likelihood \(y_i \sim \operatorname{Normal}(\mu, v_i + \tau^2)\). This is a
useful sanity check when you first adopt the `drmTMB` spelling.

```{r}
if (requireNamespace("metafor", quietly = TRUE)) {
  rma_fit <- metafor::rma(yi = yi, vi = vi, method = "ML", data = dat)

  comparison <- data.frame(
    quantity = c("pooled mu", "tau^2", "I^2 (%)"),
    drmTMB = c(mu_hat, tau_hat^2, 100 * I2),
    metafor = c(as.numeric(rma_fit$beta), rma_fit$tau2, rma_fit$I2)
  )
  print(comparison, row.names = FALSE, digits = 4)
}
```

The two engines return the same pooled effect and the same heterogeneity
variance. `drmTMB` is doing ML random-effects meta-analysis; it simply spells
the known sampling variance as `meta_V()` and the heterogeneity as `sigma`.

## REML for the heterogeneity

The ML estimate of \(\tau^2\) is known to be biased downward, because it does
not account for the degrees of freedom spent estimating \(\mu\). When the mean
model is fixed and you only want a better heterogeneity estimate, restricted
maximum likelihood (`REML = TRUE`) is the standard remedy. Keep ML
(`REML = FALSE`, the default) whenever you intend to compare different
fixed-effect mean models with AIC or BIC, since restricted likelihoods are not
comparable across different mean structures.

```{r}
fit_reml <- drmTMB(
  bf(yi ~ 1 + meta_V(V = vi), sigma ~ 1),
  family = gaussian(),
  data = dat,
  REML = TRUE
)

data.frame(
  estimator = c("ML", "REML"),
  pooled_mu = c(coef(fit, "mu")[[1]], coef(fit_reml, "mu")[[1]]),
  tau = c(sigma(fit)[1], sigma(fit_reml)[1]),
  tau_squared = c(sigma(fit)[1]^2, sigma(fit_reml)[1]^2)
)
```

The REML between-study variance is slightly larger than the ML one, as expected.

## Meta-regression: moderators on the mean

If a study-level covariate might explain part of the variation in effect sizes,
add it to the `mu` formula. This is a random-effects meta-regression: the known
sampling variances stay in `meta_V()`, `sigma` becomes the *residual* (after
moderators) between-study SD, and the new coefficient measures how the effect
size changes with the moderator.

```{r}
set.seed(202)
dat$dose <- scale(runif(K, 1, 10))[, 1] # a study-level moderator
# Give the effect size a genuine dependence on the moderator.
dat$yi <- dat$yi + 0.25 * dat$dose

fit_mr <- drmTMB(
  bf(yi ~ 1 + dose + meta_V(V = vi), sigma ~ 1),
  family = gaussian(),
  data = dat
)

coef(fit_mr, "mu")
```

The `dose` coefficient is the change in the pooled effect per one-SD change in
the moderator. After fitting a moderator, the residual `sigma` is the
between-study heterogeneity that the moderator did *not* explain; comparing it
with the no-moderator `sigma` shows how much heterogeneity the moderator
absorbed.

```{r}
c(
  residual_tau_no_moderator = unname(sigma(fit)[1]),
  residual_tau_with_moderator = unname(sigma(fit_mr)[1])
)
```

## Multiple effect sizes per study

The worked example above has one effect size per study, so `sigma` carries the
whole between-study story. When a study contributes several effect sizes, two
levels of variation appear: a study-level random effect for the studies, and a
residual for the effect sizes within a study. Those are different questions, and
`drmTMB` keeps them in different places:

```r
# Schematic: several effect sizes per study (not evaluated here).
drmTMB(
  bf(yi ~ 1 + moderator + (1 | study) + meta_V(V = vi), sigma ~ 1),
  family = gaussian(),
  data = dat_repeated
)
```

Here `(1 | study)` is the between-study random effect and `sigma` is the
within-study residual heterogeneity, while `meta_V(V = vi)` still supplies the
known sampling variances. A grouping factor used in `(1 | study)` must have at
least one study with repeated rows; a data set with exactly one row per study is
the single-level model shown above, where `sigma` alone represents
between-study heterogeneity.

## Layered heterogeneity: an experimental contract

When repeated effects are nested in studies, three SD layers answer different
questions: known sampling uncertainty in `meta_V(V = V)`, residual
heterogeneity in `sigma`, and the SD of a study- or effect-level location
random effect in `sd(group)`. The following formulas are accepted by the
current Gaussian implementation:

```{r layered-meta-syntax, eval=FALSE}
# Study-level location-SD regression (LSS): z_study is constant within study.
drmTMB(
  bf(
    yi ~ x + (1 | study) + meta_V(V = V),
    sigma ~ z,
    sd(study) ~ z_study
  ),
  family = gaussian(), data = dat
)

# Nested effect-level location-SD regression (LSSS): effect is nested in study
# and has repeated rows.
drmTMB(
  bf(
    yi ~ x + (1 | study) + (1 | effect) + meta_V(V = V),
    sigma ~ z,
    sd(study) ~ z_study,
    sd(effect) ~ z_effect
  ),
  family = gaussian(), data = dat
)
```

These formulas are an **experimental local-contract surface**, not a claim of
calibrated inference. In the Arc 7B local sentinel, the dense-`V` LSS study-SD
profiles had non-finite endpoints despite a positive-definite Hessian. Before
using these layered fits for interval-based conclusions, inspect every profile
target and consult the current development evidence. A random term in `sigma`,
such as `sigma ~ z + (1 | study)`, is a different double-hierarchical model: it
models variation in residual SD, not variation in the SD of the location random
effect.

## Known sampling variance is not a weight

Inverse-variance weights and known sampling variances answer different
questions, and `meta_V()` is not the same as the top-level `weights` argument.

A likelihood weight multiplies a study's contribution to the log-likelihood:

\[
\ell(\theta) = \sum_i w_i \, \ell_i(\theta).
\]

A known sampling variance enters the *covariance* of the response:

\[
y_i \sim \operatorname{Normal}(\mu,\; v_i + \tau^2).
\]

So `weights = 1 / vi` is **not** the random-effects meta-analysis model. It
rescales how much each row counts toward the likelihood; it does not put `vi`
into the modelled sampling variance, and it does not let `tau^2` be estimated on
top of the known variances. For meta-analysis with known sampling variances,
use `meta_V(V = vi)`. Reserve `weights =` for genuine likelihood weights such as
externally defined case weights.

## Correlated effect sizes: a small dense covariance example

Use a dense `V` when the *sampling errors* of different effect sizes are known
to be correlated, for example because they share participants or a control
group. This remains an observation-level measurement-error input. It is not a
weight, and it does not create a latent study, phylogenetic, spatial, or
relatedness effect.

For \(n\) effect sizes, the dense route fits

\[
\mathbf y \sim \operatorname{Normal}\!\left(X\boldsymbol\beta,
  V + \sigma^2 I_n\right),
\]

where the supplied \(V\) is the known sampling covariance and `sigma` is the
estimated residual heterogeneity. The following small example is executable;
the off-diagonal entries say that nearby effect sizes have correlated sampling
errors.

```{r}
set.seed(303)
n_dense <- 8
dat_dense <- data.frame(
  yi = 0.25 + 0.10 * seq_len(n_dense) + stats::rnorm(n_dense, sd = 0.04),
  x = seq_len(n_dense)
)
V_dense <- 0.012 * outer(
  seq_len(n_dense), seq_len(n_dense),
  function(i, j) 0.55^abs(i - j)
)

# A useful preflight: V is numeric, n by n, symmetric, and PSD.
stopifnot(
  is.numeric(V_dense),
  identical(dim(V_dense), c(nrow(dat_dense), nrow(dat_dense))),
  isTRUE(all.equal(V_dense, t(V_dense))),
  min(eigen(V_dense, symmetric = TRUE, only.values = TRUE)$values) >= 0
)

fit_dense <- drmTMB(
  bf(yi ~ x + meta_V(V = V_dense), sigma ~ 1),
  family = gaussian(),
  data = dat_dense
)
check_drm(fit_dense)
```

The order is part of the data contract: row and column `i` of `V_dense` must
describe the sampling covariance for row `i` of `dat_dense`. Supply one row and
one column per original data row; do not reorder the data after building `V`.
For a dense matrix, `drmTMB` first applies the ordinary model-frame exclusions
(for example a missing response or predictor), then removes the corresponding
rows *and columns* of `V` together. Let the model frame make that exclusion and
keep the original ordering rather than filtering the data and covariance matrix
by separate rules. The retained matrix must have finite entries, be symmetric,
have a non-negative diagonal, and be positive semidefinite.

## Notes on the function names

* `meta_V(V = V)` is the current marker for known sampling variance or
  covariance. The argument may be a column of **variances** (not standard
  errors), a vector, a diagonal matrix, or a dense covariance matrix when the
  effect sizes are correlated. A dense matrix must follow the row order and
  validity contract in the preceding example.
* `meta_known_V(V = V)` is a deprecated alias kept only for backward
  compatibility. It routes to the same additive known-variance likelihood but
  emits a deprecation warning; prefer `meta_V()` in new code.
* There is intentionally no `meta_gaussian()` family and no `tau ~` syntax.
  Meta-analysis reuses `family = gaussian()` and `sigma ~ ...` so that it shares
  the distributional-regression grammar with the rest of `drmTMB`.

## Reference

Higgins, J. P. T., and Thompson, S. G. (2002). Quantifying heterogeneity in a
meta-analysis. *Statistics in Medicine*, 21(11), 1539-1558.