Hi, I need help with metafor package, please.
I have read this:
http://www.metafor-project.org/doku.php/plots:forest_plot_with_subgroups
and I adopted a code to my needs from this website.
and this:
Under the first link, there is very nice description how to do metanalysis with sub-groups - this is my question is all about:
How to set a rows argument of order argument of forest() funtion from metafor package in order to plot forest plot, that is the main plot in metaanalysis ?
I have read this as well:
https://wviechtb.github.io/metafor/reference/forest.default.html
and this:
https://cran.r-project.org/web/packages/metafor/metafor.pdf
and this:
I tried to use the code from first link and I adopted it to my needs, meaning I want to recreate this metaanalysis in metafor in R:
https://www.mdpi.com/2079-6382/6/4/21
Here is what I have done so far:
library(metafor)
library(tidyverse)
Diarr_data <- structure(list(Study_and_year = c("Vanderhoof_1999", "Arvola_1999",
"Erdeye_2004", "Erdeye_2_2004", "Duman_2005", "Cindoruk_2007",
"Zojaji_2013", "De_Vrese_2011", "Chatterjee_2013"), Events_treatment = c(7,
3, 7, 7, 14, 9, 11, 4, 19), Total_treatment = c(93, 61, 127,
117, 204, 62, 80, 30, 198), Events_Control = c(25, 9, 12, 30,
28, 19, 24, 3, 26), Total_Control = c(95, 58, 105, 117, 185,
62, 80, 29, 198), Probiotics_kind = c("L_rhamnosus_GG", "L_rhamnosus_GG",
"S_boulardii", "S_boulardii", "S_boulardii", "S_boulardii", "S_boulardii",
"L_acidophilus_L", "L_acidophilus_L"), Relative_Risk = c(0.286,
0.317, 0.482, 0.233, 0.453, 0.474, 0.458, 1.289, 0.731), lowerRR = c(0.13,
0.09, 0.197, 0.107, 0.246, 0.233, 0.241, 0.316, 0.418), upperRR = c(0.629,
1.113, 1.181, 0.51, 0.834, 0.964, 0.872, 5.265, 1.277), SE = c(0.127,
0.261, 0.251, 0.103, 0.15, 0.186, 0.161, 1.262, 0.219), Categorical_moderator = c(1,
1, 2, 2, 2, 2, 2, 3, 3), ai = c(7, 3, 7, 7, 14, 9, 11, 4, 19),
bi = c(86, 58, 120, 110, 190, 53, 69, 26, 179), ci = c(25,
9, 12, 30, 28, 19, 24, 3, 26), di = c(70, 49, 93, 87, 157,
43, 56, 26, 172)), class = c("spec_tbl_df", "tbl_df", "tbl",
"data.frame"), row.names = c(NA, -9L), spec = structure(list(
cols = list(Study_and_year = structure(list(), class = c("collector_character",
"collector")), Events_treatment = structure(list(), class = c("collector_double",
"collector")), Total_treatment = structure(list(), class = c("collector_double",
"collector")), Events_Control = structure(list(), class = c("collector_double",
"collector")), Total_Control = structure(list(), class = c("collector_double",
"collector")), Probiotics_kind = structure(list(), class = c("collector_character",
"collector")), Relative_Risk = structure(list(), class = c("collector_double",
"collector")), lowerRR = structure(list(), class = c("collector_double",
"collector")), upperRR = structure(list(), class = c("collector_double",
"collector")), SE = structure(list(), class = c("collector_double",
"collector")), Categorical_moderator = structure(list(), class = c("collector_double",
"collector")), ai = structure(list(), class = c("collector_double",
"collector")), bi = structure(list(), class = c("collector_double",
"collector")), ci = structure(list(), class = c("collector_double",
"collector")), di = structure(list(), class = c("collector_double",
"collector"))), default = structure(list(), class = c("collector_guess",
"collector")), skip = 1L), class = "col_spec"))
Diarr_data$Probiotics_kind <- factor(Diarr_data$Probiotics_kind, levels = c("L_rhamnosus_GG", "S_boulardii", "L_acidophilus_L"))
class(Diarr_data$Probiotics_kind)
levels(Diarr_data$Probiotics_kind)
res <- metafor::rma(yi, vi, method = "DL" , data = Diarr_data)
### a little helper function to add Q-test, I^2, and tau^2 estimate info
mlabfun <- function(text, res) {
list(bquote(paste(.(text),
" (Q = ", .(formatC(res$QE, digits=2, format="f")),
", df = ", .(res$k - res$p),
", p ", .(metafor:::.pval(res$QEp, digits=2, showeq=TRUE, sep=" ")), "; ",
I^2, " = ", .(formatC(res$I2, digits=1, format="f")), "%, ",
tau^2, " = ", .(formatC(res$tau2, digits=2, format="f")), ")")))}
Now main function to create forest plot:
metafor::forest(res, xlim=c(-16, 4.6), at=log(c(0.05, 0.25, 1, 4)), atransf=exp,
ilab=cbind(dat$ai, dat$bi, dat$ci, dat$di),
ilab.xpos=c(-9.5,-8,-6,-4.5), cex=0.95, ylim=c(-1, 27),
order=Diarr_data$Probiotics_kind, rows = c(**I dont know what to write down here !!!!!!**),
mlab=mlabfun("DerSimonian-Laird Model for All Studies", res),
psize=1, header="Author(s) and Year")
### I don't know what to place in rows argument just above
op <- par(cex=0.95, font=2)
text(c(-9.5,-8,-6,-4.5), 26, c("Treatment+", "Treatment-", "Control+", "Control-"))
text(c(-8.75,-5.25), 27, c("Treatment", "Control"))
text(-16, c(24,16,5), pos=4, c("L_rhamnosus_GG", "S_boulardii", "L_acidophilus_L"))
res.s <- rma(yi, vi, subset=(Diarr_data$Probiotics_kind =="L_rhamnosus_GG"), data=Diarr_data)
res.r <- rma(yi, vi, subset=(Diarr_data$Probiotics_kind == "S_boulardii"), data=Diarr_data)
res.a <- rma(yi, vi, subset=(Diarr_data$Probiotics_kind =="L_acidophilus_L"), data=Diarr_data)
metafor::addpoly(res.s, row=18.5, cex=0.75, atransf=exp, mlab=mlabfun("DL Model for Subgroup", res.s))
metafor::addpoly(res.r, row= 7.5, cex=0.75, atransf=exp, mlab=mlabfun("DL Model for Subgroup", res.r))
metafor::addpoly(res.a, row= 1.5, cex=0.75, atransf=exp, mlab=mlabfun("DL Model for Subgroup", res.a))
### fit meta-regression model to test for subgroup differences
res <- rma(yi, vi, mods = ~ Probiotics_kind, data = Diarr_data)
### add text for the test of subgroup differences
text(-16, -1.8, pos=4, cex=0.75, bquote(paste("Test for Subgroup Differences: ",
Q[M], " = ", .(formatC(res$QM, digits=2, format="f")), ", df = ", .(res$p - 1),
", p = ", .(formatC(res$QMp, digits=2, format="f")))))
My desired result should be something like this:
or this:
I was able to do it a bit, but I have been making mistakes trying to properly group studies (Study_and_year - variable in Diarr_data dataframe:
Here below studies are messed up.
So again my question is, how to correctly allocate studies to relevant groups according to "Probiotics_kind" variable ?
The "rows" argument to "order" argument of metafor::forest() is responsible for that, I think.
How does Metafor count rows ? Where is a baseline ?
I would be very grateful for any help, ideas.
Thank you