Skip to content

[R-meta] Question on Forest Plot by subgroups

14 messages · Wolfgang Viechtbauer, Gabriel Cotlier

Original
Gabriel Cotlier
# 
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()
Gabriel Cotlier
# 
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 at gmail.com> wrote:

            

      
      
      
    

        

  
  


  


      

    

    

      
      

        


  

        

      Gabriel Cotlier
    

    

      
      #
    

  


  

      
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 at gmail.com> wrote:

        

            

      
      
      
    

        

  
  


  


      

    

    

      
      

        


  

        

      Wolfgang Viechtbauer
    

    

      
      #
    

  


  

      
Dear Gabriel,

Glad to hear you figured things out and nice to report this back to the list.

Best,
Wolfgang

            

      
      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Gabriel Cotlier
    

    

      
      #
    

  


  

      
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 at maastrichtuniversity.nl> wrote:

        

            

      
      
      
    

        

  
  


  


      

    

    
2 days later

    

      
      

        


  

        

      Gabriel Cotlier
    

    

      
      #
    

  


  

      
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 at gmail.com> wrote:

        

            

      
      
      
    

        

  
  


  


      

    

    

      
      

        


  

        

      Gabriel Cotlier
    

    

      
      #
    

  


  

      
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 at gmail.com>
wrote:

            

      
      
      
    

        

  
  


  


      

    

    

      
      

        


  

        

      Gabriel Cotlier
    

    

      
      #
    

  


  

      
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 at gmail.com>
wrote:

            

      
      
      
    

        

  
  


  


      

    

    

      
      

        


  

        

      Gabriel Cotlier
    

    

      
      #
    

  


  

      
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 at gmail.com> wrote:

        

            

      
      
      
    

        

  
  


  


      

    

    

      
      

        


  

        

      Wolfgang Viechtbauer
    

    

      
      #
    

  


  

      
Dear Gabriel,

I did not look at the details of your code as it is a bit much to digest on a Friday afternoon.

However, if you have a model that contains moderators and you want to add polygons to a forest plot based on such a model, you need to use predict() with 'newmods' set in such a way that you get the desired predicted values. The results can then fed into addpoly(). Here is a simpler example to illustrate this:

library(metafor)

dat <- escalc(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg, 
              data=dat.bcg, slab=paste0(author, ", ", year))

# standard random-effects model
res <- rma(yi, vi, data=dat)
res

# set up the forest plot
forest(res, header=TRUE, top=2, shade="zebra", ylim=c(-5,15))

# meta-regression model using 'alloc' as a moderator
res <- rma(yi, vi, mods = ~ 0 + alloc, data=dat)

# obtain the predicted effects for each level of 'alloc'
pred <- predict(res, newmods=diag(3))

# add these predicted effects as polygons to the plot
addpoly(pred, rows=-3, mlab=c("Alternate", "Random", "Systematic"))

# add a horizontal line
abline(h=-2)

Hopefully, this gives you an idea how you can do something analogous in your example.

Best,
Wolfgang

            

      
      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Gabriel Cotlier
    

    

      
      #
    

  


  

      
Hello Wolfgang,
Thanks a lot for the explanation and the example code.
Indeed gives me a good idea, I will give it a try to the example code with
my data now.
Have a nice weekend!!
Kind regards,
Gabriel



On Fri, Jun 23, 2023 at 5:13?PM Viechtbauer, Wolfgang (NP) <

        
wolfgang.viechtbauer at maastrichtuniversity.nl> wrote:

        

            

      
      
      
    

        

  
  


  


      

    

    

      
      

        


  

        

      Wolfgang Viechtbauer
    

    

      
      #
    

  


  

      
I only looked at the highlighted part (where you use addpoly()) and this is also possible (assuming that the coefficients you pull out from resSep are the estimates you are interested in).

Best,
Wolfgang

            

      
      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Gabriel Cotlier
    

    

      
      #
    

  


  

      
Hello Wolfgang,

Sorry for the delay in  replay. With regards to the use :

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')

I wanted to pool out the effect size estimated for each moderator in the
model together with their CI for upper and lower bounds:

#                      estimate      se     zval    pval    ci.lb    ci.ub

# alloc_charNegative   -0.8346  0.1491  -5.5969  <.0001  -1.1269  -0.5424
 ***
# alloc_charPositive    4.9600  0.2166  22.9006  <.0001   4.5355   5.3845
 ***

Which I could do with the code above.

However, I do not understand why, using  sei =.se$resSep (in this
particular case) I got the CI for upper and lower bounds (two last columns
in the upper list above), all correctly pooled out, although --if do not
understand incorrectly--  the "se" value  is supposed to be the Standard
Error (se) and not the "ci.lb" and "ci.ub", is this correct?

Yes, indeed the use of pred(), seems a good and very practical option to
pool out same data, so I tried it as :

 addpoly(resSep$pred[1], sei = c(sepSpe$ci.lib[1] , resSep$ci.ub), row= 16,
mlab="" , efac = 0.5, col='blue')

# for the size effect:
resSep$pred[1]

# for the CI upper and lower bound :
c(sepSpe$ci.lib[1] , resSep$ci.ub)

with the following:

pred <- predict(resSep, newmods=diag(2)) # in the model I used there are
two effect size(s)
pred
#      pred     se   ci.lb   ci.ub   pi.lb  pi.ub
# 1 -0.8346 0.1491 -1.1269 -0.5424 -1.6830 0.0137
# 2  4.9600 0.2166  4.5355  5.3845  4.0575 5.8625

The I used pred in addpoly() function:

addpoly(pred$pred[1], sei= c(pred$ci.lb[1], pred$ci.ub[1]), row= 16,
mlab="" , efac = 0.5, col='blue')
Error in addpoly.default(pred$pred[1], sei = c(pred$ci.lb[1],
pred$ci.ub[1]),  :
  Length of 'vi' (or 'sei') does not match length of 'x'.

However something went wrong since I got the error above. What could I have
done wrong ?

In addition, I wanted to ask something maybe a bit more "theoretical" (if
it is possible to define it in this way...?), with regards to my
understanding of the model's argument. I ask this to make sure I am
choosing the correct model in agreement to the hypothesis I want to test.

*-   for rma.mv <http://rma.mv>() "random" argument :*

As far I have understood, the use of rma.mv() argument "random = " will
enable me to influence the model by setting the underlying assumption that
it can be added a random effects structure in agreement  with either:

a.  influencing particular separate on single  id as:  random = ~ 1 | id (for
one id )  or  for several Ids with random = list(~ 1 | id1, ~ 1 | id2)) where
the outcome obtained with different values of the id variable are assumed
to be independent and will split the variance (eg. sigma1, sigma2 etc.)
obtaining one sigma per id.

b.  Whereas for assuming dependency the nested random effect is used which
will produce "a kind of" branching effect with formulas like for
instance random
= ~ 1 | id1/id2. For example, give an id = Article (i.e., research studies)
and another id = Sample_ID (i.e., outputs of each research study --in my
case the "*r* correlations"), then a dependency is assumed since each
records in id  = Article (study) can have many id = Sample_ID ( *r*
correlations).  In addition and as a consequence this influence is observed
in the obtained effect size as well as in the fact that the variance
(sigma) is splitted to sigma1 and sigma2. Using for instance the model
notation with the formula:

res <- rma.mv(yi,
              vi,
              random = ~ 1 | Article / Sample_ID,
              data=dat)

 due to assumed dependency "a kind of" branching effect occurs as I attempt
to show in the diagram below:

[image: dep_indep.jpg]

*-   for rma.mv <http://rma.mv>() "mods" argument :*

Now, If I would like also to model the influences of different effect size
in agreement to a give vector ( i.g a column in an input data table ) then
using moderators with for instance a notation formula such as  mods = ~
mod1 - 1 ,  the influence of "grouping" by moderators' attributes (e.g.
when a moderator is chosen as vector named "Method" that has as attributes:
"method1", "method2", and "method3'') in agreement to the moderator vector
attributes ( as I understood always a categorical variable, type = char )
it is assumed  that outcome of a effect size will be influence by the
moderator and observable in a effect size displayed in the model's
summary().The model's output will therefore display efect sizes per moderator
according to moderators' attributes namely  "method1", "method2", and
"method3" ( 3 effect sizes).  In this fashion, broadly speaking, it can be
said that "mods" argument will control the influence of different
moderators in the models, giving as output an efect size for
each moderator,a kind of" grouping effect of the effect sizes. This can be
exemplified for instance with the model 's notation formula such as:

res1 <- rma.mv(yi,
                 vi,
                 mods = ~ Method -1,
                 random = ~ 1 | Article / Sample_ID,
                 data=dat)

Now, my question is:
If I would be interested in modelling the influence of a moderator which
has the binary attributes of either being "positive" or "negative", then I
should consider using a vector as categorical variable (type = char) for
instance one named as "alloc_sign"  (which take positive for "positive" for
positive correlations and "negative" for negative correlations )  then to
model the effect of *r *correlation's sign I would possibly might use the
formula: mods = ~ alloc_sign -1, in the model above, could be this correct?

Thanks a lot for your guidance and help.
Kind regards,
Gabriel



On Fri, Jun 23, 2023 at 5:21?PM Viechtbauer, Wolfgang (NP) <

        
wolfgang.viechtbauer at maastrichtuniversity.nl> wrote:

        

            

      
      
      
    

        
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://stat.ethz.ch/pipermail/r-sig-meta-analysis/attachments/20230624/40e9788f/attachment-0001.html>

-------------- next part --------------
A non-text attachment was scrubbed...
Name: dep_indep.jpg
Type: image/jpeg
Size: 224771 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-meta-analysis/attachments/20230624/40e9788f/attachment-0001.jpg>

  
  


  


      

    

    

      
      

        


  

        

      Gabriel Cotlier
    

    

      
      #
    

  


  

      
Hello

I have previously mistakenly mentioned that lines  :

addpoly(coef(resSep)[[3]], sei = resSep$se[3], row= 88, mlab="",efac = 0.2,
col='red')
addpoly(coef(resSep)[[2]], sei = resSep$se[2], row= 18, mlab="",efac = 0.2,
col='green')
addpoly(coef(resSep)[[1]], sei = resSep$se[1] ,row= 1,  mlab="", efac =
0.5, col='blue')

worked properly however I realized  that they do not plot the diamond
polygon in its correct size with respect to its height (= estimate) and
its length ( I suppose that is in function of CI length).

The correct code is, which I am using now is:

p <- predict(resSep, newmods=diag(3))
p
addpoly(p, row = c(1, 18, 88), mlab="",efac = 0.2, col = c("red", "green",
"blue"))

Regards,
Gabriel

On Sat, Jun 24, 2023 at 11:10?AM Gabriel Cotlier <gabiklm01 at gmail.com>
wrote:

            

      
      
      
    

        
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://stat.ethz.ch/pipermail/r-sig-meta-analysis/attachments/20230624/f98d0135/attachment-0001.html>

-------------- next part --------------
A non-text attachment was scrubbed...
Name: dep_indep.jpg
Type: image/jpeg
Size: 224771 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-meta-analysis/attachments/20230624/f98d0135/attachment-0001.jpg>