Random-effects meta-analysis with known sampling variances • drmTMB

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 the gate for trusting the Wald standard errors and intervals.

fit$opt$convergence # 0 means the optimiser converged
#> [1] 0
fit$sdr$pdHess # TRUE means the Hessian is positive definite
#> [1] TRUE

check_drm(fit)
#> <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.

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.

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.

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 (as above), a vector, a diagonal matrix, or a dense covariance matrix when the effect sizes are correlated.
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.