Median bias reduction in random-effects meta-analysis and meta-regression
Sophia Kyriakou Ioannis Kosmidis Nicola Sartori
stat.AP stat.ME 62F03, 62F12, 62P10








Median bias reduction in random-effects meta-analysis and meta-regression
    
S. Kyriakou 

Department of Statistical Science 

University College London 

London, WC1E 6BT, UK 




I. Kosmidis 

Department of Statistics 

University of Warwick 

Coventry, CV4 7AL, UK 

The Alan Turing Institute 

96 Euston Road, London NW1 2DB, UK 




N. Sartori 

Department of Statistical Sciences 

University of Padova 

Via Cesare Battisti, 35121 Padova, Italy 



    December 30, 2023
===============================================================================================================================================================================================================================================================================================================================================================================================




 The reduction of the mean or median bias of the maximum likelihood estimator in regular parametric models can be
achieved through the additive adjustment of the score equations. In this paper, we derive the adjusted score equations
for median bias reduction in random-effects meta-analysis and meta-regression models and derive efficient estimation
algorithms. The median bias-reducing adjusted score functions are found to be the derivatives of a penalised likelihood.
The penalised likelihood is used to form a penalised likelihood ratio statistic which has known limiting distribution and
can be used for carrying out hypothesis tests or for constructing confidence intervals for either the fixed-effect
parameters or the variance component. Simulation studies and real data applications are used to assess the
performance of estimation and inference based on the median bias-reducing penalised likelihood and compare it to
recently proposed alternatives. The results provide evidence on the effectiveness of median bias reduction in improving
estimation and likelihood-based inference.


Keywords:  Adjusted score equations; Heterogeneity; Mean bias reduction; Penalized likelihood; Random effects




§ INTRODUCTION


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, <cit.> 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 <cit.>.
Examples of recently proposed methods that attempt to improve
inference are the resampling <cit.> and double resampling
approaches<cit.>, and the mean bias-reducing penalized
likelihood (BRPL) approach in
<cit.>. Specifically,
<cit.> 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.

<cit.> show that under suitable conditions
third-order median unbiased estimators can be obtained by the solution
of a suitably adjusted score equation. The components of such median bias-reduced
estimators have, to third-order, the same probability of over- and under-estimating the true parameter. 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 <cit.>. 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
<cit.> 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 O(1),
the same arguments as in <cit.>
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.





§ COCOA INTAKE AND BLOOD PRESSURE REDUCTION DATA



Consider the setting in <cit.> who carry out a meta-analysis of five randomized controlled trials from <cit.> on the efficacy of two weeks of cocoa consumption on lowering diastolic blood
pressure. Figure <ref> is a forest plot with the
estimated mean difference in diastolic blood pressure before and after cocoa intake from each study, and the associated 95% 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 Y_i be the random variable representing the mean difference in the diastolic blood pressure after two weeks of cocoa intake in the ith study. We assume that Y_1, …, Y_5 are independent random variables where Y_i has a Normal distribution with mean the overall effect β and variance σ̂^2_i + ψ, with σ̂^2_i the estimated standard error of the effect from the ith study and ψ the heterogeneity parameter.




The forest plot in Figure <ref> depicts nominally 95%
confidence intervals for β using various alternative methods. As
is apparent, the conclusions when testing the hypothesis β = 0
can vary depending on the method used. More specifically, the Wald test using the ML estimates,
the DerSimonian and Laird method <cit.>, double
resampling <cit.>, and the likelihood ratio (LR) test
give evidence that there is a relationship between cocoa consumption
and diastolic blood pressure, with p-values 0.005, 0.006,
0.016, 0.030, respectively. On the other hand, Knapp and Hartung's
method <cit.>, the mean BRPL
ratio <cit.>, the Bartlett-corrected
LR <cit.>, and Skovgaard's test suggest that the
evidence that cocoa consumption affects diastolic blood pressure is
weaker, with p-values 0.050, 0.053, 0.058, 0.067,
respectively.



§ RANDOM-EFFECTS META-REGRESSION MODEL


Let y_i and σ̂_i^2 denote the estimate of the effect
from the ith study (i = 1, …, K) and the associated
within-study variance, respectively, and
x_i = (x_i1,…,x_ip)^ denote a p-vector of
study-specific covariates that can be used to account for the
heterogeneity across studies.

The within-study variances σ̂_i^2 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 y_1,…,y_K
are assumed to be realizations of the random variables
Y_1,…,Y_K, which are independent conditionally on independent
random effects U_1, …, U_K. The conditional distribution of
Y_i given U_i = u_i is
N(u_i + x^_iβ, σ̂_i^2), where β is an
unknown p-dimensional vector of fixed effects. The random effects
U_1,…,U_K are typically assumed to be independent with U_i
having a N(0, ψ) distribution, where ψ is a parameter that
attempts to capture the unexplained between-study heterogeneity. In
matrix notation, the random-effects meta-regression model has

    Y = X β + U + ϵ ,

where Y=(Y_1,…,Y_K)^, X is the K× p model matrix
with x_i^ in its ith row, and
ϵ=(ϵ_1,…,ϵ_K)^ is a vector of
independent errors each with a N(0,σ̂_i^2) distribution
and independent of U=(U_1,…,U_K)^. Under this
specification, the marginal distribution of Y is multivariate normal
with mean Xβ and variance-covariance matrix
Σ̂ + ψ I_K, where I_K is the K × K
identity matrix and
Σ̂ =
diag(σ̂_1^2,…,σ̂_K^2). The
random-effects meta-analysis results as a special case of
meta-regression, by setting X to be a column of ones.

The log-likelihood function for θ=(β^,ψ)^ is l(θ) = {log |W(ψ)| - R(β)^W(ψ)R(β)}/2, where |W(ψ)| denotes the determinant of W(ψ) = (Σ̂ + ψ I_K)^-1 and R(β)=y-Xβ. The gradient of the log-likelihood (score function) is

    s(θ) = [                    X^W(ψ)R(β); 1/2{R(β)^W(ψ)^2R(β) - [W(ψ)]} ]

and the ML estimator θ̂=(β̂^,ψ̂)^
is obtained as the solution of s(θ) = 0_p + 1, where 0_p
denotes a p-dimensional vector of zeros.



§ MEDIAN BIAS REDUCTION




 §.§ 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 <cit.>. <cit.> 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. <cit.> 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 j(θ)= -∂^2 l(θ)/∂θ∂θ^ be the observed information matrix (see, Appendix for its
expression), and i(θ) be the expected information matrix

    i(θ) =E_θ(j(θ)) =  [     X^W(ψ)X         0_p;        0^_p 1/2[W(ψ)^2] ] ,

with tth column i_t(θ). Let also i^t(θ) and
i^tt(θ) be the tth column and the tth diagonal element of
{i(θ)}^-1, with t ∈{1,…,p+1}. <cit.> show that a median
bias-reduced estimator θ̂^† can be obtained by
solving an adjusted score equation of the form
s^†(θ) = s(θ) + A^†(θ) = 0, where
the additive adjustment to A^†(θ) is O(1), in the
sense that A^†(θ) is bounded in absolute value by a
fixed constant after a sufficiently large value of K. The median
bias-reducing adjustment A^†(θ) has tth element

    A^†_t(θ) = 1/2[{i(θ)}^-1(P_t(θ) + Q_t(θ))] - {i_t(θ)}^K^†(θ) .

The quantities
P_t(θ) = E_θ[s(θ)s^(θ)s_t(θ)] and
Q_t(θ)=E_θ[-j(θ)s_t(θ)]
in (<ref>) are those introduced by <cit.> for mean bias-reduction, and
K^†(θ) is a (p+1)-vector with tth element K^†_t(θ) = {i^t(θ)}^K_t(θ), where K_t(θ) is another (p+1)-vector with uth element

    K_tu(θ) = [i^t(θ){i^t(θ)}^/i^tt(θ)(1/3P_u(θ)
          + 1/2Q_u(θ))]   .


In the context of meta-regression values of t and u in {1,…,p} correspond to the elements of parameter β, and t,u = p+1 correspond to parameter ψ. Given that A^†(θ) is of order O(1), θ̂^† has the same asymptotic distribution as θ̂ <cit.>, i.e. multivariate normal with mean θ and variance-covariance matrix {i(θ)}^-1, which can be consistently estimated with {i(θ̂^†)}^-1.

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

    A^†(θ)  = [                                  0_p; 1/2[W(ψ)H(ψ)] + 1/3[W(ψ)^3]/[W(ψ)^2] ] ,

where H(ψ) = X(X^W(ψ)X)^-1X^W(ψ). Substituting
(<ref>) in the expression for s^†(θ) gives
that the median bias-reducing adjusted score functions for β and
ψ are s^†_β(θ) = s_β(θ) and

    s^†_ψ(θ) = s_ψ(θ) + 1/2[W(ψ)H(ψ)] + 1/3[W(ψ)^3]/[W(ψ)^2] ,

respectively.



 §.§ Computation of median bias-reduced estimator



A direct approach for computing the estimator
θ̂^† =
(β̂^†,ψ̂^†)^ is through a
modification of the two-step iterative process in
<cit.>. At the jth iteration
(j = 1,2,…)


  * calculate β^(j) by weighted least squares as β^(j) = (X^W(ψ^(j-1))X)^-1X^W(ψ^(j-1))y

  *  solve s^†_ψ(θ^(j)(ψ)) = 0 with respect to ψ, where θ^(j)(ψ)=(β^(j),ψ)^.

In the above steps, β^(j) is the candidate value for
β̂^† at the jth iteration and ψ^(j-1) is
the candidate value for ψ̂^† at the (j-1)th
iteration. The equation in step 2 is solved numerically, by searching
for the root of the function s^†_ψ(β^(j),ψ) in
a predefined positive interval. For the computations in this
manuscript we use the <cit.> estimate of ψ
as starting value ψ^(0). The iterative process is then repeated
until the components of the score function s^†(θ) are
all less than ϵ = 1 × 10^-6 in absolute value at the
current estimates.



 §.§ Median bias-reducing penalized likelihood


Although it is not generally true that s^†(θ) is the
gradient of a suitable penalized log-likelihood, in this case
s^†(θ) is the gradient of the median BRPL

    l^†(θ) = l(θ) - 1/2log |X^W(ψ)X| - 1/6log [(W(ψ)^2)] .

Hence, θ̂^† is also the maximum median BRPL estimator. The median BRPL in
(<ref>) differs from the mean BRPL derived
in <cit.> by the term
-log [(W(ψ)^2)]/6.

An advantage of the median BRPL estimators over mean BRPL ones is that
the former are equivariant under monotone component-wise parameter
transformations <cit.>. 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 {i(θ)}^-1 in (<ref>) at
ψ^†. This is because i(θ) is a function of ψ only, and moreover the square
roots of the diagonal elements of {i(θ)}^-1 are monotone functions of ψ.



 §.§ 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
l^†(θ) for fixed τ, then the same arguments as in
<cit.> can be used to show that the logarithm of the
median BRPL ratio statistic

    2{l^†(τ̂^†,λ̂^†) - l^†(τ,λ̂^†_τ)} .

has a χ^2_(τ) asymptotic distribution, as K goes to infinity. Specifically, the
adjustment to the score function is additive and of order O(1). 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 <cit.>.



§ 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 <ref> are ψ̂ = 4.199,
ψ̂^∗ = 5.546, and ψ̂^† = 6.897,
respectively. The estimates of the common effect are
β̂ = -2.799, β̂^∗ = -2.811, and
β̂^† = -2.818, with standard errors 1.000,
1.128, and 1.242, 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 1.1 × 10^-2, 1.8 × 10^-2, and
1.1 × 10^-2 seconds, respectively.

Figure <ref> shows the value of LR, mean BRPL and median BRPL ratio statistic in (<ref>) 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 <cit.>. The horizontal line in Figure <ref> is the 95% quantile of the limiting χ^2_1 distribution,
and its intersection with the values of the statistics results in the endpoints of the corresponding 95% 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 95%
confidence intervals for β are (-6.21,0.52), (-5.73,0.05), and (-5.26,-0.40) for the median BRPL ratio statistic, mean BRPL ratio statistic, and LR statistic, respectively. The corresponding 95% confidence intervals for ψ are (1.4,58.0), (1.0,38.5), and (1.1,23.5), 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 p-values of 0.053 and 0.077.




In order to further investigate the performance of the three
approaches to estimation and inference, we performed a simulation
study where we simulated 10 000 independent samples from the
random-effects meta-analysis model with parameter values set to the ML
estimates reported earlier, i.e. β_0 = -2.799 and
ψ_0 = 4.199. Figure <ref> shows boxplots of the
estimates of β and ψ calculated from each of the 10 000
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
ψ, P_ψ_0(ψ̂≤ψ_0),
P_ψ_0(ψ̂^∗≤ψ_0) and
P_ψ_0(ψ̂^†≤ψ_0) 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
(<ref>) is in reducing the median bias of
the maximum likelihood estimator of ψ.

The simulated samples were also used to calculate the empirical
p-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 <ref> shows that the empirical
p-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 95%
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 88%, 93%, and 96% for LR, mean
BRPL ratio, and median BRPL ratio respectively, and the corresponding
coverage probabilities for ψ are 88%, 94%, and 96%,
respectively.










§ SIMULATION STUDY


More extensive simulations under the random-effects meta-analysis
model (<ref>) are performed here using the design in <cit.>. Specifically, the
data y_i, i ∈{1,…,K}, are simulated from model
(<ref>) with true fixed-effect parameter β=0.5. The
within-study variances σ̂^2_i are independently generated
from a χ^2_1 distribution and are multiplied by 0.25 before
restricted to the interval (0.009, 0.6).  Eleven values of the
between-study variance ψ ranging from 0 to 0.1 are chosen,
and the number of studies K ranges from 5 to 200. For each
combination of ψ and K considered, we simulated 10 000 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.

<cit.> compared the performance of their proposed
double resampling method with the <cit.> method,
the profile likelihood method in <cit.>, and the
resampling method in <cit.> and showed that the double
resampling method improves the accuracy of statistical
inference. Based on these results <cit.>
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 <cit.>
and <cit.> 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 <ref>. 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 l(θ),
l^∗(θ), and l^†(θ) with respect to ψ,
in the interval (0,3).

First, we compare the performance of the ML, maximum mean BRPL and
maximum median BRPL estimators in terms of percentage of
underestimation. Figure <ref> shows that the
median bias-reducing adjustment is the most effective in reducing
median bias even for small values of K. As expected, the ML and
maximum mean BRPL estimators also approach the limit of 50%
underestimation as K grows, with the latter being closer to 50%
than the former. Figure <ref> 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 K, while maximum
mean BRPL results in estimators with the smallest bias.

Figures <ref> and
<ref> 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 95% nominal level. Figure <ref> 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 95% nominal level. For small values of ψ or small and moderate number of studies K the empirical coverage of the intervals is larger than the nominal 95% 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 K increases.












Figures <ref> and
<ref> give the power of the LR,
the mean BRPL ratio, and the median BRPL ratio tests for testing the
null hypothesis β = 0.5 against various
alternatives. Specifically, we simulated 10 000 data sets under the
alternative hypothesis that parameter β is equal to
b = 0.5 + δ K^-1/2, where δ ranges from 0 to
2.25. In Figure <ref> the power
is calculated using critical values of the the asymptotic null
χ^2_1 distribution of the statistics. In
Figure <ref> 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
5%.

Figure <ref> shows that the three
tests have monotone power and for small values of K 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 K increases the three tests approach
the nominal size and provide similar power. The use of the exact
critical values in Figure <ref>
allows us to compare the performance of the tests without letting the
oversizing or the conservativeness of a test skew the power
results. Figure <ref> 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 K 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 K 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.

  

§ MEAT CONSUMPTION DATA


<cit.> investigate the association between meat
consumption and relative risk of all-cause mortality. The data
consists of 16 prospective studies, eight of which are about
unprocessed red meat consumption and eight about processed meat
consumption. Figure <ref> 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
Y_i has a N(β_0 + β_1x_i,σ̂^2_i + ψ), where
Y_i is the random variable representing the logarithm of the
relative risk reported in the ith study, and x_i takes value 1
if the consumption in the ith study is about processed red meat and
0 if it is about unprocessed meat (i = 1,…,16).
Table <ref> gives the ML estimates, the mean BRPL
estimates, and the median BRPL estimates of β_0, β_1 and
ψ, along with the corresponding standard errors for β_0 and
β_1. 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
1.2 × 10^-2, 2.4 × 10^-2, and 1.5 × 10^-2
seconds, respectively.

The LR test indicates some evidence for a higher risk associated to
the consumption of red processed meat with a p-value of 0.047. On
the other hand, the mean and median BRPL ratio tests suggest that
there is weaker evidence for higher risk, with p-values of 0.066
and 0.074, respectively. Skovgaard's test also gives weak evidence
for higher risk with p-value 0.073.



Similar to Section <ref>, we performed a simulation
study in order to further investigate the performance of the three
methods in a meta-regression context. We simulated 10 000
independent samples from the meta-regression model at the ML estimates
reported in Table <ref>. Figure <ref>
shows boxplots of the estimates of β_0, β_1, 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 72.6%, 56.6%, and
49.9% for the ML, maximum mean BRPL, and maximum median BRPL
estimators, respectively.



The simulated samples were also used to calculate the empirical
p-value distribution for the tests based on the likelihood, mean
BRPL and median BRPL ratio statistics. Table <ref>
shows that the empirical p-value distribution for the median BRPL
ratio statistic is the one closest to uniformity.





§ 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 <cit.> 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
<cit.>.

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.



§ SUPPLEMENTARY MATERIAL

The  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 bias-reducing adjustment in matrix form.



§ FUNDING


  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).


99

[Bellio & Guolo(2016)]bellio2016integrated
Bellio R and Guolo A.
Integrated likelihood inference in small sample meta-analysis for continuous outcomes.
Scand. J. Stat. 2016; 43: 191–201.

[Brockwell & Gordon(2001)]brockwell2001comparison
Brockwell SE and Gordon IR.
A comparison of statistical methods for meta-analysis.
Stat. Med. 2001; 20: 825–840.

[DerSimonian & Laird(1986)]dersimonian1986meta
DerSimonian R and Laird N.
Meta-analysis in clinical trials.
Control. Clin. Trials 1986; 7: 177–188.

[Firth(1993)]firth1993bias
Firth D.
Bias reduction of maximum likelihood estimates.
Biometrika 1993; 80: 27–38.

[Guolo & Varin(2015)]guolo2015random
Guolo A and Varin C.
Random-effects meta-analysis: The number of studies matters.
Stat. Methods Med. Res. 2017; 26: 1500–1518.

[Hardy & Thompson(1996)]hardy1996likelihood
Hardy RJ and Thompson SG.
A likelihood approach to meta-analysis with random effects.
Stat. Med. 1996; 15: 619–629.

[Huizenga et al.(2011)Huizenga et al.]huizenga2011testing
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 & Bowden(2009)]jackson2009re
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 et al.(2016)Kenne Pagui et al.]pagui2016median
Kenne Pagui EC, Salvan A and Sartori N.
Median bias reduction of maximum likelihood estimates.
Biometrika 2017; 104: 923–938.

[Knapp & Hartung(2003)]knapp2003improved
Knapp G and Hartung J.
Improved tests for a random effects meta-regression with a single covariate.
Stat. Med. 2003; 22: 2693–2710.

[Kosmidis & Firth(2009)]kosmidis2009bias
Kosmidis I and Firth D.
Bias reduction in exponential family nonlinear models.
Biometrika 2009; 96: 793–804.

[Kosmidis et al.(2017)Kosmidis et al.]kosmidis2017improving
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 & Orsini(2014)]larsson2014red
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 & Salvan(1997)]pace1997principles
Pace L and Salvan A.
Principles of statistical inference: From a neo-Fisherian perspective. 1997; London: World Scientific.

[Taubert et al.(2007)Taubert, Roesen &
  Schömig]taubert2007effect
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 et al.(2002)Van Houwelingen, Arends & Stijnen]van2002advanced
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 & Lin(2015)]zeng2015random
Zeng D and Lin DY.
On random-effects meta-analysis.
Biometrika 2015; 102: 281–294.





§ APPENDIX


The observed information matrix for the random-effects meta-regression model (<ref>) is

    j(θ) = [                       X^W(ψ)X                 X^W(ψ)^2 R(β);                 X^W(ψ)^2 R(β) R(β)^W(ψ)^3R(β) - 1/2[W(ψ)^2] ].

For this model

    P_t(θ) = -Q_t(θ) =
    [     0_p × p X^W(ψ)^2X_t; X^W(ψ)^2X_t           0 ]   (t=1,…,p) ,

and

    P_p+1(θ) =
    [ X^W(ψ)^2X       0_p;      0_p^  (W(ψ)^3) ]    and     Q_p+1(θ)=
    [   0_p × p       0_p;      0_p^ -(W(ψ)^3) ] ,

where 0_p × p is the p × p zero matrix and X_t is the tth column of X. The median bias-reducing adjustment for θ is obtained by plugging the above expressions into (<ref>). The sum P_t(θ) + Q_t(θ) (t = 1,…,p+1) is also given in the Appendix of <cit.>.



[pages=1-8]Suppl_KKS.pdf


