[R-meta] Confidence intervals for I^2 levels with rma.mv
Andreas Voldstad
@ndre@@@vo|d@t@d @end|ng |rom ke||ogg@ox@@c@uk
Mon May 6 13:59:54 CEST 2024
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.
[[alternative HTML version deleted]]
More information about the R-sig-meta-analysis
mailing list