Unusual problem with metafor:rma output

Question

I have a rather unusual problem with the output of my metaregression using metafor:rma.

When I attempt to calculate the confidence intervals for the tau, i, and H values, all of the estimates are lower than the lower 95% confidence interval:

       estimate   ci.lb   ci.ub 
tau^2    0.8804  1.8735  7.4653 
tau      0.9383  1.3688  2.7323 
I^2(%)  92.0197 96.0843 98.9876 
H^2     12.5309 25.5383 98.7765 

The code i used for the random effects meta-regression is:

res <- rma(yi = LogOdds, sei = SE, data = data, method = 'DL')

Has anyone encountered this before?

Edit 1 After playing around with the code, I discovered that if I change the method to method = "SJ", this problem is eliminated.

I'm guessing this has to do with how these parameters are being estimated. But why was this happening in the first place?

Edit 2: According to help(confint.rma.uni):

"Usually, the estimate of τ² from the random/mixed-effects model will fall within the confidence interval provided by the Q-profile method. However, this is not guaranteed. Depending on the method used to estimate τ² and the width of the confidence interval, it can happen that the confidence interval does not actually contain the estimate (trying to explain this to reviewers can be tricky). However, using the empirical Bayes or Paule-Mandel estimator of τ² when fitting the model (i.e., using method="EB" or method="PM") guarantees that the estimate of τ² falls within the confidence interval. When method="GENQ" was used to fit the model, the corresponding CI obtained via the generalized Q-statistic method is also guaranteed to contain the estimate τ²."

While this text certainly explains why this is happening, it still states that it is acceptable to publish an estimate that falls outside the confidence interval without changing the method. As such, the help section states that is it difficult to explain this to reviewers.

I seem to also have difficulty understanding why such a result is acceptable. Can anyone point me in a direction to further understand this concept?

Answer 1

The estimate of $\tau^2$ may not fall inside of the CI when the method used to estimate $\tau^2$ is not based on the same statistical principle as the method used to construct the CI.

For example, ML/REML estimation of $\tau^2$ is based on a different principle than the Q-profile method (Viechtbauer, 2010) to construct the CI. Usually, this doesn't cause any issues, but can occasionally lead to the estimate falling outside of the CI. Here is an example:

library(metafor)
dat <- structure(list(yi = c(-0.05, -1.86, -0.48, -1.15, -0.57, -0.64, 0.46, -0.53, -0.96, 2.71), vi = c(0.256, 0.928, 0.273, 0.287, 0.249, 0.103, 0.369, 0.061, 0.22, 0.824)), class = "data.frame", row.names = c(NA, -10L))
res <- rma(yi, vi, data=dat, method="REML")
confint(res)

gives the following results:

       estimate   ci.lb   ci.ub 
tau^2    0.0000  0.0399  4.2715 
tau      0.0013  0.1997  2.0668 
I^2(%)   0.0008 15.6439 95.2050 
H^2      1.0000  1.1855 20.8552 

The Q-profile method is more aligned with the PM/EB estimators, so switching to those should solve this issue.

res <- rma(yi, vi, data=dat, method="PM")
confint(res)

which yields:

       estimate   ci.lb   ci.ub 
tau^2    0.7790  0.0399  4.2715 
tau      0.8826  0.1997  2.0668 
I^2(%)  78.3599 15.6439 95.2050 
H^2      4.6211  1.1855 20.8552 

If you want to stick to ML/REML estimation, then a profile likelihood CI would guarantee that the estimate of $\tau^2$ falls inside of the CI. To obtain such a CI, you can use the rma.mv() function to fit the same model and then confint() gives you the profile likelihood CI:

dat$id <- 1:10
res <- rma.mv(yi, vi, random = ~ 1 | id, data=dat, method="REML")
confint(res)

In the output, the variance component is called $\sigma^2$, but this is the same as $\tau^2$ above:

        estimate  ci.lb  ci.ub 
sigma^2   0.0000 0.0000 2.0744 
sigma     0.0000 0.0000 1.4403 

If you install the 'devel' version of metafor (https://wviechtb.github.io/metafor/#installation), then there is also an undocumented feature of confint() that gives you the profile likelihood CI when the model was fitted with rma():

res <- rma(yi, vi, data=dat, method="REML")
confint(res, type="PL")

       estimate  ci.lb   ci.ub 
tau^2    0.0000 0.0000  2.0744 
tau      0.0013 0.0000  1.4403 
I^2(%)   0.0008 0.0000 90.6034 
H^2      1.0000 1.0000 10.6422 

References

Viechtbauer, W. (2007). Confidence intervals for the amount of heterogeneity in meta-analysis. Statistics in Medicine, 26(1), 37–52.