[R-meta] Question on Forest Plot by subgroups
Gabriel Cotlier
g@b|k|m01 @end|ng |rom gm@||@com
Sat Jun 24 10:10:50 CEST 2023
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 using maastrichtuniversity.nl> wrote:
> 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
>
> >-----Original Message-----
> >From: Gabriel Cotlier [mailto:gabiklm01 using gmail.com]
> >Sent: Friday, 23 June, 2023 16:12
> >To: Viechtbauer, Wolfgang (NP)
> >Cc: R Special Interest Group for Meta-Analysis
> >Subject: Re: [R-meta] Question on Forest Plot by subgroups
> >
> >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()
>
-------------- 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>
More information about the R-sig-meta-analysis
mailing list