Differences

This shows you the differences between two versions of the page.

--- plots:forest_plot_with_subgroups [2020/03/19 13:31] – Wolfgang Viechtbauer
+++ plots:forest_plot_with_subgroups [2021/05/16 16:06] – Wolfgang Viechtbauer
@@ Line 14: / Line 14: @@
 library(metafor)
-### decrease margins so the full space is used
+### copy BCG vaccine meta-analysis data into 'dat'
-par(mar=c(4,4,1,2))
+dat <- dat.bcg
-### fit random-effects model (use slab argument to define study labels)
+### calculate log risk ratios and corresponding sampling variances (and use
-res <- rma(ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg, measure="RR",
+### the 'slab' argument to store study labels as part of the data frame)
-           slab=paste(author, year, sep=", "), method="REML")
+dat <- escalc(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat,
+              slab=paste(author, year, sep=", "))
-### set up forest plot (with 2x2 table counts added; rows argument is used
+### fit random-effects model
-### to specify exactly in which rows the outcomes will be plotted)
+res <- rma(yi, vi, data=dat)
-forest(res, xlim=c(-16, 6), at=log(c(0.05, 0.25, 1, 4)), atransf=exp,
-       ilab=cbind(dat.bcg$tpos, dat.bcg$tneg, dat.bcg$cpos, dat.bcg$cneg),
+### 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")), ")")))}
+### set up forest plot (with 2x2 table counts added; the 'rows' argument is
+### used to specify in which rows the outcomes will be plotted)
+forest(res, xlim=c(-16, 4.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, ylim=c(-1, 27),
-       order=order(dat.bcg$alloc), rows=c(3:4,9:15,20:23),
+       order=dat$alloc, rows=c(3:4,9:15,20:23),
-       xlab="Risk Ratio", mlab="", psize=1, header="Author(s) and Year")
+       mlab=mlabfun("RE Model for All Studies", res),
+       psize=1, header="Author(s) and Year")
-### add text with Q-value, dfs, p-value, and I^2 statistic
+### set font expansion factor (as in forest() above) and use a bold font
-text(-16, -1, pos=4, cex=0.75, bquote(paste("RE Model for All Studies (Q = ",
+op <- par(cex=0.75, font=2)
-     .(formatC(res$QE, digits=2, format="f")), ", df = ", .(res$k - res$p),
-     ", p = ", .(formatC(res$QEp, digits=2, format="f")), "; ", I^2, " = ",
-     .(formatC(res$I2, digits=1, format="f")), "%)")))
-### set font expansion factor (as in forest() above) and use bold italic
+### add additional column headings to the plot
-### font and save original settings in object 'op'
+text(c(-9.5,-8,-6,-4.5), 26, c("TB+", "TB-", "TB+", "TB-"))
-op <- par(cex=0.75, font=4)
+text(c(-8.75,-5.25),     27, c("Vaccinated", "Control"))
+### switch to bold italic font
+par(font=4)
 ### add text for the subgroups
@@ Line 43: / Line 57: @@
                                "Random Allocation",
                                "Alternate Allocation"))
-### switch to bold font
-par(font=2)
-### add column headings to the plot
-text(c(-9.5,-8,-6,-4.5), 26, c("TB+", "TB-", "TB+", "TB-"))
-text(c(-8.75,-5.25),     27, c("Vaccinated", "Control"))
 ### set par back to the original settings
@@ Line 55: / Line 62: @@
 ### fit random-effects model in the three subgroups
-res.s <- rma(ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg, measure="RR",
+res.s <- rma(yi, vi, subset=(alloc=="systematic"), data=dat)
-             subset=(alloc=="systematic"), method="REML")
+res.r <- rma(yi, vi, subset=(alloc=="random"),     data=dat)
-res.r <- rma(ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg, measure="RR",
+res.a <- rma(yi, vi, subset=(alloc=="alternate"),  data=dat)
-             subset=(alloc=="random"), method="REML")
-res.a <- rma(ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg, measure="RR",
-             subset=(alloc=="alternate"), method="REML")
 ### add summary polygons for the three subgroups
-addpoly(res.s, row=18.5, cex=0.75, atransf=exp, mlab="")
+addpoly(res.s, row=18.5, cex=0.75, atransf=exp, mlab=mlabfun("RE Model for Subgroup", res.s))
-addpoly(res.r, row= 7.5, cex=0.75, atransf=exp, mlab="")
+addpoly(res.r, row= 7.5, cex=0.75, atransf=exp, mlab=mlabfun("RE Model for Subgroup", res.r))
-addpoly(res.a, row= 1.5, cex=0.75, atransf=exp, mlab="")
+addpoly(res.a, row= 1.5, cex=0.75, atransf=exp, mlab=mlabfun("RE Model for Subgroup", res.a))
+### fit meta-regression model to test for subgroup differences
+res <- rma(yi, vi, mods = ~ alloc, data=dat)
-### add text with Q-value, dfs, p-value, and I^2 statistic for subgroups
+### add text for the test of subgroup differences
-text(-16, 18.5, pos=4, cex=0.75, bquote(paste("RE Model for Subgroup (Q = ",
+text(-16, -1.8, pos=4, cex=0.75, bquote(paste("Test for Subgroup Differences: ",
-     .(formatC(res.s$QE, digits=2, format="f")), ", df = ", .(res.s$k - res.s$p),
+     Q[M], " = ", .(formatC(res$QM, digits=2, format="f")), ", df = ", .(res$p - 1),
-     ", p = ", .(formatC(res.s$QEp, digits=2, format="f")), "; ", I^2, " = ",
+     ", p = ", .(formatC(res$QMp, digits=2, format="f")))))
-     .(formatC(res.s$I2, digits=1, format="f")), "%)")))
-text(-16, 7.5, pos=4, cex=0.75, bquote(paste("RE Model for Subgroup (Q = ",
-     .(formatC(res.r$QE, digits=2, format="f")), ", df = ", .(res.r$k - res.r$p),
-     ", p = ", .(formatC(res.r$QEp, digits=2, format="f")), "; ", I^2, " = ",
-     .(formatC(res.r$I2, digits=1, format="f")), "%)")))
-text(-16, 1.5, pos=4, cex=0.75, bquote(paste("RE Model for Subgroup (Q = ",
-     .(formatC(res.a$QE, digits=2, format="f")), ", df = ", .(res.a$k - res.a$p),
-     ", p = ", .(formatC(res.a$QEp, digits=2, format="f")), "; ", I^2, " = ",
-     .(formatC(res.a$I2, digits=1, format="f")), "%)")))
 </code>