Skip to content

[R-meta] Appending a risk-of-bias traffic-light plot to a 'three-level' forest plot

5 messages · Joshua Bernal, Wolfgang Viechtbauer

Joshua Bernal
# 
Greetings Dr. Viechtbauer and members,

I ran a three-level meta-analysis with rma.mv (due to multiple
measures of a construct within-studies). I seek help on how to append
a risk-of-bias traffic-light plot to an extension of the three-level
forest plot, please find the query details below...

Goal A (metafor):
Create the three-level forest plot with each study's "aggregated"
effect estimate and CI.
(Source A: stat.ethz.ch/pipermail/r-sig-meta-analysis/2019-February/001423.html)

Goal A Issues:
1) Source A presents 2 ways to obtain "aggregated" effect sizes, with
the second way getting "exactly the same results as those from the
full model". Because the present rma.mv analysis uses the
t-distribution (test = "t"), the p-values and summary CIs obtained are
slightly different/larger. Is it possible to get the exact same
results as when test = "t" for the rma.mv model? To see what this
looks like, please find the metafor reprex below.
2) After running the reprex, you will see that Study 6 contributes one
effect size only, its CI from the "standard" rma.mv forest plot (where
all effect sizes in each study are shown) vs its CI from the
"aggregated" forest plot are different. If the argument
"tau2=res$sigma2[2]" is omitted from the code of the "aggregated"
forest plot, then the obtained summary CI for Study 6 then matches
with that of the "standard" plot. Which would be considered the
correct approach in this case?

metafor reprex:
###
library(metafor)

dat <- structure(list(author = c("Study 1", "Study 1", "Study 1", "Study 1",
"Study 1", "Study 2", "Study 2", "Study 2", "Study 2", "Study 3",
"Study 3", "Study 4", "Study 4", "Study 5", "Study 5", "Study 6"
), study = c(1, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 4, 4, 5, 5, 6),
    outcome = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14,
    15, 16), yi = c(0.52, 0.256, 0.5003, 0.3131, -0.1019, 0.3591,
    0.0341, 0.0357, 0.3942, -0.1503, 0.2671, 0.131, 0.1664, 1.1761,
    0.5799, 0.406), vi = c(0.00926, 0.00199, 0.01171, 0.0171,
    0.00524, 0.0712, 0.0971, 0.0823, 0.0912, 0.1204, 0.3536,
    0.0655, 0.0607, 0.9571, 0.6617, 0.5139)), row.names = c(NA,
16L), class = "data.frame")

res <- rma.mv(yi, vi, random = ~ 1 | author/outcome, data = dat, test
= "t", method = "REML", slab = author)
res
forest(res)

sav <- lapply(split(dat, dat$author), function(x) {
    res <- rma(yi, vi, data=x, tau2=res$sigma2[2])
    return(c(coef(res), vcov(res)))
})

sav <- data.frame(do.call(rbind, sav))
names(sav) <- c("yi", "vi")

forest(sav$yi, sav$vi, slab=rownames(sav))

rma(yi, vi, data=sav)
res
###

Goal B (robvis, development version 0.3.0.900):
Create a risk-of-bias traffic-light plot using the "Generic" (or aka
"ROB1") template.
(Source B: https://github.com/mcguinlu/robvis  ("Generic")  OR
cran.r-project.org/web/packages/robvis/vignettes/Introduction_to_robvis.html#rob1
 ("ROB1"))

Goal B Issues:
After running the robvis reprex below, you will see that an extra
judgment label ("Critical" or dark red circle with exclamation mark)
is automatically added/cannot be removed (although only 5 judgment
labels are relevant in the dataset, I am forced to specify 6 judgement
labels). I am aware there is an argument "judgement_levels" under
development, but do not know if it is relevant and if so how to take
it from here...

robvis reprex:
###
library(robvis)

rp <- structure(list(Study = c("Study 1", "Study 2", "Study 3", "Study 4",
"Study 5", "Study 6"), Bias.due.to.confounding. = c("Not applicable",
"Not applicable", "Not applicable", "Some concerns", "Low", "No information"
), Bias.in.selection.of.participants.into.the.study. = c("Not applicable",
"Not applicable", "Not applicable", "Some concerns", "Low", "Low"
), Bias.in.classification.of.interventions. = c("Not applicable",
"Not applicable", "Not applicable", "Some concerns", "Low", "Low"
), Bias.arising.from.the.randomization.process. = c("Some concerns",
"Low", "Low", "Not applicable", "Not applicable", "Not applicable"
), Bias.due.to.deviations.from.intended.intervention. = c("Low",
"Some concerns", "Low", "Low", "Some concerns", "Low"),
Bias.due.to.missing.outcome.data. = c("Some concerns",
"Low", "No information", "Some concerns", "Low", "No information"
), Bias.in.measurement.of.the.outcome. = c("Some concerns", "Low",
"Some concerns", "Some concerns", "Low", "Some concerns"),
Bias.due.to.selection.of.the.reported.result. = c("Some concerns",
"Low", "High", "Some concerns", "Low", "High"), Overall = c("Some concerns",
"Low", "High", "Some concerns", "Low", "High")), row.names = c(NA,
6L), class = "data.frame")

bor <- rob_traffic_light(rp, tool="Generic", colour="cochrane",
psize=10, overall=TRUE, judgement_labels=c("","High","Some
concerns","Low","No information","Not applicable"))

bor
###

Goal C:
Append the risk-of-bias traffic-light plot (from Goal B) to the forest
plot (from Goal A).
(Source C: mcguinlu.github.io/robvis/articles/metafor.html#appending-risk-of-bias-plots-to-forest-plots-1)

Goal C Issues:
Because the present analysis considers the bias domains of both
randomized and non-randomized controlled intervention studies and risk
of bias was respectively assessed by ROB2 and ROBINS-I, the "Generic"
template is used because it allows modifying the number of domains so
as to include the relevant domains of both tools in the traffic-light
plot. However, according to "rob_tools(forest = TRUE)", robvis
currently supports only "ROB2" and "ROBINS-I" templates for the
"rob_append_to_forest()" function.

Other sub-goals:
- Modify text of header from "Estimate [95% CI]" to "SMD [95% CI]".
- Modify text of x-axis label from "Observed Outcome" to "Standardized
Mean Difference".
- To the right of "SMD [95% CI]", add a column with the header "No. of
effect sizes" to display the number of effect sizes within each study.

Many thanks in advance for your help...

Sincerely,
Josh
Wolfgang Viechtbauer
# 
Dear Josh,

See below for my responses.

Best,
Wolfgang
I would suggest to do the following:

1) Use the 'contain' method for computing the dfs:

res <- rma.mv(yi, vi, random = ~ 1 | author/outcome, data = dat, test = "t", dfs = "contain",
              method = "REML", slab = author)

2) Then use test="t" for rma():

rma(yi, vi, data=sav, test="t")

Then the results match up exactly.
The first. If you want to show a forest plot of the aggregated data, then the aggregation involves adding the within-study variance (as you have done) and hence the CI will be wider also for study 6.
Can't help with this.
Again, can't really help with robvis.

As an alternative, you could also consider creating the traffic light part of the plot yourself. Then you have full control over what is shown and how. An example where I did this can be found here:

https://www.metafor-project.org/doku.php/plots:forest_plot_revman
In forest(), see argument 'header'.
In forest(), see argument 'xlab'.
In forest(), you can use the 'textpos' argument to create space to the right of the annotations and then add this column either via the 'ilab' argument or with text(). In fact, the link above demonstrates this.
Joshua Bernal
# 
Dear Dr. Viechtbauer,

Thanks very much for your timely help! I'm delighted to learn about a
simple alternative approach for creating the traffic light part. I
would like to follow-up on the first part of the query on rma.mv /
aggregated data...

As the help page of the rma.mv function mentions: one can set
dfs="contain" (which automatically also sets test="t")

I checked to see whether the results would be the same as test="t" if
I specify dfs="contain", and whether specifying test="t" with versus
without dfs="contain" would yield the same results.

# dfs="contain"
res <- rma.mv(yi, vi, random = ~ 1 | author/outcome, data = dat, dfs =
"contain", method = "REML", slab = author)

estimate      se    zval    pval   ci.lb   ci.ub
  0.2495  0.0738  3.3828  0.0007  0.1049  0.3940  ***

# test="t"
res <- rma.mv(yi, vi, random = ~ 1 | author/outcome, data = dat, test
= "t", method = "REML", slab = author)

estimate      se    tval  df    pval   ci.lb   ci.ub
  0.2495  0.0738  3.3828  15  0.0041  0.0923  0.4067  **

# test="t" and dfs="contain"
res <- rma.mv(yi, vi, random = ~ 1 | author/outcome, data = dat, test
= "t", dfs = "contain", method = "REML", slab = author)

estimate      se    tval  df    pval   ci.lb   ci.ub
  0.2495  0.0738  3.3828   5  0.0196  0.0599  0.4391  *

I was wondering why the results differ, and importantly when would it
be appropriate or sensible to use each of these approaches?
Ultimately, I would like to know how I should determine the
appropriate method for this part of the meta-analysis and generally
what to consider when doing so (e.g., study sample size, number of
studies, number of effect estimates per study)? For example, is it
acceptable to use test="z" for a meta-analysis of eight studies with
sample sizes of 66, 50, 38, 23, 23, 18, 12, 10 versus five studies
with sample sizes of 50, 35, 28, 12, 10; or is it more sensible to be
'conservative' and use test="t" in either one (or both) cases?

Best regards,
Josh

On Fri, Feb 18, 2022 at 5:51 PM Viechtbauer, Wolfgang (SP)
<wolfgang.viechtbauer at maastrichtuniversity.nl> wrote:
Wolfgang Viechtbauer
# 
I would suggest to use test="t" with dfs="contain". There is a bit of discussion here:

https://wviechtb.github.io/metafor/reference/rma.mv.html#tests-and-confidence-intervals

Roughly, z-tests tend to be too liberal, which we can counteract by using t-tests. But the default df calculation is quite simplistic and the "contain" method is an improvement on that. Still far from perfect, but a bit better.

Best,
Wolfgang
Joshua Bernal
# 
Dear Dr. Viechtbauer,

Thank you again for your attention and help on the matter.

Best regards,
Josh

On Sat, Feb 19, 2022 at 12:59 AM Viechtbauer, Wolfgang (SP)
<wolfgang.viechtbauer at maastrichtuniversity.nl> wrote: