Random-effects meta-analysis with known sampling variances
Source:vignettes/meta-analysis.Rmd
meta-analysis.RmdA 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 for the observed effect size from study and for its known sampling variance. The random-effects model is
Each study sees the same pooled mean but its own total variance : the known sampling variance that the primary study already quantified, plus a shared between-study variance that the meta-analysis estimates. When every study is just a noisy measurement of one common effect (a fixed-effect, or common-effect, meta-analysis); when the true study effects themselves scatter around .
The two unknowns map onto the two drmTMB distributional
parameters:
| Meta-analysis quantity | Symbol |
drmTMB parameter |
How to read it |
|---|---|---|---|
| pooled effect |
mu intercept |
the average effect across studies | |
| known sampling variance | meta_V(V = vi) |
supplied, not estimated | |
| between-study SD | sigma |
how much true effects differ across studies | |
| between-study variance | sigma^2 |
the heterogeneity variance |
The matching R syntax is
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
.
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.04400439The 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.07662336
#> sigma:(Intercept) -1.0770967 0.19102210
#> Distributional, random-effect, scale, and correlation parameters:
#> component dpar term estimate std_error scale
#> sigma distributional-scale sigma (constant) 0.3405829 0.06505886 response
#> Derived summaries: This fit has a known-dispersion override, so a single scalar
#> residual variance is not defined.
#> logLik: -16.39
#> convergence: 0Before 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: 16 checks>
#> ok: 16; notes: 0; warnings: 0; errors: 0
#> check status
#> optimizer_convergence ok
#> convergence_status ok
#> optimizer_budget ok
#> finite_objective ok
#> logsigma_clamp_active ok
#> fixed_gradient ok
#> sdreport_status ok
#> hessian_positive_definite ok
#> hessian_conditioning ok
#> standard_errors_finite ok
#> standard_errors_inflated ok
#> observations_per_parameter ok
#> dropped_rows ok
#> positive_scale ok
#> known_sampling_covariance ok
#> fixed_effect_design_size ok
#> value
#> 0
#> converged
#> iterations=10; function=16; gradient=11
#> 16.39
#> <NA>
#> max=0.000000000001502; component=beta_sigma
#> ok
#> TRUE
#> min_eig=27.40; cond=6.222
#> range=[0.07662,0.1910]
#> n_inflated=0; max_se=0.1910; reference_median=0.1338
#> n_obs=30; n_par=2; ratio=15.00
#> 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 convergence and uncertainty diagnostics are consistent with a proper interior optimum.
#> 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.
#> Minimum eigenvalue and condition number of TMB's sdreport() fixed-effect covariance (sdr$cov.fixed), inverted. These are a genuinely different read of the fit's conditioning than TMB's internal pdHess flag -- comparable across fits, not claimed to be numerically identical to any raw TMB gradient or Hessian quantity. This fit's Hessian conditioning is within the requested threshold.
#> All fixed-effect standard errors are finite.
#> No fixed-effect standard error is inflated relative to the others.
#> Observations per estimated parameter are at or above the small-samplenote threshold (10).
#> 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.4691903In 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
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
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.1159967Heterogeneity is easier to communicate as a proportion. is the share of the total variation that is between-study rather than sampling noise. With the usual “typical” within-study variance of Higgins and Thompson (2002),
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.3119931An 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
.
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.312The 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
is known to be biased downward, because it does not account for the
degrees of freedom spent estimating
.
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.3187410 0.3490881 0.1218625The 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.3201353 0.1979836The 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.
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:
A known sampling variance enters the covariance of the response:
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 effect sizes, the dense route fits
where the supplied
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: 17 checks>
#> ok: 10; notes: 6; warnings: 1; errors: 0
#> check status
#> optimizer_convergence ok
#> convergence_status note
#> optimizer_budget ok
#> finite_objective ok
#> logsigma_clamp_active note
#> fixed_gradient ok
#> sdreport_status ok
#> hessian_positive_definite ok
#> hessian_conditioning note
#> standard_errors_finite ok
#> standard_errors_inflated note
#> weakly_identified_fixed_effects warning
#> observations_per_parameter note
#> dropped_rows ok
#> positive_scale ok
#> known_sampling_covariance note
#> fixed_effect_design_size ok
#> value
#> 0
#> degenerate
#> iterations=31; function=40; gradient=32
#> -10.90
#> -13.13
#> max=0.000000002886; component=beta_sigma
#> ok
#> TRUE
#> min_eig=0.000000005680; cond=1319517833883.
#> range=[0.02132,13268.]
#> n_inflated=1; max_se=13268.; reference_median=0.06802; example=sigma:(Intercept)
#> n_inflated=1; max_se=13268.; reference_median=0.06802; example=sigma:(Intercept)
#> n_obs=8; n_par=3; ratio=2.667
#> nobs=8; dropped=0
#> min=0.000001993
#> 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.
#> The optimizer reported success but the likelihood geometry or uncertainty is degenerate (no reliable maximum or non-finite standard errors). multi_start cannot clear this state; is_converged() is FALSE.
#> Optimizer evaluation counts recorded; no eval.max or iter.max control was supplied.
#> Objective and log-likelihood are finite.
#> The fitted log(sigma) reached -13.13, below the clamp band lower bound -12; the scale ran to the lower clamp. This is expected when the residual scale is genuinely zero (for example a meta-analysis with tau = 0), but it also happens when the response is on a very small numeric scale and the scale coefficients are wrong: compare sigma() with the spread of the response and consider rescaling the response before refitting.
#> Maximum absolute fixed gradient is <= 0.001; largest component is beta_sigma.
#> TMB::sdreport() completed successfully.
#> sdreport reports a positive-definite Hessian.
#> Minimum eigenvalue and condition number of TMB's sdreport() fixed-effect covariance (sdr$cov.fixed), inverted. These are a genuinely different read of the fit's conditioning than TMB's internal pdHess flag -- comparable across fits, not claimed to be numerically identical to any raw TMB gradient or Hessian quantity. The condition number exceeds 1e+08, signalling a near-flat, weakly identified direction; a clean pdHess is necessary, not sufficient. Treat standard errors and correlations among the affected parameters with caution.
#> 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.
#> At least one fixed-effect standard error is pathologically large relative to the reference scale; the fit is not identified enough for plain is_converged() to read TRUE and multi_start cannot clear this.
#> Fewer than 10 observations per estimated parameter; the diagnostic board does not become more cautious automatically at small n, so inspect effect sizes and uncertainty before interpreting tests.
#> 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
-
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; prefermeta_V()in new code. - There is intentionally no
meta_gaussian()family and notau ~syntax. Meta-analysis reusesfamily = gaussian()andsigma ~ ...so that it shares the distributional-regression grammar with the rest ofdrmTMB.