[R-meta] Confidence intervals for I^2 levels with rma.mv

[Viechtbauer, Wolfgang (NP)]

Dear Andreas,

Please see below.

Best,
Wolfgang

> -----Original Message-----
> From: R-sig-meta-analysis <r-sig-meta-analysis-bounces using r-project.org> On Behalf
> Of Andreas Voldstad via R-sig-meta-analysis
> Sent: Monday, May 6, 2024 14:00
> To: Andreas Voldstad via R-sig-meta-analysis <r-sig-meta-analysis using r-project.org>
> Cc: Andreas Voldstad <andreas.voldstad using kellogg.ox.ac.uk>
> Subject: [R-meta] Confidence intervals for I^2 levels with rma.mv
>
> Hi all,
>
> I would like to follow up an older thread on getting confidence intervals for
> the I^2:
> https://stat.ethz.ch/pipermail/r-sig-meta-analysis/2017-September/000195.html
>
> I am conducting a meta-analysis of SMD on a dataset with many sources of
> dependency (romantic couples, timepoints, multiple groups, multiple outcome
> measures). I'd like to report the 95% CI for the I^2.
>
> I have tried applying the method outlined in the previous thread. It works with
> the konstantopoulos2011 example data. However, when applying it to my own data,
> the result is that the confidence intervals are equal to the point estimate.
>
> I am wondering if the bootstrapped confidence interval method does not work in
> the same way when inputting a V-matrix in rma.mv?
>
> To obtain CIs for the total I^2, I have used confint() on a reparametrised
> model, as suggested in the thread. This seems to work out fine. But I am still
> missing the CIs for the separate levels.
>
> Unfortunately, I am not skilled in simulating data to create a reproducible
> example. I am enclosing the code I used for model fitting below in case my error
> is in any way apparent from that:
>
> library(metafor)
>
> # mod_allES <- read.csv("EffectSizeData/mod_allES.csv")
> rhoall7 <- .7
>
> V_7 <- with(mod_allES,
>             clubSandwich::impute_covariance_matrix(vi = v,
>                                                    cluster = StudyID,
>                                                    r = rhoall7,
>                                                    return_list = T))
>
> res <- rma.mv(yi = delta,
>               V = V_7,
>               slab = StudyID,
>               data = mod_allES,
>               random = ~ 1 | StudyID/esID,
>               test = "t",
>               method = "REML") #
>
> W <- diag(1/res$vi)
> X <- model.matrix(res)
> P <- W - W %*% X %*% solve(t(X) %*% W %*% X) %*% t(X) %*% W
>
> sim <- simulate(res, nsim=10000)
>
> sav <- lapply(sim, function(x) {
>   tmp <- try(
>     rma.mv(yi = delta,

Here lies the issue. You keep fitting the model to the original data, not the simulated data. This should be:

rma.mv(x,

>            V = V_7,
>            slab = StudyID,
>            data = mod_allES,
>            random = ~ 1 | StudyID/esID,
>            test = "t",
>            method = "REML"),
>     silent=TRUE)
>
>   if (inherits(tmp, "try-error")) {
>     next
>   } else {
>     tmp
>   }})
>
> I2s1 <- sapply(sav, function(x) 100 * x$sigma2[1] / (sum(x$sigma2) + (res$k-
> res$p)/sum(diag(P))))
> I2s2 <- sapply(sav, function(x) 100 * x$sigma2[2] / (sum(x$sigma2) + (res$k-
> res$p)/sum(diag(P))))
>
> quantile(I2s1, c(.025, .975))
>
> # With the result:
> #   2.5%    97.5%
> #   82.73102 82.73102
>
> quantile(I2s2, c(.025, .975))
> #   2.5%    97.5%
> #   14.90151 14.90151
>
> Best wishes,
>
> Andreas Voldstad (he/him)
> PhD student in Psychiatry
> University of Oxford