[R-meta] confidence interval for I^2 in multilevel model
Michael Dewey
lists at dewey.myzen.co.uk
Sat Sep 9 14:35:07 CEST 2017
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
More information about the R-sig-meta-analysis
mailing list