Hi Everyone,
I am conducting a meta-analysis using the metafor package. I am interested
in obtaining an estimate by subgroup only without showing an overall
effect. This is directly from the metafor website. How would i modify this
code to only show subgroup effects? Further, I want to show weights by
subgroup. The option showweights=TRUE does not display weights by subgroup
but by the weight of each study in comparison to all studies (and not the
subgroup). You help would be appreciated.
library
<http://stat.ethz.ch/R-manual/R-devel/library/base/html/library.html>(metafor)
### to save as png filepng
<http://stat.ethz.ch/R-manual/R-devel/library/grDevices/html/png.html>(filename="forest_plot_with_subgroups.png",
res=95, width=680, height=680, type="cairo")
### decrease margins so the full space is usedpar
<http://stat.ethz.ch/R-manual/R-devel/library/graphics/html/par.html>(mar=c
<http://stat.ethz.ch/R-manual/R-devel/library/base/html/c.html>(4,4,1,2))
### load BCG vaccine datadata
<http://stat.ethz.ch/R-manual/R-devel/library/utils/html/data.html>(dat.bcg)
### fit random-effects model (use slab argument to define study labels)
res <- rma(ai=tpos, bi=tneg, ci=cpos, di=cneg, data
<http://stat.ethz.ch/R-manual/R-devel/library/utils/html/data.html>=dat.bcg,
measure="RR",
slab=paste
<http://stat.ethz.ch/R-manual/R-devel/library/base/html/paste.html>(author,
year, sep=", "), method="REML")
### set up forest plot (with 2x2 table counts added; rows argument is
used### to specify exactly in which rows the outcomes will be plotted)
forest(res, xlim=c
<http://stat.ethz.ch/R-manual/R-devel/library/base/html/c.html>(-16,
6), at=log
<http://stat.ethz.ch/R-manual/R-devel/library/base/html/log.html>(c
<http://stat.ethz.ch/R-manual/R-devel/library/base/html/c.html>(.05,
.25, 1, 4)), atransf=exp
<http://stat.ethz.ch/R-manual/R-devel/library/base/html/exp.html>,
ilab=cbind
<http://stat.ethz.ch/R-manual/R-devel/library/base/html/cbind.html>(dat.bcg$tpos,
dat.bcg$tneg, dat.bcg$cpos, dat.bcg$cneg),
ilab.xpos=c
<http://stat.ethz.ch/R-manual/R-devel/library/base/html/c.html>(-9.5,-8,-6,-4.5),
cex=.75, ylim=c
<http://stat.ethz.ch/R-manual/R-devel/library/base/html/c.html>(-1,
27),
order
<http://stat.ethz.ch/R-manual/R-devel/library/base/html/order.html>=order
<http://stat.ethz.ch/R-manual/R-devel/library/base/html/order.html>(dat.bcg$alloc),
rows=c
<http://stat.ethz.ch/R-manual/R-devel/library/base/html/c.html>(3:4,9:15,20:23),
xlab="Relative Risk", mlab="RE Model for All
Studies", psize=1)
### set font expansion factor (as in forest() above) and use bold
italic### font and save original settings in object 'op'
op <- par
<http://stat.ethz.ch/R-manual/R-devel/library/graphics/html/par.html>(cex=.75,
font=4)
### add text for the subgroupstext
<http://stat.ethz.ch/R-manual/R-devel/library/graphics/html/text.html>(-16,
c
<http://stat.ethz.ch/R-manual/R-devel/library/base/html/c.html>(24,16,5),
pos=4, c
<http://stat.ethz.ch/R-manual/R-devel/library/base/html/c.html>("Systematic
Allocation",
"Random Allocation",
"Alternate Allocation"))
### switch to bold fontpar
<http://stat.ethz.ch/R-manual/R-devel/library/graphics/html/par.html>(font=2)
### add column headings to the plottext
<http://stat.ethz.ch/R-manual/R-devel/library/graphics/html/text.html>(c
<http://stat.ethz.ch/R-manual/R-devel/library/base/html/c.html>(-9.5,-8,-6,-4.5),
26, c
<http://stat.ethz.ch/R-manual/R-devel/library/base/html/c.html>("TB+",
"TB-", "TB+", "TB-"))text
<http://stat.ethz.ch/R-manual/R-devel/library/graphics/html/text.html>(c
<http://stat.ethz.ch/R-manual/R-devel/library/base/html/c.html>(-8.75,-5.25),
27, c
<http://stat.ethz.ch/R-manual/R-devel/library/base/html/c.html>("Vaccinated",
"Control"))text
<http://stat.ethz.ch/R-manual/R-devel/library/graphics/html/text.html>(-16,
26, "Author(s) and Year", pos=4)text
<http://stat.ethz.ch/R-manual/R-devel/library/graphics/html/text.html>(6,
26, "Relative Risk [95% CI]", pos=2)
### set par back to the original settingspar
<http://stat.ethz.ch/R-manual/R-devel/library/graphics/html/par.html>(op)
### fit random-effects model in the three subgroups
res.s <- rma(ai=tpos, bi=tneg, ci=cpos, di=cneg, data
<http://stat.ethz.ch/R-manual/R-devel/library/utils/html/data.html>=dat.bcg,
measure="RR",
subset
<http://stat.ethz.ch/R-manual/R-devel/library/base/html/subset.html>=(alloc=="systematic"),
method="REML")
res.r <- rma(ai=tpos, bi=tneg, ci=cpos, di=cneg, data
<http://stat.ethz.ch/R-manual/R-devel/library/utils/html/data.html>=dat.bcg,
measure="RR",
subset
<http://stat.ethz.ch/R-manual/R-devel/library/base/html/subset.html>=(alloc=="random"),
method="REML")
res.a <- rma(ai=tpos, bi=tneg, ci=cpos, di=cneg, data
<http://stat.ethz.ch/R-manual/R-devel/library/utils/html/data.html>=dat.bcg,
measure="RR",
subset
<http://stat.ethz.ch/R-manual/R-devel/library/base/html/subset.html>=(alloc=="alternate"),
method="REML")
### add summary polygons for the three subgroups
addpoly(res.s, row
<http://stat.ethz.ch/R-manual/R-devel/library/base/html/row.html>=18.5,
cex=.75, atransf=exp
<http://stat.ethz.ch/R-manual/R-devel/library/base/html/exp.html>,
mlab="RE Model for Subgroup")
addpoly(res.r, row
<http://stat.ethz.ch/R-manual/R-devel/library/base/html/row.html>7.5,
cex=.75, atransf=exp
<http://stat.ethz.ch/R-manual/R-devel/library/base/html/exp.html>,
mlab="RE Model for Subgroup")
addpoly(res.a, row
<http://stat.ethz.ch/R-manual/R-devel/library/base/html/row.html>1.5,
cex=.75, atransf=exp
<http://stat.ethz.ch/R-manual/R-devel/library/base/html/exp.html>,
mlab="RE Model for Subgroup")
dev.off
<http://stat.ethz.ch/R-manual/R-devel/library/grDevices/html/dev.off.html>()
[[alternative HTML version deleted]]
The code you posted is totally mangled up, but it's just what can be found here: http://www.metafor-project.org/doku.php/plots:forest_plot_with_subgroups If you don't want an overall estimate, just pass the estimates and corresponding sampling variances to the forest() function (and not the model object). Use the 'rows' argument to specify where the estimates will be placed and adjust 'ylim' so give you enough space to leave gaps for headings and the subgroup estimates. Then fit models within the subgroups (the 'subset' argument is useful here) and use addpoly() to add the subgroup estimates in the appropriate rows. With text(), you can add headings as needed. If you use weights() on each subgroup model object, you can get the subgroup weights (that add up to 100% within each subgroup). It's probably easiest to just add those values with text() in an appropriate place to the plot. Best, Wolfgang -- Wolfgang Viechtbauer, Ph.D., Statistician | Department of Psychiatry and Neuropsychology | Maastricht University | P.O. Box 616 (VIJV1) | 6200 MD Maastricht, The Netherlands | +31 (43) 388-4170 | http://www.wvbauer.com> -----Original Message----- > From: R-help [mailto:r-help-bounces at r-project.org] On Behalf Of John > Peterson > Sent: Monday, December 07, 2015 01:39 > To: r-help at r-project.org > Subject: [R] metafor package > > Hi Everyone, > > I am conducting a meta-analysis using the metafor package. I am > interested > in obtaining an estimate by subgroup only without showing an overall > effect. This is directly from the metafor website. How would i modify > this > code to only show subgroup effects? Further, I want to show weights by > subgroup. The option showweights=TRUE does not display weights by > subgroup > but by the weight of each study in comparison to all studies (and not the > subgroup). You help would be appreciated. > > library <http://stat.ethz.ch/R-manual/R- > devel/library/base/html/library.html>(metafor) > ### to save as png filepng > <http://stat.ethz.ch/R-manual/R- > devel/library/grDevices/html/png.html>(filename="forest_plot_with_subgrou > ps.png", > res=95, width=680, height=680, type="cairo") > ### decrease margins so the full space is usedpar > <http://stat.ethz.ch/R-manual/R- > devel/library/graphics/html/par.html>(mar=c > <http://stat.ethz.ch/R-manual/R-devel/library/base/html/c.html>(4,4,1,2)) > ### load BCG vaccine datadata > <http://stat.ethz.ch/R-manual/R- > devel/library/utils/html/data.html>(dat.bcg) > ### fit random-effects model (use slab argument to define study labels) > res <- rma(ai=tpos, bi=tneg, ci=cpos, di=cneg, data > <http://stat.ethz.ch/R-manual/R- > devel/library/utils/html/data.html>=dat.bcg, > measure="RR", > slab=paste > <http://stat.ethz.ch/R-manual/R- > devel/library/base/html/paste.html>(author, > year, sep=", "), method="REML") > ### set up forest plot (with 2x2 table counts added; rows argument is > used### to specify exactly in which rows the outcomes will be plotted) > forest(res, xlim=c > <http://stat.ethz.ch/R-manual/R-devel/library/base/html/c.html>(-16, > 6), at=log <http://stat.ethz.ch/R-manual/R- > devel/library/base/html/log.html>(c > <http://stat.ethz.ch/R-manual/R-devel/library/base/html/c.html>(.05, > .25, 1, 4)), atransf=exp > <http://stat.ethz.ch/R-manual/R-devel/library/base/html/exp.html>, > ilab=cbind > <http://stat.ethz.ch/R-manual/R- > devel/library/base/html/cbind.html>(dat.bcg$tpos, > dat.bcg$tneg, dat.bcg$cpos, dat.bcg$cneg), > ilab.xpos=c > <http://stat.ethz.ch/R-manual/R-devel/library/base/html/c.html>(-9.5,-8,- > 6,-4.5), > cex=.75, ylim=c > <http://stat.ethz.ch/R-manual/R-devel/library/base/html/c.html>(-1, > 27), > order <http://stat.ethz.ch/R-manual/R- > devel/library/base/html/order.html>=order > <http://stat.ethz.ch/R-manual/R- > devel/library/base/html/order.html>(dat.bcg$alloc), > rows=c <http://stat.ethz.ch/R-manual/R- > devel/library/base/html/c.html>(3:4,9:15,20:23), > xlab="Relative Risk", mlab="RE Model for All Studies", psize=1) > ### set font expansion factor (as in forest() above) and use bold > italic### font and save original settings in object 'op' > op <- par <http://stat.ethz.ch/R-manual/R- > devel/library/graphics/html/par.html>(cex=.75, > font=4) > ### add text for the subgroupstext > <http://stat.ethz.ch/R-manual/R-devel/library/graphics/html/text.html>(- > 16, > c <http://stat.ethz.ch/R-manual/R- > devel/library/base/html/c.html>(24,16,5), > pos=4, c <http://stat.ethz.ch/R-manual/R- > devel/library/base/html/c.html>("Systematic > Allocation", > "Random Allocation", > "Alternate Allocation")) > ### switch to bold fontpar > <http://stat.ethz.ch/R-manual/R- > devel/library/graphics/html/par.html>(font=2) > ### add column headings to the plottext > <http://stat.ethz.ch/R-manual/R-devel/library/graphics/html/text.html>(c > <http://stat.ethz.ch/R-manual/R-devel/library/base/html/c.html>(-9.5,-8,- > 6,-4.5), > 26, c <http://stat.ethz.ch/R-manual/R- > devel/library/base/html/c.html>("TB+", > "TB-", "TB+", "TB-"))text > <http://stat.ethz.ch/R-manual/R-devel/library/graphics/html/text.html>(c > <http://stat.ethz.ch/R-manual/R-devel/library/base/html/c.html>(-8.75,- > 5.25), > 27, c <http://stat.ethz.ch/R-manual/R- > devel/library/base/html/c.html>("Vaccinated", > "Control"))text > <http://stat.ethz.ch/R-manual/R-devel/library/graphics/html/text.html>(- > 16, > 26, "Author(s) and Year", pos=4)text > <http://stat.ethz.ch/R-manual/R-devel/library/graphics/html/text.html>(6, > 26, "Relative Risk [95% CI]", pos=2) > ### set par back to the original settingspar > <http://stat.ethz.ch/R-manual/R-devel/library/graphics/html/par.html>(op) > ### fit random-effects model in the three subgroups > res.s <- rma(ai=tpos, bi=tneg, ci=cpos, di=cneg, data > <http://stat.ethz.ch/R-manual/R- > devel/library/utils/html/data.html>=dat.bcg, > measure="RR", > subset > <http://stat.ethz.ch/R-manual/R- > devel/library/base/html/subset.html>=(alloc=="systematic"), > method="REML") > res.r <- rma(ai=tpos, bi=tneg, ci=cpos, di=cneg, data > <http://stat.ethz.ch/R-manual/R- > devel/library/utils/html/data.html>=dat.bcg, > measure="RR", > subset > <http://stat.ethz.ch/R-manual/R- > devel/library/base/html/subset.html>=(alloc=="random"), > method="REML") > res.a <- rma(ai=tpos, bi=tneg, ci=cpos, di=cneg, data > <http://stat.ethz.ch/R-manual/R- > devel/library/utils/html/data.html>=dat.bcg, > measure="RR", > subset > <http://stat.ethz.ch/R-manual/R- > devel/library/base/html/subset.html>=(alloc=="alternate"), > method="REML") > ### add summary polygons for the three subgroups > addpoly(res.s, row > <http://stat.ethz.ch/R-manual/R-devel/library/base/html/row.html>=18.5, > cex=.75, atransf=exp > <http://stat.ethz.ch/R-manual/R-devel/library/base/html/exp.html>, > mlab="RE Model for Subgroup") > addpoly(res.r, row > <http://stat.ethz.ch/R-manual/R-devel/library/base/html/row.html>> 7.5, cex=.75, atransf=exp > <http://stat.ethz.ch/R-manual/R-devel/library/base/html/exp.html>, > mlab="RE Model for Subgroup") > addpoly(res.a, row > <http://stat.ethz.ch/R-manual/R-devel/library/base/html/row.html>> 1.5, cex=.75, atransf=exp > <http://stat.ethz.ch/R-manual/R-devel/library/base/html/exp.html>, > mlab="RE Model for Subgroup") > dev.off <http://stat.ethz.ch/R-manual/R- > devel/library/grDevices/html/dev.off.html>() > > [[alternative HTML version deleted]] > > ______________________________________________ > R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting- > guide.html > and provide commented, minimal, self-contained, reproducible code.
Hi John, Please keep r-help copied on the reply. What's the 'previous model'? How do you get estimates within subgroups that 'includes the overall effect'? I really cannot follow you here. Best, Wolfgang -- Wolfgang Viechtbauer, Ph.D., Statistician | Department of Psychiatry and Neuropsychology | Maastricht University | P.O. Box 616 (VIJV1) | 6200 MD Maastricht, The Netherlands | +31 (43) 388-4170 | http://www.wvbauer.com> -----Original Message----- > From: John Peterson [mailto:john.peterson.can at gmail.com] > Sent: Monday, December 07, 2015 22:14 > To: Viechtbauer Wolfgang (STAT) > Subject: Re: [R] metafor package > > Hi Dr. Viechtbauer, > Thank you very much for your reply. I tried your advice and was able to > make a forest plot with only the estimates for each subgroup.For the > estiamte for each subgroup, similar to the previous model, I random > effects model within each subgroup. However, I now find the result for > the estiamte within subgroup to be different thant the result for the > previous model. I have tried analyzing this in STATA and I get the same > result as the model which includes the overall effect. Any advice on what > may be wrong here? Thanks greatly, > John > > On 7 December 2015 at 04:02, Viechtbauer Wolfgang (STAT) > <wolfgang.viechtbauer at maastrichtuniversity.nl> wrote: > The code you posted is totally mangled up, but it's just what can be > found here: > > http://www.metafor-project.org/doku.php/plots:forest_plot_with_subgroups > > If you don't want an overall estimate, just pass the estimates and > corresponding sampling variances to the forest() function (and not the > model object). Use the 'rows' argument to specify where the estimates > will be placed and adjust 'ylim' so give you enough space to leave gaps > for headings and the subgroup estimates. Then fit models within the > subgroups (the 'subset' argument is useful here) and use addpoly() to add > the subgroup estimates in the appropriate rows. With text(), you can add > headings as needed. > > If you use weights() on each subgroup model object, you can get the > subgroup weights (that add up to 100% within each subgroup). It's > probably easiest to just add those values with text() in an appropriate > place to the plot. > > Best, > Wolfgang > > -- > Wolfgang Viechtbauer, Ph.D., Statistician | Department of Psychiatry and > Neuropsychology | Maastricht University | P.O. Box 616 (VIJV1) | 6200 MD > Maastricht, The Netherlands | +31 (43) 388-4170 | http://www.wvbauer.com > > > -----Original Message----- > > From: R-help [mailto:r-help-bounces at r-project.org] On Behalf Of John > > Peterson > > Sent: Monday, December 07, 2015 01:39 > > To: r-help at r-project.org > > Subject: [R] metafor package > > > > Hi Everyone, > > > > I am conducting a meta-analysis using the metafor package. I am > > interested > > in obtaining an estimate by subgroup only without showing an overall > > effect. This is directly from the metafor website. How would i modify > > this > > code to only show subgroup effects? Further, I want to show weights by > > subgroup. The option showweights=TRUE does not display weights by > > subgroup > > but by the weight of each study in comparison to all studies (and not > the > > subgroup). You help would be appreciated.[snip garbled code]
Hi Dr. Viechtbauer, The code provided in the metafor projects website for subgroup includes fitting a random effects model on the entire dataset and fitting a random effects model within subgroups. When I exactly follow this code, my estimates and confidence intervals for estimate within each subgroup matches with what I get in STATA so it seems to be the correct estimate (and CI). However, I don't want to present an overall effect and I want to present only the effect within each subgroup. In my second attempt, I did not run a random effects model within the entire dataset and only ran the models within each subgroup. I generated a variable corresponding to the estimate(risk ratio) and standard error which I plugged in the forest() function (i.e. forest(rr, se, .....). When I run this code, the estimate I get for each subgroup is slightly different than the estimate I get for each subgroup in comparison to when I also included the random effects model for the overall effect. (i.e. my first attempt). Is this the right approach? I was not clear when you said passing the estimates and sampling variances to the forest() function. I created a variable corresponding to the estimate and standard error and plugged those in the forest() function. I am not sure if this is the right approach. Thanks, John On 8 December 2015 at 03:47, Viechtbauer Wolfgang (STAT) < wolfgang.viechtbauer at maastrichtuniversity.nl> wrote:> Hi John, > > Please keep r-help copied on the reply. > > What's the 'previous model'? How do you get estimates within subgroups > that 'includes the overall effect'? I really cannot follow you here. > > Best, > Wolfgang > > -- > Wolfgang Viechtbauer, Ph.D., Statistician | Department of Psychiatry and > Neuropsychology | Maastricht University | P.O. Box 616 (VIJV1) | 6200 MD > Maastricht, The Netherlands | +31 (43) 388-4170 | http://www.wvbauer.com > > > -----Original Message----- > > From: John Peterson [mailto:john.peterson.can at gmail.com] > > Sent: Monday, December 07, 2015 22:14 > > To: Viechtbauer Wolfgang (STAT) > > Subject: Re: [R] metafor package > > > > Hi Dr. Viechtbauer, > > Thank you very much for your reply. I tried your advice and was able to > > make a forest plot with only the estimates for each subgroup.For the > > estiamte for each subgroup, similar to the previous model, I random > > effects model within each subgroup. However, I now find the result for > > the estiamte within subgroup to be different thant the result for the > > previous model. I have tried analyzing this in STATA and I get the same > > result as the model which includes the overall effect. Any advice on what > > may be wrong here? Thanks greatly, > > John > > > > On 7 December 2015 at 04:02, Viechtbauer Wolfgang (STAT) > > <wolfgang.viechtbauer at maastrichtuniversity.nl> wrote: > > The code you posted is totally mangled up, but it's just what can be > > found here: > > > > http://www.metafor-project.org/doku.php/plots:forest_plot_with_subgroups > > > > If you don't want an overall estimate, just pass the estimates and > > corresponding sampling variances to the forest() function (and not the > > model object). Use the 'rows' argument to specify where the estimates > > will be placed and adjust 'ylim' so give you enough space to leave gaps > > for headings and the subgroup estimates. Then fit models within the > > subgroups (the 'subset' argument is useful here) and use addpoly() to add > > the subgroup estimates in the appropriate rows. With text(), you can add > > headings as needed. > > > > If you use weights() on each subgroup model object, you can get the > > subgroup weights (that add up to 100% within each subgroup). It's > > probably easiest to just add those values with text() in an appropriate > > place to the plot. > > > > Best, > > Wolfgang > > > > -- > > Wolfgang Viechtbauer, Ph.D., Statistician | Department of Psychiatry and > > Neuropsychology | Maastricht University | P.O. Box 616 (VIJV1) | 6200 MD > > Maastricht, The Netherlands | +31 (43) 388-4170 | http://www.wvbauer.com > > > > > -----Original Message----- > > > From: R-help [mailto:r-help-bounces at r-project.org] On Behalf Of John > > > Peterson > > > Sent: Monday, December 07, 2015 01:39 > > > To: r-help at r-project.org > > > Subject: [R] metafor package > > > > > > Hi Everyone, > > > > > > I am conducting a meta-analysis using the metafor package. I am > > > interested > > > in obtaining an estimate by subgroup only without showing an overall > > > effect. This is directly from the metafor website. How would i modify > > > this > > > code to only show subgroup effects? Further, I want to show weights by > > > subgroup. The option showweights=TRUE does not display weights by > > > subgroup > > > but by the weight of each study in comparison to all studies (and not > > the > > > subgroup). You help would be appreciated. > > [snip garbled code] > >[[alternative HTML version deleted]]