Prediction intervals for random-effects meta-analysis: A confidence distribution approach

Abstract

Prediction intervals are commonly used in meta-analysis with random-effects models. One widely used method, the Higgins–Thompson–Spiegelhalter prediction interval, replaces the heterogeneity parameter with its point estimate, but its validity strongly depends on a large sample approximation. This is a weakness in meta-analyses with few studies. We propose an alternative based on bootstrap and show by simulations that its coverage is close to the nominal level, unlike the Higgins–Thompson–Spiegelhalter method and its extensions. The proposed method was applied in three meta-analyses.

Keywords

Confidence distributions coverage properties meta-analysis prediction intervals random-effects models

1 Introduction

Meta-analysis is an important tool for combining the results of a set of related studies. A common objective of meta-analysis is to estimate an overall mean effect and its confidence interval.¹ Fixed-effect models and random-effects models have been widely applied.

Fixed-effect models assume that treatment effects are equal for all studies. The estimate of the common treatment effect and its confidence interval provide valuable information for applying the results to a future study or a study not included in the meta-analysis. By contrast, random-effects models assume that the true treatment effects differ for each study. The average treatment effect across all studies and its confidence interval have been used together with heterogeneity measures that are important for generalizability. For instance, the I²-statistic^2,3 has been widely used as a heterogeneity measure. However, researchers often interpret results of fixed-effect and random-effects models in the same manner.^4,5 They tend to focus on the average treatment effect estimate and its confidence interval. It is necessary to consider that the treatment effects in each study may be different from the average treatment effect. Higgins et al.⁶ proposed a prediction interval for a treatment effect in a future study. This interval can be interpreted as the range of the predicted true treatment effect in a new study, given the realized (past) studies. A prediction interval naturally accounts for the heterogeneity and helps apply the results to a future study or a study not included in the meta-analysis. Riley et al.⁴ recommended that a prediction interval should be reported alongside a confidence interval and heterogeneity measure.

Poor coverage of the confidence intervals in random-effects meta-analysis has been studied extensively,^6–8 especially in the context of synthesis of few studies⁹ (fewer than 20). Recently, Partlett and Riley¹⁰ confirmed that prediction intervals based on established methods, including the Higgins–Thompson–Spiegelhalter (HTS) prediction interval,⁶ also could have poor coverage. No explicit solution to this problem has been found thus far.

The HTS prediction interval has a fundamental problem. It can be regarded as a plug-in estimator that replaces the heterogeneity parameter $τ^{2}$ with its point estimate ${\hat{τ}}^{2}$ . The t distribution with K–2 degrees of freedom is used to approximately account for the uncertainty of ${\hat{τ}}^{2}$ , where K is the number of studies. Replacement with the t-approximation has a detrimental impact on the coverage probability, especially when K is small, as is often the case in practice. We also confirmed in Section 3 the HTS prediction intervals suffer from severe under-coverage.

In this article, we develop a new prediction interval that is valid under more general and realistic settings of meta-analyses in medical research, including those whose K is especially small. To avoid using a plug-in estimator, we propose a parametric bootstrap approach using a confidence distribution^11–15 to account for the uncertainty of ${\hat{τ}}^{2}$ with an exact distribution estimator of $τ^{2}$ . A confidence distribution, like a Bayesian posterior, is considered as a distribution function to estimate the parameter of interest in frequentist inference.

This article is organized as follows. In Section 2, we review the random-effects meta-analysis and HTS prediction interval, and then present the new method. In Section 3, we assess the performance of the HTS prediction interval and proposed prediction interval in simulation studies. In Section 4, we apply the developed method to three meta-analysis data sets. We conclude with a brief discussion.

2 Method

2.1 The random-effects model and the exact distribution of Cochran’s Q statistic

We consider the random-effects model.^6,16–19

Let the random variable Y_k ( $k = 1, 2, \dots, K$ ) be an effect size estimate from the kth study. The random-effects model can be defined as

\begin{matrix} Y_{k} = θ_{k} + ε_{k}, \\ θ_{k} = μ + u_{k} \end{matrix}

(1)

where θ_k is the true effect size of the kth study, μ is the grand mean parameter of the average treatment effect,

ε_{k}

is the random error within a study, and u_k is a random variable reflecting study-specific deviation from the average treatment effect. It is assumed that

ε_{k}

and u_k are independent, with

ε_{k} \sim N (0, σ_{k}^{2})

and

u_{k} \sim N (0, τ^{2})

, where the within-study variances

σ_{k}^{2}

are known and replaced by their efficient estimates,^20,21 and the across-studies variance

τ^{2}

is an unknown parameter that reflects the treatment effects heterogeneity.

Under the model in equation (1), the marginal distribution of Y_k is a normal distribution with the mean μ and the variance $σ_{k}^{2} + τ^{2}$ .

Random-effects meta-analyses generally estimate μ to evaluate the average treatment effect and $τ^{2}$ to evaluate the treatment effects heterogeneity. The average treatment effect μ is estimated by

\frac{\sum_{k = 1}^{K} (σ_{k}^{2} + {\hat{τ}}^{2})^{- 1} Y_{k}}{\sum_{k = 1}^{K} (σ_{k}^{2} + {\hat{τ}}^{2})^{- 1}}

where

{\hat{τ}}^{2}

is an estimator of the heterogeneity parameter

τ^{2}

. Estimators of

τ^{2}

, such as the DerSimonian and Laird estimator,¹⁸ have been applied.²² In this paper, we discuss prediction intervals using the DerSimonian and Laird estimator

{\hat{τ}}_{DL}^{2} = max [0, {\hat{τ}}_{UDL}^{2}]

and its untruncated version

{\hat{τ}}_{UDL}^{2} = \frac{Q - (K - 1)}{S_{1} + S_{2} / S_{1}}

where

Q = \sum_{k = 1}^{K} v_{k} (Y_{k} - \bar{Y})^{2}

is Cochran’s Q statistic,

v_{k} = σ_{k}^{- 2}, \bar{Y} = \sum_{k = 1}^{K} v_{k} Y_{k} / \sum_{k = 1}^{K} v_{k}

, and

S_{r} = \sum_{k = 1}^{K} v_{k}^{r}

for r = 1, 2. Under the model in equation (1), Biggerstaff and Jackson²¹ derived the exact distribution function of Q,

F_{Q} (q; τ^{2})

, to obtain confidence intervals for

τ^{2}

. Cochran’s Q is a quadratic form that can be written as

Y^{T} AY

, where

Y = (Y_{1}, Y_{2}, \dots, Y_{K})^{T}, A = V - {vv}^{T} / v_{+}, V = diag (v_{1}, v_{2}, \dots, v_{K}), v = (v_{1}, v_{2}, \dots, v_{K})^{T}, v_{+} = \sum_{k = 1}^{K} v_{k}

, and the superscript ‘T’ denotes matrix transposition. Here and subsequently,

Z = Σ^{- 1 / 2} (Y - μ) \sim N (0, I), S = Σ^{1 / 2} A Σ^{1 / 2}, μ = (μ, μ, \dots, μ)^{T}, Σ = diag (σ_{1}^{2} + τ^{2}, σ_{2}^{2} + τ^{2}, \dots, σ_{K}^{2} + τ^{2}), 0 = (0, 0, \dots, 0)^{T}

, and

I = diag (1, 1, \dots, 1)

Lemma 1

Under the model in equation (1), Q can be expressed as $Z^{T} SZ$ ; then Q has the same distribution as the random variable $\sum_{k = 1}^{K} λ_{k} χ_{k}^{2} (1)$ , where $λ_{k} \geq 0$ are the eigenvalues of matrix S, and $χ_{1}^{2} (1), χ_{2}^{2} (1), \dots, χ_{K}^{2} (1)$ are K independent central chi-square random variables each with one degree of freedom.

Lemma 1 was proven by Biggerstaff and Jackson²¹ using the location invariance of Q (e.g. Q can ben decomposed as

\sum_{k = 1}^{K} v_{k} (Y_{k} - μ)^{2} - v_{+} (\bar{Y} - μ)^{2}

), and distribution theory of quadratic forms in normal variables.^23–25

2.2 The HTS prediction interval

Suppose $τ^{2}$ is known, $\hat{μ} \sim N (μ, SE [\hat{μ}]^{2})$ and the observation in a future study $θ_{new} \sim N (μ, τ^{2})$ , where $SE [\hat{μ}] = \sqrt{1 / \sum_{k = 1}^{K} w_{k}}$ is a standard error of $\hat{μ}$ given $τ^{2}$ , and $w_{k} = (σ_{k}^{2} + τ^{2})^{- 1}$ . Assuming independence of θ_new and $\hat{μ}$ given μ, $θ_{new} - μ \sim N (0, τ^{2} + SE [\hat{μ}]^{2})$ . To replace the unknown $τ^{2}$ by its estimator ${\hat{τ}}_{DL}^{2}$ , the following approximation is used. If $(K - 2) ({\hat{τ}}_{DL}^{2} + \hat{SE} [\hat{μ}]^{2}) / (τ^{2} + SE [\hat{μ}]^{2})$ is approximately distributed as $χ^{2} (K - 2)$ , then $(θ_{new} - \hat{μ}) / \sqrt{{\hat{τ}}_{DL}^{2} + \hat{SE} [\hat{μ}]^{2}} \sim t (K - 2)$ , where $\hat{SE} [\hat{μ}] = \sqrt{1 / \sum_{k = 1}^{K} {\hat{w}}_{k}}$ is the standard error estimator of $\hat{μ}$ , and ${\hat{w}}_{k} = (σ_{k}^{2} + {\hat{τ}}_{DL}^{2})^{- 1}$ . By this approximation, the HTS prediction interval is

[\hat{μ} - t_{K - 2}^{α} \sqrt{{\hat{τ}}_{DL}^{2} + \hat{SE} [\hat{μ}]^{2}}, \hat{μ} + t_{K - 2}^{α} \sqrt{{\hat{τ}}_{DL}^{2} + \hat{SE} [\hat{μ}]^{2}}]

where

t_{K - 2}^{α}

is the

100 (1 - α / 2)

percentile of the t distribution with K–2 degrees of freedom. The t-approximation is appropriate only when both the number of studies and heterogeneity variance are large.

Several HTS-type prediction intervals following restricted maximum likelihood (REML) estimation of $τ^{2}$ have been proposed by Partlett and Riley.¹⁰ For example, they discussed an HTS-type prediction interval following REML with the Hartung–Knapp variance estimator²⁶ (HTS-HK) that replaces $\hat{μ}, {\hat{τ}}_{DL}^{2}$ , and $\hat{SE} [\hat{μ}]^{2}$ in the HTS prediction interval with ${\hat{μ}}_{R}, {\hat{τ}}_{R}^{2}$ , and ${\hat{SE}}_{HK} [{\hat{μ}}_{R}]^{2}$ , and a HTS-type prediction interval following REML with the Sidik–Jonkman bias-corrected variance estimator²⁷ (HTS-SJ) that replaces $\hat{μ}, {\hat{τ}}_{DL}^{2}$ , and $\hat{SE} [\hat{μ}]^{2}$ in the HTS prediction interval with ${\hat{μ}}_{R}, {\hat{τ}}_{R}^{2}$ , and ${\hat{SE}}_{SJ} [{\hat{μ}}_{R}]^{2}$ , where ${\hat{τ}}_{R}^{2}$ is the REML estimator for the heterogeneity variance^22,28,29 which is an iterative solution of the equation

{\hat{τ}}_{R}^{2} = \frac{\sum_{k = 1}^{K} {\hat{w}}_{R, k}^{2} {(Y_{k} - {\hat{μ}}_{R})^{2} + 1 / \sum_{l = 1}^{K} {\hat{w}}_{R, l} - σ_{k}^{2}}}{\sum_{k = 1}^{K} {\hat{w}}_{R, k}^{2}}

{\hat{w}}_{R, k} = (σ_{k}^{2} + {\hat{τ}}_{R}^{2})^{- 1}, {\hat{μ}}_{R} = \sum_{k = 1}^{K} {\hat{w}}_{R, k} Y_{k} / \sum_{k = 1}^{K} {\hat{w}}_{R, k}

, the Hartung–Knapp variance estimator is defined as

{\hat{SE}}_{HK} [{\hat{μ}}_{R}]^{2} = \frac{1}{K - 1} \sum_{k = 1}^{K} \frac{{\hat{w}}_{R, k} (Y_{k} - {\hat{μ}}_{R})^{2}}{\sum_{l = 1}^{K} {\hat{w}}_{R, l}}

the Sidik–Jonkman bias-corrected variance estimator

{\hat{SE}}_{SJ} [{\hat{μ}}_{R}]^{2} = \frac{\sum_{k = 1}^{K} {\hat{w}}_{R, k}^{2} (1 - {\hat{h}}_{k})^{- 1} (Y_{k} - {\hat{μ}}_{R})^{2}}{(\sum_{k = 1}^{K} {\hat{w}}_{R, k})^{2}}

and

{\hat{h}}_{k} = \frac{2 {\hat{w}}_{R, k}}{\sum_{l = 1}^{K} {\hat{w}}_{R, l}} - \frac{\sum_{l = 1}^{K} {\hat{w}}_{R, l}^{2} (σ_{l}^{2} + {\hat{τ}}_{R}^{2})}{(σ_{k}^{2} + {\hat{τ}}_{R}^{2}) (\sum_{l = 1}^{K} {\hat{w}}_{R, l}^{2})^{2}}

The empirical coverage of the HTS-HK and HTS-SJ prediction intervals is close to the nominal level under large heterogeneity variance and $K \geq 5$ .¹⁰

The HTS prediction intervals show severe under-coverage under small heterogeneity variance or for few studies (see Partlett and Riley¹⁰ and Section 3). We introduce a prediction interval in which uncertainty about $τ^{2}$ is accounted for and show that it is valid under a small number of studies.

2.3 The proposed prediction interval

We address the issue discussed in Section 2.2 by constructing a new prediction interval via a parametric bootstrap with the exact distribution of ${\hat{τ}}_{UDL}^{2}$ by using a confidence distribution (see Section 2.4). The proposed method uses an approximation that differs from those used by Higgins et al.⁶ The HTS prediction interval essentially combines the following two approximations: $(\hat{μ} - μ) / \sqrt{\hat{SE} [\hat{μ}]}$ approximately distributed as N(0, 1), which is often not satisfactory,³⁰ and $(K - 2) ({\hat{τ}}_{DL}^{2} + \hat{SE} [\hat{μ}]) / (τ^{2} + SE [\hat{μ}])$ is approximately distributed as $χ^{2} (K - 2)$ .

From now on, we make the following assumptions: Let the observation in a future study $θ_{new} \sim N (μ, τ^{2}), Y_{k} \sim N (μ, σ_{k}^{2} + τ^{2})$ given $σ_{k}^{2}$ and $τ^{2}$ , and θ_new and $\bar{μ} = \sum_{k = 1}^{K} w_{k} Y_{k} / \sum_{k = 1}^{K} w_{k}$ are independent. In Hartung³⁰ and Hartung and Knapp,²⁶ it was shown that assuming normality of Y_k, $(μ - \bar{μ}) / {SE}_{H} [\bar{μ}]$ is t-distributed with K–1 degrees of freedom, and ${SE}_{H} [\bar{μ}]$ is stochastically independent of $\bar{μ}$ , where ${SE}_{H} [\bar{μ}]^{2} = \frac{1}{K - 1} \sum_{k = 1}^{K} \frac{w_{k}}{w_{+}} (Y_{k} - \bar{μ})^{2}$ , and $w_{+} = \sum_{k = 1}^{K} w_{k}$ . By replacing $τ^{2}$ in $(\bar{μ} - μ) / {SE}_{H} [\bar{μ}]$ with an appropriate estimate ${\hat{τ}}^{2}, (\hat{μ} - μ) / {\hat{SE}}_{H} [\hat{μ}]$ is approximately t-distributed with K–1 degrees of freedom, where ${\hat{SE}}_{H} [\hat{μ}]^{2} = \frac{1}{K - 1} \sum_{k = 1}^{K} \frac{{\hat{w}}_{k}}{{\hat{w}}_{+}} (Y_{k} - \hat{μ})^{2}$ , and ${\hat{w}}_{+} = \sum_{k = 1}^{K} {\hat{w}}_{k}$ . This approximation exhibits better performance than $(\hat{μ} - μ) / \sqrt{\hat{SE} [\hat{μ}]} \sim^{approx .} N (0, 1)$ , even for a few studies (see Theorem 4.4 of Hartung³⁰).

The above assumptions and results lead to a system of equations

\begin{matrix} {\begin{matrix} \frac{θ_{new} - μ}{τ} = Z \\ \frac{\bar{μ} - μ}{{SE}_{H} [\bar{μ}]} = t_{K - 1} \end{matrix} \end{matrix}

(2)

where

Z \sim N (0, 1)

and

t_{K - 1} \sim t (K - 1)

. Solving for θ_new in equation (2) yields

θ_{new} = \bar{μ} + Z τ - t_{K - 1} {SE}_{H} [\bar{μ}]

(3)

and the prediction distribution has the same distribution as θ_new (even with

τ^{2}

unknown). By replacing

τ^{2}

in equation (3) with an appropriate estimator (not an estimate), we have

{\hat{θ}}_{new} = \hat{μ} + Z {\hat{τ}}_{UDL} - t_{K - 1} {\hat{SE}}_{H} [\hat{μ}]

and an approximate prediction distribution can be given by the distribution of

{\hat{θ}}_{new}

. We use the untruncated estimator

{\hat{τ}}_{UDL}^{2}

here, because we do not need the truncation to consider the distribution of an estimator of

τ^{2}

. Hence,

Pr (c_{l} < θ_{new} < c_{u}) = 1 - α

can be approximately evaluated by the distribution of

{\hat{θ}}_{new}

. Since

{\hat{θ}}_{new}

includes three random components,

{\hat{τ}}_{UDL}^{2}

, Z, and

t_{K - 1}

, this gives the following algorithm for the proposed prediction interval.

Algorithm 1

An algorithm for the proposed prediction interval.

Generate B bootstrap samples ${\tilde{τ}}_{b}^{2}$ ( $b = 1, \dots, B$ ) that are drawn from the exact distribution of ${\hat{τ}}_{UDL}^{2}$ , z_b that are drawn from N(0, 1), and t_b that are drawn from $t (K - 1)$ .

Calculate ${\tilde{μ}}_{b} = \sum_{k = 1}^{K} {\tilde{w}}_{bk} y_{k} / \sum_{k = 1}^{K} {\tilde{w}}_{bk}$ , and ${\tilde{θ}}_{new, b} = {\tilde{μ}}_{b} + z_{b} {\tilde{τ}}_{b} - t_{b} {\tilde{SE}}_{H, b} [{\tilde{μ}}_{b}]$ , where ${\tilde{w}}_{bk} = (σ_{k}^{2} + {\tilde{τ}}_{b}^{2})^{- 1}, {\tilde{SE}}_{H, b} [{\tilde{μ}}_{b}]^{2} = \frac{1}{K - 1} \sum_{k = 1}^{K} \frac{{\tilde{w}}_{bk}}{{\tilde{w}}_{b +}} (y_{k} - {\tilde{μ}}_{b})^{2}$ , and ${\tilde{w}}_{b +} = \sum_{k = 1}^{K} {\tilde{w}}_{bk}$ .

Calculate the prediction limits c_l and c_u that are $100 \times α / 2$ and $100 \times (1 - α / 2)$ percentage points of ${\tilde{θ}}_{new, b}$ , respectively.

An R package implementing the new method with the three data sets (see Section 4) and a documentation is available at the publisher’s website, the CRAN website (https://cran.r-project.org/package=pimeta) and GitHub (https://github.com/nshi-stat/pimeta/).

2.4 Sampling from the exact distribution of the estimator of τ²

Confidence distribution is a distribution estimator that can be defined and interpreted in a frequentist framework in which the parameter is a non-random quantity. A confidence distribution for the parameter of interest φ, as described below, can be easily defined as the cumulative distribution function of a statistic. The following definition of a confidence distribution was presented in Xie and Singh.¹⁵ In the definition, Φ is the parameter space of the unknown $φ, Y$ is a random vector, and $Y$ is the sample space corresponding to sample data $y = (y_{1}, y_{2}, \dots, y_{K})^{T}$ .

Definition 1

(R1) A function $H (\cdot) = H (y, φ)$ on $Y \times Φ \to [0, 1]$ is called a confidence distribution for a parameter φ; (R2) If for each given $y \in Y, H (\cdot)$ is a cumulative distribution function on φ; (R3) At the true parameter value $φ = φ_{0}, H (φ_{0}) \equiv H (y, φ_{0})$ , as a function of the sample y, follows the uniform distribution U(0, 1).

Confidence distribution has a theoretical relationship to the fiducial approach,³¹ and recent developments^11–15 have provided useful statistical tools that are more widely applicable than classical frequentist methods. For example, Efron’s bootstrap distribution³² is a confidence distribution and a distribution estimator of φ. In meta-analysis, the Q-profile method³³ for an approximate confidence interval for

τ^{2}

can be considered as an application of confidence distribution.¹² In this section, we propose the exact distribution of

{\hat{τ}}_{UDL}^{2}

, which is a distribution function for estimating the parameter

τ^{2}

using a confidence distribution, and then develop a method of sampling from the exact distribution. A useful theorem (Theorem 1) is introduced that provides a condition for confidence distribution.

Theorem 1

If a cumulative distribution function of a statistic, $T (Y)$ , is $F_{T} (T (y); φ) \equiv F_{T} (T (Y) \leq T (y); φ)$ , and F_T is a strictly monotone (without loss of generality, assume that it is decreasing) function in φ with the parameter space $Φ = {φ : φ_{min} \leq φ \leq φ_{max}}$ for each sample y, then $H (φ) = 1 - F_{T} (T (y); φ)$ is a confidence distribution for φ that satisfies Definition 1.

Lemma 2

Under the model in equation (1), $H (τ^{2}) = 1 - F_{Q} (q; τ^{2})$ is a confidence distribution for $τ^{2}$ . The proof of Theorem 1 is easy and hence is omitted. Lemma 2 can be easily proved by using Theorem 1, because $F_{Q} (q; τ^{2})$ is a strictly decreasing function in $τ^{2}$ .³⁴ Note that we use the untruncated version of an estimator of $τ^{2}$ with the parameter space $Φ = [τ_{min}^{2}, \infty]$ , and $τ_{min}^{2}$ can be negative.

The proposed algorithm samples from the confidence distribution, $H (τ^{2}) = 1 - F_{Q} (q_{obs}; τ^{2})$ , where q_obs is the observed value of Q. By applying Lemma 2 and the inverse transformation method, if U is distributed as U(0, 1), then $H^{- 1} (U)$ follows the distribution $H (τ^{2})$ . A sample ${\tilde{τ}}^{2} = H^{- 1} (u)$ can be computed by numerical inversion³⁵ of $H ({\tilde{τ}}^{2}) = u$ , where u is an observed value of the random variable U. If $H (0) > u$ , then the sample is truncated to zero ( ${\tilde{τ}}^{2} = 0$ ). It follows from Lemma 1 that $F_{Q} (q; τ^{2})$ is the distribution function of a positive linear combination of $χ^{2}$ random variables. It can be calculated with the Farebrother’s algorithm.³⁶

3 Simulations

We assessed the properties of the HTS and proposed prediction intervals through simulations.

Simulation data were generated by the random-effects model in equation (1), assuming independent normal errors $ε_{k} \sim N (0, σ_{k}^{2})$ and $u_{k} \sim N (0, τ^{2})$ . We conducted three sets of simulations described below.

By reference to Brockwell and Gordon^7,37 and Jackson,³⁴ parameter settings that mimic meta-analyses for estimating an overall mean log odds-ratio were determined. The average treatment effect μ was fixed at 0, as no generality is lost by setting μ to zero. The across-studies variance was set to $τ^{2} = 0.01, 0.05, 0.1, 0.2, 0.3, 0.4$ , or 0.5;^38,39 mean I² values were 29.8%, 66.0%, 79.1%, 88.2%, 91.8%, 93.7%, or 94.9%, respectively. The within-study variances $σ_{k}^{2}$ were generated from a scaled $χ^{2}$ distribution with one degree of freedom, multiplied by 0.25, and then truncated to lie within $[0.009, 0.6]$ . The number of studies was set to $K = 3, 5, 10, 15, 20$ , or 25.

In reference to Partlett and Riley,¹⁰ parameter settings were determined to evaluate the empirical performance of prediction intervals under various relative degrees of heterogeneity scenarios. The within-study variances $σ_{k}^{2}$ were generated from $σ^{2} χ^{2} (n - 1) / (n - 1)$ , an average within-study variance was set to $σ^{2} = 0.1$ , and the study sample size was set to n = 30, where $χ^{2} (n - 1)$ is a random number from a $χ^{2}$ distribution with n–1 degrees of freedom. The degree of heterogeneity is controlled using the ratio $τ^{2} / σ^{2}$ . The heterogeneity parameter was set to $τ^{2} = 0.01, 0.05, 0.1$ , or 1, which corresponds to $τ^{2} / σ^{2} = 0.1, 0.5, 1$ , or 10; mean I² values were 9.1%, 33.3%, 50.0%, or 90.9%, respectively. In addition to the above situation where (ii-a) with all studies of similar size, we consider situations (ii-b) with one large study or (ii-c) with one small study (i.e. a within study variance 10 times smaller or 10-fold higher than the others; one study was randomly selected and the within study variance was set to $σ_{k}^{2} / 10$ or $10 σ_{k}^{2}$ ). The average treatment effect μ was fixed at 1. The number of studies was set to $K = 3, 5, 7, 10, 25$ , or 100.

We generated data for $K \times 2 \times 2$ tables using a method similar to that used by Sidik and Jonkman²² and considered meta-analyses based on log odds-ratios. The heterogeneity variance was set to $τ^{2} = 0.01, 0.1, 0.2, 0.4$ , or 0.6. The average treatment effect was set to $μ = 0, - 0.5$ , or 0.5 to assess the impact of standard errors for odds ratios. The number of studies was set to $K = 3, 6, 12, 24, 48$ , or 96. For each $τ^{2}$ and K, we first generated θ_k from $N (μ, τ^{2})$ . For each study, the sample sizes were set to be equal $n_{0 k} = n_{1 k}$ and were randomly selected from integers between 20 and 200. The responses of the control group, $X_{0 k}$ , were generated from a binomial distribution $Bin (n_{0 k}, p_{0 k})$ , and the probability $p_{0 k}$ was randomly drawn from a uniform distribution $U (0.05, 0.65)$ . The responses of the treatment group, $X_{1 k}$ , were generated from a binomial distribution $Bin (n_{1 k}, p_{1 k})$ and probability $p_{1 k} = p_{0 k} exp {θ_{k}} (1 - p_{0 k} + p_{0 k} exp {θ_{k}})$ . Finally, we constructed an estimator of θ_k as $Y_{k} = log [X_{1 k} (n_{0 k} - X_{0 k}) / {X_{0 k} (n_{1 k} - X_{1 k})}]$ , its variance estimator as ${\hat{σ}}_{k}^{2} = 1 / X_{1 k} + 1 / (n_{1 k} - X_{1 k}) + 1 / X_{0 k} + 1 / (n_{0 k} - X_{0 k})$ , and we used ${\hat{σ}}_{k}^{2}$ rather than $σ_{k}^{2}$ . If any cells are empty, we added 0.5 to each cell for all K tables. Mean I² values were 7.1%, 42.3%, 58.7%, 73.1%, or 79.7%, which correspond to $τ^{2} = 0.01, 0.1, 0.2, 0.4$ , or 0.6, respectively.

For each setting, we simulated 25,000 replications. For each method, two-tailed 95% prediction intervals were calculated. The number of bootstrap samples B was set to 5000. The coverage probability was estimated by the proportion of simulated prediction intervals containing the result of a future study θ_new that was generated from a normal distribution $N (μ, τ^{2})$ .

The results of simulation (i) are presented in Figure 1. The coverage probabilities of the HTS prediction interval were approximately 90%, far short of the nominal level of 95%. The under-coverage of the HTS prediction interval reflects the rough t-approximation; thus, the source of the problem is substitution of an estimate for $τ^{2}$ or ignoring uncertainty in $τ^{2}$ . The results show that the HTS-HK and HTS-SJ prediction intervals are also deficient. The coverage probabilities for the HTS-HK and HTS-SJ prediction intervals almost retained the nominal level except in situations where the relative degree of heterogeneity is small or moderate. For example, the coverage probabilities of the HTS-HK prediction interval were 82.8%–98.5% for $τ^{2} = 0.01$ , 83.8%–96.5% for $τ^{2} = 0.05$ , and 85.3%–95.6% for $τ^{2} = 0.1$ ; the coverage probabilities of the HTS-SJ prediction interval were 79.7%–96.7% for $τ^{2} = 0.01$ , 80.9%–93.9% for $τ^{2} = 0.05$ , and 83.9%–92.9% for $τ^{2} = 0.1$ . By contrast, the coverage probabilities for the proposed prediction interval almost always retained the nominal level. The only exception was when K = 25 and $τ^{2} = 0.01$ , where the coverage probability for the proposed prediction interval was 93.0%, slightly below the nominal level. In this case, the coverage probability for the HTS, HTS-HJ, and HTS-SJ prediction intervals were even smaller, at 82.4%, 82.8%, and 82.0%, respectively. Analyses using very few studies (K < 5) pose problems in random-effects models, as discussed by Higgins et al.⁶ Nevertheless, the proposed method performed well even when K = 3. The nominal level was attained for nearly all values of the heterogeneity parameter in the proposed prediction interval.

Figure 1.

Simulation results (i): the performance of the HTS and proposed prediction intervals. The number of studies $K = 3, 5, 10, 15, 20$ , or 25. The number of simulations was 25,000. Methods: Proposed, the proposed prediction interval; HTS, the HTS prediction interval; HTS-HK, the HTS-type prediction interval following REML with the Hartung–Knapp variance estimator; HTS-SJ, the HTS-type prediction interval following REML with the Sidik–Jonkman bias-corrected variance estimator.

The results of simulation (ii-a), with all studies of similar size, are presented in Figure 2. The results show that all HTS prediction intervals are also deficient except for $τ^{2} = 0.001, 0.05$ . The coverage probabilities for all HTS prediction intervals almost retained the nominal level for $τ^{2} = 1$ and $K \geq 5$ . The coverage probabilities were too large for K = 3 and too small for K = 5–25 and $τ^{2} = 0.1$ . In the case of $τ^{2} = 0.01, 0.05$ , the coverage probabilities of the HTS-HJ and HTS-SJ prediction intervals were too small for K = 5–100, and the coverage probability of the HTS prediction interval was too small for K = 10–100. By contrast, the coverage probabilities for the proposed prediction interval almost always retained the nominal level. The only exception was when $τ^{2} = 0.01$ , where the coverage probabilities for the proposed prediction interval were 92.8%–96.7%, slightly below the nominal level for K = 100. The results of simulation (ii-b), with one large study, are presented in Figure 2. The coverage probabilities appear to be relatively poor compared to the balanced case (ii-a) even in large heterogeneity variance not only for the HTS-HK prediction interval, but also for the HTS-SJ prediction interval. Moreover, the performance of the HTS-SJ prediction interval was somewhat poorer (showed under coverage) compared to the HTS-HK prediction interval. In contrast, the coverage probabilities for the proposed prediction interval nearly always retained the nominal level. The results of simulation (ii-c), with one small study, are presented in Supplementary Figure S1. The coverage probabilities were similar to those of the balanced case (ii-a).

Figure 2.

Simulation results (ii): the performance of the HTS and proposed prediction intervals (ii-a) with all studies of similar size and (ii-b) with one large study. The heterogeneity parameters $τ^{2} = 0.01, 0.05, 0.1$ , or 1. The number of simulations was 25,000. Methods: Proposed, the proposed prediction interval; HTS, the HTS prediction interval; HTS-HK, the HTS-type prediction interval following REML with the Hartung–Knapp variance estimator; HTS-SJ, the HTS-type prediction interval following REML with the Sidik–Jonkman bias-corrected variance estimator.

The results of simulation (iii) for μ = 0 are presented in Figure 3. The coverage probabilities for all HTS prediction intervals were too large for K = 3, too small for $K = 6, 12, 24, 48, 96$ and $τ^{2} \leq 0.1$ , and nearly retained the nominal level for $τ^{2} = 0.6$ or K = 96. In contrast, the coverage probabilities for the proposed prediction interval nearly retained the nominal level, except for K = 96 and $τ = 0.01$ . The results of simulation (iii) for $μ = - 0.5, 0.5$ are presented in Supplementary Figures S2 and S3. The coverage probabilities for $μ = - 0.5, 0.5$ were similar to those of the case for μ = 0.

Figure 3.

Simulation results (iii): the performance of the HTS and proposed prediction intervals for μ = 0. The number of studies $K = 3, 6, 12, 24, 48$ , or 96. The number of simulations was 25,000. Methods: Proposed, the proposed prediction interval; HTS, the HTS prediction interval; HTS-HK, the HTS-type prediction interval following REML with the Hartung–Knapp variance estimator; HTS-SJ, the HTS-type prediction interval following REML with the Sidik–Jonkman bias-corrected variance estimator.

In summary, the HTS prediction intervals have insufficient coverage, except when the relative degree of heterogeneity is large and may show severe under-coverage under realistic meta-analysis settings, possibly providing misleading results and interpretation. In contrast, the proposed prediction interval achieves the nominal level of coverage.

4 Applications

We applied the methods to the following three published random-effects meta-analyses.

Set-shifting data: Higgins et al.⁶ re-analyzed data⁴⁰ that included 14 studies evaluating the set-shifting ability in people with eating disorders by using a prediction interval. Standardized mean differences in the time taken to complete Trail Making Test between subjects with eating disorders and healthy controls were collected. Positive estimates indicate impairment in set shifting ability in people with eating disorders.

Pain data: The pain data^4,41 included 22 studies comparing the treatment effect of antidepressants on reducing pain in patients with fibromyalgia syndrome. The treatment effects were summarized using standardized mean differences on a visual analog scale for pain between the antidepressant group and control group. Negative estimates indicate the reduction of pain in the antidepressant group.

Systolic blood pressure (SBP) data: Riley et al.⁴ analyzed a hypothetical meta-analysis. They generated a data set of 10 studies examining the same antihypertensive drug. Negative estimates suggested reduced blood pressure in the treatment group.

These data sets are reproduced in Figure 4. The number of bootstrap samples B was set to 50,000.

Figure 4.

The three data sets and summary results: (a) Set-shifting data⁴⁰ (K = 14), (b) Pain data⁴¹ (K = 22), and (c) SBP data⁴ (K = 10). CI: confidence interval; PI: prediction interval.

Table 1 presents estimates of the average treatment effect and its confidence interval, heterogeneity measures, the P-value for the test of heterogeneity, the proposed prediction interval, and the HTS prediction intervals. None of the confidence intervals for the average treatment effect included 0 (set-shifting data:

[0.19, 0.53]

; pain data:

[- 0.55, - 0.30]

; SBP data:

[- 0.48, - 0.18]

). This means that on average the interventions are significantly effective. However, small, moderate, and large heterogeneity were observed in the three data sets (set-shifting data:

{\hat{τ}}_{DL}^{2} = 0.023, {\hat{I}}^{2} = 22.5 %

; pain data:

{\hat{τ}}_{DL}^{2} = 0.185, {\hat{I}}^{2} = 44.9 %

; SBP data:

{\hat{τ}}_{DL}^{2} = 0.023, {\hat{I}}^{2} = 70.5 %

). Accounting for heterogeneity, prediction intervals would provide additional relevant statistical information. There were large differences between the 95% confidence interval and prediction intervals, even in the case of small heterogeneity.

Table 1.

Results from the three data sets: the average treatment effect ( $\hat{μ}$ ) and its 95% confidence interval, heterogeneity measures ( ${\hat{τ}}_{DL}^{2}, {\hat{τ}}_{R}^{2}$ and ${\hat{I}}^{2}$ ), the P-value for the test of heterogeneity, the proposed prediction interval, and the HTS prediction intervals.

		Set-shifting	Pain	SBP
Data		(K = 14)	(K = 22)	(K = 10)
$\hat{μ}$ (DL)		0.36	–0.43	–0.33
95%CI (DL)		[0.19, 0.53]	[–0.55, –0.30]	[–0.48, –0.18]
${\hat{τ}}_{DL}^{2}$		0.023	0.034	0.023
${\hat{τ}}_{R}^{2}$		0.013	0.025	0.070
${\hat{I}}^{2}$ (DL)		22.5%	44.9%	70.5%
P-value for heterogeneity		0.209	0.012	<0.001
95%PI	Proposed	[–0.13, 0.85]	[–0.89, 0.02]	[–0.88, 0.23]
	HTS	[–0.02, 0.74]	[–0.84, –0.02]	[–0.76, 0.09]
	HTS-HK	[0.05, 0.67]	[–0.78, –0.06]	[–0.99, 0.33]
	HTS-SJ	[0.06, 0.67]	[–0.77, –0.07]	[–0.98, 0.33]
length of	Proposed	0.98	0.91	1.10
95%PI	HTS	0.76	0.82	0.85
	HTS-HK	0.62	0.72	1.32
	HTS-SJ	0.61	0.72	1.31

HTS: Higgins–Thompson–Spiegelhalter; SBP: systolic blood pressure.

As shown in Figure 4 and summarized in Table 1, the proposed prediction intervals were substantially wider than the HTS prediction intervals in all three analyses. The proposed prediction intervals were 29%, 11%, and 31% wider than the HTS prediction intervals for the set-shifting data, pain data, and SBP data, respectively. As observed in Section 3, the HTS-HK and HTS-SJ prediction intervals showed similar results. The proposed prediction intervals were 58%, 27% wider and 17% narrower than the HTS-HK (or HTS-SJ) prediction intervals for the set-shifting data, pain data, and SBP data, respectively. Only for the SBP data, the proposed prediction interval was narrower than the HTS-HK prediction interval; this is because the two intervals were based on different heterogeneity variance estimators and ${\hat{τ}}_{DL}^{2} < {\hat{τ}}_{R}^{2}$ .

The prediction intervals may lead to different interpretations of the results. In the set-shifting data, the HTS-HK and HTS-SJ prediction intervals did not include 0, but the proposed prediction interval included 0. For the pain data, the HTS, HTS-HK, and HTS-SJ prediction intervals did not include 0, in a frequentist sense, suggesting that the intervention may be beneficial in most subpopulations. In contrast, the proposed prediction interval included 0, indicating that the intervention may not be beneficial in some subpopulations. However, taking a Bayesian perspective, all the prediction intervals suggest that there is a large probability and the treatment will be effective in a new population. The simulation results in Section 3 suggest that the HTS prediction intervals could have under-coverage in situations where the relative degree of heterogeneity is small or moderate. Since ${\hat{τ}}_{DL}^{2}$ of three data sets and ${\hat{τ}}_{R}^{2}$ of set-shifting and pain data were small ( $\approx 0.02$ ), it may be too narrow under realistic situations and may provide misleading results.

5 Discussion and conclusion

For the random-effects model in meta-analysis, the average treatment effect and its confidence interval have been used with heterogeneity measures such as the I²-statistic and $τ^{2}$ . However, the results from random-effects models have sometimes been misinterpreted. Thus, the new concept “prediction interval” was proposed, which is useful in applying the results to other subpopulations and in decision making. The HTS prediction intervals have a theoretical problem, namely that its rough t-approximation could have a detrimental impact on the coverage probability. We have presented an appropriate prediction interval to account for the uncertainty in $τ^{2}$ by using a confidence distribution.

Simulation studies showed that the HTS prediction intervals could have severe under-coverage for realistic meta-analysis settings and might lead to misleading results and interpretation. The simulation results suggested that the HTS prediction interval may be too narrow when analyzing a small number of studies. This interval is valid when $K \geq 25$ , but in many meta-analyses K is much smaller than 25. The HTS-HK and HTS-SJ prediction intervals may be too narrow when the relative degree of heterogeneity is small. By contrast, the coverage probabilities for the proposed prediction interval satisfactorily retained the nominal level. Although Higgins et al.⁶ cautioned that the random-effects model may not work well under very small numbers of studies (K < 5), the proposed method performed well even when K = 3. Since the heterogeneity parameter had very little effect on the performance of the proposed prediction interval, the method would be valid regardless of the value of the heterogeneity parameter. Moreover, all prediction intervals (i.e. the random-effects model in equation (1)) assume normality of the between-study distribution of true effects, $u_{k} \sim N (0, τ^{2})$ , but the assumption may not be true in practice. A full Bayesian approach may be useful for constructing a suitable prediction interval.

Applications to the three published random-effects meta-analyses concluded that substantially different results and interpretation might be obtained from the prediction intervals. Since the HTS prediction interval is always narrower and the HTS-HK and HTS-SJ prediction intervals are narrower when the heterogeneity parameter is small or moderate, we should be cautious in using and interpreting these approaches.

In conclusion, we showed that the proposed prediction interval works well and is suitable for random-effects meta-analysis. As shown in the three illustrative examples, quite different results and interpretations are obtained using our new method. Extensions of these results to other complicated models such as network meta-analysis are now warranted.

Supplemental Material

Supplemental material for Prediction intervals for random-effects meta-analysis: A confidence distribution approach

Supplemental material for Prediction intervals for random-effects meta-analysis: A confidence distribution approach by Kengo Nagashima, Hisashi Noma and Toshi A Furukawa in Statistical Methods in Medical Research

Footnotes

Acknowledgements

The authors would like to thank two reviewers for very insightful and constructive comments that substantially improved this article.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This study was supported by CREST from the Japan Science and Technology Agency (Grant number: JPMJCR1412).

Supplemental material

Supplemental material for this article is available online.

References

Borenstein

Hedges

Higgins

JPT

, et al. Introduction to meta-analysis, Chichester: Wiley, 2009.

Higgins

JPT

Thompson

. Quantifying heterogeneity in a meta-analysis. Stat Med 2002; 21: 1539–1558.

Higgins

JPT

Thompson

Deeks

, et al. Measuring inconsistency in meta-analyses. BMJ 2003; 327: 557–560.

Riley

Higgins

JPT

Deeks

. Interpretation of random effects meta-analyses. BMJ 2011; 342: d549–d549.

Riley

Gates

Neilson

, et al. Statistical methods can be improved within Cochrane pregnancy and childbirth reviews. J Clin Epidemiol 2011; 64: 608–618.

Higgins

JPT

Thompson

Spiegelhalter

. A re-evaluation of random-effects meta-analysis. J R Stat Soc Ser A Stat Soc 2009; 172: 137–159.

Brockwell

Gordon

. A comparison of statistical methods for meta-analysis. Stat Med 2001; 20: 825–840.

Noma

. Confidence intervals for a random-effects meta-analysis based on Bartlett-type corrections. Stat Med 2011; 30: 3304–3312.

Kontopantelis

Springate

Reeves

. A re-analysis of the Cochrane library data: the dangers of unobserved heterogeneity in meta-analyses. PLoS One 2013; 8: e69930–e69930.

10.

Partlett

Riley

. Random effects meta-analysis: coverage performance of 95% confidence and prediction intervals following REML estimation. Stat Med 2017; 36: 301–317.

11.

Schweder

Hjort

. Confidence and likelihood. Scand J Stat 2002; 29: 309–332.

12.

Schweder

Hjort

. Confidence, likelihood, probability: statistical inference with confidence distributions, New York: Cambridge University Press, 2016.

13.

Singh

Xie

Strawderman

. Combining information from independent sources through confidence distributions. Ann Stat 2005; 33: 159–183.

14.

Singh

Xie

Strawderman

Confidence distribution (CD) – distribution estimator of a parameter. In: Liu

Strawderman

Zhang

(eds). Complex Datasets and inverse problems: tomography, networks, and beyond 2007; vol. 54, Beachwood: Institute of Mathematical Statistics, pp. 132–150.

15.

Xie

Singh

. Confidence distribution, the frequentist distribution estimator of a parameter: a review. Int Stat Rev 2013; 81: 3–39.

16.

Cochran

. Problems arising in the analysis of a series of similar experiments. J R Stat Soc 1937; 4(Supplment): 102–118.

17.

Cochran

. The combination of estimates from different experiments. Biometrics 1954; 10: 101–129.

18.

DerSimonian

Laird

. Meta-analysis in clinical trials. Control Clin Trials 1986; 7: 177–188.

19.

Whitehead

. A general parametric approach to the meta-analysis of randomized clinical trials. Stat Med 1991; 10: 1665–1677.

20.

Biggerstaff

Tweedie

. Incorporating variability of estimates of heterogeneity in the random effects model in meta-analysis. Stat Med 1997; 16: 753–768.

21.

Biggerstaff

Jackson

. The exact distribution of Cochran’s heterogeneity statistic in one-way random effects meta-analysis. Stat Med 2008; 27: 6093–6110.

22.

Sidik

Jonkman

. A comparison of heterogeneity variance estimators in combining results of studies. Stat Med 2007; 26: 1964–1981.

23.

Scheffé

. The analysis of variance, New York: Wiley, 1959.

24.

Graybill

. Theory and application of the linear model, North Scituate: Duxbury Press, 1976.

25.

Mathai

Provost

. Quadratic forms in random variables: theory and applications, New York: Marcel Dekker, 1992.

26.

Hartung

Knapp

. On tests of the overall treatment effect in meta-analysis with normally distributed responses. Stat Med 2001; 20: 1771–1782.

27.

Sidik

Jonkman

. Robust variance estimation for random effects meta-analysis. Comput Stat Data Anal 2006; 50: 3681–3701.

28.

Harville

. Maximum likelihood approaches to variance component estimation and to related problems. J Am Stat Assoc 1977; 72: 320–339.

29.

Raudenbush

Bryk

. Empirical Bayes meta-analysis. J Educ Stat 1985; 10: 75–98.

30.

Hartung

. An alternative method for meta-analysis. Biom J 1999; 41: 901–916.

31.

Fisher

. The fiducial argument in statistical inference. Ann Eugen 1935; 6: 391–398.

32.

Efron

B. R. A

. Fisher in the 21st century. Stat Sci 1998; 13: 95–122.

33.

Viechtbauer

. Confidence intervals for the amount of heterogeneity in meta-analysis. Stat Med 2007; 26: 37–52.

34.

Jackson

. Confidence intervals for the between-study variance in random effects meta-analysis using generalised Cochran heterogeneity statistics. Res Synth Meth 2013; 4: 220–229.

35.

Forsythe

Malcolm

Moler

. Computer methods for mathematical computations, Englewood Cliffs: Prentice-Hall, 1977.

36.

Farebrother

. Algorithm AS 204: the distribution of a positive linear combination of

χ^{2}

random variables. J R Stat Soc Ser C Appl Stat 1984; 33: 332–339.

37.

Brockwell

Gordon

. A simple method for inference on an overall effect in meta-analysis. Stat Med 2007; 26: 4531–4543.

38.

Turner

Davey

Clarke

, et al. Predicting the extent of heterogeneity in meta-analysis, using empirical data from the Cochrane Database of Systematic Reviews. Int J Epidemiol 2012; 41: 818–827.

39.

Rhodes

Turner

Higgins

JPT

. Predictive distributions were developed for the extent of heterogeneity in meta-analyses of continuous outcome data. J Clin Epidemiol 2015; 68: 52–60.

40.

Roberts

Tchanturia

Stahl

, et al. A systematic review and meta-analysis of set-shifting ability in eating disorders. Psychol Med 2007; 37: 1075–1084.

41.

Häuser

Bernardy

Üçeyler

, et al. Treatment of fibromyalgia syndrome with antidepressants: a meta-analysis. JAMA 2009; 301: 198–209.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

0.12 MB

Prediction intervals for random-effects meta-analysis: A confidence distribution approach

Abstract

Keywords

1 Introduction

2 Method

2.1 The random-effects model and the exact distribution of Cochran’s Q statistic

Lemma 1

2.2 The HTS prediction interval

2.3 The proposed prediction interval

Algorithm 1

2.4 Sampling from the exact distribution of the estimator of τ2

Definition 1

Theorem 1

Lemma 2

3 Simulations

4 Applications

5 Discussion and conclusion

Supplemental Material

Supplemental material for Prediction intervals for random-effects meta-analysis: A confidence distribution approach

Footnotes

Acknowledgements

Declaration of conflicting interests

Funding

Supplemental material

References

Supplementary Material

2.4 Sampling from the exact distribution of the estimator of τ²