What is already known?

• • The routine report of a random-effects meta-analysis presents the estimated mean parameter and its confidence interval; however, these do not reveal the heterogeneity, that is, the differences in the underlying effects of individual studies.

• • The prediction interval (PI) in random-effects meta-analysis is a useful tool to assess the heterogeneity; it intends to contain the true effect of a new similar study with a pre-specified probability.

• • The Higgins–Thompson–Spiegelhalter PI method cannot keep the nominal mean coverage for cases when the number of studies is small, but the parametric bootstrap method can maintain appropriate mean coverage even if there are only three involved studies.

What is new?

• • If the number of involved studies is not large enough, the distribution of coverage probabilities is skewed. It means that with high probability, the coverage probability is close to 1; however, the mean coverage can still be close to the nominal level.

• • The distribution of coverages reveals that even if the nominal mean coverage is maintained, the interpretation of a single published PI: “an interval that covers 95% of the true effects distribution” is not valid, because by definition, the PI covers this distribution only on average. In other words, if one considers a published PI with a nominal level of 0.95, it will not be true that approximately 95% of the upcoming individual studies will be within the published interval.

Potential impact for RSM readers

• • Even if the nominal mean coverage is approximately achieved, researchers using the PI should be cautious with its interpretation.

• • Further research is needed regarding how new methods can be developed or how the currently existing PI methods can be altered to fulfill more strict conditions than the simple mean coverage criterion.

2.1 Parameter estimation in the random-effects model

The conventional aim of meta-analysis is to estimate μ and τ². An unbiased estimate of μ can be given by a weighted mean: $\widehat{\unicode{x3bc}}=\frac{\sum \nolimits_{\mathrm{k}}{\widehat{\mathrm{\varTheta}}}_{\mathrm{k}}\;{\mathrm{w}}_{\mathrm{k}}}{\sum \nolimits_{\mathrm{k}}{\mathrm{w}}_{\mathrm{k}}}$ . Weights are conventionally chosen as ${\mathrm{w}}_{\mathrm{k}}=\mathsf{1}/\left({\unicode{x3c3}^2}_{\mathrm{k}}+{\mathsf{\tau}}^2\right)$ , because the variance of this estimator, $\mathsf{Var}\left(\widehat{\unicode{x3bc}}\right)=\mathsf{1}/\sum \nolimits_{\mathrm{k}}{\mathrm{w}}_{\mathrm{k}}$ provides the uniformly minimum variance estimator for μ.Reference Viechtbauer ¹⁹ Note that in the formula for weights, neither ${\unicode{x3c3}^2}_{\mathrm{k}}$ , nor ${\unicode{x3c4}}^2$ are known quantities. The standard approach treats ${\unicode{x3c3}^2}_{\mathrm{k}}$ -s as fixed and known constants and replaces them with ${\widehat{\unicode{x3c3}}}_{\mathrm{k}}^2$ , the within-study variance estimates. Also, ${\unicode{x3c4}}^2$ is replaced by its estimate, ${\widehat{\unicode{x3c4}}}^2$ . Hence, in practice, weights are estimated as ${\widehat{\mathrm{w}}}_{\mathrm{k}}=\mathsf{1}/\left({\widehat{\unicode{x3c3}}}_{\mathrm{k}}^2+{\widehat{\unicode{x3c4}}}^2\right)$ and the estimate for $\mathsf{Var}\left(\widehat{\unicode{x3bc}}\right)$ is calculated with these weights as $\widehat{\mathrm{Var}}\left(\widehat{\unicode{x3bc}}\right)=\mathsf{1}/\sum \nolimits_{\mathrm{k}}{\widehat{\mathrm{w}}}_{\mathrm{k}}$ . This approach, like the within-study normality, is justifiable when the sample sizes are appropriately large.Reference Veroniki, Jackson and Bender ⁶ There are numerous methods in the literature for the estimation of the ${\unicode{x3c4}}^2$ parameter, a good description of these can be found in Veroniki et al.Reference Veroniki, Jackson and Viechtbauer ²⁰ We used two frequently applied ${\unicode{x3c4}}^2$ estimators to construct PI, the method of moments-based estimator suggested by DerSimonian and Laird (DL),Reference DerSimonian and Laird ¹⁶ and the restricted maximum likelihood (REML) estimator.Reference Raudenbush, Cooper, Hedges and Valentine ²¹ Hartung and KnappReference Hartung and Knapp ²² and Sidik and JonkmanReference Sidik and Jonkman ²³ independently proposed an alternative variance estimator for the μ parameter, denoted in the following by ${\widehat{\mathrm{Var}}}_{\mathrm{HKSJ}}\left(\widehat{\unicode{x3bc}}\right)$ . This estimator takes into account that the ${\mathrm{w}}_{\mathrm{k}}$ weights are not known but estimated quantities, and they use it with t distribution to construct a confidence interval for μ. ${\widehat{\mathrm{Var}}}_{\mathrm{HKSJ}}\left(\widehat{\unicode{x3bc}}\right)$ is usually larger than $\widehat{\mathrm{Var}}\left(\widehat{\unicode{x3bc}}\right)$ , but in some rare cases it can be the opposite way.

Heterogeneity is frequently assessed by the Q test introduced by Cochran,Reference Cochran ²⁴ testing the null hypothesis of homogeneity ( ${\mathrm{H}}_0:{\unicode{x3c4}}^2=\mathsf{0}$ ), and by the I² statistic or the magnitude of ${\widehat{\unicode{x3c4}}}^2$ . The I² statistic was proposed by Higgins and Thompson,Reference Higgins and Thompson ²⁵ and it can be interpreted as the proportion of total variance in the effect estimates that can be attributed to the between-study variance.

3.1 Formal definition of the prediction interval

We proceed to use the notation introduced in Section 2, stating that we have K studies that can be modelled by the data-generating process represented by (1). Let F denote the cumulative distribution function of the true effects distribution. The inputs of the meta-analysis are realizations of the random variables ${\widehat{\mathrm{\varTheta}}}_{\mathrm{k}}$ along with the corresponding within-study variances ${\widehat{\unicode{x3c3}}}_{\mathrm{k}}^2$ . Let D denote the input variables, that is, the collection of the random variables ${\widehat{\mathrm{\varTheta}}}_{\mathrm{k}}$ ${\widehat{\unicode{x3c3}}}_{\mathrm{k}}^2$ . Based on a realization of D, the applied PI method outputs a prediction interval; let L(D) and U(D) denote the lower and upper PI boundaries. Note that if a realization of D is given, the PI is a deterministic function of this realization. However, L(D) and U(D) are random variables, as the input data D can be considered as a 2K dimensional vector of random variables. The length of the interval can be defined as

(2) $$\begin{align}\ell =\mathsf{U}\left(\mathrm{D}\right)- \mathsf{L}\left(\mathsf{D}\right)\end{align}$$

and the covered probability as

(3) $$\begin{align}\mathsf{C}=\mathsf{F}\left(\mathsf{U}\left(\mathrm{D}\right)\right)-\mathsf{F}\left(\mathsf{L}\left(\mathrm{D}\right)\right).\end{align}$$

Definition 1. In random-effects meta-analysis, the random interval (L(D), U(D)) with E[C] = 1 – α is called a prediction interval with expected coverage probability 1 – α.

Let ${\mathrm{\varTheta}}_{\mathrm{New}}$ denote the true effect in a new study, that is, a random variable having the same distribution as ${\mathrm{\varTheta}}_{\mathrm{k}}$ and which is independent of D. The tower rule implies that

(4) $$\begin{align}\mathsf{E}\left[\mathrm{C}\right]=\mathsf{P}\left[\;{\mathrm{\varTheta}}_{\mathrm{New}}\in \left(\;\mathsf{L}\left(\mathrm{D}\right),\mathsf{U}\left(\mathrm{D}\right)\;\right)\;\right].\end{align}$$

(For proof, see Appendix Proof 1 of the Supplementary Material). The above equality makes it clear that E[C] can be approximated by performing several times the following simulation procedure: simulating the input data D from model (1), calculating the interval (L(D), U(D)), simulating independently an additional true effect ${\mathrm{\varTheta}}_{\mathrm{New}}$ , then checking if ${\mathrm{\varTheta}}_{\mathrm{New}}\in \left(\;\mathsf{L}\left(\mathrm{D}\right),\mathsf{U}\left(\mathrm{D}\right)\;\right)$ . If the number of simulations is large, the relative frequency of the cases when the calculated interval contains the newly generated true effect will be close to E[C]. All previous studies that investigated any PI method in the random-effects meta-analysis assessed the coverage performance only by approximating E[C] based on equation (4) using the above-described simulation approach.Reference Brannick, French, Rothstein, Kiselica and Apostoloski ^{1 ,} Reference Nagashima, Noma and Furukawa ^{8 –} Reference Gnambs and Schroeders ¹¹

3.2 Methods to construct the prediction interval

3.2.1 The Higgins–Thompson–Spiegelhalter (HTS) method

Higgins, Thompson, and Spiegelhalter were the first in their 2009 article to introduce the prediction interval in the context of the frequentist random-effects meta-analysis.Reference Higgins, Thompson and Spiegelhalter ⁵ They construct their PI estimator assuming $\widehat{\unicode{x3bc}}\sim \mathsf{N}\left(\mathsf{\mu},\mathsf{Var}\left(\widehat{\unicode{x3bc}}\right)\;\right)$ and ${\mathrm{\varTheta}}_{\mathrm{New}}\sim \mathsf{N}\left(\mathsf{\mu},{\unicode{x3c4}}^2\right)$ and the independence of ${\mathrm{\varTheta}}_{\mathrm{New}}$ and $\widehat{\unicode{x3bc}}$ . These infer ${\mathrm{\varTheta}}_{\mathrm{New}}-\widehat{\unicode{x3bc}}\sim \mathrm{N}\left(0,{\unicode{x3c4}}^2+\mathsf{Var}\left(\widehat{\unicode{x3bc}}\right)\;\right)$ . They argue that as ${\unicode{x3c4}}^2$ has to be estimated, and it impacts both parts of the variance of ${\mathrm{\varTheta}}_{\mathrm{New}}-\widehat{\unicode{x3bc}}$ , the t-distribution with K − 2 degrees of freedom should be used to construct the PI. Assuming $\frac{{\mathrm{\varTheta}}_{\mathrm{New}}-\widehat{\unicode{x3bc}}}{\sqrt{{\widehat{\unicode{x3c4}}}^2+\widehat{\mathrm{Var}}\left(\widehat{\unicode{x3bc}}\right)}}\sim {\mathrm{t}}_{\mathrm{K}-2}$ , they propose the (1 – α) expected coverage level PI as

(5) $$\begin{align}\widehat{\unicode{x3bc}}\pm {\mathrm{t}}_{\mathrm{K}-2}^{\unicode{x3b1}}\sqrt{{\widehat{\unicode{x3c4}}}^2+\widehat{\mathrm{Var}}\left(\widehat{\unicode{x3bc}}\right),}\end{align}$$

where ${\mathrm{t}}_{\mathrm{K}-2}^{\unicode{x3b1}}$ is the (1 − α/2) quantile of the t-distribution with K − 2 degrees of freedom.

Note that using the ${\mathrm{t}}_{\mathrm{K}-2}$ distribution has no convincing theoretical basis; Viechtbauer, for example, offers the standard normal distribution to construct the PI as the default method in his software.Reference Viechtbauer ²⁶ Partlett and RileyReference Partlett and Riley ⁹ also used modified versions of (5) in their work.

3.2.2 The parametric bootstrap method

Nagashima, Noma, and Furukawa published a bootstrap PI construction method,Reference Nagashima, Noma and Furukawa ⁸ which is a parametric method because it utilizes heavily the initial assumption of the normality of random effects. They show that based on this assumption, ${\mathrm{\varTheta}}_{\mathrm{New}}\sim \overline{\unicode{x3bc}}+\mathsf{Z}\mathsf{\tau }-{\mathrm{t}}_{\mathrm{K}-1}\sqrt{{\mathrm{Var}}_{\mathrm{HKSJ}}\left(\overline{\unicode{x3bc}}\right)}$ , where $\overline{\unicode{x3bc}}$ is an estimate for $\mathsf{\mu}$ using the theoretical weights, assuming ${\unicode{x3c3}^2}_{\mathrm{k}}$ and ${\unicode{x3c4}}^2$ are known, Z is a random variable with standard normal distribution, t_K-1 is a random variable with t distribution with K-1 degrees of freedom, and $\sqrt{{\mathrm{Var}}_{\mathrm{HKSJ}}\left(\overline{\unicode{x3bc}}\right)}$ is the theoretical standard error of $\overline{\unicode{x3bc}}$ defined by Hartung and KnappReference Hartung and Knapp ²² and Sidik and Jonkman.Reference Sidik and Jonkman ²³ As a next step, they define a plug-in estimator for ${\mathrm{\varTheta}}_{\mathrm{New}}$ and give an approximate predictive distribution for ${\widehat{\mathrm{\varTheta}}}_{\mathrm{New}}$ as

(6) $$\begin{align}{\widehat{\mathrm{\varTheta}}}_{\mathrm{New}}\sim \widehat{\unicode{x3bc}}+\mathsf{Z}{\widehat{\unicode{x3c4}}}_{\mathrm{UDL}}-{\mathrm{t}}_{\mathrm{K}-1}\sqrt{{\widehat{\mathrm{Var}}}_{\mathrm{HKSJ}}\left(\widehat{\unicode{x3bc}}\right)},\end{align}$$

where $\widehat{\unicode{x3bc}}$ and $\sqrt{{\widehat{\mathrm{Var}}}_{\mathrm{HKSJ}}\left(\widehat{\unicode{x3bc}}\right)}$ are computed as ${\unicode{x3c3}^2}_{\mathrm{k}}$ and ${\unicode{x3c4}}^2$ are estimated quantities, and ${\widehat{\unicode{x3c4}}}_{\mathrm{UDL}}$ is the untruncated DerSimonian and Laird estimator for τ.Reference DerSimonian and Laird ¹⁶ As the formula for ${\widehat{\mathrm{\varTheta}}}_{\mathrm{New}}$ contains three random variables (Z, t_K-1, and ${\widehat{\unicode{x3c4}}}_{\mathrm{UDL}}$ ), they generate B independent random realizations from these distributions and compute ${\widehat{\mathrm{\varTheta}}}_{\mathrm{New}}$ for each realization based on (6). Generating realizations from the distribution of ${\widehat{\unicode{x3c4}}}_{\mathrm{UDL}}$ is difficult; it is based on the distribution of the Q statistics described by Biggerstaff and Jackson.Reference Biggerstaff and Jackson ²⁷ For the boundaries of the (1 – α) expected coverage level PI, they simply take the α/2 and (1 − α/2) quantiles of the B realizations of ${\widehat{\mathrm{\varTheta}}}_{\mathrm{New}}$ .

3.3 PI estimators in common MA software

Since its introduction for the random-effects meta-analysis in the frequentist framework in 2009,Reference Higgins, Thompson and Spiegelhalter ⁵ the PI became an optional part in most of the commonly used MA software and program packages. The available and the default methods show a notable variety. We investigated the newest versions of the MA software that we think are commonly used and made a table showing the default and the other available PI methods (Table 1). The most common ones are the different HTS-type intervals. The parametric bootstrap method is only available in the meta R package.Reference Balduzzi, Rücker and Schwarzer ²⁹ The ensemble method is not implemented in any of these software, but is available freely in a supplementary spreadsheet file published along with their article.Reference Wang and Lee ¹²

Table 1 Default and other available prediction interval estimators in common meta-analysis software

Abbreviations: HTS, Higgins–Thompson–Spiegelhalter method; REML, Restricted maximum likelihood estimation of ${\boldsymbol{\unicode{x3c4}}}^{\mathbf{2}}$ parameter; DL, DerSimonian and Laird estimation of ${\boldsymbol{\unicode{x3c4}}}^{\mathbf{2}}$ parameter; HKSJ, Hartung–Knapp–Sidik–Jonkman variance estimation for the μ parameter; t_K-2, method calculated with t distribution with K-2 degrees of freedom; z, method calculated with standard normal distribution; “-,” denotes that the given option is not available.

* In metafor v. 4.8–0, the PI and CI are calculated based on the critical values from the standard normal distribution as the default setting. Otherwise, the PI is calculated based on the CI method selected, so that the PI is identical to the CI in the absence of observed heterogeneity. For more details, see the reference manual (https://cran.r-project.org/web/packages/metafor/metafor.pdf).

4.1 Generating study numbers and sample sizes

First, we fix K, the number of studies involved in the meta-analysis, as K in {3, 4, 5, 7, 10, 15, 20, 30, 100}; k represents the index of a given study in the meta-analysis, k = (1, 2, …,K). Then we determine the total sample size in a given study, N_k, by one of two ways. The sample sizes are either equal in all included studies, in this case N₁ = N₂ = … = N_K in {30, 50, 100, 200, 500, 1000, 2000}, or if they are mixed, the vector (50, 100, 500) is repeated periodically to N_K, namely N₁ = 50, N₂ = 100, N₃ = 500, N₄ = 50 … N_K. For the sample sizes in the two groups, N_{E, k} and N_{C, k}, we simply halve the total sample size: N_{E, k} = N_{C, k} = N_k/2.

4.2 Simulation of standard errors

We fix the population variance (V) of the continuous variable in the two groups as V_{E, k} = V_{C, k} = 10. This does not lead to loss of generality, as we will control the standard error of the mean difference by controlling the sample size. As a next step, we calculate the theoretical variance of the mean difference, ${\unicode{x3c3}^2}_{\mathrm{k}}$ , and simulate its observed value in the sample, ${\widehat{\unicode{x3c3}}}_{\mathrm{k}}^2$ . For these, we first generate the sample variance, $\widehat{\mathrm{V}}$ , in the two groups as a scaled random realization from a ${\unicode{x3c7}}^2$ distribution with appropriate degrees of freedom: ${\widehat{\mathrm{V}}}_{\mathrm{E},\mathrm{k}}\sim {\unicode{x3c7}}^2\left(\mathsf{df}={\mathrm{N}}_{\mathrm{E},\mathrm{k}}-\mathsf{1}\right)\ast \frac{{\mathrm{V}}_{\mathrm{E},\mathrm{k}}}{{\mathrm{N}}_{\mathrm{E},\mathrm{k}}-\mathsf{1}}$ and ${\widehat{\mathrm{V}}}_{\mathrm{C},\mathrm{k}}\sim {\unicode{x3c7}}^2\left(\mathsf{df}={\mathrm{N}}_{\mathrm{C},\mathrm{k}}-\mathsf{1}\right)\ast \frac{{\mathrm{V}}_{\mathrm{C},\mathrm{k}}}{{\mathrm{N}}_{\mathrm{C},\mathrm{k}}-\mathsf{1}}$ and compute the observed squared standard error as ${\widehat{\unicode{x3c3}}}_{\mathrm{k}}^2=\frac{{\widehat{\mathrm{V}}}_{\mathrm{E},\mathrm{k}}}{{\mathrm{N}}_{\mathrm{E},\mathrm{k}}}+\frac{{\widehat{\mathrm{V}}}_{\mathrm{C},\mathrm{k}}}{{\mathrm{N}}_{\mathrm{C},\mathrm{k}}}$ and its theoretical version as ${\unicode{x3c3}}_{\mathrm{k}}^2=\frac{{\mathrm{V}}_{\mathrm{E},\mathrm{k}}}{{\mathrm{N}}_{\mathrm{E},\mathrm{k}}}+\frac{{\mathrm{V}}_{\mathrm{C},\mathrm{k}}}{{\mathrm{N}}_{\mathrm{C},\mathrm{k}}}$ .

4.3 Simulation of true effects and observed effects

To simulate the true effects, ${\mathrm{\varTheta}}_{\mathrm{k}}$ , and their observed value, ${\widehat{\mathrm{\varTheta}}}_{\mathrm{k}}$ , we first fix τ², the variance of the true effects as τ² in {0.1, 0.2, 0.3, 0.5, 1, 2, 5}, and μ, the expected value of the true effects, as μ = 0. To simulate the true effects, we used the normal distribution and five other distributions to investigate how robust the methods are for the deviation from normality. We simulated ${\mathrm{\varTheta}}_{\mathrm{k}}$ from slightly, moderately, and highly skewed skew-normal distributions with skewness coefficients of 0.5, 0.75, and 0.99, respectively. We also simulated true effects from a bimodal distribution mimicking the case where there is an underlying factor that divides the studies into two subgroups. We generated random numbers from a bimodal distribution as a mixture of two normals using the R package EnvStats v. 2.8.1.,Reference Millard ³¹ and from the skew-normal distributions using the sn package v. 2.1.1.Reference Azzalini ³² As an extremity, we also investigated the case where the true effects follow the uniform distribution. The density functions of these distributions are visualized in the Online Material file.Reference Matrai, Koi, Sipos and Farkas ³³ As each investigated simulation factors are fully crossed, the number of the total investigated scenarios is as follows: (N) * (K) * (τ²) * (true effects distributions) = 8 * 9 * 7 * 6 = 3024.

4.4 Tested PI methods

We tested the parametric bootstrap method, the ensemble method, and six different versions of the HTS-type PI-s. The HTS-type methods differ in the τ² estimator, the $\mathsf{Var}\left(\widehat{\unicode{x3bc}}\right)$ estimator, and the used distribution. Table 2 summarizes these versions with their notation used in the following. To all the HTS-type methods that are constructed with t distribution, we will refer simply as HTS (t).

Table 2 HTS-type prediction interval estimators tested in the simulation

Abbreviations: HTS, Higgins–Thompson–Spiegelhalter method; REML, Restricted maximum likelihood estimation of ${\unicode{x3c4}}^2$ parameter; DL, DerSimonian and Laird estimation of ${\unicode{x3c4}}^2$ parameter; HKSJ, Hartung–Knapp–Sidik–Jonkman variance estimation for the μ parameter; t (df = K–2), method calculated with t distribution with K–2 degrees of freedom.

For each scenario, we simulated 5000 meta-analysis realizations and computed the PI with each investigated method, with the exception of the parametric bootstrap method, which we evaluated for 1000 realizations for each scenario. The parametric bootstrap method is a highly computationally intensive method requiring many bootstrap samples. This parameter can be set by the B parameter of the pima function in the pimeta R package. The default setting of the function is B = 25000, Nagashima et al. prepared their simulation with B = 5000 repetitions.Reference Nagashima, Noma and Furukawa ⁸ We also used the B = 5000 setting in our simulation. We used the high-performance computing system of the Hungarian Governmental Agency for IT Development (KIFÜ) to evaluate the parametric bootstrap method and a simple personal computer for the other methods. The simulation was programmed in R v.4.1.3. ³⁴

4.5 Heterogeneity measures

We measured the heterogeneity of the scenarios in two ways. We calculated the ratio $\mathsf{v}=\frac{\unicode{x3c4}^2}{\overline{\unicode{x3c3}^2}}$ following Partlett and Riley,Reference Partlett and Riley ⁹ where $\overline{\unicode{x3c3}^2}$ is the mean within-study variance. We also considered the mean of the observed I² measures and used the DL method to estimate τ² for the I² computation. We chose our simulation parameters so that both v and I² cover a wide range of values. Our investigated scenarios cover the range of 0.08–250 of v values and the 7%–99% range of the mean I² values.

4.6 Investigated performance measures

We assessed the coverage performance of the PI methods based on the empirical distribution of the coverage probabilities defined by equation (3). We visualized and examined these distributions on histograms. We computed the mean of the coverage probabilities and the median coverage probability, approximating E[C] and Median[C], respectively. To gain more insight into the behavior of the investigated intervals, we also analyzed the length of the PI. To be able to compare the observed interval length (ℓ) to a reasonable ideal length, we defined T, the theoretical length, as the difference between the 0.975 and 0.025 quantiles of the true effects distribution. We assessed ℓ by comparing it to this theoretical length and investigated the mean of the observed length relative to the theoretical length, approximating $\mathsf{E}\left[\frac{\ell }{\mathrm{T}}\right]$ .

5.1 Histogram of coverage probabilities

The histograms of coverage probabilities show the same pattern for each method: for any given sample size and τ² combination, the histograms are very left-skewed unless the number of involved studies is extremely large (K = 100). As there are more and more studies in the meta-analysis, the coverage distributions get less and less left-skewed and start cumulating on the nominal level of 95%, showing an ideal, fairly symmetric shape. This pattern is true in general for all scenarios, but for smaller heterogeneity (Figures 1 and 2), this process is slower, and even for a large number of studies, the distributions remain left-skewed. For higher heterogeneity scenarios (Figures 4 and 5), as the number of studies increases, the coverage distributions reach a concentrated, less skewed shape faster. For low heterogeneity scenarios, it can be observed that the distributions of the HTS type PIs are multi-modal; for example, the HTS-DL (t_K-2) method with K = 30 studies produces a mode of about 0.4 coverage (Figure 1c). The cause of this phenomenon is that the τ² parameter estimate is often 0 in these cases, which frequently results in low coverage. The ${\widehat{\unicode{x3c4}}}^2=\mathsf{0}$ cases are also the reason why the ensemble method gives a mode coverage probability of 0 for scenarios with low heterogeneity and low study number (Figure 1i).

5.2 Mean coverage probability

The HTS (t) methods give a very high, nearly 100% mean coverage for K = 3 studies, and then as the number of studies increases, their mean coverage drops below the nominal level. For low heterogeneity, these methods are unable to produce the nominal mean coverage, and when the number of involved studies is extremely large (K = 100), the mean coverage of these methods is still about 2–3% below 95% (Figure 3a). When the heterogeneity is larger (Figure 6a), this drop is less dramatic, and they yield a mean coverage close to the nominal level even for K = 15 studies. The HTS-DL (z) and the ensemble methods give a very low mean coverage for the lower heterogeneity scenarios, and they only retain mean coverages close to the nominal level when both heterogeneity and the number of studies are large. The mean coverage of the parametric bootstrap method remains close to the nominal level irrespective of heterogeneity and study number. This method is able to produce mean coverages close to 95% even when the included number of studies is very low (K = 3–5).

Figure 3 Mean coverage probability (a), Median coverage probability (b), and Mean observed length relative to the theoretical length (c) of the investigated prediction interval methods as a function of the number of involved studies (horizontal axis) for a low heterogeneity simulation scenario (N = 100, τ² = 0.2, I² = 33%, v = 0.5).

Abbreviations: HTS, Higgins–Thompson–Spiegelhalter method; DL, DerSimonian and Laird estimation of ${\unicode{x3c4}}^2$ parameter; HKSJ, Hartung–Knapp–Sidik–Jonkman variance estimation for the μ parameter; t_K-2, method calculated with t distribution with K-2 degrees of freedom; z, method calculated with standard normal distribution.

Figure 4 Histograms showing the coverage probability distribution of the HTS-DL (t_K-2), the parametric bootstrap, and the ensemble prediction interval methods for a high heterogeneity simulation scenario (N = 100, τ² = 1, I² = 71%, v = 2.5). The vertical green lines on the histograms indicate 95% coverage probability. The number of involved studies is constant in each column, showing the distribution for 5, 10, 30, and 100 studies. Results for the HTS-DL (t_K-2) method are represented in the first row of histograms with letters (a), (b), (c), and (d), the parametric bootstrap method is represented in the second row of histograms with letters (e), (f), (g), and (h), and the ensemble method is represented in the third row of histograms with letters (i), (j), (k), and (l).

Abbreviations: HTS-DL (t_K-2), Higgins–Thompson–Spiegelhalter method with the DerSimonian and Laird estimation of ${\unicode{x3c4}}^2$ parameter and using the t distribution with K-2 degrees of freedom.

Figure 5 Histograms showing the coverage probability distribution of the HTS-HKSJ (t_K-2), the HTS-DL (t_K-1), and the HTS-DL (z) prediction interval methods for a high heterogeneity simulation scenario (N = 100, τ² = 1, I² = 71%, v = 2.5). The vertical green lines on the histograms indicate 95% coverage probability. The number of involved studies is constant in each column, showing the distribution for 5, 10, 30, and 100 studies. Results for the HTS-HKSJ (t_K-2) method are represented in the first row of histograms with letters (a), (b), (c), and (d), the HTS-DL (t_K-1) method is represented in the second row of histograms with letters (e), (f), (g), and (h), and the HTS-DL (z) method is represented in the third row of histograms with letters (i), (j), (k), and (l).

Abbreviations: HTS-HKSJ (t_K-2), Higgins–Thompson–Spiegelhalter method with the Hartung–Knapp–Sidik–Jonkman variance estimation for the μ parameter and using the t distribution with K-2 degrees of freedom; HTS-DL (t_K-1), Higgins–Thompson–Spiegelhalter method with the DerSimonian and Laird estimation of ${\unicode{x3c4}}^2$ parameter and using the t distribution with K-1 degrees of freedom; HTS-DL (z), Higgins–Thompson–Spiegelhalter method with the DerSimonian and Laird estimation of ${\unicode{x3c4}}^2$ parameter and using the standard normal distribution.

Figure 6 Mean coverage probability (a), Median coverage probability (b), and Mean observed length relative to the theoretical length (c) of the investigated prediction interval methods as a function of the number of involved studies (horizontal axis) for a high heterogeneity simulation scenario (N = 100, τ² = 1, I² = 71%, v = 2.5).

Abbreviations: HTS, Higgins–Thompson–Spiegelhalter method; DL, DerSimonian and Laird estimation of ${\unicode{x3c4}}^2$ parameter; HKSJ, Hartung–Knapp–Sidik–Jonkman variance estimation for the μ parameter; t_K-2, method calculated with t distribution with K-2 degrees of freedom; z, method calculated with standard normal distribution.

5.3 Median coverage

The HTS (t) and the parametric bootstrap methods yield a higher than 95% median coverage in general, and they only retain 95% median coverage for scenarios with very large study numbers (Figures 3b and 6b). The parametric bootstrap method tends to give higher median coverages than the HTS (t) methods, especially for the lower heterogeneity scenarios. The HTS-DL (z) and the ensemble methods give a median coverage below 95% especially for the lower heterogeneity, lower study number scenarios, and as the study number increases, they approach the 95% level from below 95%.

5.4 Mean observed length relative to the theoretical length

The HTS (t) and the parametric bootstrap methods produce, on average, much longer intervals than the theoretical length, especially when the number of studies is low. As the study number increases, their mean observed length gets closer and closer to the theoretical length (Figures 3c and 6c). The parametric bootstrap method tends to give even longer intervals on average than the HTS (t) methods, especially when the heterogeneity is lower. The HTS-DL (z) and the ensemble methods yield a shorter mean observed interval compared to the HTS (t) and bootstrap methods when the number of involved studies is low (K < 15).

5.5 Impact of non-normality

When the true effects distributions depart from normal, the histograms are similar to the case when they follow the normal distribution. The only exception is where the random effects are drawn from the uniform distribution. In this case, a spike remains at >99% coverage for all methods irrespective of other circumstances, and the histograms do not approach a concentrated, symmetric form as the number of studies increases.Reference Matrai, Koi, Sipos and Farkas ³³

The mean and median coverages are very similar to the normal distribution case described above when the true effects distribution is skewed or bimodal. For the uniform scenario, methods give higher mean coverage compared to the normal case, and the median coverage is nearly 100% regardless of any other parameters. Methods tend to give slightly higher mean observed length compared to the theoretical length when the true effects have a non-normal distribution, especially when their distribution is uniform.