[R-meta] Confidence Intervals for I² values for three-level model

[Sarah Roesch]

Dear all, 

I would like to calculate confidence Intervals for I² in a three-level model.
Here, I closely followed the instructions provided by Wolfgang on this post: https://stat.ethz.ch/pipermail/r-sig-meta-analysis/2017-August/000150.html,

yet with a slightly different code because I prefer the list notation.
Thereby, I perfectly figured out I² on Level-2 and Level-3 and the according confidence intervals as well as overall I² and the corresponding confidence intervals. 
In our paper, however, we would like to report I² on Level 1. Figuring out I² on Level 1 is not the problem, but I am currently having difficulties figuring out confidence intervals
for the I², since the confint function of the metafor-package only provides confidence intervals for Level 2 and Level 3.

My code is the following:
### to calculate I² and CIs on Level-2 and Level-3
overall <- rma.mv(yi, vi, random = list(~ 1 | esID, ~ 1 | searchID), tdist= TRUE, method="REML", data=alldata) 
W <- diag(1/alldata$vi) # brings sample variance into a diganoal structure
X <- model.matrix(overall)
P <- W - W %*% X %*% solve(t(X) %*% W %*% X) %*% t(X) %*% W
I2 <- 100 * sum(overall$sigma2) / (sum(overall$sigma2) + (overall$k-overall$p)/sum(diag(P))) # total I²
I2 # This value is the Higgins I2 for the three-level model
100 * res$tau2 / (res$tau2 + (res$k-res$p)/sum(diag(P))) ### total I^2
# look at variance on both levels separately
100 * overall$sigma2 / (sum(overall$sigma2) + (overall$k-overall$p)/sum(diag(P)))
100 * confint(overall)[[1]]$random[1,2:3] /(sum(overall$sigma2) + (overall$k-overall$p)/sum(diag(P))) ### CI for es-level I^2
100 * confint(overall)[[2]]$random[1,2:3] /(sum(overall$sigma2) + (overall$k-overall$p)/sum(diag(P))) ### CI for the report-level I^2

### to calculate I² and CI for total variance
overall.innerouter <- rma.mv(yi, vi, random = ~ esID | searchID, data=alldata)
sav <- confint(overall.innerouter)
X <- model.matrix(overall.innerouter)
P <- W - W %*% X %*% solve(t(X) %*% W %*% X) %*% t(X) %*% W
100 * overall.innerouter$tau2 / (overall.innerouter$tau2 + (overall.innerouter$k-overall.innerouter$p)/sum(diag(P))) ### total I^2
100 * sav[[1]]$random[1,2:3] / (sav[[1]]$random[1,2:3] + (overall.innerouter$k-overall.innerouter$p)/sum(diag(P))) ### CI for the total I^2

Could anyone help out with that issue?

Thanks so much!

Sarah