Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Change confidence interval format in package metafor forest graph?

Assume the code below (as given in Viechtbauer, 2010):

library(metafor)
data("dat.bcg", package = "metafor")
dat <- escalc(measure = "RR", ai = tpos, bi = tneg, ci = cpos, di = cneg, data = dat.bcg, append = TRUE)
res <- rma(ai = tpos, bi = tneg, ci = cpos, di = cneg, data = dat, measure = "RR")
forest(res, slab = paste(dat$author, dat$year, sep = ", "), xlim = c(-16, 6), at = log(c(0.05, 0.25, 1, 4)), atransf = exp, ilab = cbind(dat$tpos, dat$tneg, dat$cpos, dat$cneg), ilab.xpos = c(-9.5, -8, -6, -4.5), cex = 0.75)
op <- par(cex = 0.75, font = 2)
text(c(-9.5, -8, -6, -4.5), 15, c("TB+", "TB-", "TB+", "TB-"))
text(c(-8.75, -5.25), 16, c("Vaccinated", "Control"))
text(-16, 15, "Author(s) and Year", pos = 4)
text(6, 15, "Relative Risk [95% CI]", pos = 2)
par(op)

This gives a forest graph as below:

Forest plot

So how can I change the format of confidence intervals in the graph? Is it possible to replace brackets with parentheses and use "to" instead of ","? How about using "-" or long hypen instead of ","? This should change i.e. [0.13, 1.26] to (0.13 to 1.26) or (0.13 – 1.26)

like image 755
M. Er Avatar asked Aug 10 '16 04:08

M. Er


People also ask

What is Yi in Metafor?

yi|θi ∼ N(θi,vi). In other words, the observed effect sizes or outcomes are assumed to be unbiased and normally dis- tributed estimates of the corresponding true effects/outcomes with sampling variances equal to vi. (where vi is just the square of the standard errors of the estimates).

What is the effect size in forest plot?

1) of the effect size: studies with a larger weight are given a larger square, while studies with a lower weight have a smaller square. Conventionally, a forest plot should also contain the effect size data that was used to perform the meta-analysis. This provides others with the data needed to replicate our results.

What data do you need for a forest plot?

Additional data to the forest plot are (1) n, (2) mean and standard deviation, (3) P of each primary study, (4) meta-analysis overall 95% confidence interval, and (5) heterogeneity test result (with statistical significance level).


2 Answers

You need to do some hacking of the code for forest.rma. Several steps:

After displaying the current version of the code by typing the function name:

 forest.rma   # Copy the name and the code and paste into the console 
              #  Add an assignment operator `<-`
              # leave off the bytecode and environment notations at the bottom

Or you can do this in an editor, which would probably be the preferred method since you might then want to save this code to a .Rprofile file.

1) Add parameters to the argument list:

forest.rma <- 
function (x, annotate = TRUE, addfit = TRUE, addcred = FALSE, 
    showweights = FALSE, xlim, alim, clim, ylim, at, steps = 5, 
    level = x$level, digits = 2, refline = 0, xlab, slab, mlab, 
    ilab, ilab.xpos, ilab.pos, order, transf, atransf, targs, 
    rows, efac = 1, pch = 15, psize, col, border, lty, cex, cex.lab, 
    cex.axis, annosep = " , ", bkt=c("[", "]"), ...) 
{  #  ....not showing all the _long_ function body 
   # Scroll down to almost the bottom of the function body

2) Find and change arguments to the annotext cbind-assignment. There are several places where annotext might get constructed, but only one of them matches your "format target". Find the one that looks like this:

# annotext <- cbind(annotext[, 1], " [ ", annotext[, 
#                2], " , ", annotext[, 3], " ]")

Change to this:

annotext <- cbind(annotext[, 1], bkt[1], annotext[, 
                 2], annosep, annotext[, 3], bkt[2] )
# hit enter to get the modification to hold in your workspace

3) Now assign the correct environment to the function so it can play well with its siblings:

environment(forest.rma) <- environment(forest.default)
# if you forget this step you see this error:

Error in forest.rma(res, slab = paste(dat$author, dat$year, sep = ", "), : could not find function ".setlab"

And call it with the new arguments of your choosing:

png(); forest(res, slab = paste(dat$author, dat$year, sep = ", "), xlim = c(-16, 6), at = log(c(0.05, 0.25, 1, 4)), atransf = exp, ilab = cbind(dat$tpos, dat$tneg, dat$cpos, dat$cneg), ilab.xpos = c(-9.5, -8, -6, -4.5), cex = 0.75, annosep=" to ", bkt = c( "(", ")" ) )
op <- par(cex = 0.75, font = 2)
text(c(-9.5, -8, -6, -4.5), 15, c("TB+", "TB-", "TB+", "TB-"))
text(c(-8.75, -5.25), 16, c("Vaccinated", "Control"))
text(-16, 15, "Author(s) and Year", pos = 4)
text(6, 15, "Relative Risk [95% CI]", pos = 2)
dev.off()

enter image description here

like image 152
IRTFM Avatar answered Sep 28 '22 18:09

IRTFM


Here is a solution that does not require changing the code of the forest.rma() function. I use annotate=FALSE so that the function does not annotate the forest plot and instead add those annotations myself.

library(metafor)
data("dat.bcg", package = "metafor")
dat <- escalc(measure = "RR", ai = tpos, bi = tneg, ci = cpos, di = cneg, data = dat.bcg, append = TRUE)

res <- rma(ai = tpos, bi = tneg, ci = cpos, di = cneg, data = dat, measure = "RR")

### note the use of annotate=FALSE in forest()
forest(res, slab = paste(dat$author, dat$year, sep = ", "), xlim = c(-16, 6), 
       at = log(c(0.05, 0.25, 1, 4)), atransf = exp, 
       ilab = cbind(dat$tpos, dat$tneg, dat$cpos, dat$cneg), ilab.xpos = c(-9.5, -8, -6, -4.5), 
       cex = 0.75, annotate=FALSE) 

op <- par(cex = 0.75, font = 2)
text(c(-9.5, -8, -6, -4.5), 15, c("TB+", "TB-", "TB+", "TB-"))
text(c(-8.75, -5.25), 16, c("Vaccinated", "Control"))
text(-16, 15, "Author(s) and Year", pos = 4)
text(6, 15, "Relative Risk [95% CI]", pos = 2)

### add annotations manually
tmp <- summary(dat, transf=exp)[,c("yi","ci.lb","ci.ub")] ### for the individual studies
tmp <- rbind(tmp, with(predict(res, transf=exp), c(pred, ci.lb, ci.ub))) ### add model estimate and CI bounds
sav <- apply(tmp, 2, formatC, format="f", digits=2)
annotext <- apply(sav, 1, function(x) {paste0(x[1], " (", x[2], " to ", x[3], ")")})
text(6, c(res$k:1, -1), annotext, pos=2, font=1)

par(op)
like image 43
Wolfgang Avatar answered Sep 28 '22 18:09

Wolfgang