# [R-meta] Question on Forest Plot by subgroups

Gabriel Cotlier g@b|k|m01 @end|ng |rom gm@||@com
Fri Jun 23 15:14:33 CEST 2023

```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)
)
# 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)

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)
>>
>>
>> ############################################################
>>
>> # 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)
>>
>> 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)
>>>> >>
>>>> >> 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]]

```