Random-effects meta-analysis with known sampling variances

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 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

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.

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)
#>   study          yi         vi
#> 1     1  0.43189395 0.04977587
#> 2     2  0.23230076 0.02350599
#> 3     3  0.22189632 0.07677472
#> 4     4  0.07652755 0.07261523
#> 5     5 -0.51182514 0.03998846
#> 6     6  0.10958169 0.04400439

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

Fitting the model

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

summary(fit)
#> <summary.drmTMB>
#> estimator: ML
#>                     estimate  std_error
#> mu:(Intercept)     0.3190113 0.07662335
#> sigma:(Intercept) -1.0770968 0.19102211
#> Distributional, random-effect, scale, and correlation parameters:
#>                  component  dpar       term  estimate  std_error    scale
#> sigma distributional-scale sigma (constant) 0.3405829 0.06505886 response
#> logLik: -16.39
#> convergence: 0

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.

is_converged(fit)                 # optimizer convergence
#> [1] TRUE
is_converged(fit, include_hessian = TRUE) # also requires a positive-definite Hessian
#> [1] TRUE
diagnostics <- check_drm(fit)
diagnostics[, c("check", "status", "value", "message")]
#> <drm_check: 13 checks>
#> ok: 13; notes: 0; warnings: 0; errors: 0
#>                      check status
#>      optimizer_convergence     ok
#>           optimizer_budget     ok
#>           finite_objective     ok
#>      logsigma_clamp_active     ok
#>             fixed_gradient     ok
#>            sdreport_status     ok
#>  hessian_positive_definite     ok
#>     standard_errors_finite     ok
#>   standard_errors_inflated     ok
#>               dropped_rows     ok
#>             positive_scale     ok
#>  known_sampling_covariance     ok
#>   fixed_effect_design_size     ok
#>                                                                                   value
#>                                                                                       0
#>                                                 iterations=10; function=16; gradient=11
#>                                                                                   16.39
#>                                                                                    <NA>
#>                                                   max=0.000002714; component=beta_sigma
#>                                                                                      ok
#>                                                                                    TRUE
#>                                                                  range=[0.07662,0.1910]
#>                                           n_inflated=0; max_se=0.1910; median_se=0.1338
#>                                                                      nobs=30; dropped=0
#>                                                                              min=0.3406
#>                                            type=diagonal; n=30; range=[0.02309,0.09655]
#>  total_mb=0.005508; max_cols=1; largest=mu; largest_class=matrix; largest_density=1.000
#>                                                                              message
#>                                                        nlminb convergence code is 0.
#>  Optimizer evaluation counts recorded; no eval.max or iter.max control was supplied.
#>                                             Objective and log-likelihood are finite.
#>                                   The log(sigma) clamp is not active at the optimum.
#>        Maximum absolute fixed gradient is <= 0.001; largest component is beta_sigma.
#>                                              TMB::sdreport() completed successfully.
#>                                        sdreport reports a positive-definite Hessian.
#>                                         All fixed-effect standard errors are finite.
#>                   No fixed-effect standard error is inflated relative to the others.
#>                   No rows were dropped by model-frame or known-covariance filtering.
#>                                     All fitted scale values are finite and positive.
#>                         Known sampling covariance is recorded through meta_V(V = V).
#>                          Dense fixed-effect design matrices are modest for this fit.

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().

mu_hat <- coef(fit, "mu")[["(Intercept)"]]
mu_hat
#> [1] 0.3190113

confint(fit, parm = "mu:(Intercept)")[, c("parm", "lower", "upper")]
#>                   parm     lower     upper
#> 1 fixef:mu:(Intercept) 0.1688323 0.4691903

In this run the pooled estimate is about 0.319, 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.

tau_hat <- sigma(fit)[1]
c(tau = unname(tau_hat), tau_squared = unname(tau_hat^2))
#>         tau tau_squared 
#>   0.3405829   0.1159967

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}. \]

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)
)
#> tau_squared   typical_v  I2_percent 
#>   0.1159967   0.0538076  68.3119888

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.

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)
}
#>   quantity drmTMB metafor
#>  pooled mu  0.319   0.319
#>      tau^2  0.116   0.116
#>    I^2 (%) 68.312  68.312

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.

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)
)
#>   estimator pooled_mu       tau tau_squared
#> 1        ML 0.3190113 0.3405829   0.1159967
#> 2      REML 0.3187518 0.3490881   0.1218625

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.

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")
#> (Intercept)        dose 
#>    0.320136    0.197983

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.

c(
  residual_tau_no_moderator = unname(sigma(fit)[1]),
  residual_tau_with_moderator = unname(sigma(fit_mr)[1])
)
#>   residual_tau_no_moderator residual_tau_with_moderator 
#>                   0.3405829                   0.3363327

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:

# 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:

# 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.

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)
#> <drm_check: 13 checks>
#> ok: 11; notes: 2; warnings: 0; errors: 0
#>                      check status
#>      optimizer_convergence     ok
#>           optimizer_budget     ok
#>           finite_objective     ok
#>      logsigma_clamp_active     ok
#>             fixed_gradient     ok
#>            sdreport_status     ok
#>  hessian_positive_definite     ok
#>     standard_errors_finite     ok
#>   standard_errors_inflated   note
#>               dropped_rows     ok
#>             positive_scale     ok
#>  known_sampling_covariance   note
#>   fixed_effect_design_size     ok
#>                                                                                   value
#>                                                                                       0
#>                                                 iterations=31; function=40; gradient=32
#>                                                                                  -10.90
#>                                                                                    <NA>
#>                                                max=0.000000002886; component=beta_sigma
#>                                                                                      ok
#>                                                                                    TRUE
#>                                                                  range=[0.02132,13268.]
#>                n_inflated=1; max_se=13268.; median_se=0.1147; example=sigma:(Intercept)
#>                                                                       nobs=8; dropped=0
#>                                                                         min=0.000001881
#>   type=matrix; n=8; storage=dense; density=1.000; size_mb=0.0006943; rank=8; cond=9.477
#>  total_mb=0.002609; max_cols=2; largest=mu; largest_class=matrix; largest_density=1.000
#>                                                                                                                                                                                                                                                                                                                                                                                                                     message
#>                                                                                                                                                                                                                                                                                                                                                                                               nlminb convergence code is 0.
#>                                                                                                                                                                                                                                                                                                                                         Optimizer evaluation counts recorded; no eval.max or iter.max control was supplied.
#>                                                                                                                                                                                                                                                                                                                                                                                    Objective and log-likelihood are finite.
#>                                                                                                                                                                                                                                                                                                                                                                          The log(sigma) clamp is not active at the optimum.
#>                                                                                                                                                                                                                                                                                                                                               Maximum absolute fixed gradient is <= 0.001; largest component is beta_sigma.
#>                                                                                                                                                                                                                                                                                                                                                                                     TMB::sdreport() completed successfully.
#>                                                                                                                                                                                                                                                                                                                                                                               sdreport reports a positive-definite Hessian.
#>                                                                                                                                                                                                                                                                                                                                                                                All fixed-effect standard errors are finite.
#>  At least one standard error is finite but extremely large despite a positive-definite Hessian (pdHess = TRUE), which signals a near-flat, weakly identified direction the Hessian did not resolve -- often a correlation or SD running to a boundary. Confirm with a likelihood profile (profile()) and consider a simpler model before interpreting the affected parameter; a clean Hessian is necessary, not sufficient.
#>                                                                                                                                                                                                                                                                                                                                                          No rows were dropped by model-frame or known-covariance filtering.
#>                                                                                                                                                                                                                                                                                                                                                                            All fitted scale values are finite and positive.
#>                                                                                                                                                                                                                              Known sampling covariance is recorded as a dense matrix with finite non-negative diagonal; treat this as a small-to-moderate path until sparse or block-sparse storage has benchmark evidence.
#>                                                                                                                                                                                                                                                                                                                                                                 Dense fixed-effect design matrices are modest for this fit.

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

Reference

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