Trace plots in metafor

[SamCH93]

Classification

Feature Request

Summary

Hi Wolfgang,

First, thank you very much for developing and maintaining this awesome package!

I recently read the preprint "How trace plots help interpret meta-analysis
results" from Christian Röver, David Rindskopf, and Tim Friede
(https://doi.org/10.48550/arXiv.2306.17043), perhaps you have also seen it.
The idea is to show the overall effect size estimate + the study specific effect
size estimates as a function of the heterogeneity tau and also show the
prior/posterior of tau below that plot, below one figure from the paper

I like the idea and think it is really useful for simultaneously assessing the
uncertainty of effect sizes and heterogeneity and give analysts an indication
for the sensitivity of their results. While these authors propose a Bayesian
version of the trace plot (available in bayesmeta) they also show possible
frequentist extensions. Instead of the prior/posterior, they show the Q-test
statistic as a function of the heterogeneity.

IMO instead of the Q-test statistic it would perhaps be more useful to show the
p-value function or confidence curve (i.e., all possible confidence intervals
stacked on top of each other) as this is more analogue to the posterior
distribution. Below is a hacky prototype for the most simple case of
random-effects meta-analysis without moderators. I used confint.rma.uni to
compute the confidence curve, but this doesn't work so well for getting the
confidence curve for the whole tau range to zero, for sure there is a better and
more efficient way. I also added confidence bands for the overall effect, which
the authors of the preprint didn't do.

Is there any possibility that trace plots could become a feature in metafor?

Best,
Samuel

library(metafor)
traceplot <- function(object, upperQuantile = 99.9) {
    levelSeq <- seq(0, upperQuantile, length.out = 100)
    ## confidence curve
    cc <- sapply(X = levelSeq, FUN = function(level) {
        confint.rma.uni(object = object, level = level)$random[2,2:3]
    })
    lower <- cc[1,]
    upper <- cc[2,]
    ## effect sizes
    tauseq <- seq(0, max(upper), length.out = 100)
    estimates <- t(sapply(X = tauseq, FUN = function(tau) {
        fit <- rma.uni(yi = object$yi, vi = object$vi, tau2 = tau^2)
        est <- fit$beta[1,1]
        se <- fit$se
        ci <- est + c(-1, 1)*qnorm(p = 1 - fit$level)*se
        c("estimate" = est, "lower" = ci[1], "upper" = ci[2],
          "blup" = blup(fit)$pred)
    }))

    oldpar <- par(no.readonly = TRUE)
    on.exit(par(oldpar))
    layout(mat = matrix(c(1, 1, 2), ncol = 1))
    ## upper plot
    mar1 <- oldpar$mar + c(0, 0, 0, 2)
    mar1[1] <- 0
    par(mar = mar1)
    blupInd <- seq(4, ncol(estimates))
    blupcols <- adjustcolor(col = seq(1, length(blupInd)), alpha = 0.7)
    plot(x = tauseq, y = estimates[,1], type = "n", bty = "n", xaxt = "n", xlab = NA,
         ylab = "Effect size", las = 1, panel.first = grid(ny = NA, lty = 2),
         ylim = c(min(estimates), max(estimates)))
    polygon(x = c(tauseq, rev(tauseq)), y = c(estimates[,2], rev(estimates[,3])),
            col = adjustcolor(col = 1, alpha.f = 0.1), border = NA)
    lines(x = tauseq, y = estimates[,1], lwd = 2, lty = 1)
    matlines(x = tauseq, y = estimates[,blupInd], type = "l", lty = 1, col = blupcols)
    mtext(text = blupInd - 3, side = 4, at = estimates[nrow(estimates),blupInd],
          line = -1, col = blupcols, las = 1, cex = 0.75)
    mtext(text = bquote(bold("overall")), side = 4, at = estimates[nrow(estimates),1],
          line = -1, col = 1, las = 1, cex = 0.75)

    ## lower plot
    mar2 <- oldpar$mar + c(0, 0, 0, 2)
    mar2[3] <- 0
    par(mar = mar2)
    pSeq <- 100 - levelSeq # p-value function
    plot(x = c(rev(lower), upper), y = c(rev(pSeq), pSeq), type = "l",
         lwd = 1.5, xlab = bquote("Heterogeneity" ~ tau),
         xlim = c(0, max(upper)), ylab = NA, yaxt = "n", ylim = c(0, 120),
         panel.first = grid(ny = NA, lty = 2), bty = "n")
    pbks <- seq(0, 100, 10)
    axis(side = 2, at = pbks, las = 1)
    mtext(text = bquote(italic(P) * "-value (%)"), side = 2, line = 2.8, cex = 0.7)
    axis(side = 4, at = pbks, labels = 100 - pbks, las = 1)
    mtext("Confidence (%)", side = 4, line = 2.8, cex = 0.7)
    CItau <- confint(object = object)$random[2,2:3]
    arrows(x0 = CItau[1], x1 = CItau[2], y0 = object$level*100, angle = 90,
           code = 3, length = 0.05)
}

### example
dat <- escalc(measure = "RR", ai = tpos, bi = tneg, ci = cpos, di = cneg,
              data = dat.bcg)
object <- rma(yi, vi, data = dat, method = "REML", level = 95)
traceplot(object = object)

sessionInfo()

#> R version 4.3.1 (2023-06-16)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Debian GNU/Linux 12 (bookworm)
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.21.so;  LAPACK version 3.11.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> time zone: Europe/Zurich
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] metafor_4.2-0       numDeriv_2016.8-1.1 metadat_1.2-0      
#> [4] Matrix_1.6-0       
#> 
#> loaded via a namespace (and not attached):
#> [1] compiler_4.3.1  tools_4.3.1     mathjaxr_1.6-0  nlme_3.1-162   
#> [5] grid_4.3.1      lattice_0.20-45

[wviechtb]

Hi Samuel. Yes, I am aware of the paper. I think it would make even more sense to add the profile likelihood for tau^2 at the bottom (as profile() provides). In any case, I will add the addition of trace plots to my to-do list and report back here when I have gotten around to this (could take a while).