Highlights
What is already known?
-
Random-effects meta-analyses with few studies often struggle to accurately estimate between-study heterogeneity, resulting in biased effect estimates and poorly covered confidence intervals. For normally distributed outcomes, the Hartung–Knapp–Sidik–Jonkman (HKSJ) method adjusts the CI for the global mean. Rover et al. proposed a modified Knapp–Hartung (mKH) method to improve performance with a small number of studies.
What is new?
-
We proposed two new approaches to construct confidence intervals for the global mean: one based on fiducial inference and another using modified signed log-likelihood ratio test statistics for improved small-sample performance. Simulations showed that the new methods yielded coverage probabilities closer to nominal levels and shorter intervals than the existing methods.
Potential impact for RSM readers
-
Random-effects meta-analyses with a small number of studies are a common challenge, particularly in areas such as rare disease research. The proposed approaches offer readers effective alternatives with improved performance for meta-analyses with few studies, helping researchers make more informed decisions.
1 Introduction
Conducting a meta-analysis often involves a random-effects model that can account for both within-study and between-study variability by using a normal–normal hierarchical model.Reference Rover, Knapp and Friede1 For outcomes that follow a normal distribution, the data are represented as k estimates
$Y_i \in \mathbb {R}$
with known within-trial variances
$s_i^2>0$
:
$$ \begin{align} Y_{i} \sim N\left(\theta_i,s_{i}^{2}\right), \text { for } i=1, \ldots, k, \end{align} $$
where the
$\theta _i$
s vary from trial to trial around a global mean
$\mu $
according to the normal distribution
$$ \begin{align} \theta_{i} \sim N\left(\mu,\tau^{2}\right), \text { for } i=1, \ldots, k. \end{align} $$
Here,
$\tau ^2$
is a heterogeneity variance component, representing the variability between trials or studies. The integrated and simplified model under consideration, as discussed in,Reference Rover, Knapp and Friede1 is
$$ \begin{align} Y_{i} \sim N\left(\mu, \tau^{2}+s_{i}^{2}\right), \text { for } i=1, \ldots, k, \end{align} $$
where
$\mu $
is the unknown parameter of interest. The within-trial variances
$s_i^2$
, assumed to be known, depend on the group sizes in each trial and on the within-trial variances of the outcome variables in the ith trial. The latter are typically estimated, even though we proceed under the standard assumption that these are known. In many applications, only a limited number of studies is available for a meta-analysis, often fewer than five. Estimating the between-study heterogeneity
$\tau ^2$
is particularly challenging in these situations. Inaccurate heterogeneity estimates can result in biased effect estimates and overly narrow confidence intervals.Reference Mathes and Kuss2 Deriving reliable estimates and confidence intervals from only a small number of studies is challenging, especially in the context of rare diseases.
The problem of interest in our work is the computation of confidence limits for the underlying mean
$\mu $
when the number of studies k is small. Mathes and KussReference Mathes and Kuss2 demonstrated that for binary outcomes and in the presence of random effects, the beta-binomial model was the most appropriate model for meta-analysis involving a small number of studies, offering a balanced trade-off between coverage probability and power. For normally distributed outcomes, there are several methods to construct confidence intervals for
$\mu $
. However, there is no universally accepted satisfactory method, especially for small k. In the present work, the focus is on the outcomes following a normal distribution.
A natural point estimate of
$\mu $
can be obtained as a weighted average of the
$Y_i$
s in (3), weighted by the inverse of the variance in the model (3). The resulting point estimate, say
$\hat {\mu }_{\tau }$
, and the associated normal distribution, are given by
$$ \begin{align} \hat{\mu}_{\tau}=\frac{\sum_{i=1}^{k}\left(\tau^{2}+s_{i}^{2}\right)^{-1} Y_{i}}{\sum_{i=1}^{k}\left(\tau^{2}+s_{i}^{2}\right)^{-1}} \sim N\left(\mu,\hat{\sigma}_{\mu}^2 = \left \{\sum_{i=1}^{k}\left(\tau^{2}+s_{i}^{2}\right)^{-1}\right\}^{-1}\right). \end{align} $$
Using a normal approximation for
$\hat {\mu }_{\tau }$
, an approximate
$(1-\alpha )\%$
confidence interval for
$\mu $
is given by
$$ \begin{align} \hat{\mu}_{\tau} \pm \hat{\sigma}_{\mu} z_{(1-\alpha/2)}, \end{align} $$
where
$ z_{(1-\alpha /2)} $
is the
$(1-\alpha /2)$
-quantile of the standard normal distribution. In practice,
$\tau ^2$
is replaced with a suitable estimate, and the normal approximation is satisfactory when the number of studies k is large, with small within-trial variances or when heterogeneity is minimal.Reference Rover, Knapp and Friede1
An adjusted CI introduced by Hartung and Knapp,Reference Hartung and Knapp3 and Sidik and Jonkman,Reference Sidik and Jonkman4 usually referred to as the Hartung–Knapp–Sidik–Jonkman (HKSJ) method is given by
$$ \begin{align} \hat{\mu}_{\tau} \pm \sqrt{q} \hat{\sigma}_{\mu} t_{(k-1);(1-\alpha/2)}, \end{align} $$
where
$q = \frac {1}{k-1} \sum _i \frac {(y_i - \hat {\mu }_{\tau })^2}{\left (\tau ^{2}+s_{i}^{2}\right )}$
, and
$t_{(k-1);(1-\alpha /2)}$
is the
$(1-\alpha /2)$
-quantile of the Student-t distribution with
$(k-1)$
degrees of freedom. For a small number of studies k, Rover et al.Reference Rover, Knapp and Friede1 proposed the modified Knapp–Hartung (mKH) method, replacing q by
$$ \begin{align} q^* = \max\{1, q\}. \end{align} $$
The CI of
$\mu $
using the HKSJ method in (6) is generally wider than the CI based on the normal approximation in (5) since the student-t quantile is larger than the corresponding normal quantile. However, q can become arbitrarily small, and if
$q < \frac {z_{(1-\alpha /2)}}{t_{(k-1);(1-\alpha /2)}}$
, the HKSJ method’s adjusted interval would be shorter than the normal approximation-based interval, which is counter-intuitive. The mKH method provides a more conservative approach, with error rates that are more closely aligned with the nominal level, particularly when dealing with a small number of studies with varying sizes or precision levels.Reference Rover, Knapp and Friede1
In this work, two novel approaches are proposed to construct confidence intervals for
$\mu $
, and their performance is evaluated numerically through a simulation study, and compared with the existing methods. The first approach is based on the idea of fiducial inference, and the second approach is based on two modifications to the signed log-likelihood ratio test statistic to obtain accurate small sample performance.Reference DiCiccio, Martin and Stern5 For details on fiducial inference, we refer to the works of Hannig et al.Reference Hannig6,
Reference Hannig, Iyer, Lai and Lee7 Earlier, fiducial quantities were investigated by Weerahandi,Reference Weerahandi8 referring to it as generalized pivotal quantities.
The following sections first describe the proposed methods and algorithms for obtaining confidence limits of
$\mu $
in (3). Next, the performance of these methods is evaluated through the results of a simulation study and two applied examples. The paper concludes with a discussion of the findings.
2 Methods
To obtain the confidence limits for
$\mu $
in the model (3), especially when the number of studies k is small, two methods are proposed. One method is based on the fiducial approach, while the other involves modifications of the signed log-likelihood ratio test statistic. The fiducial approach is primarily based on the work of Iyer et al.Reference Iyer, Wang and Mathew9 The second method, which focuses on small sample asymptotics, includes two modifications to the signed log-likelihood ratio test to improve performance in small sample scenarios, as detailed by DiCiccio et al.Reference DiCiccio, Martin and Stern5
2.1 The fiducial approach
Consider the point estimate of
$\mu $
, denoted by
$\hat {\mu }_{\tau }$
, given in equation (4). Let the corresponding residual sum of squares be denoted by
$R\left (\tau ^{2}\right )$
. Defining
$$ \begin{align} \boldsymbol{y} = \left(Y_{1}, Y_{2}, \ldots, Y_{k}\right)^{\prime}, \quad \text{and} \quad \boldsymbol{D} = \operatorname{diag}\left(s_{1}^{2}, s_{2}^{2}, \ldots, s_{k}^{2}\right) + \tau^{2} I_{k}, \end{align} $$
the quantity
$R\left (\tau ^{2}\right )$
can be expressed as
$$ \begin{align} \begin{aligned} R\left(\tau^{2}\right) &= \left(\boldsymbol{y} - \hat{\mu}_{\tau} \mathbf{1}_{k}\right)^{\prime} \boldsymbol{D}^{-1}\left(\boldsymbol{y} - \hat{\mu}_{\tau} \mathbf{1}_{k}\right) \\ &= \boldsymbol{y}^{\prime}\left[\boldsymbol{D}^{-1} - \frac{\boldsymbol{D}^{-1} \mathbf{1}_{k} \mathbf{1}_{k}^{\prime} \boldsymbol{D}^{-1}}{\mathbf{1}_{k}^{\prime} \boldsymbol{D}^{-1} \mathbf{1}_{k}}\right] \boldsymbol{y} \sim \chi_{k-1}^{2}, \end{aligned} \end{align} $$
where
$\mathbf {1}_{k}$
is a
$k \times 1$
vector of ones, and
$\chi _{r}^{2}$
denotes the central chi-square distribution with r degrees of freedom. The quantities
$\hat {\mu }_{\tau }$
and
$R\left (\tau ^{2}\right )$
are independent, as shown by the standard results on linear models under the normality assumption.
To understand the behavior of
$R\left (\tau ^{2}\right )$
as a function of
$\tau ^{2}$
, consider another representation of
$R\left (\tau ^{2}\right )$
. Let
$\mathbf {Z}$
be a
$k \times (k-1)$
matrix whose columns are orthonormal and orthogonal to the vector
$\mathbf {1}_{k}$
(i.e.,
$\mathbf {Z}^{\prime } \mathbf {Z} = I_{k-1}$
and
$\mathbf {Z}^{\prime } \mathbf {1}_{k}$
is a vector of zeros). Then, the following result holds:
$$\begin{align*}\boldsymbol{D}^{-1} - \frac{\boldsymbol{D}^{-1} \mathbf{1}_{k} \mathbf{1}_{k}^{\prime} \boldsymbol{D}^{-1}}{\mathbf{1}_{k}^{\prime} \boldsymbol{D}^{-1} \mathbf{1}_{k}} = \mathbf{Z}\left(\mathbf{Z}^{\prime} \boldsymbol{D} \mathbf{Z}\right)^{-1} \mathbf{Z}^{\prime}, \end{align*}$$
see Appendix M, Section M.4.f in the work by Searle et al.Reference Searle, Casella and McCulloch10 Thus,
$R\left (\tau ^{2}\right )$
can be represented as
$$ \begin{align} R\left(\tau^{2}\right) = \boldsymbol{y}^{\prime} \mathbf{Z}\left(\mathbf{Z}^{\prime} \boldsymbol{D} \mathbf{Z}\right)^{-1} \mathbf{Z}^{\prime} \boldsymbol{y}, \end{align} $$
where
$\boldsymbol {D}$
is the diagonal matrix given in Equation (8). From this representation, it follows that
$R\left (\tau ^{2}\right )$
is a decreasing function of
$\tau ^{2}$
, reaching its maximum value when
$\tau ^{2} = 0$
, and is also a convex function of
$\tau ^{2}$
. These properties were also proved by Iyer et al.Reference Iyer, Wang and Mathew9 based on the representation (9). Let
$$ \begin{align} \boldsymbol{D}_{0} = \operatorname{diag}\left(s_{1}^{2}, s_{2}^{2}, \ldots, s_{k}^{2}\right) \end{align} $$
which is the value of
$\boldsymbol {D}$
in (8) when
$\tau ^{2} = 0$
. Thus,
$R\left (\tau ^{2}\right ) \leq R(0)$
. Note that
$$ \begin{align} \begin{aligned} R(0) &= \left(\boldsymbol{y} - \hat{\mu}_{0} \mathbf{1}_{k}\right)^{\prime} \boldsymbol{D}_{0}^{-1}\left(\boldsymbol{y} - \hat{\mu}_{0} \mathbf{1}_{k}\right) \\ &= \boldsymbol{y}^{\prime}\left[\boldsymbol{D}_{0}^{-1} - \frac{\boldsymbol{D}_{0}^{-1} \mathbf{1}_{k} \mathbf{1}_{k}^{\prime} \boldsymbol{D}_{0}^{-1}}{\mathbf{1}_{k}^{\prime} \boldsymbol{D}_{0}^{-1} \mathbf{1}_{k}}\right] \boldsymbol{y} \\ &= \boldsymbol{y}^{\prime} \mathbf{Z}\left(\mathbf{Z}^{\prime} \boldsymbol{D}_{0} \mathbf{Z}\right)^{-1} \mathbf{Z}^{\prime} \boldsymbol{y}, \end{aligned} \end{align} $$
where
$\hat {\mu }_{0} = \frac {\sum _{j=1}^{k} \left (s_{j}^{2}\right )^{-1} Y_{j}}{\sum _{j=1}^{k}\left (s_{j}^{2}\right )^{-1}}$
. We note that
$R(0)$
is the usual
$Q-$
statistic to test if
$\tau ^2 = 0$
, and
$R(\tau ^2)$
is the generalized
$Q-$
statistic.Reference Viechtbauer11
We shall now explain the derivation of a fiducial quantity for
$\mu $
, whose percentiles can be used to compute confidence limits for
$\mu $
. Since the fiducial quantity for
$\mu $
depends on the fiducial quantity for
$\tau ^{2}$
, we shall first obtain a fiducial quantity for the latter. Let
$\boldsymbol {y}_{o}$
denote the observed value of
$\boldsymbol {y}$
. Furthermore, let
$$ \begin{align} \begin{aligned} R_{o}\left(\tau^{2}\right) &= \boldsymbol{y}_{o}^{\prime}\left[\boldsymbol{D}^{-1} - \frac{\boldsymbol{D}^{-1} \mathbf{1}_{k} \mathbf{1}_{k}^{\prime} \boldsymbol{D}^{-1}}{\mathbf{1}_{k}^{\prime} \boldsymbol{D}^{-1} \mathbf{1}_{k}}\right] \boldsymbol{y}_{o} \\ \text{and} \quad R_{o}(0) &= \boldsymbol{y}_{o}^{\prime}\left[\boldsymbol{D}_{0}^{-1} - \frac{\boldsymbol{D}_{0}^{-1} \mathbf{1}_{k} \mathbf{1}_{k}^{\prime} \boldsymbol{D}_{0}^{-1}}{\mathbf{1}_{k}^{\prime} \boldsymbol{D}_{0}^{-1} \mathbf{1}_{k}}\right] \boldsymbol{y}_{o}. \end{aligned} \end{align} $$
Let
$U^{2} \sim \chi _{k-1}^{2}$
. Following the work of Iyer et al.,Reference Iyer, Wang and Mathew9 a fiducial quantity for
$\tau ^{2}$
is the solution to the equation
$R_{o}\left (\tau ^{2}\right ) = U^{2}$
, provided
$U^{2}$
is between 0 and
$R_{o}(0)$
. Thus
$\tilde {\tau }^{2}$
, the fiducial quantity for
$\tau ^{2}$
, is given by
$$ \begin{align} \tilde{\tau}^{2} = \left\{ \begin{array}{l} \text{Solution to } R_{o}\left(\tau^{2}\right) = U^{2}, \text{ if } 0 < U^{2} < R_{o}(0), \\ 0, \text{ otherwise }. \end{array} \right. \end{align} $$
Here, we would like to note that Paule and MandelReference Paule and Mandel12 suggested an estimator of
$\tau ^2$
as the solution to
$R_o(\tau ^2) = k-1$
.
To get a fiducial quantity for
$\mu $
, note that
$\hat {\mu }_{\tau }$
in (4) can be written as
$$\begin{align*}\hat{\mu}_{\tau} = \left(\mathbf{1}_{k}^{\prime} \boldsymbol{D}^{-1} \mathbf{1}_{k}\right)^{-1} \mathbf{1}_{k}^{\prime} \boldsymbol{D}^{-1} \boldsymbol{y} \sim N\left(\mu, \left(\mathbf{1}_{k}^{\prime} \boldsymbol{D}^{-1} \mathbf{1}_{k}\right)^{-1}\right). \end{align*}$$
Thus, we can write
$$ \begin{align} \mu = \hat{\mu}_{\tau} - \left(\mathbf{1}_{k}^{\prime} \boldsymbol{D}^{-1} \mathbf{1}_{k}\right)^{-\frac{1}{2}} Z, \end{align} $$
where
$Z = \left (\mathbf {1}_{k}^{\prime } \boldsymbol {D}^{-1} \mathbf {1}_{k}\right )^{\frac {1}{2}}\left (\hat {\mu }_{\tau } - \mu \right ) \sim N(0,1)$
. The expression on the right-hand side of the above equation depends on the unknown parameter
$\tau ^{2}$
and the observable random variable
$\boldsymbol {y}$
(in the expression for
$\hat {\mu }_{\tau }$
). A fiducial quantity for
$\mu $
is obtained from the above equation by replacing
$\tau ^{2}$
with the fiducial quantity
$\tilde {\tau }^{2}$
in (14) and by replacing
$\boldsymbol {y}$
with its observed value
$\boldsymbol {y}_{o}$
. Thus, let
$$ \begin{align} \tilde{\boldsymbol{D}} = \operatorname{diag}\left(s_{1}^{2}, s_{2}^{2}, \ldots, s_{k}^{2}\right) + \tilde{\tau}^{2} I_{k}, \quad \text{and} \quad \hat{\mu}_{\tilde{\tau}} = \left(\mathbf{1}_{k}^{\prime} \tilde{\boldsymbol{D}}^{-1} \mathbf{1}_{k}\right)^{-1} \mathbf{1}_{k}^{\prime} \tilde{\boldsymbol{D}}^{-1} \boldsymbol{y}_{o}. \end{align} $$
A fiducial quantity for
$\mu $
, denoted by
$\tilde {\mu }$
, is then given by
$$ \begin{align} \tilde{\mu} = \hat{\mu}_{\tilde{\tau}} - \left(\mathbf{1}_{k}^{\prime} \tilde{\boldsymbol{D}}^{-1} \mathbf{1}_{k}\right)^{-\frac{1}{2}} Z. \end{align} $$
We note that
$\tilde {\mu }$
is a function of two independent random variables: the chi-square random variable
$U^{2}$
with
$\mathrm {df} = k-1$
used in the computation of
$\tilde {\tau }^{2}$
in (14), and a standard normal random variable Z. Once the observed data vector
$\boldsymbol {y}_{o}$
is available, a large number of values, say m, of the pair
$\left (U^{2}, Z\right )$
can be generated. These m pairs of values, along with the observed data
$\boldsymbol {y}_{o}$
, can be used to compute m realizations of the fiducial quantities
$\tilde {\tau }^{2}$
and
$\tilde {\mu }$
. The percentiles of the m realizations of
$\tilde {\mu }$
provide confidence limits for
$\mu $
.
In the above approach, point estimates of
$\mu $
and
$\tau ^{2}$
are not explicitly used. An algorithm for implementing the fiducial approach to obtain confidence limits for
$\mu $
is outlined below.
2.2 Small sample asymptotics
It is well known that under regularity conditions, the signed log-likelihood ratio test statistic is asymptotically standard normal, providing first-order accuracy. Small-sample asymptotics involve modifying the signed log-likelihood ratio test statistic to achieve accurate performance in small samples. In order to introduce the basic idea, consider an
$n \times 1$
vector of observations
$\boldsymbol {y}$
whose distribution depends on a
$d \times 1$
parameter vector
$\boldsymbol {\theta }$
. Let
$\ell (\boldsymbol {\theta }) = \ell (\boldsymbol {\theta }; \boldsymbol {y})$
be the log-likelihood function. Suppose that the parameter
$\boldsymbol {\theta }$
is partitioned as
$\left (\psi , \boldsymbol {\lambda }^{\prime }\right )^{\prime }$
, where
$\psi $
is a scalar parameter of interest, and
$\boldsymbol {\lambda }$
is a nuisance parameter. Let
$\hat {\boldsymbol {\theta }} = \left (\hat {\psi }, \hat {\boldsymbol {\lambda }}^{\prime }\right )^{\prime }$
denote the MLE,
$\hat {\boldsymbol {\lambda }}_{\psi }$
denote the constrained MLE of
$\boldsymbol {\lambda }$
for a fixed
$\psi $
, and write
$\hat {\boldsymbol {\theta }}_{\psi } = \left (\psi , \hat {\boldsymbol {\lambda }}_{\psi }^{\prime }\right )^{\prime }$
. Large sample inference for
$\psi $
can be based on the signed log-likelihood root
$$ \begin{align} r(\psi) = \operatorname{sign}(\hat{\psi} - \psi) \left[2\left\{\ell(\hat{\boldsymbol{\theta}}) - \ell(\hat{\boldsymbol{\theta}}_{\psi})\right\}\right]^{1 / 2}, \end{align} $$
where
$\operatorname {sign}(x)$
is +1 if
$x> 0$
and -1 if
$x < 0$
. Under regular conditions,
$r(\psi )$
is asymptotically standard normal with an error of the order
$O\!\left (n^{-1 / 2}\right )$
. Several modifications of
$r(\psi )$
are available in the literature to improve small sample performance.Reference Brazzale, Davison and Reid13 Here, we consider two modifications proposed by DiCiccio et al.Reference DiCiccio, Martin and Stern5
2.2.1 Modification I
To describe the first modification, let
$E(r(\psi )) = m(\psi , \boldsymbol {\lambda })$
and
$\operatorname {Var}(r(\psi )) = 1 + v(\psi , \boldsymbol {\lambda })$
. Now define
$$ \begin{align} r_{1}(\psi) = \frac{r(\psi) - m\left(\psi, \hat{\boldsymbol{\lambda}}_{\psi}\right)}{\left\{1 + v\left(\psi, \hat{\boldsymbol{\lambda}}_{\psi}\right)\right\}^{1 / 2}}. \end{align} $$
For a fixed value of
$\psi $
,
$r_{1}(\psi )$
is thus a standardization of
$r(\psi )$
using the mean and variance evaluated at
$\boldsymbol {\lambda } = \hat {\boldsymbol {\lambda }}_{\psi }$
. DiCiccio et al.Reference DiCiccio, Martin and Stern5 have shown that under regular conditions,
$r_{1}(\psi )$
has an asymptotic standard normal distribution, and
$P\left [r_{1}(\psi ) \geq t\right ] = [1 - \Phi (t)] + O\left (n^{-3 / 2}\right )$
. In other words, the tail area approximation based on the asymptotic standard normal distribution of
$r_{1}(\psi )$
is significantly more accurate than that based on the asymptotic standard normal distribution of
$r(\psi )$
.
In order to carry out inference based on the asymptotic standard normal distribution of
$r_{1}(\psi )$
, it is necessary to compute the mean and variance of
$r(\psi )$
; namely, the quantities
$m\left (\psi , \hat {\boldsymbol {\lambda }}_{\psi }\right )$
and
$1 + v\left (\psi , \hat {\boldsymbol {\lambda }}_{\psi }\right )$
. For a specified value of
$\psi $
, these can be obtained numerically, as noted in the work of DiCiccio et al.Reference DiCiccio, Martin and Stern5 After specifying
$\psi $
, compute
$\hat {\boldsymbol {\lambda }}_{\psi }$
, and generate a parametric bootstrap sample of size n when
$\psi $
takes the specified value and
$\boldsymbol {\lambda } = \hat {\boldsymbol {\lambda }}_{\psi }$
. Let
$r^{*}(\psi )$
denote the value of
$r(\psi )$
calculated from the parametric bootstrap sample. Generate several parametric bootstrap samples of size n, say B samples, and compute B values of
$r(\psi )$
, denoted by
$r_{j}^{*}(\psi )$
for
$j = 1, 2, \ldots , B$
. The mean and variance of
$r_{j}^{*}(\psi )$
for
$j = 1, 2, \ldots , B$
, provide estimates of
$m\left (\psi , \hat {\boldsymbol {\lambda }}_{\psi }\right )$
and
$1 + v\left (\psi , \hat {\boldsymbol {\lambda }}_{\psi }\right )$
, respectively, to be used in the calculation of
$r_{1}(\psi )$
. To compute confidence limits for
$\psi $
based on
$r_{1}(\psi )$
, equate
$r_{1}(\psi )$
to standard normal percentiles and solve for
$\psi $
. For example, a
$95\%$
upper confidence limit for
$\psi $
can be obtained as the solution for
$\psi $
to equation
$r_{1}(\psi ) = -1.96$
.
2.2.2 Modification II
A second modification suggested by DiCiccio et al.Reference DiCiccio, Martin and Stern5 involves estimating the tail probability
$P(r(\psi ) \leq t)$
using Monte Carlo simulation at the parameter value
$\boldsymbol {\theta } = \hat {\boldsymbol {\theta }}_{\psi }$
, rather than relying on the standard normal distribution for
$r(\psi )$
. For example, to test the null hypothesis
$\psi = \psi _{0}$
against the one-sided alternative
$\psi> \psi _{0}$
, a p-value can be calculated as
$P\left (r\left (\psi _{0}\right ) \leq r_{\text {obs}}\right )$
, evaluated at the parameter value
$\boldsymbol {\theta }_{0} = \hat {\boldsymbol {\theta }}_{\psi _{0}}$
, where
$r_{\text {obs}}$
is the observed value of
$r\left (\psi _{0}\right )$
. The confidence limits for
$\psi $
can then be obtained by inverting the test.
2.2.3 Algorithms
In the context of the random-effects meta-analysis model given in (3), the parameter of interest is
$\psi = \mu $
, and the nuisance parameter is the scalar
$\lambda = \tau ^{2}$
. Thus, in terms of the notation used for the two modifications given above, the parameter vector is
$\boldsymbol {\theta } = \left (\mu , \tau ^{2}\right )^{\prime }$
. We use
$\ell (\boldsymbol {\theta }) = \ell \left (\mu , \tau ^{2}\right )$
to denote the log-likelihood function under the model (3).
The algorithms proposed to implement the two modifications to obtain the confidence limits of
$\mu $
under the model (3) are presented below.
For the numerical implementation of the proposed methods, a grid of values for
$\mu _{0}$
can be considered for both algorithms. The grid of values
$\mu $
can be selected to cover a range around the MLE of
$\mu $
, since the upper confidence limit will be larger than the MLE of
$\mu $
, and the lower confidence limit will be smaller. For the algorithm for modification I, to calculate a two-sided confidence interval for
$\mu $
, the upper and lower confidence limits are obtained by finding the values of
$\mu _0$
for which
$r_{1}\left (\mu _{0}\right )$
is close to the standard normal quantiles according to a specified threshold. For the algorithm for modification II, the confidence limits for
$\mu $
are determined by finding the values of
$\mu _0$
that make the p values
$P\left (r\left (\mu _{0}\right )^* \leq r\left (\mu _{0}\right ) \right )$
and
$P\left (r\left (\mu _{0}\right )^* \geq r\left (\mu _{0}\right ) \right )$
close to the significance level, once again according to a specified threshold.
3 Results
3.1 Simulation settings
In order to evaluate the performance of the proposed methods, a simulation study similar to that of Rover et al.Reference Rover, Knapp and Friede1 was conducted for the modified HKSJ (mKH) method with a small number of studies (k). The “relative amount of heterogeneity” is expressed in terms of the measure
$I^2$
defined as
$$ \begin{align} I^2 = \frac{\tau^2}{\tau^2 + \tilde{s}^2}, \end{align} $$
where
$ \tilde {s} $
is an “average” standard error among the study-specific
$ s_i $
values.Reference Rover, Knapp and Friede1,
Reference Higgins and Thompson14 In the simulation study,
$ \tilde {s}^2 $
was calculated as the arithmetic mean of the squared standard errors, i.e.,
$$ \begin{align} \tilde{s}^2 = \frac{1}{k} \sum_{i=1}^k \frac{2\sigma^2}{n_i} \end{align} $$
with
$ n_i $
being the group size of trial
$ i $
(
$ i = 1, \ldots , k $
) and
$ \sigma $
is a scaling factor. Since the results depend solely on the ratio
$ \tau /\sigma $
,
$ \sigma = 1 $
was set in the simulations.Reference IntHout, Ioannidis and Borm15 The simulation setup we have used and also adopted in the work of Rover et al.,Reference Rover, Knapp and Friede1 is taken from the study by IntHout et al.Reference IntHout, Ioannidis and Borm15
Similarly to the work by Rover et al.,Reference Rover, Knapp and Friede1 four different meta-analyses scenarios were considered in the simulation: (A) trials of equal size, (B) equally sized trials but including one small trial, (C) with half large and half small trials, and (D) equally sized trials and one large trial. A large trial is 10 times the size of a small trial, which means that the associated standard errors differ by approximately a factor of 3.Reference Rover, Knapp and Friede1 The number of trials k ranged from 2 to 10, that is,
$k = 2, \ldots , 10$
. The true level of heterogeneity was set at
$I^2 = 0.5$
; see (20). The average group size
$n_i$
in each trial was set at around 100. Without loss of generality, the true value of
$\mu $
was set to zero. The
$95\%$
confidence interval (CI) of
$\mu $
was constructed using the proposed methods and compared to the HKSJ methodReference Hartung and Knapp3 and the modified HKSJ (mKH) method for small k.Reference Rover, Knapp and Friede1
The simulation for normally distributed outcomes under the model (3) involved several key steps, as detailed in Appendix 2 of the work of IntHout et al.Reference IntHout, Ioannidis and Borm15 For each scenario with
$I^2 = 0.5$
, the variance
$\tau ^2$
was calculated using equations (20) and (21). For each trial
$ i $
, the ‘true’ trial-specific effect size
$\theta _i$
in (2) was drawn from a normal distribution with mean
$\mu = 0$
and variance
$\tau ^2$
. The trial outcomes were then generated from a normal distribution with mean
$\theta _i$
and variance
$\frac {2\sigma ^2}{n_i} = 2/n_i$
(recall that the scaling factor
$\sigma $
is chosen to be equal to one). Additionally, for each trial i with group size
$n_i$
, the variance of the trial outcome, i.e., the quantity
${s}_i^2$
, was generated based on a
$\chi ^2$
distribution with
$2n_i - 2$
degrees of freedom, divided by
$(n_i - 1)n_i$
; see the work of IntHout et al.Reference IntHout, Ioannidis and Borm15 For each simulation, the generated observed trial outcomes
$\boldsymbol {y}_{o}$
and within-trial variances
${s}_i^2$
for
$i = 1, \ldots , k,$
were used as input for the proposed algorithms to construct
$95\%$
CIs of
$\mu $
. For the fiducial approach, 5,000 realizations of the fiducial quantities were used to obtain the confidence limits. For asymptotic modifications I & II, 1000 bootstrap samples were used, i.e.,
$B = 1,000$
.
Plots of the 95% CI coverage, median and 90th quantiles of the CI width ratios compared to the HKSJ method for trials with half large/small sizes and average size around 100 (with a ceiling (
$k/2$
) of small sizes for kodd), by Asymptotic Modification I (Asym_I) method, Asymptotic Modification II (Asym_II) method, Fiducial approach, HKSJ method, and modified HKSJ (mKH) method.
Plots of the 95% CI coverage, median and 90th quantiles of the CI width ratios compared to the HKSJ method for trials with equal sizes of
$n = 100$
, by Asymptotic Modification I (Asym_I) method, Asymptotic Modification II (Asym_II) method, Fiducial approach, HKSJ method, and modified HKSJ (mKH) method.
Plots of the 95% CI coverage, median and 90th quantiles of the CI width ratios compared to the HKSJ method for the case of equally-sized trials with one large trial and average size around 100, by Asymptotic Modification I (Asym_I) method, Asymptotic Modification II (Asym_II) method, Fiducial approach, HKSJ method, and modified HKSJ (mKH) method.
Plots of the 95% CI coverage, median and 90th quantiles of the CI width ratios compared to the HKSJ method for the case of equally-sized trials with one small trial and average size around 100, by Asymptotic Modification I (Asym_I) method, Asymptotic Modification II (Asym_II) method, Fiducial approach, HKSJ method, and modified HKSJ (mKH) method.
To compare the proposed methods with the existing ones, the CIs of
$\mu $
were also constructed using the HKSJ methodReference Hartung and Knapp3 and the modified HKSJ (mKH) method.Reference Rover, Knapp and Friede1 The DerSimonian–Laird (DL) estimatorReference DerSimonian and Laird16 of the between-study variance
$\tau ^2$
was used for both the HKSJ and mKH methods in the simulation. According to Rover et al.Reference Rover, Knapp and Friede1 and Sánchez-Meca and Marín-Martínez,Reference Sánchez-Meca and Marín-Martínez17 the constructed CIs for both the HKSJ and mKH methods are similar if using the DL estimator, the restricted maximum likelihood (REML) estimator,Reference Viechtbauer18 and the Paule–Mandel (PM) estimator.Reference Paule and Mandel12
A total of 5,000 meta-analyses were simulated, and all simulations were performed using R.19 In order to compare the different approaches for computing the confidence intervals for
$\mu $
, the coverage probabilities and the expected widths of the CIs were calculated, where the latter was expressed as a ratio to the expected length of the interval obtained by the HKSJ method. The median of these ratios and the 90th percentile of the ratios were also evaluated, as in the work of Rover et al.Reference Rover, Knapp and Friede1
3.2 Simulation results
For trials with half large and half small sizes, with the average size around 100, the simulation results are shown in Figure 1. The results include
$95\%$
CI coverage, median and 90th quantiles of the expected CI length ratios mentioned above, for the number of trials
$k = 2, \ldots , 10$
. A large trial was 10 times the size of a small trial, with the overall trial average being around 100, and the number of small-sized trials was chosen to be
$ (k+1)/2 $
or
$ k/2 $
depending on whether k is odd or even. For example, for
$ k = 3 $
, the group sizes for the trials would be
$ n = (25, 25, 250) $
; for
$ k = 2 $
, the group sizes would be
$ n = (18, 180) $
.
The simulation results demonstrate that the proposed new methods, including the fiducial approach and the small sample asymptotic methods, particularly Asymptotic Modification II, perform better than the modified HKSJ (mKH) method when the number of trials is small. The HKSJ method underperforms for small
$ k $
, and while its performance is indeed enhanced by the modified HKSJ (mKH) method, the proposed new methods offer superior results. The fiducial approach consistently provides more accurate confidence interval (CI) coverage. Asymptotic Modification II, in particular, shows more accurate CI coverage and shorter CI lengths for small numbers of trials (
$ k = 2, 3 $
), compared to other methods. Additionally, Asymptotic Modification I also offers a more accurate CI coverage when
$ k> 4 $
. The fiducial approach yields larger median CI lengths, but with less variability in CI lengths, especially for
$ k = 2, 3 $
. On the other hand, Asymptotic Modification I results in CI lengths similar to mKH but with a larger median, and Asymptotic Modification II shows smaller CI lengths for
$ k = 2 $
and less variability of CI lengths for
$ k = 3, 4 $
.
Figure 2 displays the simulation results for trials with equal group sizes of 100. In this scenario, the HKSJ method can handle small k effectively. It works very well when the analyzed studies are of equal size (i.e., have equal standard errors), as shown in the work of Rover et al.Reference Rover, Knapp and Friede1
Figures 3 and 4 show the corresponding simulation results for the scenarios with one large trial, and with one small trial, respectively, with an average group size around 100. The results and conclusions are similar as in the case of trials with half-large and half-small sizes. The HKSJ method underperforms for small
$ k $
, while the mKH method improves its performance. The fiducial approach consistently provides more accurate CI coverage close to the nominal 95%. Asymptotic Modification II maintains more accurate CI coverage, outperforming the mKH method for small numbers of trials. For the CI lengths, the results reveal that the fiducial approach yields larger median CI lengths with less variability, while Asymptotic Modification II shows smaller and more stable CI lengths for small
$ k $
, indicating its robustness.
In summary, for computing a CI of
$\mu $
in unbalanced settings with a small number of studies k, Asymptotic Modification I shows results comparable to those based on existing methods. However, the fiducial approach and Asymptotic Modification II show better results, with more accurate CI coverages and shorter CI lengths. The findings highlight the superior performance of the newly proposed methods in terms of both CI coverage and length, making them excellent choices for practical applications.
4 Examples
An alternative methodology, not focused on in our work, is a Bayesian approach (see the works of Bender et al.,Reference Bender, Friede and Koch20 and Friede et al.,Reference Friede, Röver, Wandel and Neuenschwander21) which introduces stability through the incorporation of prior knowledge about the heterogeneity parameter
$\tau ^2$
through (weakly) informative prior distributions. The two examples of few studies illustrated in the work of Bender et al.Reference Bender, Friede and Koch20 are used here to demonstrate the performance of the proposed methods. These examples come from dossier assessments conducted by the Institute for Quality and Efficiency in Health Care.Reference Bender, Friede and Koch20 Assuming that the logarithm of hazard ratios (HR) or relative risks (RR) follows a normal distribution, 95% confidence intervals for the data in both examples were generated using the HKSJ method, the modified HKSJ (mKH) method, the Bayesian approach, and the proposed methods, including the Fiducial approach and Asymptotic Modification I and II. For Bayesian approaches, half-normal priors for the heterogeneity parameter
$\tau $
with scales of 0.5 (Bayesian-HN(0.5)) and 1.0 (Bayesian-HN(1)) were applied.
In the first example, the added benefit of belatacept compared to ciclosporin A in combination with corticosteroids and mycophenolate mofetil was assessed as the appropriate comparator therapy for the prophylaxis of graft rejection in adults receiving a renal transplant.Reference Bender, Friede and Koch20 Only
$k = 2$
studies were available. The results for this example are illustrated in Figure 5. The HKSJ method yielded a wide 95% CI due to insufficient data to reliably estimate heterogeneity, leading to a conservative result. The mKH method produced the same result as HKSJ in this case. The Bayesian approach provides an alternative, yielding narrower intervals than HKSJ and mKH. However, prespecification of the prior distribution for between-study variation
$\tau ^2$
is crucial to interpreting the intervals in a frequentist sense, as a hypothesis test.Reference Bender, Friede and Koch20 The application of a half-normal prior with scale 0.5 for
$\tau $
results in a statistically significant pooled effect estimate, while a scale of 1 does not. Among the proposed methods, Asymptotic Modifications I and II generate similar but slightly shorter confidence intervals compared to the Bayesian method with a half-normal prior of scale 1. The Fiducial approach provides slightly wider intervals than the asymptotic modification methods but remains narrower than those from HKSJ and mKH.
A forest plot showing random-effects meta-analysis results for the first belatacept example, by HKSJ method, modified HKSJ (mKH) method, Bayesian method using half-normal priors for
$\tau $
with scales 0.5 (Bayesian-HN(0.5)) and 1 (Bayesian-HN(1)), Fiducial approach, and Asymptotic Modification I (Asym_I) and II (Asym_II) methods.
A forest plot showing random-effects meta-analysis results for the second sipuleucel-T example, by HKSJ method, modified HKSJ (mKH) method, Bayesian method using half-normal priors for
$\tau $
with scales 0.5 (Bayesian-HN(0.5)) and 1 (Bayesian-HN(1)), Fiducial approach, and Asymptotic Modification I (Asym_I) and II (Asym_II) methods.
In the second example, the added benefit of sipuleucel-T compared to the appropriate comparator therapy to treat asymptomatic or minimally symptomatic metastatic (non-visceral) castrate-resistant prostate cancer was assessed in male adults for whom chemotherapy is not yet clinically indicated.Reference Bender, Friede and Koch20 Here,
$k = 3$
relevant studies were available. The results for this example are shown in Figure 6. Bayesian approaches yield shorter confidence intervals, indicating a statistically significant difference between treatment groups, to the disadvantage of sipuleucel-T. The mKH method produces the same wide confidence interval as HKSJ, including the relative risk of 1. The three proposed methods generate wider confidence intervals, particularly on the right-hand side, covering the finding from the second study with a large upper confidence limit for the relative risk, which may suggest better coverage of the true value. Although Bayesian approaches tended to give qualitatively different results relative to the other methods, the choice of method generally had negligible impact on the conclusions drawn in this example.
5 Discussion
In this work, we have proposed two novel methods for computing confidence limits for
$\mu $
in random-effects meta-analysis, particularly suited for a small number of studies. One method is based on fiducial inference, while the other involves two small sample asymptotic procedures, referred to as Asymptotic Modifications I and II. Our simulation studies indicate that for unbalanced settings, when the number of studies k is small, the proposed methods show significant improvements compared to existing methods. The fiducial approach and Asymptotic Modification II exhibit superior performance, providing more accurate CI coverages and shorter CI lengths compared to those based on the modified HKSJ method proposed by Rover et al.Reference Rover, Knapp and Friede1 Also, Asymptotic Modification I provides results that are comparable to those based on the modified HKSJ method. We believe that these findings have important practical implications since meta-analyses with a small number of studies are frequently encountered in applications. It is noted that for
$k=2$
, none of the competing approaches yields useful information, so the differences in relative performance of the approaches do not have practical relevance.
For methods based on small sample asymptotics, other modifications of the statistic
$r(\psi )$
given in (18) are also available in the literature, with the goal of achieving improved small sample performance; see the work of Brazzale et al.Reference Brazzale, Davison and Reid13 Such modifications have in fact been applied in the context of the model and problem investigated in our work; see the work of Guolo.Reference Guolo22 Figures 1-3 in the work of GuoloReference Guolo22 show that such modifications exhibit a performance quite similar to what we have noted for Asymptotic Modifications I and II considered in our investigation. We also want to point out that in the context of inter-laboratory studies where the model (3) is assumed, but with unknown within-trial variances, the fiducial approach was developed in the work of TianReference Tian23 and small sample asymptotic procedures were investigated by Sharma and Mathew.Reference Sharma and Mathew24 Furthermore, Zejnullahi and Hedges Reference Zejnullahi and Hedges25 introduced two alternative robust variance estimators were introduced by Zejnullahi and HedgesReference Zejnullahi and Hedges25 for small meta-analyses, evaluating degrees-of-freedom adjustments for confidence intervals.
In addition to the methodological contributions, we are developing an R package that will offer a practical tool for the easy implementation of the proposed methods.
Future research could explore the impact of estimated within-trial variability on the performance of the CIs of
$\mu $
. Furthermore, the median or mean of the fiducial quantities could be used as a point estimate of
$\mu $
. The performance of such a point estimate of
$\mu $
could be evaluated and compared with the maximum likelihood estimate (MLE), the restricted MLE (REML), and the estimate resulting from the DerSimonian–Laird (DL) approach. These are currently under investigation.
Acknowledgements
Our sincere thanks to Professor Javier Cabrera of the Department of Statistics at Rutgers University for his valuable discussions and suggestions on statistical methods and simulation study design.
Author contributions
Conceptualization: T.M., D.A.; Methodology and Investigation: T.M., D.A., Y.D., G.C.; Formal analysis, Data curation, Software and Visualization: Y.D., G.C.; Supervision and Validation: T.M., D.A.; Writing—original draft: Y.D., T.M., G.C.; Writing—review and editing: T.M., D.A., Y.D., G.C.
Competing interest statement
The authors declare that no competing interests exist.
Data availability statement
Program codes, simulated data summary, and real data examples supporting the findings of this study can be found in the Supplementary Material.
Funding statement
The authors declare that no specific funding has been received for this article.
Supplementary material
The supplementary material for this article can be found at http://doi.org/10.1017/rsm.2025.10022.
Open access
