[R-meta] Question on Forest Plot by subgroups
Gabriel Cotlier
g@b|k|m01 @end|ng |rom gm@||@com
Fri Jun 23 09:36:15 CEST 2023
Hello all,
I have a question and need help on how to:
- plot a forest plot by subgroups using a model defining subgroups (each
groups corresponds to positive and negative correlations respectively by
moderator (mods = ) using a categorical variable for separating both groups
named as "alloc_char" wirth values "positive" and "negative"
respectively,as mods. = ~ alloc_char -1, also including a random effect
argument for inter and intra variability as random =Article/ Sample_ID in
the following manner:
resSep <- rma.mv(yi,
vi,
mods = ~ alloc_char-1,
random = ~ 1 | Article / Sample_ID,
data=dat)
while for the total effect size at last line of the forest plot line # -1
the overall effect size of all data that includes no moderator therefor
together "positives" and "negatives" together with the model for overall
estimate estimate as follows:
res <- rma.mv(yi,
vi,
random = ~ 1 | Article / Sample_ID,
data=dat)
For this purpose I applied the following code below, however, my problem is
that since I want in the forest plot the result of resSep which separates
size effects by "positive" and "negative" in a single model with moderators
I do not know how to include these two overall model results from resSep under
the each subgroup by moderator ("negatives" and "positives") while at the
same time at line # -1 have plotted the overall results of the model for
the whole datasets (no separation by moderators into "positive" and
"negatives").
This is because the results for the model with moderator and "no
subsetting" by positive and negative give different results to the same
model subsetting by positive and negatives :
## with moderator:
resSep <- rma.mv(yi,
vi,
mods = ~ alloc_char-1,
random = ~ 1 | Article / Sample_ID,
data=dat)
## with moderator and subsetting by positives and negatives:
# MODEL FOR POSITIVES
res.s <- rma.mv(yi,
vi,
mods = ~ alloc_char -1,
random = ~ 1 | Article / Sample_ID,
subset = c(alloc_char == "Positive"),
data=dat)
# MODEL FOR NEATIVE
res.r <- rma.mv(yi,
vi,
mods = ~ alloc_char -1,
random = ~ 1 | Article / Sample_ID,
subset = c(alloc_char == "Negative"),
data=dat)
The problem is that suseting the same model ( res.s and res.r ) gives
different size effect for each subset than if the model were estimated as
single model (resSep) and I would like to have consistency in the forest
plot plotting the efect size for negatives and positives as in the single
model with moderators resSep
Please find the complete code below.
Thanks a lot for your help and guidance.
Kind regards,
Gabrel
#################################### CODE
####################################################
## celaring all
rm(list = ls())
dev.off()
library(metafor)
##################################### ---LOAD DATA---
############################################################
## Load data
library(readxl)
dat <- read_excel(file.choose())
# View(dat)
##################################### ---MODEL---
############################################################
## Transformation of Pearson's Product-moment correlation coefficients (r)
to Fisher's Z
dat <-escalc(measure = "ZCOR", ri = ri, ni = ni, data = dat)
# View(dat)
## MODEL COMPLETE DATA
res <- rma.mv(yi,
vi,
# mods = ~ alloc_char + Computational_simulation -1 , #
mods = ~ mod1 + mod2 + mod3 - 1
random = ~ 1 | Article / Sample_ID, # ~ 1 |
studyid/sampleid,
data=dat)
res
## Calculate I^2 for all dataset together
W <- diag(1/res$vi)
X <- model.matrix(res)
P <- W - W %*% X %*% solve(t(X) %*% W %*% X) %*% t(X) %*% W
I_TOT = 100 * sum(res$sigma2) / (sum(res$sigma2) +
(res$k-res$p)/sum(diag(P)))
I_TOT
# MODEL USING MOODERARTOPR FOR POSITIVE AND NEGATIVE
resSep <- rma.mv(yi,
vi,
mods = ~ alloc_char-1,
random = ~ 1 | Article / Sample_ID, # ~ 1 |
studyid/sampleid,
data=dat)
resSep
########################################### ---GRAPH---
#####################################################
png(filename="ForestComputerSimulationPositiveAndNegatives.png",
res=95, width=1300, height=880, type="cairo")
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(I_TOT, digits=1, format="f")),
"%, ",
paste(sigma^2,sep = "_", "1"), " = ",
.(formatC(res$sigma2[1], digits=3, format="f")), " ",
paste(sigma^2,sep = "_", "2"), " = ",
.(formatC(res$sigma2[2], digits=3, format="f")),")")))}
# mlabfun1 <- function(text, resSep) {
# list(bquote(paste(.(text),
# " (Q = ", .(formatC(resSep$QE, digits=2, format="f")),
# ", df = ", .(resSep$k - resSep$p),
# ", p ", .(metafor:::.pval(resSep$QEp, digits=2,
showeq=TRUE, sep=" ")), "; ",
# # I^2, " = ", .(formatC(res$I2, digits=1,
format="f")), "%, ",
# paste(sigma^2,sep = "_", "1"), " = ",
.(formatC(resSep$sigma2[1], digits=3, format="f")), " ",
# paste(sigma^2,sep = "_", "2"), " = ",
.(formatC(resSep$sigma2[2], digits=3, format="f")),")")))}
forest(res,
top = 2,
# xlim=c(-8, 6),
# alim = c(-3.36077, 5.181815),
alim = c(-2, 6),
showweights = TRUE,
ilab=cbind(ni), # content of extra columns of information
ilab.xpos=c(-4.6), # position of extra columns of information
cex=1,
ylim = c(-1, 22.5), # rows in the y-axis from min-row to max-row
with text.
slab = paste(dat$Authors, dat$Year, sep = ", "),
order = alloc_char,
rows=c(2:13,17:18), # this divide by groups the studies (3 groups) =
aggregate 4 spaces between lines
mlab=mlabfun("RE Model", res),,
efac = 0.5, # size of the polygon (the diamond)
)
### set font expansion factor (as in forest() above) and use a bold font
op <- par(cex=0.75, font=2)
### add additional column headings to the plot
text(c(-10.8), res$k+6, c("Author(s) and Year"), cex = 1.4, font = 2)
text(c(-4.6), res$k+6, c("Samples size"), cex = 1.4, font = 2)
text(c(8.3), res$k+6, c("Weight"),cex = 1.4, font=2)
text(c(8.3), res$k+5.3, c("%"),cex = 1.4, font=2)
text(c(10.8), res$k+6, c("Fisher's zr [95% CI]"),cex = 1.4, font=2)
text(c(res$k+9), c("Computer Simulation Dataset"), cex = 2, font = 2)
### switch to bold italic font
par(font=4)
### add text for the subgroups
text(-12.3, c(14,19), pos=4, c("Positive",
"Negative"), cex = 1.4)
### set par back to the original settings
par(op)
### fit random-effects model in the three subgroups
# MODEL FOR POSITIVES
res.s <- rma.mv(yi,
vi,
mods = ~ alloc_char -1,
random = ~ 1 | Article / Sample_ID,
subset = c(alloc_char == "Positive"),
data=dat)
# calculate Positive Heterogeneity
W <- diag(1/res.s$vi)
X <- model.matrix(res.s)
P <- W - W %*% X %*% solve(t(X) %*% W %*% X) %*% t(X) %*% W
I_POS = 100 * sum(res.s$sigma2) / (sum(res.s$sigma2) +
(res.s$k-res.s$p)/sum(diag(P)))
I_POS
# MODEL FOR NEATIVE
res.r <- rma.mv(yi,
vi,
mods = ~ alloc_char -1,
random = ~ 1 | Article / Sample_ID,
subset = c(alloc_char == "Negative"),
data=dat)
# calculate Negative Heterogeneity
W <- diag(1/res.r$vi)
X <- model.matrix(res.r)
P <- W - W %*% X %*% solve(t(X) %*% W %*% X) %*% t(X) %*% W
I_NEG = 100 * sum(res.r$sigma2) / (sum(res.r$sigma2) +
(res.r$k-res.r$p)/sum(diag(P)))
I_NEG
### add summary polygons for the three subgroups
addpoly(res.s, row= 16, mlab="" ,efac = 0.5, col='blue')
addpoly(res.r, row= 1, mlab="", efac = 0.5, col='green')
# text(-12, -1, pos = 4, cex = 0.90, bquote(paste("RE Model (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(I_TOT, digits=1, format="f")), "%, ",
# paste(sigma^2,sep = "_",
"1"), " = ", .(formatC(res$sigma2[1], digits=3, format="f")), " ",
# paste(sigma^2,sep = "_",
"2"), " = ", .(formatC(res$sigma2[2], digits=3, format="f"))
# )))
text(-12.5, 1, pos = 4, cex = 0.90, bquote(paste("RE Model (Q = ",
.(formatC(res.r$QE, digits=2, format="f")),
", df = ", .(res.r$k -
res.r$p),
", p = ",
.(metafor:::.pval(res.r$QEp, digits=2, showeq=TRUE, sep=" ")), "; ",
I^2, " = ",
.(formatC(I_NEG, digits=1, format="f")), "%, ",
paste(sigma^2,sep = "_",
"1"), " = ", .(formatC(res.r$sigma2[1], digits=3, format="f")), " ",
paste(sigma^2,sep = "_",
"2"), " = ", .(formatC(res.r$sigma2[2], digits=3, format="f"))
)))
text(-12.5, 16, pos = 4, cex = 0.90, bquote(paste("RE Model (Q = ",
.(formatC(res.s$QE, digits=2, format="f")),
", df = ", .(res.s$k -
res.s$p),
", p = ",
.(metafor:::.pval(res.s$QEp, digits=2, showeq=TRUE, sep=" ")), "; ",
I^2, " = ",
.(formatC(I_POS, digits=1, format="f")), "%, ",
paste(sigma^2,sep = "_",
"1"), " = ", .(formatC(res.s$sigma2[1], digits=3, format="f")), " ",
paste(sigma^2,sep = "_",
"2"), " = ", .(formatC(res.s$sigma2[2], digits=3, format="f"))
)))
dev.off()
On Wed, Jun 21, 2023 at 7:24 AM Gabriel Cotlier <gabiklm01 using gmail.com> wrote:
> Hello Wolfgang,
> Thanks a lot!
> I am happy to participate on the mailing list and share coding experience
> in meta-analysis with metafor.
> Kind regards,
> Gabriel
>
>
> On Tue, Jun 20, 2023 at 3:07 PM Viechtbauer, Wolfgang (NP) <
> wolfgang.viechtbauer using maastrichtuniversity.nl> wrote:
>
>> Dear Gabriel,
>>
>> Glad to hear you figured things out and nice to report this back to the
>> list.
>>
>> Best,
>> Wolfgang
>>
>> >-----Original Message-----
>> >From: R-sig-meta-analysis [mailto:
>> r-sig-meta-analysis-bounces using r-project.org] On
>> >Behalf Of Gabriel Cotlier via R-sig-meta-analysis
>> >Sent: Tuesday, 20 June, 2023 9:04
>> >To: R Special Interest Group for Meta-Analysis
>> >Cc: Gabriel Cotlier
>> >Subject: Re: [R-meta] Question on Forest Plot by subgroups
>> >
>> >Hello all,
>> >
>> >Yesterday I posted a solution for the plotting forest plot by subgroups
>> for
>> >my particular case and I was actually making the mistake when choosing
>> the
>> >variable for the argument "order" which I corrected in the previously
>> >posted code (found my own mistake in the argument "order").
>> >
>> >But I also posted that was redundant information in my case to use
>> >subset=(alloc
>> >== 1) for the subgroups models: res.s, res.r and res.a and indeed it is
>> if
>> >I use "alloc" as follows:
>> >
>> >res.s <- rma(yi, subset=(alloc == 1), vi, data=dat)
>> >res.r <- rma(yi, subset=(alloc == 1), vi, data=dat)
>> >res.a <- rma(yi, vi, subset=(alloc == 1), data=dat)
>> >
>> >However, I found out that this gives the same overall result for each
>> >sub-group as for the total data.
>> >Therefore I correct it herein as follows: to get at the end of each
>> >subgroup its particular overall results the use subset= c() is needed.
>> >It worked for me as follows:
>> >
>> >### fit random-effects model in the three subgroups
>> >res.s <- rma(yi, vi, subset= c(Type == "Rainfed"), data=dat)
>> >res.r <- rma(yi, vi, subset = c(Type == "Experimental"), data=dat)
>> >res.a <- rma(yi, vi, subset = c(Type == "Computational simulation"),
>> >data=dat)
>> >
>> >In this way it will give the overall model results for each subgroup
>> >correctly together with the overall result for the complete dataset at
>> the
>> >end.
>> >
>> >Sorry for the confusion.
>> >Kind regards,
>> >Gabriel
>> >
>> >On Mon, Jun 19, 2023 at 5:49 PM Gabriel Cotlier <gabiklm01 using gmail.com>
>> wrote:
>> >
>> >> Hello all,
>> >>
>> >> I found the solution to the problem to the code I plotted previously.
>> The
>> >> problem was that I was ordering by "alloc" which is a column of sing()
>> of
>> >> correlations and not by the column that divide the groups (3 groups
>> >> expressed in a column named "Type" : Computational
>> >> simulation","Experimental","Natural Rainfed") which in my case is named
>> >> "Type",after correcting this I also remove unnecessary redundant
>> >> information subset=(alloc == 1) in each of the models res.s, res.r and
>> >> res.a. Please find the corrected working code below with highlighted
>> >> locations where modifications were made :
>> >>
>> >> ############################## CODE
>> >> ###########################################
>> >>
>> >> pdf(file="ForestPlotCompleteDataSetRahul_Test.pdf", width = 12, height
>> =
>> >> 35)
>> >>
>> >> 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")), ")")))}
>> >>
>> >> forest(res,
>> >> top = 2,
>> >> #xlim=c(-8, 6),
>> >> #alim = c(-3.36077, 5.181815),
>> >> alim = c(-3.8, 6),
>> >> #at=log(c(0.05, 0.25, 1, 4)),
>> >> #atransf=exp,
>> >> #width = 4,
>> >> showweights = TRUE,
>> >> ilab=cbind(ni), # content of extra columns of information
>> >> ilab.xpos=c(-4.6), # position of extra columns of information
>> >> cex=0.75,
>> >> ylim=c(-1, 161), # rows in the y-axis from min-row to max-row
>> with
>> >> text.
>> >> slab = paste(dat$Authors, dat$Year, sep = ", "),
>> >> order=Type,
>> >> rows=c(2:15,19:85,89:156), # this divide by groups the studies
>> (3
>> >> groups)
>> >> mlab=mlabfun("RE Model for All Studies", res),
>> >> #psize= 1 ,
>> >> efac = 0.2, # size of the polygon (the diamond)
>> >> # header="Author(s) and Year")
>> >> )
>> >> # abline(h=res$k+8, lwd=2, col="white")
>> >>
>> >> op <- par(cex=0.75, font=2)
>> >>
>> >> ### add additional column headings to the plot
>> >> text(c(-15), res$k+10, c("Author(s) and Year"), cex = 1.2, font = 2)
>> >> text(c(-4.6), res$k+10, c("Samples size"), cex = 1.2, font = 2)
>> >> text(c(9.7), res$k+10, c("Weight"),cex = 1.2, font=2)
>> >> text(c(9.7), res$k+9.3, c("%"),cex = 1.2, font=2)
>> >> text(c(12.6), res$k+10, c("Fisher's zr [95% CI]"),cex = 1.2, font=2)
>> >> text(c(res$k+13), c("Complete Dataset"), cex = 2, font = 2)
>> >>
>> >> par(font=4)
>> >>
>> >> ### add text for the subgroups
>> >> text(-17.2, c(16,86,157), pos=4, c("Computational simulation",
>> >> "Experimental",
>> >> "Natural Rainfed"), cex = 1.3)
>> >>
>> >> par(op)
>> >>
>> >> ### fit random-effects model in the three subgroups
>> >> res.s <- rma(yi, vi, data=dat)
>> >> res.r <- rma(yi, vi, data=dat)
>> >> res.a <- rma(yi, vi, data=dat)
>> >>
>> >> ### add summary polygons for the three subgroups
>> >> addpoly(res.s, row= 88, mlab=mlabfun("RE Model for Subgroup",
>> res.s),efac
>> >> = 0.2, col='red')
>> >> addpoly(res.r, row= 18, mlab=mlabfun("RE Model for Subgroup",
>> res.r),efac
>> >> = 0.2, col='green')
>> >> addpoly(res.a, row= 1, mlab=mlabfun("RE Model for Subgroup",
>> res.a),efac =
>> >> 0.2, col='blue')
>> >>
>> >>
>> >> dev.off()
>> >>
>> >>
>> >>
>> >>
>> >>
>> >> On Mon, Jun 19, 2023 at 3:46 PM Gabriel Cotlier <gabiklm01 using gmail.com>
>> >> wrote:
>> >>
>> >>> Hello all,
>> >>>
>> >>> I have taken an example from here
>> >>> <
>> https://www.metafor-project.org/doku.php/plots:forest_plot_with_subgroups
>> >
>> >and
>> >>> modified it to my case study of Fisher's z correlations meta-analysis
>> with
>> >>> the code below.
>> >>> I would like to plot a forest plot in which each of the subgroups ( 3
>> >>> groups expressed in a column named "Type" : Computational
>> >>> simulation","Experimental","Natural Rainfed") each showing the studies
>> >>> (rows) with positive correlations and those with negative correlations
>> >>> divided in agreement to the line of no effect located at zero. I have
>> used
>> >>> a column named "alloc" which has a 1 for those correlations (rows)
>> with
>> >>> positive value and -1 for those with negative value, however for some
>> >>> reason it does not work as expected. It mixes the lines from
>> >>> different groups, ploting them where they do not correspond and also
>> does
>> >>> not plot the model summary for the second and third group.
>> >>> I tried to change several arguments to solve it, however
>> unsuccessfully.
>> >>>
>> >>> What could be the error?
>> >>> Thanks a lot for your help.
>> >>>
>> >>> ############################## CODE
>> >>> ###########################################
>> >>>
>> >>> pdf(file="ForestPlotCompleteDataSetbyGroups.pdf", width = 12, height
>> = 35)
>> >>>
>> >>> 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")), ")")))}
>> >>> forest(res,
>> >>> top = 2,
>> >>> # xlim=c(-8, 6),
>> >>> # alim = c(-3.36077, 5.181815),
>> >>> alim = c(-3.8, 6),
>> >>> # at=log(c(0.05, 0.25, 1, 4)),
>> >>> # atransf=exp,
>> >>> # width = 4,
>> >>> showweights = TRUE,
>> >>> ilab=cbind(ni), # content of extra columns of information
>> >>> ilab.xpos=c(-4.6), # position of extra columns of information
>> >>> cex=0.75,
>> >>> ylim=c(-1, 161), # rows in the y-axis from min-row to max-row
>> with
>> >>> text.
>> >>> slab = paste(dat$Authors, dat$Year, sep = ", "),
>> >>> order=alloc,
>> >>> rows=c(2:15,19:85,89:156), # this divide by groups the studies
>> (3
>> >>> groups) from column named "Type"
>> >>> mlab=mlabfun("RE Model for All Studies", res),
>> >>> #psize= 1 ,
>> >>> efac = 0.2, # size of the polygon (the diamond)
>> >>> # header="Author(s) and Year")
>> >>> )
>> >>>
>> >>> op <- par(cex=0.75, font=2)
>> >>>
>> >>> text(c(-15), res$k+10, c("Author(s) and Year"), cex = 1.2, font = 2)
>> >>> text(c(-4.6), res$k+10, c("Samples size"), cex = 1.2, font = 2)
>> >>> text(c(9.7), res$k+10, c("Weight"),cex = 1.2, font=2)
>> >>> text(c(9.7), res$k+9.3, c("%"),cex = 1.2, font=2)
>> >>> text(c(12.6), res$k+10, c("Fisher's zr [95% CI]"),cex = 1.2, font=2)
>> >>> text(c(res$k+13), c("Complete Dataset"), cex = 2, font = 2)
>> >>>
>> >>> par(font=4)
>> >>>
>> >>> text(-17.2, c(16,86,157), pos=4, c("Computational simulation",
>> >>> "Experimental",
>> >>> "Natural Rainfed"), cex = 1.3)
>> >>> par(op)
>> >>>
>> >>> ### fit random-effects model in the three subgroups
>> >>> res.s <- rma(yi, vi, subset=(alloc == 1), data=dat)
>> >>> res.r <- rma(yi, vi, subset=(alloc == 2), data=dat)
>> >>> res.a <- rma(yi, vi, subset=(alloc == 3), data=dat)
>> >>>
>> >>> ### add summary polygons for the three subgroups
>> >>> addpoly(res.s, row= 88, mlab=mlabfun("RE Model for Subgroup",
>> res.s),efac
>> >>> = 0.2, col='red')
>> >>> addpoly(res.r, row= 18, mlab=mlabfun("RE Model for Subgroup",
>> res.r),efac
>> >>> = 0.2, col='green')
>> >>> addpoly(res.a, row= 1, mlab=mlabfun("RE Model for Subgroup",
>> res.a),efac
>> >>> = 0.2, col='blue')
>> >>>
>> >>> dev.off()
>>
>
[[alternative HTML version deleted]]
More information about the R-sig-meta-analysis
mailing list