[R-meta] Question on Forest Plot by subgroups
Gabriel Cotlier
g@b|k|m01 @end|ng |rom gm@||@com
Sat Jun 24 15:47:12 CEST 2023
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 using gmail.com>
wrote:
> 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/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>
More information about the R-sig-meta-analysis
mailing list