Getting Started with drmeta

Overview

drmeta implements Design-Robust Meta-Analysis (DR-Meta; Hait, 2025), a variance-function random-effects framework that embeds causal design robustness directly into the heterogeneity structure of meta-analysis.

The key idea: between-study variance \(\tau^2(\text{DR}_i)\) is modelled as a monotone-decreasing function of a study-level design robustness index \(\text{DR}_i \in [0,1]\). Studies with stronger causal identification receive greater weight through the total variance \[\sigma_i^2 = v_i + \tau^2(\text{DR}_i;\,\psi),\] rather than through ad-hoc quality multipliers.

DR-Meta nests classical fixed-effects (when \(\tau_0^2 = 0\)) and random-effects (when \(\text{DR}_i = c\) constant) meta-analysis as special cases.

Installation

# From GitHub (once released):
# devtools::install_github("subirhait/drmeta")

# Local install (development):
# devtools::install("path/to/drmeta")

Step 1: Build the Design Robustness Index

The DR index \(\text{DR}_i \in [0,1]\) summarises causal identification strength per study. The drmeta package provides two helpers.

From design type labels

library(drmeta)

designs <- c("RCT", "IV", "DiD", "matching", "OLS", "OLS", "RCT", "DiD")
dr <- dr_from_design(designs)
data.frame(Design = designs, DR = dr)
#>     Design   DR
#> 1      RCT 1.00
#> 2       IV 0.75
#> 3      DiD 0.60
#> 4 matching 0.50
#> 5      OLS 0.25
#> 6      OLS 0.25
#> 7      RCT 1.00
#> 8      DiD 0.60

From continuous sub-scores

k <- 8
balance  <- c(0.95, 0.80, 0.60, 0.70, 0.30, 0.25, 0.90, 0.65)
overlap  <- c(0.90, 0.75, 0.55, 0.80, 0.40, 0.35, 0.85, 0.70)
design_s <- dr_from_design(designs)   # reuse labels from above

dr_composite <- dr_score(
  balance  = balance,
  overlap  = overlap,
  design   = design_s,
  weights  = c(2, 2, 3)   # design quality weighted most
)
round(dr_composite, 3)
#> [1] 0.957 0.764 0.586 0.643 0.307 0.279 0.929 0.643
#> attr(,"subscores")
#>   balance overlap design        DR
#> 1    0.95    0.90   1.00 0.9571429
#> 2    0.80    0.75   0.75 0.7642857
#> 3    0.60    0.55   0.60 0.5857143
#> 4    0.70    0.80   0.50 0.6428571
#> 5    0.30    0.40   0.25 0.3071429
#> 6    0.25    0.35   0.25 0.2785714
#> 7    0.90    0.85   1.00 0.9285714
#> 8    0.65    0.70   0.60 0.6428571

Step 2: Fit the DR-Meta Model

set.seed(42)
k  <- 20
dr <- runif(k, 0.1, 0.95)
vi <- runif(k, 0.008, 0.055)

# True data-generating model (Sections 3–4, Hait 2025):
#   tau0sq = 0.04, gamma = 2, mu = 0.30
tau2_true <- 0.04 * exp(-2 * dr)
yi <- rnorm(k, mean = 0.30, sd = sqrt(vi + tau2_true))

# Fit DR-Meta (exponential variance function, REML)
fit <- drmeta(yi = yi, vi = vi, dr = dr)
print(fit)
#> 
#> --- DR-Meta: Design-Robust Random-Effects Model ---
#> 
#> Variance function : exponential
#> Estimation method : REML
#> Studies (k)       : 20
#> 
#> Pooled estimate   : 0.2502
#> Std. error        : 0.0504
#> 95% CI            : [0.1514, 0.3491]
#> z = 4.9617,  p = 0.0000
#> 
#> tau0^2 (DR=0)     : 0.0283
#> gamma             : 0.4338
summary(fit)
#> 
#> === DR-Meta Summary ===
#> 
#> Call:
#>   drmeta(yi = yi, vi = vi, dr = dr) 
#> 
#> --- Pooled Effect ---
#>   mu   = 0.2502  (SE = 0.0504)
#>   95% CI: [0.1514, 0.3491]
#>   z = 4.9617,  p = 0.0000
#> 
#> --- Variance-Function Parameters ---
#>   tau0^2 = 0.0283  (heterogeneity at DR=0)
#>   gamma  = 0.4338  (decay rate)
#>   vfun   = "exponential",  method = "REML"
#> 
#> --- Model Fit ---
#>   logLik (ML)   = 19.6192
#>   logLik (REML) = 16.6321
#>   AIC = -33.2384,  BIC = -30.2512
#> 
#> Converged: TRUE

Compare exponential vs linear variance function

fit_lin <- drmeta(yi, vi, dr, vfun = "linear")
cat("Exponential AIC:", round(AIC(fit), 2), "\n")
#> Exponential AIC: -33.24
cat("Linear      AIC:", round(AIC(fit_lin), 2), "\n")
#> Linear      AIC: -33.25

Step 3: Visualise

Forest plot

dr_forest(fit, main = "DR-Meta Forest Plot — Simulated Data")

Variance function

dr_plot_vfun(fit, main = "Estimated Variance Function")

Weight vs DR (Lemma 3)

dr_plot(fit, main = "DR-Meta Weights vs Design Robustness")

Step 4: Heterogeneity Decomposition (Proposition 6)

het <- dr_heterogeneity(fit)

cat("--- Heterogeneity Summary ---\n")
#> --- Heterogeneity Summary ---
print(het$summary)
#>    k       Q df   pval tau2_mean   I2     H2
#> 1 20 19.2252 19 0.4425    0.0217 1.17 1.0119

cat("\n--- Proposition 6 Decomposition ---\n")
#> 
#> --- Proposition 6 Decomposition ---
print(het$decomposition)
#>   tau2_design_explained tau2_residual tau2_total   R2_DR
#> 1                0.0217             0      6e-04 35.7437
cat(sprintf("\n%.1f%% of total heterogeneity is design-explained (R2_DR).\n",
            het$decomposition$R2_DR * 100))
#> 
#> 3574.4% of total heterogeneity is design-explained (R2_DR).

cat("\n--- Per-Study Q Contributions (top 5) ---\n")
#> 
#> --- Per-Study Q Contributions (top 5) ---
top5 <- head(het$contributions[order(-het$contributions$pct_Q), ], 5)
print(top5)
#>      study    DR tau2_i    q_i pct_Q
#> 19 Study19 0.504 0.0227 4.3803 22.78
#> 5   Study5 0.645 0.0214 3.4038 17.70
#> 8   Study8 0.214 0.0258 2.5166 13.09
#> 16 Study16 0.899 0.0191 1.8310  9.52
#> 4   Study4 0.806 0.0199 1.6765  8.72

Step 5: Influence Diagnostics (LOO)

loo <- dr_loo(fit)
cat("Influential studies:\n")
#> Influential studies:
print(subset(loo$summary, influential == TRUE))
#>  [1] study        DR           est_loo      ci.lb_loo    ci.ub_loo   
#>  [6] tau0sq_loo   gamma_loo    delta_mu     delta_tau0sq delta_gamma 
#> [11] influential 
#> <0 rows> (or 0-length row.names)
cat("\nFull LOO summary (first 5 rows):\n")
#> 
#> Full LOO summary (first 5 rows):
print(head(loo$summary, 5))
#>    study        DR   est_loo ci.lb_loo ci.ub_loo tau0sq_loo   gamma_loo
#> 1 Study1 0.8775851 0.2495530 0.1467965 0.3523095 0.02328814 0.001875022
#> 2 Study2 0.8965141 0.2738357 0.1747192 0.3729523 0.04867336 1.565702714
#> 3 Study3 0.3432186 0.2508325 0.1485454 0.3531196 0.04366491 0.983644766
#> 4 Study4 0.8058805 0.2380564 0.1382161 0.3378968 0.02956543 0.554658594
#> 5 Study5 0.6454837 0.2246103 0.1369916 0.3122290 0.01436547 0.711650788
#>        delta_mu delta_tau0sq delta_gamma influential
#> 1 -0.0006945935 -0.004987028  -0.4318844       FALSE
#> 2  0.0235881283  0.020398191   1.1319433       FALSE
#> 3  0.0005848714  0.015389746   0.5498853       FALSE
#> 4 -0.0121911588  0.001290268   0.1208992       FALSE
#> 5 -0.0256373003 -0.013909691   0.2778914       FALSE

Step 6: Publication Bias

pb <- dr_pub_bias(fit)
cat("--- PET ---\n");   print(pb$PET)
#> --- PET ---
#>             estimate       se      zval       pval
#> intercept  0.4189576 0.172371  2.430557 0.01507562
#> sei       -1.0296246 1.005491 -1.024001 0.30583463
cat("\n--- PEESE ---\n"); print(pb$PEESE)
#> 
#> --- PEESE ---
#>             estimate        se      zval         pval
#> intercept  0.3490149 0.1047119  3.333097 0.0008588502
#> vi        -3.3607913 3.1209101 -1.076863 0.2815416631
cat("\n--- Recommendation ---\n"); cat(pb$recommendation, "\n")
#> 
#> --- Recommendation ---
#> PET slope is not significant (p = 0.3058): no strong evidence of small-study bias.
#>   Use main DR-Meta pooled estimate: 0.2502.
dr_funnel(fit, main = "DR-Meta Funnel Plot")

Connection to the Literature

DR-Meta is closely related to the meta-analytic location-scale model (Viechtbauer & Lopez-Lopez, 2022), in which both the mean and variance of the true effect distribution can depend on study-level covariates. DR-Meta contributes a specific causal motivation for the scale component: the variance function is an index of susceptibility to confounding, grounded in design robustness theory (Rubin, 2008; Rosenbaum, 2010).

It is complementary to bias-model frameworks (Turner et al., 2009; Rhodes et al., 2015; Mathur & VanderWeele, 2020), which adjust the location (mean) for anticipated confounding. DR-Meta adjusts the scale (variance), so the two families can be combined in future extensions.

References