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

[Viechtbauer Wolfgang (SP)]

Hi Michael,

Let me address your questions (from your mail). Indeed, I was a bit hasty in my original post. So, the goal is compute I^2-type statistics for a multilevel meta-analytic model and corresponding CIs. Here is an example:

library(metafor)
dat <- get(data(dat.konstantopoulos2011))

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

W <- diag(1/dat$vi)
X <- model.matrix(res)
P <- W - W %*% X %*% solve(t(X) %*% W %*% X) %*% t(X) %*% W

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

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

So, total I^2 is computed with:

(sigma^2.1 + sigma^2.2) / (sigma^2.1 + sigma^2.2 + s^2)

and the I^2 values for each level with:

sigma^2.1 / (sigma^2.1 + sigma^2.2 + s^2)
sigma^2.2 / (sigma^2.1 + sigma^2.2 + s^2)

where s^2 is the 'typical' sampling variance.

In principle, one could also consider another set of I^2 values, namely:

sigma^2.1 / (sigma^2.1 + s^2)
sigma^2.2 / (sigma^2.2 + s^2)

and for these, it is easy to get 95% CIs as I described. That is:

### alternative I^2 values for each level
100 * res$sigma2 / (res$sigma2 + (res$k-res$p)/sum(diag(P)))

### 95% CIs
sav <- confint(res)
100 * sav[[1]]$random[1,2:3] / (sav[[1]]$random[1,2:3] + (res$k-res$p)/sum(diag(P)))
100 * sav[[2]]$random[1,2:3] / (sav[[2]]$random[1,2:3] + (res$k-res$p)/sum(diag(P)))

But these I^2 values do not add up to the total I^2 and it does make more sense to use 'sigma^2.1 + sigma^2.2 + s^2' in the denominator (so that it reflects the sum of all variance components). Unfortunately, getting CIs for:

sigma^2.1 / (sigma^2.1 + sigma^2.2 + s^2)
sigma^2.2 / (sigma^2.1 + sigma^2.2 + s^2)

is a bit more tricky, since these values are a function of more than one variance component (confint() only gives us CIs for each component, not pairs of components -- in which case we would need 2d confidence regions).

Another approach we could take here is parametric bootstrapping. There was a slight error in simulate.rma(), so if you want to try this, make sure you install the devel version as described here:

http://www.metafor-project.org/doku.php/installation#development_version

Okay, so now you can do:

sim <- simulate(res, nsim=10000)
sav <- lapply(sim, function(x) {
   tmp <- try(rma.mv(x, vi, random = ~ 1 | district/school, data=dat), silent=TRUE) 
   if (inherits(tmp, "try-error")) {
      next
   } else {
      tmp
   }})

I2s <- sapply(sav, function(x) 100 * x$sigma2[1] / (sum(x$sigma2) + (res$k-res$p)/sum(diag(P))))
quantile(I2s, c(.025, .975))

I2s <- sapply(sav, function(x) 100 * x$sigma2[2] / (sum(x$sigma2) + (res$k-res$p)/sum(diag(P))))
quantile(I2s, c(.025, .975))

You will have to be a bit patient if you use nsim=10000. This will give you 95% parametric bootstrap (percentile) CIs for the two I^2 statistics. For the total I^2, we would use:

I2s <- sapply(sav, function(x) 100 * sum(x$sigma2) / (sum(x$sigma2) + (res$k-res$p)/sum(diag(P))))
quantile(I2s, c(.025, .975))

For the total I^2, we can take another approach by fitting the same model as above, but parameterized in a different way (see: http://www.metafor-project.org/doku.php/analyses:konstantopoulos2011):

### multivariate parameterization
res <- rma.mv(yi, vi, random = ~ factor(school) | district, data=dat)
res

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

### CI for total I^2
sav <- confint(res)
100 * sav[[1]]$random[1,2:3] / (sav[[1]]$random[1,2:3] + (res$k-res$p)/sum(diag(P)))

The parametric bootstrap CI and the approach with confint() (which gives us profile likelihood CIs) yield quite similar CIs for the total I^2 in this case.

I hope this clarifies things.

Best,
Wolfgang

-----Original Message-----
From: R-sig-meta-analysis [mailto:r-sig-meta-analysis-bounces at r-project.org] On Behalf Of Michael Dewey
Sent: Saturday, 09 September, 2017 14:35
To: r-sig-meta-analysis at r-project.org
Subject: [R-meta] confidence interval for I^2 in multilevel model

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