Design-indexed heterogeneity with drmeta

Question and model

drmeta examines whether residual between-study heterogeneity follows a prespecified ordered study-design score. Its model is

\[ y_i \sim N\{x_i^\top\beta,\ v_i + \tau_0^2\exp(-\gamma d_i)\}. \]

The default restriction \(\gamma\geq0\) encodes a constant or decreasing variance function. It is a hypothesis to check. It does not make the design score a quality weight, remove bias, or adjust a design-related mean shift.

This vignette uses the 58-study CBT recidivism data in metadat. The score is declared before fitting: nonequivalent groups 0, matched groups 0.5, and randomized trials 1.

dat <- metadat::dat.landenberger2005
design <- tolower(trimws(as.character(dat$design)))
dat$dr <- c(nonequiv = 0, match = .5, rct = 1)[design]
stopifnot(!anyNA(dat$dr))
dat <- metafor::escalc("OR", ai = n.cbt.non, bi = n.cbt.rec,
                       ci = n.ctrl.non, di = n.ctrl.rec, data = dat)
table(dat$dr)
#> 
#>   0 0.5   1 
#>  16  23  19

Fit a minimum comparison set

The four fits separate constant heterogeneity, a directional scale relation, an unrestricted scale relation, and a location-plus-scale model.

constant <- drmeta(dat$yi, dat$vi, dat$dr, gamma_fixed = 0,
                   slab = dat$study)
directional <- drmeta(dat$yi, dat$vi, dat$dr, slab = dat$study)
unrestricted <- drmeta(dat$yi, dat$vi, dat$dr, constrained = FALSE,
                       slab = dat$study)
joint <- drmeta(dat$yi, dat$vi, dat$dr, mods = 1 - dat$dr,
                slab = dat$study)

data.frame(
  model = c("constant", "directional", "unrestricted", "joint"),
  beta0 = vapply(list(constant, directional, unrestricted, joint),
                 function(x) unname(x$beta[1]), numeric(1)),
  tau0sq = vapply(list(constant, directional, unrestricted, joint),
                  function(x) x$tau0sq, numeric(1)),
  gamma = vapply(list(constant, directional, unrestricted, joint),
                 function(x) x$gamma, numeric(1))
)
#>          model     beta0    tau0sq         gamma
#> 1     constant 0.4225970 0.1046084  0.0000000000
#> 2  directional 0.4225966 0.1046072  0.0000000000
#> 3 unrestricted 0.4225988 0.1045828 -0.0004723874
#> 4        joint 0.4387907 0.1084085  0.0000000000

The full-data directional estimate is at \(\gamma=0\). The unrestricted value is also effectively zero. This says the monotone decreasing pattern is not supported. It does not establish that heterogeneity is identical across the three design categories.

Check the shape before interpreting zero as flatness

A categorical scale fit can reveal a pattern that the one-parameter monotone curve cannot represent. Supply an ordered factor so the printed group order matches the design score.

ordered_design <- factor(design,
                         levels = c("nonequiv", "match", "rct"))
shape <- dr_shape_check(directional, ordered_design)
shape
#> Grouped scale-shape diagnostic
#> 
#>     group  n      tau2
#>  nonequiv 16 5.555e-08
#>     match 23 1.571e-01
#>       rct 19 1.801e-02
#> 
#> Pattern: interior peak 
#> LR = 5.991  df = 2  p = 0.05

The full data have fitted \(\tau^2\) values of approximately 0, 0.157, and 0.018 for nonequivalent, matched, and randomized studies. The categorical versus constant-scale likelihood-ratio statistic is about 5.99 (\(p\approx .05\)). The middle peak is outside the shape of a monotone exponential curve, so the boundary estimate partly reflects shape misspecification.

The likelihood-ratio reference is approximate when a grouped variance is near zero. Treat this as a diagnostic, and report that limitation.

Boundary inference and influence

When the directional estimate is zero, drmeta_bootstrap_gamma() reports the boundary and does not simulate an uninformative null distribution.

drmeta_bootstrap_gamma(directional, B = 999, seed = 20260928)
#> Parametric-bootstrap test of the scale gradient
#> 
#> H0: gamma = 0   vs   H1: gamma > 0   (boundary null)
#> 
#> The constrained estimate is at the boundary (gamma = 0).
#> The likelihood-ratio statistic is zero by construction, so the
#> bootstrap is uninformative and was not run. The data provide no
#> support for a positive scale gradient.
#> 
#> tau0^2: 0.1046 
#> 
#> Refit with constrained = FALSE to check whether the unrestricted
#> gradient is negative, which would contradict the constraint rather
#> than merely fail to support it.

Study labels flow from slab into the leave-one-out results.

loo <- dr_loo(directional)
anderson <- loo[loo$study == "Anderson (2002)", ]
anderson[, c("study", "est_loo", "tau0sq_loo", "gamma_loo")]
#>              study  est_loo tau0sq_loo gamma_loo
#> 54 Anderson (2002) 0.339541 0.04156371 0.4570625
sum(loo$gamma_loo > 0, na.rm = TRUE)
#> [1] 23

Deleting Anderson (2002) gives a constant-mean estimate near 0.340, baseline variance near 0.0416, and \(\gamma\approx0.457\), a fitted 36.7% variance decrease over scores 0 to 1. Twenty-three of the 58 single deletions give a positive gradient. These results show instability in the point estimate. They do not license deletion.

The positive Anderson-deletion fit also needs boundary-aware inference.

keep <- dat$study != "Anderson (2002)"
without_anderson <- drmeta(dat$yi[keep], dat$vi[keep], dat$dr[keep],
                           slab = dat$study[keep])
anderson_test <- drmeta_bootstrap_gamma(without_anderson, B = 999,
                                        seed = 20260928)
c(LR = anderson_test$statistic, p = anderson_test$p.value)
#>         LR          p 
#> 0.06497798 0.24700000

The verified run gives LR about 0.065 and \(p=0.247\). Deletion changes the point estimate, while the inferential conclusion remains a lack of clear support for a positive gradient. After this deletion, the categorical pattern also weakens (LR about 1.55, \(p=0.46\)).

Score sensitivity

Changing category spacing and changing the contrast answer different questions. A randomized-versus-rest score gives \(\hat\gamma\approx1.78\) in this example, but it reduces the scale predictor to two support points.

randomized_vs_rest <- as.numeric(dat$dr == 1)
binary_fit <- drmeta(dat$yi, dat$vi, randomized_vs_rest,
                     slab = dat$study)
#> Warning in drmeta(dat$yi, dat$vi, randomized_vs_rest, slab = dat$study): The
#> design-robustness index has fewer than three distinct values; gamma may be
#> weakly identified.
binary_fit$gamma
#> [1] 1.781187

matched_at_075 <- ifelse(dat$dr == .5, .75, dat$dr)
spacing_fit <- drmeta(dat$yi, dat$vi, matched_at_075,
                      slab = dat$study)
c(primary = directional$gamma,
  matched_at_075 = spacing_fit$gamma,
  randomized_vs_rest = binary_fit$gamma)
#>            primary     matched_at_075 randomized_vs_rest 
#>           0.000000           0.000000           1.781187

Prespecify the primary coding. Use sensitivity fits to show dependence on defensible alternatives, not to select the most favorable result.

Numerical validation

The package includes a regression check against a general location-scale implementation. For metafor::rma(scale = ~ dr), the parameter map is \(\alpha_0=\log(\tau_0^2)\) and \(\alpha_1=-\gamma\). With ML and no constraint, the two packages agree to numerical tolerance, including the log-likelihood. With REML, point estimates agree, while raw log-likelihood values use different additive constants and should not be compared across packages. This is an implementation check rather than a separate analysis goal.

mf <- metafor::rma(yi, vi, scale = ~ dr, data = dat, method = "ML")
dm <- drmeta(dat$yi, dat$vi, dat$dr, constrained = FALSE,
             method = "ML", slab = dat$study)
stopifnot(
  abs(unname(dm$beta[1] - mf$beta[1])) < 1e-5,
  abs(unname(dm$tau0sq - exp(mf$alpha[1]))) < 1e-5,
  abs(unname(dm$gamma + mf$alpha[2])) < 1e-5,
  abs(as.numeric(logLik(dm)) - as.numeric(logLik(mf))) < 1e-5
)

Reporting checklist

Report the score definition and support, the direction of the effect measure, the constant, directional, unrestricted, and joint fits, boundary status, fitted variance contrasts within observed support, the grouped shape check, and leave-one-out changes. A scale model changes weights. It does not remove a score-related location difference or establish a causal effect of design.

sessionInfo()
#> R version 4.6.0 (2026-04-24 ucrt)
#> Platform: x86_64-w64-mingw32/x64
#> Running under: Windows 11 x64 (build 26200)
#> 
#> Matrix products: default
#>   LAPACK version 3.12.1
#> 
#> locale:
#> [1] LC_COLLATE=C                          
#> [2] LC_CTYPE=English_United States.utf8   
#> [3] LC_MONETARY=English_United States.utf8
#> [4] LC_NUMERIC=C                          
#> [5] LC_TIME=English_United States.utf8    
#> 
#> time zone: America/Chicago
#> tzcode source: internal
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] metafor_5.0-1       numDeriv_2016.8-1.1 metadat_1.6-0      
#> [4] Matrix_1.7-5        drmeta_0.2.3       
#> 
#> loaded via a namespace (and not attached):
#>  [1] nlme_3.1-169      cli_3.6.6         knitr_1.51        rlang_1.2.0      
#>  [5] xfun_0.60         otel_0.2.0        jsonlite_2.0.0    htmltools_0.5.9  
#>  [9] sass_0.4.10       rmarkdown_2.31    grid_4.6.0        evaluate_1.0.5   
#> [13] jquerylib_0.1.4   fastmap_1.2.0     yaml_2.3.12       lifecycle_1.0.5  
#> [17] compiler_4.6.0    mathjaxr_2.0-0    rstudioapi_0.19.0 lattice_0.22-9   
#> [21] digest_0.6.39     R6_2.6.1          bslib_0.11.0      tools_4.6.0      
#> [25] cachem_1.1.0