[R-meta] Which method to compute 95% confidence intervals around individual effect sizes (standardized mean differences) ?

Thu Jun 3 15:42:44 CEST 2021

Hi Dakis,

It's good that you mention this (for the record), since this is indeed an important issue when using split() and I have fallen into this trap myself before.

I actually added an 'order' argument recently to the bldiag() function that can also help with this issue. See here:

https://wviechtb.github.io/metafor/reference/bldiag.html

In any case, glad to hear that this is resolved now.

Best,
Wolfgang

>-----Original Message-----
>From: Dakis-Yaoba OUEDRAOGO [mailto:dakis-yaoba.ouedraogo using mnhn.fr]
>Sent: Thursday, 03 June, 2021 15:25
>To: Viechtbauer, Wolfgang (SP)
>Cc: r-sig-meta-analysis
>Subject: Re: Which method to compute 95% confidence intervals around individual
>effect sizes (standardized mean differences) ?
>
>Dear Wolfang,
>
>thank you very much for your answer and your patience.
>I was finally able to fix the problem : there was an error in my computation of
>the matrix V.
>
>This was the error :
>In my dataset, the column "dat$study" was not coded 1, 2, 3, etc. but was a
>character vector, a combination of the article ID and the number of the experiment
>within the article.
>During the split, dat$study is automatically converted into a factor and the order
>of the factor was different than the initial order of the studies in the dataset
>
>> unique(dat$study)
>[1] "ScreenTA_13502_1_1" "ScreenTA_13502_2_2" "ScreenTA_13496_1_1"
>"ScreenTA_13496_2_2" # initial order of the studies in the dataset
>> names(split(dat, dat$study))
>[1] "ScreenTA_13496_1_1" "ScreenTA_13496_2_2" "ScreenTA_13502_1_1"
>"ScreenTA_13502_2_2" # order of the studies in the matrix V
>
>Thus, in my model, the matrix V was not associated to the right yi because the
>order of the studies was not the same.
>
>I corrected this with
>> dat$study <- factor(dat$study, levels=unique(dat$study))
>> names(split(dat, dat$study))
>[1] "ScreenTA_13502_1_1" "ScreenTA_13502_2_2" "ScreenTA_13496_1_1"
>"ScreenTA_13496_2_2"
>
>And now, all 3 methods to compute CI give the same results.
>
>Many thanks again, through your anwers I could discover and fix this big problem
>in my code.
>Best wishes,
>Dakis
>
>----- Mail original -----
>De: "Wolfgang Viechtbauer, SP" <wolfgang.viechtbauer using maastrichtuniversity.nl>
>À: "Dakis-Yaoba OUEDRAOGO" <dakis-yaoba.ouedraogo using mnhn.fr>
>Cc: "r-sig-meta-analysis" <r-sig-meta-analysis using r-project.org>
>Envoyé: Mardi 1 Juin 2021 19:35:10
>Objet: RE: Which method to compute 95% confidence intervals around individual
>effect sizes (standardized mean differences) ?
>
>Dear Dakis,
>
>Method 3 should give the same results. The yi and vi elements in 'mod' are just
>the original data -- these are not modified or changed by the model that is
>fitted.
>
>Here is a fully reproducible example (using one of the Gleser & Olkin, 2009,
>examples) showing that all three methods give the exact same CIs:
>
>library(metafor)
>
>dat <- data.frame(study=c(1,1,2,3,3,3), trt=c(1,2,1,1,2,3),
>                  ai=c( 40, 40, 10,150,150,150),
>n1i=c(1000,1000,200,2000,2000,2000),
>                  ci=c(100,150, 15, 40, 80, 50),
>n2i=c(4000,4000,400,1000,1000,1000))
>dat$pti <- with(dat, ci / n2i)
>dat$pci <- with(dat, ai / n1i)
>dat <- escalc(measure="OR", ai=ai, ci=ci, n1i=n1i, n2i=n2i, data=dat)
>
>calc.v <- function(x) {
>   v <- matrix(1/(x$n1i[1]*x$pci[1]*(1-x$pci[1])), nrow=nrow(x), ncol=nrow(x))
>   diag(v) <- x$vi
>   v
>}
>
>V <- bldiag(lapply(split(dat, dat$study), calc.v))
>
>res <- rma.mv(yi, V, mods = ~ factor(trt) - 1, data=dat)
>res
>
>cbind(dat$yi - sqrt(dat$vi)*qnorm(0.05/2, lower.tail=FALSE),
>      dat$yi + sqrt(dat$vi)*qnorm(0.05/2, lower.tail=FALSE))
>
>cbind(summary(dat)$ci.lb, summary(dat)$ci.ub)
>
>cbind(res$yi - sqrt(res$vi)*qnorm(0.05/2, lower.tail = FALSE),
>      res$yi + sqrt(res$vi)*qnorm(0.05/2, lower.tail = FALSE))
>
>Best,
>Wolfgang
>
>>-----Original Message-----
>>From: Dakis-Yaoba OUEDRAOGO [mailto:dakis-yaoba.ouedraogo using mnhn.fr]
>>Sent: Tuesday, 25 May, 2021 11:11
>>To: Viechtbauer, Wolfgang (SP)
>>Cc: r-sig-meta-analysis
>>Subject: Re: Which method to compute 95% confidence intervals around individual
>>effect sizes (standardized mean differences) ?
>>
>>Dear Wolfang,
>>
>>thank you for your answer, methods 1 and 2 gave indeed the same results, sorry I
>>did not ckeck correctly. Unfortunately confidence intervals computed from method
>1
>>and 3 are different.
>>Sometimes they are very different one being positive and the other negative
>(which
>>is my main problem).
>>Please find below an example from my data with problematic cases highlighted with
>>*!!!*
>>(cilow1 and ciup1 are lower CI computed with method 1, whereas cilow3 and ciup3
>>are upper CI computed from method 3 (rma.mv model))
>>
>>
>>    cilow1 cilow3   ciup1   ciup3
>> 1:  0.316  0.600   3.721   3.436
>> 2:  3.758  3.223  11.928  12.463
>> 3:  1.103  1.744   5.284   4.643
>> 4:  3.027  4.369   9.975   8.634
>> 5:  3.029  3.882   9.980   9.127
>> 6: 10.258 15.743  30.346  24.861
>> 7: 11.055 14.887  32.647  28.815
>> 8:  7.519  7.137  22.475  22.857
>> 9: 36.817 31.138 107.667 113.346
>>10: -0.804 -1.088   2.032   2.317
>>11:  4.374  4.909  13.613  13.079
>>
>>12: -0.581 -1.222   2.318   2.959
>>13:  1.175 -0.166   5.440   6.782 *!!!*
>>14:  1.921  1.068   7.166   8.019
>>15:  4.305 -1.180  13.422  18.907 *!!!*
>>16:  6.964  3.132  20.893  24.725
>>17:  7.929  8.311  23.649  23.267
>>18: 42.736 48.415 124.944 119.265
>>
>>Looking specifically at rows 12 to 18, this is the same experiment testing
>>increasing concentration of a chemical from row 12 to 18 with all compared to the
>>same control. When I look to the raw data, confidence intervals provided with
>>method 1 (cilow1) seem logical because there is an increasing effect of the
>>chemical with increasing concentration, whereas the lack of effect indicated by
>>cilow3 for line 15 seem very strange.
>>Based on this I think I will finally choose the method 1, but it would be very
>>nice to understand why the model gives different confidence intervals.
>>
>>Best wishes,
>>Dakis
>>
>>----- Mail original -----
>>De: "Viechtbauer, Wolfgang (SP)" <wolfgang.viechtbauer using maastrichtuniversity.nl>
>>À: "Dakis-Yaoba OUEDRAOGO" <dakis-yaoba.ouedraogo using mnhn.fr>, "r-sig-meta-analysis"
>><r-sig-meta-analysis using r-project.org>
>>Envoyé: Dimanche 23 Mai 2021 19:06:01
>>Objet: RE: Which method to compute 95% confidence intervals around individual
>>effect sizes (standardized mean differences) ?
>>
>>Dear Dakis,
>>
>>All three methods you listed should give you identical results.
>>
>>Best,
>>Wolfgang
>>
>>>-----Original Message-----
>>>From: R-sig-meta-analysis [mailto:r-sig-meta-analysis-bounces using r-project.org] On
>>>Behalf Of Dakis-Yaoba OUEDRAOGO
>>>Sent: Tuesday, 18 May, 2021 17:00
>>>To: r-sig-meta-analysis using r-project.org
>>>Subject: [R-meta] Which method to compute 95% confidence intervals around
>>>individual effect sizes (standardized mean differences) ?
>>>
>>>Dear all,
>>>
>>>I am gathering studies about effects of various chemicals on corals, and for
>>every
>>>concentration of a chemical I computed the standardized mean differences for
>>>several outcomes. My final goal is to know for which maximal concentration of a
>>>chemical no significant effect is observed, and for which minimal concentration
>>of
>>>a chemical a significant effect is observed.
>>>To get this I computed the 95% confidence intervals around the standardized mean
>>>differences to identify the effect sizes that are significantly/non
>significantly
>>>different from zero.
>>>
>>>I am quite confused about how to properly compute these condidence intervals. I
>>>could see 3 different types of CI :
>>>
>>>1/ 95% CI assuming a normal distribution
>>>data$cilow <- data$yi - sqrt(data$vi)*qnorm(0.05/2, lower.tail = FALSE)
>>>data$ciup <- data$yi + sqrt(data$vi)*qnorm(0.05/2, lower.tail = FALSE)
>>>
>>>2/ 95% Wald-type confidence intervals
>>>These are computed using the summary.escalc() function in metafor
>>>
>>>3/ The confidence intervals computed from a multi-level model rma.mv where a
>>>variance-covariance matrix V is specified to take into account that I have
>>several
>>>concentrations compared to the same control (Correction of Gleser & Olkin, I
>>>followed the tutorial in
>>>http://www.metafor-
>>>project.org/doku.php/analyses:gleser2009#quantitative_response_variable
>>>to compute V)
>>>
>>>mod <- rma.mv(yi=yi, V=V, mods= ~1, random= ~1 | ID_case, data=data,
>>>method="REML")
>>>
>>>With forest(mod) I can see the individual study 95%CI on the plot and I can get
>>>them with
>>>mod$yi - sqrt(mod$vi)*qnorm(0.05/2, lower.tail = FALSE)
>>>mod$yi + sqrt(mod$vi)*qnorm(0.05/2, lower.tail = FALSE)
>>>
>>>Because the method chosen to compute the 95% CI around the individual
>>standardized
>>>mean differences will greatly influence the conclusions about the problematic
>>>chemical concentrations, I will greatly appreciate any help, comment or advise
>on
>>>my issue.
>>>
>>>Best wishes,
>>>Dakis