[R-meta] Question on Forest Plot by subgroups
Gabriel Cotlier
g@b|k|m01 @end|ng |rom gm@||@com
Fri Jun 23 16:11:57 CEST 2023
Hello all,
Sorry for the many emails but here is the solution that finally works
nicely !!
The issue was that I just needed to define the addpoly() function arguments
properly as in the code below highlighted in yellow.
With this code I could produce a forest plot combining two different models
:
1. one general for all dataset obtaining its estimate and CI (at line # -1)
2. one by grouping using 2 moderators obtaining 2 estimates (one per
moderator which are particular variables of interest of the meta-analys)
and their CI.
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(I_TOT, 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),
# 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=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 devide by groups the studeies (3 groups)
=ageregar 4 esopacios entre lineas
mlab=mlabfun("RE Model", res),
# mlab="",
# psize= 1 ,
efac = 0.5, # size of the polygon (the diamond)
# header="Author(s) and Year")
)
# abline(h=res$k+8, lwd=2, col="white")
### 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)
### add summary polygons for the three subgroups
addpoly(coef(resSep)[[2]], sei = resSep$se[2], row= 16, mlab="" , efac =
0.5, col='blue')
addpoly(coef(resSep)[[1]], sei = resSep$se[1] ,row= 1, mlab="", efac =
0.5, col='green')
# text(x=-1, y=-2:80, labels=c(-2:80), col = "red" )
text(-12.5, 1, pos = 4, cex = 0.90, bquote(paste("RE Model (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(I_NEG, 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"))
)))
text(-12.5, 16, pos = 4, cex = 0.90, bquote(paste("RE Model (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(I_POS, 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"))
)))
text(8.5, 1, pos = 4, cex = 0.90, round(coef(resSep)[[1]],3))
text(8.5, 16, pos = 4, cex = 0.90, round(coef(resSep)[[2]],3))
dev.off()
On Fri, Jun 23, 2023 at 4:14 PM Gabriel Cotlier <gabiklm01 using gmail.com> wrote:
> Hello all,
>
> I think I solve partially the problem since I can plot the forest graph
> with the subgroup size effect in it together with the model size effect of
> the complete dataset however it is missing to plot the
> diamond corresponding to the size effect of the groups for both
> moderartor generated groups in resSep using addpoly() funcion that gives
> the following error:
>
> Error in addpoly.rma(resSep, row = 16, mlab = "", efac = 0.5, col =
> "blue") :
> Fitted model should not contain moderators.
>
> Error in addpoly.rma(resSep, row = 1, mlab = "", efac = 0.5, col =
> "green") :
> Fitted model should not contain moderators.
>
> Does anyone know how the diamonds polygon for the size effect could be
> plotted ?
> Thanks a lot.
> Regards,
> Gabriel
>
> #########################################-- code for
> 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(I_TOT, 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),
> # at=log(c(0.05, 0.25, 1, 4)),
> # atransf=exp,
> # width = 4,
> showweights = TRUE,
> ilab=cbind(ni), # content of extra columns of informaiton
> ilab.xpos=c(-4.6), # position of extra columns od 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 devide by groups the studeies (3 groups)
> =ageregar 4 esopacios entre lineas
> mlab=mlabfun("RE Model", res),
> # mlab="",
> # psize= 1 ,
> efac = 0.5, # size of the polygon (the diamond)
> # header="Author(s) and Year")
> )
> # abline(h=res$k+8, lwd=2, col="white")
>
> ### 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)
>
> ### add summary polygons for the three subgroups
> addpoly(resSep, row= 16, mlab="" , efac = 0.5, col='blue')
> addpoly(resSep, row= 1, mlab="", efac = 0.5, col='green')
>
> # text(x=-1, y=-2:80, labels=c(-2:80), col = "red" )
>
> text(-12.5, 1, pos = 4, cex = 0.90, bquote(paste("RE Model (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(I_NEG, 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"))
> )))
>
>
> text(-12.5, 16, pos = 4, cex = 0.90, bquote(paste("RE Model (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(I_POS, 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"))
> )))
>
> text(10.8, 1, pos = 4, cex = 0.90, round(coef(resSep)[[1]],3))
> text(10.8, 16, pos = 4, cex = 0.90, round(coef(resSep)[[2]],3))
>
> dev.off()
>
>
>
> On Fri, Jun 23, 2023 at 12:05 PM Gabriel Cotlier <gabiklm01 using gmail.com>
> wrote:
>
>> Hello all,
>>
>> I've been thinking and maybe a possible solution could be to have two
>> models one for all together positives and negatives (with no moderator) to
>> better capture overall effect size of the dataset :
>>
>> res <- rma.mv(yi,
>> vi,
>> random = ~ 1 | Article / Sample_ID,
>> data=dat)
>>
>> and another model using moderator for better capture effect size between
>> positives and negatives :
>>
>> resSep <- rma.mv(yi,
>> vi,
>> mods = ~ alloc_char-1,
>> random = ~ 1 | Article / Sample_ID,
>> data=dat)
>>
>> and then when reporting (plotting in the forest plot) the model's summary
>> of each models use mlab =mala" " in both forest() and addpoly() and
>> instead use text() function to text the summary for both the whole
>> dataset model res the by subgroups model based on moderators named resSep while
>> for the size effect for each model maybe possible to be reported in the
>> plot also using again text() at the correspondent location/ position
>> line for the total data set (res) and each of the subgroup ( resSep ) by
>> extracting the value of the efect size corresponding to each moderator with
>> coef(resSep)[1] and coef(resSep) [2].
>>
>> My question in this case would be how to place the extracted size effect
>> at its correspondent location using the function text()?
>>
>> Maybe could this be a possible solution?
>> Thanks a lot.
>> Kind regards,
>> Gabriel
>>
>> On Fri, Jun 23, 2023 at 10:36 AM Gabriel Cotlier <gabiklm01 using gmail.com>
>> wrote:
>>
>>> 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