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

Mon, May 6, 2024 4:59 AM #

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,

           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
Please don?t feel obliged to read or respond to my email outside your own working hours.

Wolfgang Viechtbauer

Mon, May 6, 2024 6:33 AM #

Dear Andreas,

Please see below.

Best,
Wolfgang

-----Original Message-----
From: R-sig-meta-analysis <r-sig-meta-analysis-bounces at 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 at r-project.org>
Cc: Andreas Voldstad <andreas.voldstad at 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,