# [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