[R-meta] confidence interval for I^2 in multilevel model

[Michael Dewey]

Apologies for starting a new thread (with the same subject line) as I 
deleted the old one.

In a post in the thread Wolfgang posted some code to answer the 
question. When i tried it I did not get quite what I expected and so, at 
his suggestion, I am re-posting the code here for further comment.

============= cut here ===================

With confint(), you can get CIs for variance components. You can plug 
the CI bounds for a variance component into the equation for I^2 to 
obtain a CI for I^2. Here is an example:
dat <- get(data(dat.konstantopoulos2011))

res <- rma.mv(yi, vi, random = ~ 1 | district/school, data=dat)
res

sav <- confint(res)
sav

W <- diag(1/dat$vi)
X <- model.matrix(res)
P <- W - W %*% X %*% solve(t(X) %*% W %*% X) %*% t(X) %*% W
100 * res$sigma2 / (res$sigma2 + (res$k-res$p)/sum(diag(P))) ### 
district and school level I^2
100 * sav[[1]]$random[1,2:3] / (sav[[1]]$random[1,2:3] + 
(res$k-res$p)/sum(diag(P))) ### CI for district-level I^2
100 * sav[[2]]$random[1,2:3] / (sav[[2]]$random[1,2:3] + 
(res$k-res$p)/sum(diag(P))) ### CI for the school-level I^2

100 * sum(res$sigma2) / (sum(res$sigma2) + (res$k-res$p)/sum(diag(P))) 
### total I^2

But you cannot use this approach to get a CI for the total I^2. For 
this, you have to use the 'multivariate parameterization' of this model 
(see: 
http://www.metafor-project.org/doku.php/analyses:konstantopoulos2011). So:

res <- rma.mv(yi, vi, random = ~ factor(study) | district, data=dat)
res

sav <- confint(res)
sav

W <- diag(1/dat$vi)
X <- model.matrix(res)
P <- W - W %*% X %*% solve(t(X) %*% W %*% X) %*% t(X) %*% W
100 * res$tau2 / (res$tau2 + (res$k-res$p)/sum(diag(P))) ### total I^2
100 * sav$random[1,2:3] / (sav$random[1,2:3] + 
(res$k-res$p)/sum(diag(P))) ### CI for the total I^2

-- 
Michael
http://www.dewey.myzen.co.uk/home.html