Median bias reduction in randomeffects metaanalysis and metaregression
Abstract
Randomeffects models are frequently used to synthesise information from different studies in metaanalysis. While likelihoodbased inference is attractive both in terms of limiting properties and of implementation, its application in randomeffects metaanalysis 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 randomeffects metaregression, which allows the mean effect size to vary with studyspecific covariates.
Keywords: Adjusted score equations; Heterogeneity; Mean bias reduction; Penalized likelihood; Random effects
= 1 \pdfinclusioncopyfonts=1
1 Introduction
Metaanalysis 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 randomeffects metaanalysis 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 randomeffects metaanalysis can be problematic in the usual metaanalytic 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 biasreducing 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 BRPLbased inference outperforms its competitors in terms of inferential performance.
Kenne Pagui et al. (2016) show that under suitable conditions thirdorder median unbiased estimators can be obtained by the solution of a suitably adjusted score equation. Such median biasreduced estimators have componentwise and to thirdorder the same probability of over and underestimating the true parameter value. A key property of these estimators, not shared with the mean biasreduced ones, is that any monotone componentwise transformation of the estimators results automatically in median biasreduced estimators of the transformed parameters (Kenne Pagui et al., 2016). Such equivariance property can be useful in the context of randomeffects metaanalysis where the Fisher information and, hence, the asymptotic variances of various likelihoodbased estimators depend only on the heterogeneity parameter.
In this paper, we derive the median biasreducing adjusted score functions for randomeffects metaanalysis and metaregression. 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 fixedeffect 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 likelihoodbased inference.
2 Cocoa intake and blood pressure reduction data
Consider the setting in Bellio & Guolo (2016) who carry out a metaanalysis 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 Waldtype confidence intervals. Four out of the five studies reported a reduction of diastolic blood pressure from cocoa intake.
The randomeffects metaanalysis 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 Bartlettcorrected 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 Randomeffects metaregression model
Let and denote the estimate of the effect from the th study and the associated withinstudy variance, respectively, and denote a vector of studyspecific covariates that can be used to account for the heterogeneity across studies.
The withinstudy variances are usually assumed to be estimated wellenough 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 betweenstudy heterogeneity. In matrix notation, the randomeffects metaregression model has
(1) 
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 variancecovariance matrix , where is the identity matrix and . The randomeffects metaanalysis results as a special case of metaregression, by setting to be a column of ones.
The loglikelihood function for is , where denotes the determinant of and . The gradient of the loglikelihood (score function) is
(2) 
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 biasreduced 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 CornishFisher asymptotic expansion.
Let be the observed information matrix (see, Appendix for its expression), and be the expected information matrix
(3) 
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 biasreduced 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 biasreducing adjustment has th element
(4) 
The quantities and in (4) are those introduced by Kosmidis & Firth (2009) for mean biasreduction, and is a vector with th element , where is another vector with th element
In the context of metaregression 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 variancecovariance matrix , which can be consistently estimated with .
After some algebra (see Appendix for details) the median biasreducing adjustment for the randomeffects metaanalysis and metaregression models has the form
(5) 
where . Substituting (5) in the expression for gives that the median biasreducing adjusted score functions for and are and
respectively.
4.2 Computation of median biasreduced estimator
A direct approach for computing the estimator is through a modification of the twostep 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 biasreducing penalized likelihood
Although it is not generally true that is the gradient of a suitable penalized loglikelihood, in this case is the gradient of the median BRPL
(6) 
Hence, is also the maximum median BRPL estimator. The median BRPL in (6) differs from the mean BRPL derived in Kosmidis et al. (2017) by the term .
An advantage of the median BRPL estimators over mean BRPL ones is that the former are equivariant under monotone componentwise parameter transformations (Kenne Pagui et al., 2016). In the context of randomeffects metaanalysis and metaregression, this equivariance implies that not only we get a median biasreduced estimator of , but we also get median biasreduced 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 likelihoodbased inference
For inference about either the components of the fixedeffect parameters or the betweenstudy 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
(7) 
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 loglikelihood (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 metaanalysis model in Section 2 are , , and , respectively. The estimates of the common effect are , and , with standard errors , , and , respectively. The biasreduced estimates of and, as a consequence, the corresponding estimated standard errors for are larger than their ML counterparts, which is typical in randomeffects metaanalysis. 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 runtime for the twostep 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 randomeffects metaanalysis 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 overcorrect 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 simulationbased 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 twosided 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.
LR  5.7  8.4  11.8  18.2  34.6  57.8  79.2  91.8  96.1  98.0  99.3 

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 randomeffects metaanalysis model (1) are performed here using the design in Brockwell & Gordon (2001). Specifically, the data , , are simulated from model (1) with true fixedeffect parameter . The withinstudy variances are independently generated from a distribution and are multiplied by before restricted to the interval . Eleven values of the betweenstudy 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 withinstudy 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 randomeffect parameters obtained from the three methods are calculated using variants of the twostep 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 betweenstudy 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 biasreducing 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 onesided and twosided 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 twosided 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 twostep 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 runtimes 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 runtime taken per fit for the twostep 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 allcause 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 metaanalysis. The results from the studies point towards the conclusion that high consumption of red meat, in particular processed red meat, is associated with higher allcause mortality.
We consider the randomeffects metaregression 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 fixedeffect 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.
Method  

ML  0.099 (0.044)  0.106 (0.061)  0.009 
[0.004,0.189]  [0.022,0.244]  [0.003,0.030]  
Mean BRPL  0.095 (0.050)  0.110 (0.069)  0.012 
[0.020,0.199]  [0.040,0.264]  [0.003,0.042]  
Median BRPL  0.093 (0.052)  0.111 (0.072)  0.013 
[0.027,0.203]  [0.048,0.271]  [0.004,0.048] 
Similar to Section 5, we performed a simulation study in order to further investigate the performance of the three methods in a metaregression context. We simulated independent samples from the metaregression 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.
LR  2.2  4.5  7.7  13.1  28.0  50.0  71.7  86.6  92.1  95.3  97.7 

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 randomeffects metaanalysis and metaregression models and describe the associated inferential procedures.
We show that the solution of the median biasreducing adjusted score equations is equivalent to maximizing a penalized loglikelihood. 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 twostep 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 componentwise parameter transformations, which, in the case of randomeffects metaregression, leads to median biasreduced standard errors.
As randomeffect models are widely used in practice, the median BRPL method is likely to be useful in models with more complex randomeffect structures, such as linear mixed models.
Supplementary material
The R
code for replication of examples is provided as online supplementary material.
Acknowledgements
The authors thank E.C. Kenne Pagui for providing the formula for the median biasreducing adjustment in matrix form.
Funding
The first author was partially supported by the Erasmus+ programme of the European Union which funded her for a 2month 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).
Appendix
The observed information matrix for the randomeffects metaregression model (1) is
For this model
and
where is the zero matrix and is the th column of . The median biasreducing adjustment for is obtained by plugging the above expressions into (4). The sum () is also given in the Appendix of Kosmidis et al. (2017).
References
 Bellio R and Guolo A. Integrated likelihood inference in small sample metaanalysis for continuous outcomes. Scand. J. Stat. 2016; 43: 191–201.
 Brockwell SE and Gordon IR. A comparison of statistical methods for metaanalysis. Stat. Med. 2001; 20: 825–840.
 DerSimonian R and Laird N. Metaanalysis 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. Randomeffects metaanalysis: The number of studies matters. Stat. Methods Med. Res. 2017; 26: 1500–1518.
 Hardy RJ and Thompson SG. A likelihood approach to metaanalysis with random effects. Stat. Med. 1996; 15: 619–629.
 Huizenga HM, Visser i and Dolan CV. Testing overall and moderator effects in random effects metaregression. Br. J. Math. Stat. Psychol. 2011; 64: 1–19.
 Jackson D and Bowden J. A reevaluation of the “quantile approximation method” for random effects metaanalysis. 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 metaregression 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 likelihoodbased inference in metaanalysis and metaregression. Biometrika 2017; 104: 489–496.
 Larsson SC and Orsini N. Red meat and processed meat consumption and allcause mortality: A metaanalysis. Am. J. Epidemiol. 2014; 179: 282–289.
 Pace L and Salvan A. Principles of statistical inference: From a neoFisherian perspective. 1997; London: World Scientific.
 Taubert D, Roesen R and Schömig E. Effect of cocoa and tea intake on blood pressure: a metaanalysis. Arch. Intern. Med. 2007; 167: 626–634.
 Van Houwelingen HC, Arends LR and Stijnen T. Advanced methods in metaanalysis: multivariate approach and metaregression. Stat. Med. 2002; 21: 589–624.
 Zeng D and Lin DY. On randomeffects metaanalysis. Biometrika 2015; 102: 281–294.