Median bias reduction in random-effects meta-analysis and meta-regression
Random-effects models are frequently used to synthesise information from different studies in meta-analysis. While likelihood-based inference is attractive both in terms of limiting properties and of implementation, its application in random-effects meta-analysis may result in misleading conclusions, especially when the number of studies is small to moderate. The current paper shows how methodology that reduces the asymptotic bias of the maximum likelihood estimator of the variance component can also substantially improve inference about the mean effect size. The results are derived for the more general framework of random-effects meta-regression, which allows the mean effect size to vary with study-specific covariates.
Keywords: Adjusted score equations; Heterogeneity; Mean bias reduction; Penalized likelihood; Random effects
= 1 \pdfinclusioncopyfonts=1
Meta-analysis is a core tool for synthesizing the results from independent studies investigating a common effect of interest. One of the main challenges when combining results from multiple studies is the variability or heterogeneity in the design and the methods employed in each study. Accounting for and quantifying that heterogeneity is critical when drawing inferences about the common effect. In this direction, DerSimonian & Laird (1986) introduced the random-effects meta-analysis model, which expresses the heterogeneity between studies in terms of a variance component that can be estimated through standard estimation techniques.
Nevertheless, there is ample evidence that frequentist inference based on random-effects meta-analysis can be problematic in the usual meta-analytic scenario where the number of studies is small or moderate. Specifically, the estimation of the heterogeneity parameter can be highly imprecise, which in turn results in misleading conclusions (Van Houwelingen et al., 2002; Guolo & Varin, 2015; Kosmidis et al., 2017). Examples of recently proposed methods that attempt to improve inference are the resampling (Jackson & Bowden, 2009) and double resampling approaches(Zeng & Lin, 2015), and the mean bias-reducing penalized likelihood (BRPL) approach in Kosmidis et al. (2017). Specifically, Kosmidis et al. (2017) show that maximization of the BRPL results in an estimator of the heterogeneity parameter that has notably smaller bias than maximum likelihood (ML) with small loss in efficiency, and illustrate that BRPL-based inference outperforms its competitors in terms of inferential performance.
Kenne Pagui et al. (2016) show that under suitable conditions third-order median unbiased estimators can be obtained by the solution of a suitably adjusted score equation. Such median bias-reduced estimators have component-wise and to third-order the same probability of over- and under-estimating the true parameter value. A key property of these estimators, not shared with the mean bias-reduced ones, is that any monotone component-wise transformation of the estimators results automatically in median bias-reduced estimators of the transformed parameters (Kenne Pagui et al., 2016). Such equivariance property can be useful in the context of random-effects meta-analysis where the Fisher information and, hence, the asymptotic variances of various likelihood-based estimators depend only on the heterogeneity parameter.
In this paper, we derive the median bias-reducing adjusted score functions for random-effects meta-analysis and meta-regression. The adjusted score functions are found to correspond to a median BRPL, whose logarithm differs from the logarithm of the mean BRPL in Kosmidis et al. (2017) by a simple additive term that depends on the heterogeneity parameter. Since the adjustments to the score function for mean and median bias reduction are both of order , the same arguments as in Kosmidis et al. (2017) are used to obtain a median BRPL ratio statistic with known asymptotic null distribution that can be used for carrying out hypothesis tests and constructing confidence regions or intervals for either the fixed-effect or the heterogeneity parameter. Simulation studies and real data applications are used to assess the performance of estimation and inference based on the median BRPL, and compare it to recently proposed alternatives, including the mean BRPL. The results provide evidence on the effectiveness of median bias reduction in improving estimation and likelihood-based inference.
2 Cocoa intake and blood pressure reduction data
Consider the setting in Bellio & Guolo (2016) who carry out a meta-analysis of five randomized controlled trials from Taubert et al. (2007) on the efficacy of two weeks of cocoa consumption on lowering diastolic blood pressure. Figure 1 is a forest plot with the estimated mean difference in diastolic blood pressure before and after cocoa intake from each study, and the associated Wald-type confidence intervals. Four out of the five studies reported a reduction of diastolic blood pressure from cocoa intake.
The random-effects meta-analysis model is used to synthesize the evidence from the five studies. In particular, let be the random variable representing the mean difference in the diastolic blood pressure after two weeks of cocoa intake in the th study. We assume that are independent random variables where has a Normal distribution with mean the overall effect and variance , where is the estimated standard error of the effect from the th study and is the heterogeneity parameter.
The forest plot in Figure 1 depicts nominally confidence intervals for using various alternative methods. As is apparent, the conclusions when testing the hypothesis can vary depending on the method used. More specifically, the Wald test using the ML estimates, the DerSimonian and Laird method (DerSimonian & Laird, 1986), double resampling (Zeng & Lin, 2015), and the likelihood ratio (LR) test give evidence that there is a relationship between cocoa consumption and diastolic blood pressure, with -values , , , , respectively. On the other hand, Knapp and Hartung’s method (Knapp & Hartung, 2003), the mean BRPL ratio (Kosmidis et al., 2017), the Bartlett-corrected LR (Huizenga et al., 2011), and Skovgaard’s test suggest that the evidence that cocoa consumption affects diastolic blood pressure is weaker, with -values , , , , respectively.
3 Random-effects meta-regression model
Let and denote the estimate of the effect from the th study and the associated within-study variance, respectively, and denote a -vector of study-specific covariates that can be used to account for the heterogeneity across studies.
The within-study variances are usually assumed to be estimated well-enough to be considered as known and equal to the values reported in each study. Then the observations are assumed to be realizations of the random variables , which are independent conditionally on independent random effects . The conditional distribution of given is , where is an unknown -dimensional vector of fixed effects. The random effects are typically assumed to be independent with having a distribution, where is a parameter that attempts to capture the unexplained between-study heterogeneity. In matrix notation, the random-effects meta-regression model has
where , is the model matrix with in its th row, and is a vector of independent errors each with a distribution and independent of . Under this specification, the marginal distribution of is multivariate normal with mean and variance-covariance matrix , where is the identity matrix and . The random-effects meta-analysis results as a special case of meta-regression, by setting to be a column of ones.
The log-likelihood function for is , where denotes the determinant of and . The gradient of the log-likelihood (score function) is
and the ML estimator is obtained as the solution of , where denotes a -dimensional vector of zeros.
4 Median bias reduction
4.1 The method
A popular method for reducing the mean bias of ML estimates in regular statistical models is through the adjustment of the score equation (Firth, 1993; Kosmidis & Firth, 2009). Kenne Pagui et al. (2016) propose an extension of the adjusted score equation approach which can be used to obtain median bias-reduced estimators. Specifically, under the model, the new estimator has a distribution with median closer to the “true” parameter value than the ML estimator. Kenne Pagui et al. (2016) consider the median as a centering index for the score, and the adjusted score function for median bias reduction then results by subtracting from the score its approximate median, obtained using a Cornish-Fisher asymptotic expansion.
Let be the observed information matrix (see, Appendix for its expression), and be the expected information matrix
with th column . Let also and be the th column and the th diagonal element of , with . Kenne Pagui et al. (2016) show that a median bias-reduced estimator can be obtained by solving an adjusted score equation of the form , where the additive adjustment to is , in the sense that is bounded in absolute value by a fixed constant after a sufficiently large value of . The median bias-reducing adjustment has th element
In the context of meta-regression values of and in correspond to the elements of parameter , and correspond to parameter . Given that is of order , has the same asymptotic distribution as (Kenne Pagui et al., 2016), i.e. multivariate normal with mean and variance-covariance matrix , which can be consistently estimated with .
After some algebra (see Appendix for details) the median bias-reducing adjustment for the random-effects meta-analysis and meta-regression models has the form
where . Substituting (5) in the expression for gives that the median bias-reducing adjusted score functions for and are and
4.2 Computation of median bias-reduced estimator
A direct approach for computing the estimator is through a modification of the two-step iterative process in Kosmidis et al. (2017). At the th iteration ()
calculate by weighted least squares as
solve with respect to , where .
In the above steps, is the candidate value for at the th iteration and is the candidate value for at the th iteration. The equation in step 2 is solved numerically, by searching for the root of the function in a predefined positive interval. For the computations in this manuscript we use the DerSimonian & Laird (1986) estimate of as starting value . The iterative process is then repeated until the components of the score function are all less than in absolute value at the current estimates.
4.3 Median bias-reducing penalized likelihood
Although it is not generally true that is the gradient of a suitable penalized log-likelihood, in this case is the gradient of the median BRPL
An advantage of the median BRPL estimators over mean BRPL ones is that the former are equivariant under monotone component-wise parameter transformations (Kenne Pagui et al., 2016). In the context of random-effects meta-analysis and meta-regression, this equivariance implies that not only we get a median bias-reduced estimator of , but we also get median bias-reduced estimates of the standard errors for by calculating the square roots of the diagonal elements of in (3) at . This is because is a function of only, and moreover the square roots of the diagonal elements of are monotone functions of .
4.4 Penalized likelihood-based inference
For inference about either the components of the fixed-effect parameters or the between-study heterogeneity we propose the use of the median BRPL ratio. If and is the maximizer of for fixed , then the same arguments as in Kosmidis et al. (2017) can be used to show that the logarithm of the median BRPL ratio statistic
has a asymptotic distribution, as goes to infinity. Specifically, the adjustment to the score function is additive and of order . As a result, the extra terms in the asymptotic expansion of the logarithm of the median BRPL that depend on the penalty and its derivatives disappear as information increases, and the expansion has the same leading term as that of the log-likelihood (see, for example, Pace & Salvan, 1997, Section 9.4).
5 Cocoa intake and blood pressure reduction data (revisited)
The ML estimate, the maximum mean BRPL estimate and the maximum median BRPL estimate of the heterogeneity parameter in the meta-analysis model in Section 2 are , , and , respectively. The estimates of the common effect are , and , with standard errors , , and , respectively. The bias-reduced estimates of and, as a consequence, the corresponding estimated standard errors for are larger than their ML counterparts, which is typical in random-effects meta-analysis. The iterative process used for computing the ML, maximum mean BRPL, and maximum median BRPL estimates converged in 4, 5, and 11 iterations, respectively. The computational run-time for the two-step iterative process which computes the ML, maximum mean BRPL, and maximum median BRPL estimates is , , and seconds, respectively.
Figure 2 shows the value of LR, mean BRPL and median BRPL ratio statistic in (7) for a range of values of , when is either or . Here and in the following simulation studies we compare median BRPL ratio statistic with only LR and mean BRPL ratio statistics because the mean BRPL ratio statistic is a strong competitor against other alternatives in terms of inferential performance (Kosmidis et al., 2017). The horizontal line in Figure 2 is the quantile of the limiting distribution, and its intersection with the values of the statistics results in the endpoints of the corresponding confidence intervals. For both and , the confidence intervals based on the LR statistic are the narrowest and the confidence intervals based on the median BRPL ratio statistic are the widest. Specifically, the confidence intervals for are , , and for the median BRPL ratio statistic, mean BRPL ratio statistic, and LR statistic, respectively. The corresponding confidence intervals for are , , and , respectively. Contrary to the LR test, the mean BRPL and median BRPL ratio tests suggest that there is only weak evidence that cocoa consumption affects diastolic blood pressure with -values of and .
In order to further investigate the performance of the three approaches to estimation and inference, we performed a simulation study where we simulated independent samples from the random-effects meta-analysis model with parameter values set to the ML estimates reported earlier, i.e. and . Figure 3 shows boxplots of the estimates of and calculated from each of the simulated samples. The distributions of the three alternative estimators for are similar. On the other hand, the ML estimator of has a large negative mean bias, maximum median BRPL tends to over-correct for that bias, while maximum mean BRPL almost fully corrects for the bias of ML estimator. The distribution of the median BRPL estimates has the heaviest right tail. The simulation-based estimates of the probabilities of underestimation for , , and are 0.708, 0.591 and 0.493 for the ML, maximum mean BRPL, and maximum median BRPL, respectively, illustrating how effective maximizing the median BRPL in (6) is in reducing the median bias of the maximum likelihood estimator of .
The simulated samples were also used to calculate the empirical -value distribution for the two-sided tests that each parameter is equal to the true values based on the LR statistic, the mean BRPL ratio statistic, and the median BRPL ratio statistic. Table 1 shows that the empirical -value distribution for the mean and median BRPL ratio statistics are closest to uniformity, with the latter being slightly more conservative than the former. The coverage probability of the confidence intervals of based on the mean BRPL ratio and the median BRPL ratio are notably closer to the nominal level than those based on the likelihood ratio. Specifically, the coverage probabilities for are , , and for LR, mean BRPL ratio, and median BRPL ratio respectively, and the corresponding coverage probabilities for are , , and , respectively.
|Mean BRPL ratio||1.6||3.7||6.7||12.2||28.4||52.8||76.6||90.9||95.5||97.9||99.1|
|Median BRPL ratio||0.6||1.8||4.1||8.6||23.1||48.5||74.2||90.0||95.0||97.5||99.1|
6 Simulation study
More extensive simulations under the random-effects meta-analysis model (1) are performed here using the design in Brockwell & Gordon (2001). Specifically, the data , , are simulated from model (1) with true fixed-effect parameter . The within-study variances are independently generated from a distribution and are multiplied by before restricted to the interval . Eleven values of the between-study variance ranging from to are chosen, and the number of studies ranges from to . For each combination of and considered, we simulated data sets initializing the random number generator at a common state. The within-study variances where generated only once and kept fixed while generating the samples.
Zeng & Lin (2015) compared the performance of their proposed double resampling method with the DerSimonian & Laird (1986) method, the profile likelihood method in Hardy & Thompson (1996), and the resampling method in Jackson & Bowden (2009) and showed that the double resampling method improves the accuracy of statistical inference. Based on these results Kosmidis et al. (2017) compared the performance of their mean BRPL approach with the double resampling method and illustrated that the former results in confidence intervals with coverage probabilities closer to the nominal level that the alternative methods.
We take advantage of the results reported in Zeng & Lin (2015) and Kosmidis et al. (2017) and evaluate the performance of estimation and inference based only on the median BRPL with that based on the likelihood and the mean BRPL. The estimators of the fixed and random-effect parameters obtained from the three methods are calculated using variants of the two-step algorithm described in Section 4.2. In the second step of the algorithm the candidate values for the ML, and maximum mean and median BRPL estimators of the between-study variance are calculated by searching for the root of the partial derivatives of , , and with respect to , in the interval .
First, we compare the performance of the ML, maximum mean BRPL and maximum median BRPL estimators in terms of percentage of underestimation. Figure 4 shows that the median bias-reducing adjustment is the most effective in reducing median bias even for small values of . As expected, the ML and maximum mean BRPL estimators also approach the limit of underestimation as grows, with the latter being closer to than the former. Figure 5 shows that maximum median BRPL is also effective in reducing the mean bias of the ML estimator of but only for moderate to large values of , while maximum mean BRPL results in estimators with the smallest bias.
Figures 6 and 7 show the estimated coverage probability for the one-sided and two-sided confidence intervals for based on the LR, mean BRPL ratio and median BRPL ratio statistics at the nominal level. Figure 8 shows the estimated coverage probability for the two-sided confidence intervals for based on the LR, mean BRPL ratio and median BRPL ratio statistics at the nominal level. For small values of or small and moderate number of studies the empirical coverage of the intervals is larger than the nominal level. In general, the confidence intervals based on mean and median BRPL ratio have empirical coverage that is closer to the nominal level with the latter having generally better coverage. The differences between the three methods diminish as the number of studies increases.
Figures 9 and 10 give the power of the LR, the mean BRPL ratio, and the median BRPL ratio tests for testing the null hypothesis against various alternatives. Specifically, we simulated data sets under the alternative hypothesis that parameter is equal to , where ranges from to . In Figure 9 the power is calculated using critical values of the the asymptotic null distribution of the statistics. In Figure 10 the power is calculated using critical values based on the exact null distribution of each statistic, obtained by simulation under the null hypothesis. In this way, the three tests are calibrated to have size .
Figure 9 shows that the three tests have monotone power and for small values of the LR test yields the largest power. This is because the LR test is oversized, while the mean and median BRPL ratio tests are slightly more conservative and this conservativeness comes at the cost of lower power. As the number of studies increases the three tests approach the nominal size and provide similar power. The use of the exact critical values in Figure 10 allows us to compare the performance of the tests without letting the oversizing or the conservativeness of a test skew the power results. Figure 10 shows that the power of the median BRPL ratio test is almost identical to that of the mean BRPL ratio test, and both tests have larger power than the LR test. Again, inference based on either of the two penalized likelihoods becomes indistinguishable from classical likelihood inference as the number of studies increases.
Across all and values considered, the average number of iterations taken per fit for the two-step iterative process to converge is 6.20, 5.75 and 5.86 iterations for ML, maximum mean BRPL, and maximum median BRPL, respectively. The average computational run-times for ML, maximum mean BRPL, and maximum median BRPL are 0.005, 0.021, and 0.017 seconds, respectively. Figures 1 and 2 in the Supplementary material show the average number of iterations and the average computational run-time taken per fit for the two-step iterative process to converge for each value of and used in the simulation study. The results show that in all cases estimation is achieved rapidly and after a small number of iterations for all three methods, with only negligible overhead with the two bias reducing methods.
7 Meat consumption data
Larsson & Orsini (2014) investigate the association between meat consumption and relative risk of all-cause mortality. The data consists of prospective studies, eight of which are about unprocessed red meat consumption and eight about processed meat consumption. Figure 11 displays the information provided by each study in the meta-analysis. The results from the studies point towards the conclusion that high consumption of red meat, in particular processed red meat, is associated with higher all-cause mortality.
We consider the random-effects meta-regression model, assuming that has a , where is the random variable representing the logarithm of the relative risk reported in the th study, and takes value if the consumption in the th study is about processed red meat and if it is about unprocessed meat . Table 2 gives the ML estimates, the mean BRPL estimates, and the median BRPL estimates of , and , along with the corresponding standard errors for and . The median BRPL estimate of and the standard errors of the fixed-effect parameters are the largest. The iterative process for computing the ML, maximum mean BRPL, and maximum median BRPL estimates converged in eight, nine, and twelve iterations, in , , and seconds, respectively.
The LR test indicates some evidence for a higher risk associated to the consumption of red processed meat with a -value of . On the other hand, the mean and median BRPL ratio tests suggest that there is weaker evidence for higher risk, with -values of and , respectively. Skovgaard’s test also gives weak evidence for higher risk with -value 0.073.
|ML||0.099 (0.044)||0.106 (0.061)||0.009|
|Mean BRPL||0.095 (0.050)||0.110 (0.069)||0.012|
|Median BRPL||0.093 (0.052)||0.111 (0.072)||0.013|
Similar to Section 5, we performed a simulation study in order to further investigate the performance of the three methods in a meta-regression context. We simulated independent samples from the meta-regression model at the ML estimates reported in Table 2. Figure 12 shows boxplots of the estimates of , , and . Maximum likelihood underestimates the parameter , while mean BRPL and median BRPL almost fully compensate for the negative bias of ML estimates, with the latter having a slightly heavier right tail. The percentages of underestimation are , , and for the ML, maximum mean BRPL, and maximum median BRPL estimators, respectively.
The simulated samples were also used to calculate the empirical -value distribution for the tests based on the likelihood, mean BRPL and median BRPL ratio statistics. Table 3 shows that the empirical -value distribution for the median BRPL ratio statistic is the one closest to uniformity.
|Mean BRPL ratio||1.3||3.0||5.6||11.1||25.9||49.8||73.8||89.0||94.2||96.9||98.6|
|Median BRPL ratio||1.0||2.5||4.9||9.9||25.1||49.7||74.7||89.8||94.8||97.5||98.9|
8 Concluding remarks
In this paper we derive the adjusted score equations for the median bias reduction of the ML estimator for random-effects meta-analysis and meta-regression models and describe the associated inferential procedures.
We show that the solution of the median bias-reducing adjusted score equations is equivalent to maximizing a penalized log-likelihood. The logarithm of that penalized likelihood differs from the logarithm of the mean BRPL in Kosmidis et al. (2017) by a simple additive term. The computation of the maximum median BRPL estimators can be performed through a two-step iteration that involves a weighted least squares update and the solution of a nonlinear equation with respect to a scalar parameter, and which converges rapidly, as illustrated by the computational times and number of iterations reported in the paper. The reported times and number of iterations were computed using a workstation with 24 cores at 2.90GHz and 80GB memory running under the CentOS 7 operating system, using one core per data set.
Using various settings we were able to retrieve enough information on the performance of the maximum median BRPL estimators. All our simulation studies illustrate that use of the median BRPL succeeds in achieving median centering in estimation, and results in confidence intervals with good coverage properties. Furthermore, while tests based on the LR suffer from size distortions, the median BRPL ratio statistic results in tests with size and power properties, sometimes better to those of the mean BRPL ratio statistic in Kosmidis et al. (2017).
The main advantage of the maximum median BRPL estimators from the maximum mean BRPL ones is their equivariance under monotone component-wise parameter transformations, which, in the case of random-effects meta-regression, leads to median bias-reduced standard errors.
As random-effect models are widely used in practice, the median BRPL method is likely to be useful in models with more complex random-effect structures, such as linear mixed models.
R code for replication of examples is provided as online supplementary material.
The authors thank E.C. Kenne Pagui for providing the formula for the median bias-reducing adjustment in matrix form.
The first author was partially supported by the Erasmus+ programme of the European Union which funded her for a 2-month traineeship at the University of Padova. The work of the second author was supported by the Alan Turing Institute under the EPSRC grant EP/N510129/1 (Turing award number TU/B/000082) and part of it was completed when he was a Senior Lecturer at University College London. The third author was supported by the Italian Ministry of Education under the PRIN 2015 grant 2015EASZFS_003, and by the University of Padova (PRAT 2015).
The observed information matrix for the random-effects meta-regression model (1) is
For this model
where is the zero matrix and is the th column of . The median bias-reducing adjustment for is obtained by plugging the above expressions into (4). The sum () is also given in the Appendix of Kosmidis et al. (2017).
- Bellio R and Guolo A. Integrated likelihood inference in small sample meta-analysis for continuous outcomes. Scand. J. Stat. 2016; 43: 191–201.
- Brockwell SE and Gordon IR. A comparison of statistical methods for meta-analysis. Stat. Med. 2001; 20: 825–840.
- DerSimonian R and Laird N. Meta-analysis in clinical trials. Control. Clin. Trials 1986; 7: 177–188.
- Firth D. Bias reduction of maximum likelihood estimates. Biometrika 1993; 80: 27–38.
- Guolo A and Varin C. Random-effects meta-analysis: The number of studies matters. Stat. Methods Med. Res. 2017; 26: 1500–1518.
- Hardy RJ and Thompson SG. A likelihood approach to meta-analysis with random effects. Stat. Med. 1996; 15: 619–629.
- Huizenga HM, Visser i and Dolan CV. Testing overall and moderator effects in random effects meta-regression. Br. J. Math. Stat. Psychol. 2011; 64: 1–19.
- Jackson D and Bowden J. A re-evaluation of the “quantile approximation method” for random effects meta-analysis. Stat. Med. 2009; 28: 338–348.
- Kenne Pagui EC, Salvan A and Sartori N. Median bias reduction of maximum likelihood estimates. Biometrika 2017; 104: 923–938.
- Knapp G and Hartung J. Improved tests for a random effects meta-regression with a single covariate. Stat. Med. 2003; 22: 2693–2710.
- Kosmidis I and Firth D. Bias reduction in exponential family nonlinear models. Biometrika 2009; 96: 793–804.
- Kosmidis I, Guolo A and Varin C. Improving the accuracy of likelihood-based inference in meta-analysis and meta-regression. Biometrika 2017; 104: 489–496.
- Larsson SC and Orsini N. Red meat and processed meat consumption and all-cause mortality: A meta-analysis. Am. J. Epidemiol. 2014; 179: 282–289.
- Pace L and Salvan A. Principles of statistical inference: From a neo-Fisherian perspective. 1997; London: World Scientific.
- Taubert D, Roesen R and Schömig E. Effect of cocoa and tea intake on blood pressure: a meta-analysis. Arch. Intern. Med. 2007; 167: 626–634.
- Van Houwelingen HC, Arends LR and Stijnen T. Advanced methods in meta-analysis: multivariate approach and meta-regression. Stat. Med. 2002; 21: 589–624.
- Zeng D and Lin DY. On random-effects meta-analysis. Biometrika 2015; 102: 281–294.