Methods to calculate uncertainty in the estimated overall effect size from a random-effects meta-analysis

3.1.1 WT methods

The WTz approach is the most popular technique for calculating a CI for μ,11 and a 95% CI is given by

where z_0.975 is the 0.975 quantile of the standard normal distribution. Any τ² estimator can be used when computing a WTz CI.50, 113, 114 This method often has coverage probability considerably below nominal 0.95 level14, 15, 45, 115-121 when k is small and/or τ² is large.13, 45, 52, 79, 118, 119, 122-125 Brockwell and Gordon13 stated that the greatest source of error in the method is the use of a normal approximation for . Despite the widespread use of the WTz method, it ignores uncertainty of the estimates of τ² and v_i in .

2. Wald-type t-distribution (WTt) CIs (method 2)

A slight modification of the WTz CI is the WTt approach, where the t distribution with k − 1 degrees of freedom is used, as opposed to the normal distribution. Although the two distributions converge asymptotically, the t quantile is larger than the z quantile associated with the WTz method. Hence, the WTt approach results in a wider CI and was proposed in order to increase the coverage probability, especially when the number of studies is small.117, 126 A 95% CI can be obtained by

where t_k−1,0.975 is the 0.975 quantile of the t distribution with k − 1 degrees of freedom. Any τ² estimator can be used to compute a WTt CI.50, 113

3. Quantile approximation (WTqa) CIs (method 3)

Brockwell and Gordon122 proposed the WTqa method as an alternative to WTz method in an attempt to achieve better coverage. The method resembles WTz and WTt, but instead of using normal or t distributions, it approximates the 0.025 and 0.975 quantiles of the distribution of the statistic

that are required for the 95% CI for μ. Hence, the WTqa uses different quantiles than the ones used in the WTz and WTt approaches. Let b_k be the quantile approximation function, which monotonically decreases as a nonlinear function of k; then, a 95% CI is calculated as

The quantiles b_k are estimated via a Monte Carlo simulation process of samples of the M statistic with b_k equal to the average of 0.025 and 0.975 absolute quantiles of the distribution, thus accounting for any small asymmetry in the distribution of M around zero.122 To obtain the function b_k, Brockwell and Gordon fit a regression equation for the quantiles as a function of k. The resulting regression equation (for k = 1, 2, …, 30) is

However, both number of studies k and the magnitude of τ² may impact on the performance of the WTqa method,122 and changes in the distribution of the within-study variances can importantly impact on b_k.127 Although WTqa approach has been criticized on the grounds that it is, at best, very difficult to obtain suitable critical values b_k that apply to all meta-analyses,127 we include it in this paper for completeness. As a conservative approach, Jackson and Bowden127 suggested the use of the standard normal quantile instead, and to assess the robustness of the findings via a sensitivity analysis of alternative quantiles. Brockwell and Gordon122 developed the WTqa method using the DerSimonian and Laird11 estimator of τ², but WTqa could, in principle, be implemented for any alternative τ² estimator. However, it is not advised to develop the WTqa CI further.127

3.1.2 Hartung-Knapp/Sidik-Jonkman CIs (methods 4 and 5)

Hartung and Knapp14 and Sidik and Jonkman15 independently introduced the Hartung-Knapp/Sidik-Jonkman (HKSJ) CI (method 4) to handle meta-analyses that include a small number of studies. This method is based on the S statistic, which follows a t distribution with k − 1 degrees of freedom,

with

and where Q_gen is the generalized Q statistic:

A 95% CI for μ is given by

Although the method does not take into account the uncertainty in τ², the use of a different statistical approximation to the usual WT CI may improve accuracy. Also, the HKSJ method can be applied with any τ² estimator and is exact for known variance components.14 For meta-analysis software where this method is not available yet, IntHout et al16 suggested an approach to convert WTz CIs easily to HKSJ CIs. The extension to meta-regression was investigated by Knapp and Hartung,125 and a generalization of the method to multivariate meta-analysis was explored by Jackson and Riley.128

When all variance components, including the between-study variances, are fixed and known, the expected value of Q_gen is k − 1, which equals the degrees of freedom of the associated χ² distribution.129-131 Hence, the small-sample adjustment q will tend to be close to 1. However, q may in fact turn out to be much smaller than 1, such as in cases where the effect sizes are very homogeneous or when the number and/or size of studies is small. This leads to a narrower CI than the WTt approach and can also lead to a narrower CI compared with the WTz method.132, 133 Although the t quantile is always larger than the z quantile associated with the WTz method, explaining in part why the HKSJ CI performs better than the WTz CI, in the case of , the HKSJ CI will be narrower than the WTz CI.

Wiksten et al132 show that if , then we estimate , and further that the variance of the estimate of μ simplifies to , with and Q ≤ (k − 1) when the DerSimonian and Laird11 estimator of τ² is used. Therefore, the variance of the estimated effect size is always smaller or equal for the HKSJ method than the WTz method when this estimator of the between-study variance is zero. The possibility that the variance of the estimated effect size from the HKSJ method can be smaller than the variance of the WTz method was discussed by Knapp and Hartung,125 who proposed a simple modification to the procedure. The authors125 suggested using q^* instead of q (method 5)

to ensure more conservative results. However, this practice may be overly conservative, leading to loss of power.134, 135

Sidik and Jonkman15 recommend using the HKSJ CI, but instead of , they suggest applying the sandwich variance estimator:

This is a robust estimator of , where the inverse of the study weights ( ) is estimated through the squared sample residuals ( ) from the data, rather than assuming .

However, Sidik and Jonkman120 state that is biased when k is small, and hence, they suggest a bias corrected estimator of (see Sidik and Jonkman120 for details). An alternative approach based on the expected information and on appropriately modified degrees of freedom of the t distribution was suggested by Kenward and Roger.136 These alternative expressions for could also be used in WT CIs but have not been adopted in practice, so we do not explore their use further here.

3.1.3 Likelihood-based methods

The profile likelihood (PL) method has been established in meta-analysis by Hardy and Thompson43 and is based on the likelihood ratio statistic, which unlike the WTz approach allows for asymmetric intervals. For y_i~Ν(μ, v_i + τ²), the log-likelihood function of the parameter vector (μ, τ²) is given by

Maximum likelihood estimates of (μ, τ²) can be found by maximizing lnL(μ, τ²) under the restriction τ² ≥ 0. The PL function is based on the log-likelihood function and uses an iterative process that provides CIs for μ that allow for the fact that τ² needs to be estimated as well. Since the PL approach profiles over τ², it accounts for, but does not fully allow for, the uncertainty in τ². This is because asymptotic results are required when using this method. However, the PL method is anticipated to be more accurate than the WT methods in smaller samples.

The profile log-likelihood for μ is defined as

where (μ) is the maximum likelihood estimator for τ²as μ varies.43 A 95% CI for μ can be obtained as the values that satisfy (see Hardy and Thompson43 equation 11)

where is the 0.05 quantile of the χ² distribution with 1 degree of freedom. It has been shown that for small k and τ², iterative algorithms are less likely to converge to a single value.13

2. Higher order likelihood inference methods (methods 7 and 8)

As Reid137 explains, the main asymptotic properties of likelihood-based inference include (a) consistent, asymptotically normal, and efficient maximum likelihood estimators; (b) an asymptotically normally distributed score statistic with mean zero; and (c) an asymptotic chi-squared distributed likelihood ratio statistic. For example, the PL CI in the previous section relies upon the third of these standard results. As Reid137 also explains, higher order asymptotic results for likelihood-based inference are also available. Some higher order likelihood inference methods have recently been applied to meta-analysis, which is a situation where they may be thought to be especially valuable. This is because the number of studies is often small, so that the commonly used “lower-order” asymptotic approximations to the likelihood function will be inadequate. Higher order likelihood–based methods therefore have the potential to produce more accurate results in meta-analysis, and several proposals for this have been made. We briefly summarize the methods here, but the details are technical, and so we refer the reader to the articles cited below for more information.

The Bartlett-type correction of the likelihood ratio statistic was first introduced by Bartlett (method 7).138 Noma17 explains how to apply this to random-effects meta-analysis and so use a higher-order approximation than the PL method above. Noma17 also explains how to use the score statistic to compute CIs and subsequently derives a higher order Bartlett-type adjustment to this score. Skovgaard139 proposed an alternative higher order approximation to the profile log-likelihood (method 8), and Guolo79 explains how to apply this to random-effects meta-analysis. For details on the method, we direct the reader elsewhere.79, 140 The higher order asymptotic methods have higher degree of accuracy, but in some cases (eg, when the between-study variance is close to zero), they may produce numerically unstable maximum likelihood estimates due to the discontinuity of the statistic.79, 140-142 In such cases, a bias reduction approach is suggested.143 Hence, the Bartlett-type correction (method 7) and the Skovgaard statistic (method 8) are the two main proposals for higher order approximations when using methods based on the PL.79

3.1.4 Henmi and Copas CIs (method 9)

Henmi and Copas (HC)30 propose an alternative strategy for obtaining intervals for μ that are less sensitive to publication bias than the widely used WTz method. Since the fixed-effect estimates assign larger weight to bigger studies, and study size is one component among others that is associated with the overall effect size in the presence of publication bias, this method centres the CI on a fixed-effect estimate. This is because the fixed-effect estimates are less sensitive to publication bias than the random-effects estimates. To allow for heterogeneity, they first estimate the variance of the fixed-effect estimate under the random-effects model as

Henmi and Copas30 then derive an approximation to the resulting pivot G that is used for making inferences about μ

assuming that the DerSimonian and Laird11 estimator of the between-study variance is used. Hence, approximate CIs can be computed. This can be thought of as a hybrid approach, where the fixed-effect estimate is accompanied by a CI that allows for between-study heterogeneity under the assumptions made in the random-effects model. A limitation of the approach is that the fixed-effect estimate is not fully efficient under the random-effects model, but Henmi and Copas30 argue that it is “better to use a method that is more robust to publication bias, even if this means sacrificing some efficiency under the standard setting.” A much simpler, but less accurate, way to implement Henmi and Copas' idea would be to assume that the pivot G approximately follows a standard normal distribution, but this would ignore all uncertainty in the between-study variance. This simpler approach could also be used with alternative estimators of the between-study variance. Alternatively, one could apply the IVher model suggested by Doi et al,144 which uses quasi-likelihood approaches and is performed under the fixed-effect assumption. Doi et al144 show that the IVher model favors larger trials, retains the nominal coverage probability, and exhibits lower variance of the overall effect size as opposed to the random-effects model irrespective of the degree of the between-study heterogeneity.

3.1.5 Biggerstaff and Tweedie CIs (method 10)

Biggerstaff and Tweedie (BT)53 proposed the use of different study-specific weights to those more conventionally used in the random-effects model, and estimated μ along with its variance using the weight , so as to acknowledge for τ² variability in the computation of CI for μ. Acknowledging the uncertainty in the weights allows greater uncertainty in the estimation of μ.53 The weights are the expected value of the random-effects weights (calculated using the estimated τ²) rather than the usual random-effects observed weights:

The weights depend on the density form of τ² and were derived using the DerSimonian and Laird11 estimator for the between-study variance. Alternative estimators could also be used, in principle, when using this method, provided that their distribution, and so the expected weights used by the method, can be evaluated. The variance of is estimated as

Assuming normality, a 95% CI can be obtained as

Biggerstaff and Tweedie53 use an approximate distribution to obtain the expected weights, but this has been improved upon by Preuß and Ziegler54 who used the exact weights through the exact cumulative distribution function of Q, where 145 Biggerstaff and Tweedie53 provided the algorithm to implement the method in SAS.

3.1.6 Resampling methods

Zeng and Lin (ZL)146 examine the distribution of the estimated overall effect size under the random-effects model and find that it is not asymptotically normally distributed for a finite number of studies k. This makes intuitive sense, because the textbook result that a linear combination of normal random variables is normally distributed requires that the coefficients in this linear combination are constants. When estimating the overall effect size however, these coefficients are proportional to the weights and so are functions of the estimated between-study variance. We require a large number of studies in order to estimate this variance accurately enough to treat the weights as fixed constants.

Recognizing that the estimated overall effect size is not asymptotically normally distributed for small k, Zeng and Lin146 suggest a resampling procedure to obtain the distribution of this estimate, assuming that the DerSimonian and Laird11 estimator of τ² is to be used. Briefly, they simulate values of τ² using the DerSimonian and Laird11 estimating equation (where the individual study results used in this estimation are simulated from the fitted random-effects model). They then simulate estimated average effect sizes using the sampled τ² to calculate the weights in the estimating equation for (as given in section 2.1, where the individual study results used in this estimation are simulated from the random-effects model centred at the estimated overall effect and where the between study variance is taken to be the sampled value used to compute the weights). By repeating both aspects of this sampling process B times, B² estimates provide an empirical distribution of estimated overall effects that can be used to compute CIs and make inferences.146 This resampling procedure could be modified to accommodate alternative estimators of τ², by instead calculating alternative estimates at the first stage, but this idea would need to be critically evaluated before it could be accepted.

2. Bootstrap CIs (methods 12 and 13)

Nonparametric bootstrapping is a way to approximate the sampling distribution of a statistic by resampling, from the sample itself, with replacement. Parametric bootstrapping instead samples from a fitted model. Both forms of bootstrapping can be used to make a variety of inferences but are most usually used to quantify the uncertainty in point estimates through the computation of standard errors and CIs. Briefly, bootstrap datasets are sampled (either nonparametrically, method 12,147, 148 or parametrically, method 13),149 from which the bootstrap statistics (the statistic of interest calculated using the bootstrap datasets) are calculated. Then, the empirical distribution of the bootstrap statistics is taken to approximate the distribution of the statistic of interest. Hence, measures of the uncertainty in the statistic, such as standard errors and CIs, can be calculated from this empirical distribution. In our context, this statistic is the estimated overall effect size.

Permutation tests have been suggested primarily to assess the true statistical significance of an observed finding under the null hypothesis of the absence of effect, especially in meta-analyses with a small number of studies.152 This method can be extended and used for calculating CIs for the effect size. These tests are especially appropriate when the included studies in a meta-analysis may not be considered randomly sampled from a larger population of studies. Confidence intervals can be constructed by inverting hypothesis tests, where parameter values that are not rejected by the hypothesis test lie within the corresponding CI.

Follmann and Proschan (FP)117 begin by considering a permutation method for testing the null hypothesis H₀ : μ = 0. Their argument assumes that the distributions of the outcome data y_i are symmetric117 and this is implied by the random-effects model. Under the null hypothesis, the sign of y_i is equally likely to be positive or negative for a symmetric y_i. There are 2^k possible permutations of the signs of the values of the outcome data y_i. We take Z^p = ( , ,..., ) to be the pth of these 2^k permutations; for example, with k = 5 studies, (+1, +1, −1, −1, +1) is one of the 32 possible permutations. We define to be the pth permutation of the outcome data corresponding to Z^p. The central idea is that under the null hypothesis, H₀ : μ = 0 and because the distributions of the y_i are symmetric, all 2^k permutations X^p are equally likely. Hence, Follman and Proschan117 propose the null distribution where all values of

are equally likely, where are the normalized (sum to one) random-effects weights described in section 2.1, where the between-study variance is estimated using X^p as outcome data. Hence, is the estimated average effect size using the pth permutation of signs of the outcome data. The two-sided P value based on the group permutation method proposed by Follmann and Proschan117 is simply the proportion of the absolute values of that are more than or equal to the absolute value of the estimated average effect under the random-effects model using the observed data. Follman and Proschan117 describe their procedure in terms of the DerSimonian and Laird11 estimator of the between-study variance, but, in principle, alterative estimators could be used. If k is too large for all permutations to be evaluated, then the permutation distribution can be approximated by instead simulating a large number of permutations.117

The procedure described above can be extended so that we instead test the hypothesis H₀ : μ = c, and invert this hypothesis test to give the bounds of CIs. For further details on the FP method, we refer the reader to Follman and Proschan.117 This method has been suggested as an alternative approach to the WT and likelihood-based approaches that assume normality of the observed effects, but it can be computationally demanding. The discrete nature of the permutation distribution will ensure that the CI maintains the desired coverage probability, but in general, this coverage probability will be larger than the nominal level.

3.1.7 Bayesian credible intervals (method 15)

Bayesian credible intervals (CrIs) for the overall effect size can be obtained within a Bayesian framework using specialized software and the Markov chain Monte Carlo (MCMC) technique, such as WinBUGS153 or SAS PROC MCMC. Some advantages of the Bayesian approach include (1) incorporation of uncertainty in model parameters (μ, τ²); (2) derivation of CrIs from the posterior distribution; and (3) use of informative prior distributions on the model parameters. However, the use of informative priors for the effect size parameters has been discouraged by some researchers due to potential inclusion of bias.154 The use of vague priors allows the analysis to remain data driven. On the contrary, the use of informative priors for the between-study variance has been suggested to increase confidence in the overall effect size, especially when few studies are included in a meta-analysis.20, 21 Informative priors for τ² under several treatment comparison types and outcome settings are available for dichotomous data20 and for continuous data.21 Friede et al155 suggest Bayesian CrIs perform well even in rare diseases with a small number of studies when the appropriate prior for τ² is applied. In rare diseases and small populations, the use of half-normal priors, with expectation 0 and variance 0.25 or 1 for τ², has been recommended when log-odds ratios are used to measure the effect size.155 Vague priors can also be applied for τ², but caution is needed, as results are sensitive to the prior specification, especially when the number of studies is small.156 This is because the choice of prior may impact on the estimation of the between-study variance and consequently on the estimated overall effect size and the width of its CI. Other difficulties that have been associated with the derivation of Bayesian CrIs include the complication of determining whether convergence is achieved, the need to burn-in when using MCMC, and the impact of MC error. Alternative methods to implement a Bayesian meta-analysis are available by using numerical integration, importance sampling, and data augmentation as described by Turner et al57 and Rhodes et al.157 For practical application, the R package bayesmeta is available.58

3.2.1 WT methods (methods 1, 2, and 3)

The performance of the popular WTz method has been assessed in several studies, and it is poor when compared with other methods. Simulations suggest that the WTz performs worse in terms of coverage for small numbers of studies (k < 16) compared with the PL and the WTt methods, whereas for large k, all three methods perform well.158 The performance of the WTz method though does not only depend on the number of studies but also depend on the τ² estimator employed and its magnitude.45 The WTz coverage has been found to differ by up to 0.05 between different τ² estimators, up to 0.30 between meta-analysis samples, and up to 0.20 across between-study variance values ranging from small to large τ². Coverage has been found to be as low as 0.65 (at 0.95 nominal level) when I² (defined as the percentage of the total variability in a set of effect sizes that is due to between-study variability beyond what is expected by within-study random error) is 90% and two or three studies are included in a meta-analysis, but it tends towards the nominal level as the number of studies increases.113

To increase coverage, the t distribution can be used, which produces wider CIs than those obtained by the standard normal distribution, especially when τ² and k are small.15, 117 The coverage probability is therefore higher with the WTt approach, but it depends on the estimator and the magnitude of τ², as well as on the number of studies.45 Simulation showed that the WTt CI is less affected by the number of studies compared with the WTz CI.52 Although WTt coverage may be more robust to changes in the τ² magnitude compared with WTz when few studies are included in a meta-analysis, it has been found to differ by up to 0.05 depending on the τ² estimator used and the number of studies.113 For large meta-analysis samples (eg, k ≥ 20), the coverage of the 95% WTt CI may be below the nominal level, but it becomes conservative (close to 1) when k is small.113, 122, 158

Alternatively, the WTqa method is easy to implement and produces intervals with better coverage in comparison with the WTz method.122 A simulation study45 showed that different estimators of the between-study variance may impact on coverage and that the WTqa method is associated with higher coverage than WTz and HKSJ CIs, but the HKSJ method produced values closer to the nominal level. The same study showed that the WTqa method has similar coverage to the WTt method. For small k, coverage of the WTt method is well above the nominal level and higher than that for the WTqa method, but as k increases, the differences in coverage are not so important.122 Simulations have also shown that WTqa outperforms WTz, PL, and ZL approaches but it is very conservative.146

3.2.2 HKSJ methods (method 4 and 5)

The HKSJ approach (method 4) is often preferred, as in the case of small k, it is conservative and on average produces wider CIs with more adequate type I error compared with the WTz method.16, 119, 141, 152, 159 The HKSJ method provides exact inference when all study sizes are equal and the random-effects model is true, resulting in better inference than WTz,160 but also provides more accurate inference in small meta-analyses with different study sizes.14, 15, 116 Several studies suggested that the HKSJ method has coverage close to the nominal level and that it is not influenced by the magnitude or estimator of τ².16, 45, 113, 115, 119, 121, 125, 132, 135, 155 Nevertheless, Knapp and Hartung125 recommend using the PM161 estimator for the between-study variance along with the HKSJ method to obtain CIs for μ so as to get a cohesive approach based on Q_gen.161, 162 Sanchez-Meca and Marin Martinez45 recommend using the HKSJ method, as it is additionally insensitive to the number of trials. Simulation studies suggest that HKSJ has good coverage when the effect measure is the log-odds ratio,15, 119 the standardized mean difference,121 the mean difference, and the risk difference.118 The coverage of the 95% HKSJ CI is generally better than the WTz and WTt coverages, but it is suboptimal in meta-analyses with binary outcomes and rare events, as shown in simulated meta-analyses where the odds-ratio was used as the measure of effect.113

A real-life data study of 920 Cochrane meta-analyses with k ≥ 3 showed that the WTz method yielded more often statistically significant results compared with the HKSJ method (45% vs 35% of meta-analyses).163 IntHout et al16 found similar results in their real-life data study with 434 Cochrane meta-analyses with dichotomous data (43% vs 34%) and 255 Cochrane meta-analyses with continuous data (51% vs 40%). It is recommended that caution is needed when fewer than five studies of unequal sizes are included in the meta-analysis.16 Wiksten et al132 in their empirical evaluation including 157 meta-analyses with dichotomous data and k ≥ 4 found that in the presence of heterogeneity (using the DerSimonian and Laird11 estimator [ ] or the Cochran Q statistic [P < 0.10]),11, 164 the P value for the overall effect size was typically greater when using the HKSJ than the WTz method. However, they comment that the HKSJ method is not always more conservative when .

It has been shown that in the absence of heterogeneity, the coverage of HKSJ may be smaller than the WTz coverage providing narrower CIs.15, 115, 118, 121, 132, 133, 165 This was more prevalent in cases with rare events.132 Jackson et al133 raise a variety of concerns about the use of the HKSJ method, including (1) the implications of the modification for any given meta-analysis are hard to predict, (2) HKSJ can result in shorter CIs for the overall effect size than the WTz method, and (3) the coverage of the HKSJ CI might be anticipated to be low when . However, in simulation studies conducted by Röver et al,135 Viechtbauer et al,134 and Sanchez-Meca and Marin Martinez,45 HKSJ worked well even in the absence of heterogeneity. This is in line with the simulations by Gonnermann et al,166 but in the presence of only two studies for τ² = 0, HKSJ is associated with very low power compared with WTz (15% vs 60%), which may be due to the wider CI, whereas for mild to moderate τ², both methods have poor control of type I error. A simulation study compared HKSJ with the small sample modification suggested by Knapp and Hartung125 and indicated that the use of the modified HKSJ (method 5) is preferable when few studies of varying size and precision are available.135 Another simulation study suggested the use of the modified HKSJ approach instead of the common HKSJ and WTz approaches when dichotomous data are considered.167 However, for few studies (and particularly for k = 2) and as the between-study variance decreases, the modified HKSJ tends to be overconservative, and selection between the methods is a matter of power vs type I error.133-135

3.2.3 Likelihood-based methods (methods 6, 7, and 8)

The PL method is often preferred to the WTz method, as it is associated with a higher coverage closer to the nominal level, even when k is relatively small.122, 143 Jackson et al158 showed that the PL method performed well and better than the WTz and WTt methods in meta-analyses with few studies (k ≤ 8) with coverage close to the nominal level. However, coverage decreases as τ² increases and/or k decreases.43 Simulations suggest that the PL CI is less affected by the number of studies in a meta-analysis compared with the WTz CI, but both WTz and PL have poor coverage control, as they yield values below the nominal level.52

Simulations found that the Bartlett-type correction CI (method 7) improves coverage properties over the WTz, WTt, and PL methods that their coverage deviates the nominal level as τ² increases and/or k decreases.17, 52 Although the Bartlett-type correction CI has a satisfactory power compared with the WTz, WTt, and PL CIs52 and performs well when τ² = 0,142 caution is needed for k ≤ 5, as it tends to be overconservative.17 The Skovgaard statistic CI (method 8) is associated with coverage closer to the nominal level compared with the WTz and PL CIs, which is remarkable for small k.17, 79, 141 Both the Skovgaard statistic CI and the Bartlett-type correction CI perform very satisfactorily regarding coverage and yield similar results.141

3.2.4 HC and BT (methods 9 and 10)

Simulations showed that in the absence of publication bias and for k > 10, the HC method yields better coverage than WTz, HKSJ, PL, and BT methods, whereas for k < 10, the HKSJ and PL methods perform best.30 The same study showed that when publication bias is present and for k > 10, HC improved coverage compared with WTz, HKSJ, PL, and BT methods and showed less bias than the fixed-effect model. Also, the WTz and BT methods have comparable coverage probabilities with coverage below the nominal level,54, 122 but coverage is increased for the exact weights.54

3.2.5 Resampling methods (methods 11, 12, 13, and 14)

Zeng and Lin146 showed that the ZL CI outperforms both WTz and PL CIs for small k in terms of coverage. Another simulation study showed that the FP CI controls coverage better than WTz, WTt, and PL and is closely followed by the Bartlett-type correction CI, but the latter is slightly more powerful especially for small k.52 The same study showed that the FP CI and the Bartlett-type correction CI were less affected by the number of studies than WTz, PL, and WTt methods.

3.2.6 Bayesian CrIs (method 15)

Simulation studies showed that Bayesian intervals produce intervals with coverage closer to the nominal level compared with the HKSJ, modified HKSJ, and PL CIs,155, 168 and they tend to be smaller than the HKSJ CI even in situations with similar or larger coverage.155 However, the performance of the Bayesian CrIs may vary depending on the prior assigned to the between-study variance.156