# Random-Effect Meta-Analysis with Robust Between-Study Variance

Keisuke Hanada and Tomoyuki Sugimoto Graduate School of Engineering Science,
Osaka University. Address:1-3 Machikaneyama-cho, Toyonaka, Osaka, Japan. Tel:
+81-(0)80-1882-4281 E-Mail<EMAIL_ADDRESS>School of
Engineering Science, Osaka University.

###### Abstract

Meta-analyses are widely employed to demonstrate strong evidence across
numerous studies. On the other hand, in the context of rare diseases, meta-
analyses are often conducted with a limited number of studies in which the
analysis methods are based on theoretical frameworks assuming that the
between-study variance is known. That is, the estimate of between-study
variance is substituted for the true value, neglecting the randomness with the
between-study variance estimated from the data. Consequently, excessively
narrow confidence intervals for the overall treatment effect for meta-analyses
have been constructed in only a few studies. In the present study, we propose
overcoming this problem by estimating the distribution of between-study
variance using the maximum likelihood-like estimator. We also suggest an
approach for estimating the overall treatment effect via the distribution of
the between-study variance. Our proposed method can extend many existing
approaches to allow more adequate estimation under a few studies. Through
simulation and analysis of real data, we demonstrate that our method remains
consistently conservative compared to existing methods, which enables meta-
analyses to consider the randomness of the between-study variance.

Keywords: Bias-robust estimator; Meta-analysis; Quantile estimation; Random-
effects model, Small sample.

##  1 Introduction

Meta-analyses can provide strong evidence by integrating multiple studies with
a similar research question (DerSimonian and Laird, 1986). Conducting numerous
large-scale studies can be challenging in fields such as rare diseases because
the number of patients with each rare disease is small, and the number of
studies is limited. Random-effects meta-analyses allow for the incorporation
of between-study heterogeneity as the between-study variance. This
heterogeneity encompasses various factors, such as differences in study
conduct, outcome variable measurement methods, protocol discrepancies, and
unobserved confounders not accounted for in the data.

However, it is known that the expected variance of the overall treatment
effect in random-effects meta-analyses is less than the Cramér-Rao lower
bounds (Li et al., 1994; Viechtbauer, 2005). In other words, an
underestimation of the between-study variance may be possible under a few
studies. This underestimation can lead to the erroneous conclusion that
ineffective treatments are effective, potentially resulting in harmful
outcomes for patients. To address this issue in scenarios with few studies,
conservative estimation methods have been proposed (Hartung, 1999; Röver et
al., 2015; Michael et al., 2019). Among these methods, the Hartung-Knapp-
Sidik-Jonkman (HKSJ) approach stands out for its applicability regardless of
the methods chosen for estimating between-study variance (Hartung, 1999).
Although the HKSJ approach may be underestimated when including studies with
large or small sample sizes, the modified Knapp-Hartung method gives a
correction in the situation where the HKSJ approach underestimates (Röver et
al., 2015).

Conservative estimation methods in meta-analyses, such as the HKSJ approach,
address the issue of underestimating between-study variance with few studies.
On the other hand, these methods disregard the inherent randomness of the
estimated between-study variance. However, ignoring this randomness may lead
to an inaccurate estimation of the overall treatment effect. Typical interval
estimation methods for the overall treatment effect like DerSimonian–Laird and
restricted maximum likelihood assume the knowledge of the true between-study
variance (DerSimonian and Laird, 1986; Raudenbush, 2009; Sidik and Jonkman,
2007). In actual data analyses, however, we usually only rely on an estimate
of the between-study variance obtained from the data. The estimated value is
then substituted for the true value within these methods. Similarly, the HKSJ
approach utilizes the estimated between-study variance, but neglects the
inherent randomness associated with this variance estimate.

The inherent randomness associated with estimating between-study variance has
not been considered traditionally in meta-analyses for two primary reasons.
First, the distribution of the between-study variance often includes unknown
parameters. Second, a universally accepted method to determine the exact
distribution of the between-study variance is lacking. In a single meta-
analysis, we can only obtain a single estimate of the between-study variance.
This makes it challenging to identify the unknown parameter within the
distribution of the between-study variance. For example, although the method
based on the almost-exact distribution can provide accurate estimates of the
overall treatment effect, the true value of the between-study variance is
still required to determine the unknown parameter (Hanada and Sugimoto, 2023).
While some estimation methods, such as the DerSimonian-Laird (DL), are known
for the distribution of the between-study variance, the distribution of the
between-study variance by the general estimation method remains elusive (Sidik
and Jonkman, 2022).

One potential approach to address the unknown parameters in the distribution
of the between-study variance involves assuming that the between-study
variance estimator is a maximum likelihood-like estimator, which we call the
M-estimator. This allows for selecting the optimal unknown parameters within
the distribution of the between-study variance. The M-estimator composes a
general class of estimators that extend the concept of maximum likelihood
estimation and are particularly useful for robust estimation (Huber, 1992). A
specific optimization function is chosen based on several conditions and
characteristics in constructing an M-estimator. Among the functions, the
absolute value function provides the most bias-robust estimator and is
considered to be optimal for estimating the location parameters (Rousseeuw,
1982).

This paper proposes a novel estimation procedure for the overall treatment
effect. It incorporates the previously neglected randomness associated with
between-study variance. This approach leverages the framework of the
M-estimators. We further derive the distributions for the commonly used
general moment estimator and Sidik-Jonkman (SJ) estimator of between-study
variance because our approach needs the distribution of the between-study
variance. We additionally introduce an application of the M-estimator approach
within the framework of the HKSJ approach. To evaluate the performance of our
proposed approach compared with existing methods, we conduct simulations and
analyse a practical meta-analysis, focusing specifically on scenarios with a
limited number of studies.

This paper is organized as follows. Section 2 introduces the random-effect
meta-analysis using typical methods and proposes an approach based on the
M-estimator. Section 3 describes the theoretical properties of the proposed
approach, shows the distributions of a general moment estimator and SJ
estimator of between-study variance, and proposes an extension to some
conservative approach. Section 4 describes the setting and results of computer
simulations and presents the results of an analysis of the data from patients
with juvenile idiopathic arthritis (JIA) patients. Finally, Section 5
concludes.

##  2 Methods

We discuss typical approaches to random-effect meta-analysis and propose a
novel inference approach for the overall treatment effect based on a bias-
robust estimation approach.

###  2.1 Random-Effect Meta-Analysis

We conduct a meta-analysis using $K$ studies. The $k$-th study treatment
effect $\theta_{k}$ is estimated as

$\displaystyle\hat{\theta}_{k}|\theta_{k},\sigma_{k}^{2}\sim
N(\theta_{k},\sigma_{k}^{2}),$

where $\hat{\theta}_{k}$ and $\sigma_{k}^{2}$ are an estimator of true effect
$\theta_{k}$ and a true within-study variance in the $k$-th study,
respectively. The random-effect meta-analysis assumes that the $k$-th study
effect $\theta_{k}$ is distributed as

$\displaystyle\theta_{k}|\theta,\tau^{2}\sim N(\theta,\tau^{*2}),$

where $\theta$ and $\tau^{*2}$ are the true overall treatment effect and
between-study variance, respectively. The model can be rewritten by
marginalization as

$\displaystyle\hat{\theta}_{k}|\theta,\sigma_{k}^{2},\tau^{*2}\sim
N(\theta,\sigma_{k}^{2}+\tau^{*2}).$ (1)

Model (1) aims to infer the true overall treatment effect $\theta$. The
between-study variance $\tau^{*2}$, which represents an effect difference
between studies, affects the estimation of $\theta$.

The overall treatment effect can be estimated in the following form:

$\displaystyle\hat{\theta}(\tau^{*2})$
$\displaystyle=\frac{\sum_{k=1}^{K}w_{k}(\tau^{*2})\hat{\theta}_{k}}{\sum_{k=1}^{K}w_{k}(\tau^{*2})},$
(2)

where $w_{k}(\tau^{*2})=(\sigma_{k}^{2}+\tau^{*2})^{-1}$. We can calculate the
expectation $E[\hat{\theta}(\tau^{*2})]=\theta$ and the variance
$V[\hat{\theta}(\tau^{*2})]=\left(\sum_{k=1}^{K}w_{k}(\tau^{*2})\right)^{-1}$.
Thus, the estimator $\hat{\theta}(\tau^{*2})$ is unbiased and consistent. The
test statistic
$T(\tau^{*2})=(\hat{\theta}(\tau^{*2})-\theta_{0})/\sqrt{V[\hat{\theta}]}$ is
standard normally distributed in a null hypothesis, $\theta=\theta_{0}$. The
$100(1-\alpha)\%$ confidence interval (CI) is obtained as

$\displaystyle\hat{\theta}(\tau^{*2})\pm
z_{\alpha/2}\left(\sum_{k=1}^{K}w_{k}(\tau^{*2})\right)^{-1},$ (3)

where $z_{\alpha/2}$ is the $100(1-\alpha/2)\%$ point of the standard normal
distribution. However, $\tau^{*2}$ is usually unknown in real data analysis,
so a between-study variance estimator $\hat{\tau}^{2}$ may be substituted for
$\tau^{*2}$ in equation (2). The estimation of $\tau^{*2}$ is discussed by
many authors. For example, DerSimonian and Laird (1986) and Sidik and Jonkman
(2007) propose moment-based estimation, and (Raudenbush, 2009) provides
restricted maximum likelihood estimation.

###  2.2 Proposed M-estimator of $\tau^{*2}$ from $\hat{\tau}^{2}$

Although the CI (3) with true between-study variance $\tau^{*2}$ achieves a
nominal level coverage probability, the CI of (3) in which $\tau^{*2}$ is
inserted by some estimates $\hat{\tau}^{2}$ may not reach the nominal level
because $\hat{\tau}^{2}$ is regarded as a constant despite being a random
variable. We consider the marginal distribution of the overall treatment
effect with the between-study variance estimator using the framework of the
M-estimator and propose the CI and p-value based on the marginal distribution.

Let $f_{\hat{\tau}^{2}}(\bullet|\tau^{*2})$ be the probability density
function of $\hat{\tau}^{2}$ that depends on some parameter $\tau^{*2}$.
Because the estimated between-study variance $\hat{\tau}^{2}$ is constructed
via the (empirical) expectation given the true $\tau^{*2}$ in general, the
parameter located at $\tau^{*2}$ can be formulated as a location parameter of
$f_{\hat{\tau}}(\bullet|\tau^{*2})$. The distribution function of the overall
treatment effect is given by marginalizing for between-study variance as
follows:

$\displaystyle F_{\hat{\theta}}(y|\tau^{*2})$
$\displaystyle=\int_{-\infty}^{y}\int_{-\infty}^{\infty}f_{\hat{\theta}|\hat{\tau}^{2}}(z|x)f_{\hat{\tau}^{2}}(x|\tau^{*2})dxdz$
$\displaystyle=\int_{-\infty}^{\infty}F_{\hat{\theta}|\hat{\tau}^{2}}(y|x)f_{\hat{\tau}^{2}}(x|\tau^{*2})dx,$

where $\hat{\theta}$ is normally distributed with mean $\theta$ and variance
$1/\sum_{k=1}^{K}w_{k}(x)$ given $\hat{\tau}^{2}=x$. In fact, considering that
the true parameter $\tau^{*2}$ is unknown in real data analysis, we have to
consider how we estimate $\tau^{*2}$ as well as the variability.

Let us consider constructing the between-study variance estimator
$\hat{\tau}^{2}$ by the M-estimation from the distribution
$f_{\hat{\tau}^{2}}(\bullet|\tau^{2})$. Define the estimation function
$R(\tau^{2})$ for the M-estimator (which is called by Huber (1992)) by

$\displaystyle R(\tau^{2})$ $\displaystyle=\underset{t\geq
0}{\mathrm{argmin}}~{}E[\rho(\hat{\tau}^{2}-t)|\tau^{2}]$
$\displaystyle=\underset{t\geq
0}{\mathrm{argmin}}~{}\int_{-\infty}^{\infty}\rho(x-t)f_{\hat{\tau}^{2}}(x|\tau^{2})dx.$

When $\rho(u)$ is a real-valued convex function satisfying the symmetry
condition $\rho(u)=\rho(-u)$, the estimation function $R(\tau^{2})$ can
estimate a location parameter depending on the convex function (Takeuchi,
2020). If an estimator $\hat{\tau}^{2}$ is the M-estimator, we have the
estimation equation

$\displaystyle\hat{\tau}^{2}$ $\displaystyle=R(\tau^{2}).$ (4)

Therefore, we can find a location parameter $\tau^{2}$ of
$f_{\hat{\tau}^{2}}(x|\tau^{2})$ that best fits an estimated value
$\hat{\tau}^{2}$ by solving the equation (4) for $\tau^{2}$. Thus, a new
estimator of the location parameter
$\hat{\tau}_{M}^{2}=R^{-1}(\hat{\tau}^{2})$ can be obtained. We provide two
propositions below, but the details of statements and proofs are shown in
Supplementary material A.2 because they do not align with our main objective:
i) the variance of the estimated parameter $\hat{\tau}_{M}^{2}$ can be
expressed as a function of the estimated between-study variance
$\hat{\tau}^{2}$; and ii) the estimated location parameter
$\hat{\tau}_{M}^{2}$ is consistent when $\hat{\tau}^{2}$ is consistent. We can
estimate the distribution function of the overall treatment effect
$F_{\hat{\theta}}(y|\tau^{*2})$ by

$\displaystyle F_{\hat{\theta}}(y|\hat{\tau}_{M}^{2})$
$\displaystyle=\int_{-\infty}^{\infty}F_{\hat{\theta}|\hat{\tau}^{2}}(y|x)f_{\hat{\tau}^{2}}(x|\hat{\tau}_{M}^{2})dx,$
(5)

and the $100(1-\alpha)\%$ CI of the overall treatment effect can be
approximately constructed as

$\displaystyle\left(F_{\hat{\theta}}^{-1}\left(\frac{\alpha}{2}|\hat{\tau}_{M}^{2}\right),F_{\hat{\theta}}^{-1}\left(1-\frac{\alpha}{2}|\hat{\tau}_{M}^{2}\right)\right).$

The one- and two-sided p-values for the null hypothesis
$H_{0}:\theta=\theta_{0}$ can be derived from the distribution function (5) as
$1-F_{\hat{\theta}}(|\hat{\theta}(\hat{\tau}^{2})-\theta_{0}||\hat{\tau}_{M}^{2})$
and
$2(1-F_{\hat{\theta}}(|\hat{\theta}(\hat{\tau}^{2})-\theta_{0}||\hat{\tau}_{M}^{2}))$,
respectively.

##  3 Theoretical Aspect of M-estimator Approach

We present some properties of the M-estimator approach in Section 3.1.
Similary, we investigate some estimators of the between-study variance and
discuss their exact distributions in Section 3.2. In parallel, we apply the
M-estimator approach to two estimation procedures for more accurate estimation
in Section 3.3.

###  3.1 Properties of M-estimator of Between-Study Variance

First, we show one of the main results for a property of the CI of the overall
treatment effect based on the M-estimator approach.

Theorem 1: Suppose that
$E[\hat{\tau}^{2}|\hat{\tau}_{M}^{2}]\leq\hat{\tau}^{2}$ under the random-
effect model (1). Let $\hat{\theta}_{L}$ and $\hat{\theta}_{U}$ be the
$100\alpha/2\%$ and $100(1-\alpha/2)\%$ points of
$F_{\hat{\theta}}(\bullet|\hat{\tau}_{M}^{2})$, respectively. We have

$\displaystyle F_{\hat{\theta}|\hat{\tau}^{2}}(q|\hat{\tau}^{2})$
$\displaystyle>F_{\hat{\theta}}(q|\hat{\tau}_{M}^{2})\,\,\,\,\,(q>\theta)$
$\displaystyle F_{\hat{\theta}|\hat{\tau}^{2}}(q|\hat{\tau}^{2})$
$\displaystyle<F_{\hat{\theta}}(q|\hat{\tau}_{M}^{2})\,\,\,\,\,(q<\theta)$
$\displaystyle F_{\hat{\theta}|\hat{\tau}^{2}}(q|\hat{\tau}^{2})$
$\displaystyle=F_{\hat{\theta}}(q|\hat{\tau}_{M}^{2})\,\,\,\,\,(q=\theta),$

and the following relation on confidence limits $\hat{\theta}_{L}$ and
$\hat{\theta}_{U}$

$\displaystyle
F_{\hat{\theta}|\hat{\tau}^{2}}(\hat{\theta}_{L}|\hat{\tau}^{2})\geq
F_{\hat{\theta}}(\hat{\theta}_{L}|\hat{\tau}^{2}_{M})$
$\displaystyle\to_{p}\frac{\alpha}{2}~{}~{}~{}~{}~{}~{}~{}(K\to\infty)$
$\displaystyle
F_{\hat{\theta}|\hat{\tau}^{2}}(\hat{\theta}_{U}|\hat{\tau}^{2})\leq
F_{\hat{\theta}}(\hat{\theta}_{U}|\hat{\tau}^{2}_{M})$
$\displaystyle\to_{p}1-\frac{\alpha}{2}~{}~{}(K\to\infty),$

where $\to_{p}$ denotes the convergence in probability.

The proof of Theorem 1 is shown in Supplementary material A.1. Theorem 1
implies that the CI based on the M-estimator approach is more conservative
than the CI without the M-estimator approach. Furthermore, the Type-I error
based on the M-estimator approach converges to a nominal significance level as
the number of studies increases.

We focus on the convex function $\rho(u)=|u|$ because it is the most bias-
robust estimator (Rousseeuw, 1982). The M-estimator can be rewritten as

$\displaystyle\hat{\tau}^{2}$ $\displaystyle=\underset{t\geq
0}{\mathrm{argmin}}~{}E_{\hat{\tau}^{2}}[|X-t||\hat{\tau}^{2}_{M}]$
$\displaystyle=\underset{t\geq
0}{\mathrm{argmin}}~{}\left[\int_{t}^{\infty}(x-t)f_{\hat{\tau}^{2}}(x|\hat{\tau}_{M}^{2})dx-\int_{-\infty}^{t}(x-t)f_{\hat{\tau}^{2}}(x|\hat{\tau}_{M}^{2})dx\right]$
$\displaystyle=\left\\{t\in[0,\infty)\left|\int_{t}^{\infty}f_{\hat{\tau}^{2}}(x|\hat{\tau}_{M}^{2})dx-\int_{-\infty}^{t}f_{\hat{\tau}^{2}}(x|\hat{\tau}_{M}^{2})dx=0\right.\right\\}$
$\displaystyle=F_{\hat{\tau}^{2}}^{-1}\left(\left.\frac{1}{2}\right|\hat{\tau}^{2}_{M}\right).$

Hence, the location parameter can be estimated as

$\displaystyle\hat{\tau}_{M}^{2}$
$\displaystyle=\left\\{\tau^{2}\in[0,\infty)\left|F_{\hat{\tau}^{2}}(\hat{\tau}^{2}|\tau^{2})=\frac{1}{2}\right.\right\\}.$

Considering that the estimated between-study variance $\hat{\tau}^{2}$ is the
median of the probability distribution
$f_{\hat{\tau}^{2}}(\bullet|\hat{\tau}^{2}_{M})$, if
$E[\hat{\tau}^{2}|\hat{\tau}_{M}^{2}]=\hat{\tau}_{M}^{2}$ and
$E[(\hat{\tau}^{2}-\hat{\tau}_{M}^{2})^{3}|\hat{\tau}^{2}_{M}]\geq 0$, then
the condition of Theorem 1 is satisfied,
$\hat{\tau}_{M}^{2}=E[\hat{\tau}^{2}|\hat{\tau}_{M}^{2}]\leq
F_{\hat{\tau}^{2}}^{-1}(\frac{1}{2}|\hat{\tau}_{M}^{2})=\hat{\tau}^{2}$. This
inequality is satisfied when the distribution has a long right tail. Hence, it
is found that the estimator of the between-study variance with a long right
tail distribution satisfies the conditions of Theorem 1. In Section 3.2, we
discuss the distribution of the estimated between-study variance.

###  3.2 Distribution of Between-Study Variance Estimator $\hat{\tau}^{2}$

We illustrate the distribution of the general moment estimator of the between-
study variance because our proposed M-estimator needs the probability density
function. After that, we apply the distribution to the DL and SJ methods,
which are special cases of the general moment estimator.

The general moment estimator of the between-study variance can be obtained as
(Veroniki et al., 2016):

$\displaystyle\hat{\tau}_{G}$ $\displaystyle=\max\\{0,\hat{\tau}_{G,u}\\},$

with untruncated estimator

$\displaystyle\hat{\tau}_{G,u}$
$\displaystyle=\frac{Q_{a}-\left(\sum_{k=1}^{K}a_{k}\sigma_{k}^{2}-A_{1}^{-1}\sum_{k=1}^{K}a_{k}^{2}\sigma_{k}^{2}\right)}{A_{1}-A_{2}/A_{1}},$

where $Q_{a}=\sum_{k=1}^{K}a_{k}(\hat{\theta}_{k}-\hat{\theta}_{a})$ is the
generalized Cochran between-study variance statistic, $a_{k}$ represents the
weights assigned to each study that is equal to any positive value,
$A_{r}=\sum_{k=1}^{K}a_{k}^{r}$, and
$\hat{\theta}_{a}=A_{1}^{-1}\sum_{k=1}^{K}a_{k}\hat{\theta}_{k}$.

Theorem 2: Under the random-effect model (1), the exact distribution of
$\hat{\tau}_{G,u}$ is given by

$\displaystyle\hat{\tau}_{G,u}\sim\frac{\sum_{r=1}^{R}\lambda_{r}\chi_{r(1)}^{2}-\left(\sum_{k=1}^{K}a_{k}\sigma_{k}^{2}-A_{1}^{-1}\sum_{k=1}^{K}a_{k}^{2}\sigma_{k}^{2}\right)}{A_{1}-A_{2}/A_{1}},$

where $\lambda_{1},\cdots,\lambda_{R}$ are non-negative eigenvalues of matrix
$VW$, $V$ is the variance-covariance matrix of
$(\hat{\theta}_{1}-\hat{\theta}_{a},\cdots,\hat{\theta}_{K}-\hat{\theta}_{a})$,
and the element $(k,k^{\prime})$ of $V$ is given as

$\displaystyle V_{kk^{\prime}}$
$\displaystyle=\left\\{\begin{array}[]{lc}\left(\sigma_{k}^{2}+\frac{1}{A_{1}^{2}}\sum_{i=1}^{K}a_{i}^{2}\sigma_{i}^{2}-\frac{2a_{k}\sigma_{k}^{2}}{A_{1}}\right)+\left(1+\frac{A_{2}}{A_{1}^{2}}-\frac{2a_{k}}{A_{1}}\right)\tau^{2}&(k=k^{\prime})\\\
\left(-\frac{a_{k}\sigma_{k}^{2}+a_{k^{\prime}}\sigma_{k^{\prime}}^{2}}{A_{1}}+\frac{1}{A_{1}^{2}}\sum_{i=1}^{K}a_{i}^{2}\sigma_{i}^{2}\right)+\left(\frac{A_{2}}{A_{1}^{2}}-\frac{a_{k}+a_{k^{\prime}}}{A_{1}}\right)\tau^{2}&(k\neq
k^{\prime})\end{array}\right.,$ (8)

$W$ is a diagonal matrix with $a_{1},\cdots,a_{K}$, and
$\chi_{1(1)}^{2},\cdots,\chi_{R(1)}^{2}$ are independently and identically
distributed as a chi-square distribution with $1$ degree of freedom.

The proof of Theorem 2 is shown in Supplementary material A.3. The most
popular DL estimator is a special case of a general moment estimator with
weights $a_{k}=1/\sigma_{k}^{2}$. The between-study variance estimator is
given as

$\displaystyle\hat{\tau}_{DL}$
$\displaystyle=\max(0,\hat{\tau}_{u}^{2}),\,\,\,\hat{\tau}_{u}^{2}=\frac{Q-(K-1)}{W_{1}(0)-W_{2}(0)/W_{1}(0)}$
$\displaystyle Q$
$\displaystyle=\sum_{k=1}^{K}w_{k}(0)(\hat{\theta}_{k}-\hat{\theta}(0))^{2},$

where $W_{r}(\tau^{2})=\sum_{k=1}^{K}w_{k}^{r}(\tau^{2})$. The element
$(k,k^{\prime})$ of the variance-covariance matrix of
$(\hat{\theta}_{1}-\hat{\theta}(0),\cdots,\hat{\theta}_{K}-\hat{\theta}(0))$
is given as

$\displaystyle V_{kk^{\prime}}$
$\displaystyle=\left\\{\begin{array}[]{lc}\left(\sigma_{k}^{2}-\frac{1}{W_{1}(0)}\right)+\left(1+\frac{W_{2}(0)}{W_{1}^{2}(0)}-\frac{2}{\sigma_{k}^{2}W_{1}(0)}\right)\tau^{*2}&(k=k^{\prime})\\\
-\frac{1}{W_{1}(0)}+\left(\frac{W_{2}(0)}{W_{1}^{2}(0)}-\frac{1/\sigma_{k}^{2}+1/\sigma_{k^{\prime}}^{2}}{W_{1}(0)}\right)\tau^{*2}&(k\neq
k^{\prime})\end{array}\right..$

Thus, we can construct a CI of the between-study variance from the probability
distribution function as

$\displaystyle\left(F_{\hat{\tau}_{DL}^{2}}^{-1}\left(\frac{\alpha}{2}|\tau^{*2}\right),F_{\hat{\tau}_{DL}^{2}}^{-1}\left(1-\frac{\alpha}{2}|\tau^{*2}\right)\right),$

with unknown parameter $\tau^{*2}$. M-estimator $\hat{\tau}_{M}^{2}$ is
available to estimate the between-study variance instead of the unknown
parameter $\tau^{*2}$.

The SJ estimator is another moment estimator considered in a different
framework than the general moment estimator. The SJ between-study variance
estimator is given as

$\displaystyle\hat{\tau}_{SJ}^{2}$
$\displaystyle=\frac{1}{K-1}\sum_{k=1}^{K}v_{k}^{-1}(\hat{\theta}_{k}-\bar{\theta}_{v})^{2},$

with

$\displaystyle\bar{\theta}_{v}$
$\displaystyle=\frac{\sum_{k=1}^{K}v^{-1}_{k}\hat{\theta}_{k}}{\sum_{k=1}^{K}v_{k}^{-1}}~{}\mathrm{and}~{}v_{k}=1+\frac{\sigma_{k}^{2}}{\tau^{*2}}.$

Typically,
$\hat{v}_{k}=1+K\sigma_{k}^{2}/\sum_{k=1}^{K}(\hat{\theta}_{k}-\bar{\theta})^{2}$
is instead of $v_{k}$, where
$\bar{\theta}=\frac{1}{K}\sum_{k=1}^{K}\hat{\theta}_{k}$. In our approach, the
SJ estimator with an unknown parameter $\tau^{2}$ can be estimated as the
solution of the following equation:

$\displaystyle\hat{\tau}_{SJ}(\hat{\tau}^{2}_{M})$
$\displaystyle=\underset{t\geq
0}{\mathrm{argmin}}~{}E_{\hat{\tau}_{SJ}^{2}}[\rho(X-t)|\hat{\tau}_{M}^{2}],$

where $X$ is a random variable distributed as $\hat{\tau}_{SJ}^{2}$ with the
unknown parameter $\tau^{*2}$ instead of $\hat{\tau}_{M}^{2}$. We show the
distribution of the SJ estimator of the between-study variance as the next
theorem.

Theorem 3: Under the random-effect model (1), the SJ estimator of between-
study variance is distributed as

$\displaystyle\hat{\tau}_{SJ}^{2}(\tau^{2})\sim\frac{\tau^{2}}{K-1}\chi_{K-1}^{2},$

where $\chi_{K-1}^{2}$ is a weighted chi-square distribution with $K-1$
degrees of freedom.

The proof of Theorem 3 is shown in Supplementary material A.4. In Sidik and
Jonkman (2005), $(K-1)\hat{\tau}_{SJ}^{2}/\tau^{*2}$ is approximated by
$\chi_{K-1}^{2}$. The remarkable point of Theorem 3 states that it is not an
approximation but strictly follows $\chi_{K-1}^{2}$. Thus, we can construct a
CI as well as the DL method. From Theorem 3, $\hat{\tau}_{SJ}^{2}$ satisfies
the condition of Theorem 1,
$\hat{\tau}_{SJ}^{2}=F^{-1}(\frac{1}{2}|\hat{\tau}_{M}^{2})\geq
E[\hat{\tau}_{SJ}^{2}|\hat{\tau}_{M}^{2}]$. The CI for the overall treatment
effect by the SJ method using an M-estimator is more conservative than that of
the typical SJ method.

###  3.3 Expansion to Some Estimation Procedure

We expand two estimation procedures, the HKSJ approach and the method proposed
by Michael et al. (2019), to apply the M-estimator approach for more accurate
estimation. Using several estimators for the between-study variance tends to
shorten the interval estimation of the overall treatment effect in the meta-
analysis, prompting proposals to modify the estimation method. A method
proposed to adjust the test statistic to approximate a $t$-distribution with
$K-1$ degrees of freedom has been put forward, ensuring the accuracy of both
testing and CI (Hartung, 1999; Hartung and Knapp, 2001a, b; Sidik and Jonkman,
2002). This approach is known as the HKSJ method. For instance, when adjusting
the test statistic of the DL method, the HKSJ method modifies the test
statistic as follows:

$\displaystyle T_{HKSJ}$
$\displaystyle=\frac{\sqrt{\sum_{k=1}^{K}w_{k}(\hat{\tau}_{DL}^{2})}\hat{\theta}_{DL}}{\sqrt{\left\\{\sum_{k=1}^{K}w_{k}(\hat{\tau}_{DL}^{2})(\hat{\theta}_{k}-\hat{\theta}_{DL})^{2}/(K-1)\right\\}}}$

The test statistic $T_{HKSJ}$ approximately follows a $t$-distribution with
$K-1$ degrees of freedom. Therefore, the CI for the overall treatment effect
$\theta$ can be computed as

$\displaystyle\hat{\theta}_{DL}\pm
t_{K-1,\alpha/2}\sqrt{\frac{\sum_{k=1}^{K}w_{k}(\hat{\tau}_{DL}^{2})(\hat{\theta}_{k}-\hat{\theta}_{DL})^{2}}{(K-1)\sum_{k=1}^{K}w_{k}(\hat{\tau}_{DL}^{2})}},$

where $t_{K-1,\alpha/2}$ denotes the $100(1-\alpha/2)\%$ point of the
$t$-distribution with $K-1$ degrees of freedom. In many cases, the HKSJ
approach is more accurate than the ordinary approach directly applied.
However, it is known that including studies with an extremely large or small
number of subjects under a small number of studies reduces the performance
(Röver et al., 2015).

We consider applying our method to the HKSJ approach. The probability
distribution function of $\hat{\theta}$ with HKSJ correction is given as

$\displaystyle F_{\hat{\theta}}(y)$
$\displaystyle=\int_{-\infty}^{\infty}\Pr\left\\{\hat{\theta}\leq
y|\tau^{2}=x\right\\}f_{\hat{\tau}^{2}}(x|\hat{\tau}_{M}^{2})dx$
$\displaystyle=\int_{-\infty}^{\infty}\Pr\left\\{Y\leq(y-\hat{\theta}(x))\left(\frac{\sum_{k=1}^{K}w_{k}(x)(\hat{\theta}_{k}-\hat{\theta}(x))^{2}}{(K-1)\sum_{k=1}^{K}w_{k}(x)}\right)^{-\frac{1}{2}}\right\\}f_{\hat{\tau}^{2}}(x|\hat{\tau}_{M}^{2})dx$
$\displaystyle=\int_{-\infty}^{\infty}F_{Y}\left((y-\hat{\theta}(x))\left(\frac{\sum_{k=1}^{K}w_{k}(x)(\hat{\theta}_{k}-\hat{\theta}(x))^{2}}{(K-1)\sum_{k=1}^{K}w_{k}(x)}\right)^{-\frac{1}{2}}\right)f_{\hat{\tau}^{2}}(x|\hat{\tau}_{M}^{2})dx,$

where $f_{Y}$ is the probability density function of the $t$-distribution with
$K-1$ degrees of freedom. Thus, we can construct $100(1-\alpha)\%$ CI as

$\displaystyle\left(F_{\hat{\theta}}^{-1}\left(\frac{\alpha}{2}|\hat{\theta}_{1},\dots,\hat{\theta}_{K},\hat{\tau}_{M}^{2}\right),F_{\hat{\theta}}^{-1}\left(1-\frac{\alpha}{2}|\hat{\theta}_{1},\dots,\hat{\theta}_{K},\hat{\tau}_{M}^{2}\right)\right).$

The one- and two-sided p-values for the null hypothesis
$H_{0}:\theta=\theta_{0}$ can also be derived as
$1-F_{\hat{\theta}}(|\hat{\theta}(\hat{\tau}^{2})-\theta_{0}||\hat{\theta}_{1},\dots,\hat{\theta}_{K},\hat{\tau}_{M}^{2})$
and
$2(1-F_{\hat{\theta}}(|\hat{\theta}(\hat{\tau}^{2})-\theta_{0}||\hat{\theta}_{1},\dots,\hat{\theta}_{K},\hat{\tau}_{M}^{2}))$,
respectively.

We next expand an exact inference proposed by Michael et al. (2019) using the
M-estimator approach. Michael’s method first computes the lower and upper
bound of the between-study variance $[\tau_{min}^{2},\tau_{max}^{2}]$ via
interval estimation for the between-study variance that satisfied
$\Pr\\{\tau_{min}^{2}\leq\tau^{2}\leq\tau_{max}^{2}\\}\geq 1-\beta$. The
coverage probability of the CI for $\theta$ is $1-\alpha-\beta\approx
1-\alpha$ for the small $\beta$. However, the CI of the between-study variance
cannot sometimes be calculated. We consider applying the CI of the between-
study variance based on the M-estimator to the lower and upper bounds of the
between-study variance. The lower and upper bounds of between-study variance
are given as

$\displaystyle[\tau_{min}^{2},\tau_{max}^{2}]=\left[0,F^{-1}_{\hat{\tau}^{2}}\left(1-\beta|\hat{\tau}_{M}^{2}\right)\right].$

Since the CI based on the M-estimator can be computed at any
$\\{(\hat{\theta}_{1},\sigma_{1}^{2}),\dots,(\hat{\theta}_{K},\sigma_{K}^{2})\\}$,
the CI of the overall treatment effect by Michael’s method can be extended to
calculate at any study treatment effects and within-study variances without
prior setting of the lower and upper bound of between-study variance.

##  4 Numerical Results

We evaluate existing and our proposed methods via simulation and re-analysis
of real data to compare existing methods. First, we introduce a setup for our
simulation. Next, we show the result of our simulation of the CI for the
overall treatment effect. We re-analyse a meta-analysis with a small number of
studies.

###  4.1 Simulation Setup

We compare the CIs of the overall treatment effect by the DL method, the SJ
method, and our proposed M-estimator approach for the DL method (mDL) and SJ
method (mSJ). The DL and SJ methods are applied by metafor in R (Viechtbauer,
2010). The conservative approach by Michael et al. (2019) (Mi) and the
expansion by M-estimator approach (mMi) are also compared in our simulation by
rma.exact in R. We show the result of the HKSJ approach in Supplementary
materials A.5 to compare the overall treatment effect of each method.

We set some parameters for our simulations from IntHout et al. (2014) because
they compare meta-analysis methods under a few studies. The number of studies
is set from small to medium $(K=2,3,\cdots,10,15,20)$. We apply four within-
study variance scenarios as follows:

  * 1)

equal size study $(\sigma_{1}^{2}=\cdots=\sigma_{K}^{2}=1)$,

  * 2)

one small study, 1/10th of other studies
$(\sigma_{1}^{3}=\cdots=\sigma_{K-1}^{2}=0.1,\sigma_{K}^{2}=1)$,

  * 3)

50/50 small and large studies
$(\sigma_{1}^{2}=\cdots=\sigma_{K/2}^{2}=0.1,\sigma_{K/2+1}^{2}=\cdots=\sigma_{K}^{2}=1)$,

  * 4)

various studies $(\sigma_{k}^{2}=0.1+0.9\frac{k-1}{K-1}$, $k=1,\dots,K$).

The between-study heterogeneity is considered small to large
$(I^{2}=0.25,0.50,0.90)$. The between-study variance can be calculated from
within-study variance and between-study heterogeneity,
$\tau^{2}=K^{-1}\sum_{k=1}^{K}\sigma_{k}^{2}\frac{I^{2}}{1-I^{2}}$. The
overall treatment effect is set to $\theta=0$, and 10,000 evaluations are
performed for each parameter setting.

###  4.2 Simulation Results

The coverage probabilities and lengths of the CI of $\theta$ by comparison and
our proposed methods are shown in Figures 1 and 2, respectively. The
M-estimator approach (mDL, mSJ) always had a larger coverage probability than
the DL and SJ methods. The existence of a larger between-study heterogeneity,
led to a lower coverage probability under a small number of studies. It was
confirmed that the coverage probability reached a nominal level as the number
of studies increased, but convergence was slower when the heterogeneity was
larger. In the equal size study scenario, all methods can analyse the nominal
level of the coverage probability if the between-study heterogeneity is small.
However, in the other scenarios, the CIs of the DL and SJ methods narrowed,
and the coverage probabilities decreased below the nominal level. The Mi
method is also conservative, but the CI for the overall treatment effect
cannot be calculated when the search range for the between-study variance
cannot be calculated. The mMi method is as conservative as the Mi method, and
has the further advantage that it can always compute CI without computational
failures. We illustrate the number of computational failures in Figure 3 of
Supplementary material A.6. In particular, the calculation of the Mi method
becomes infeasible when the number of studies is small. Because the CI of the
treatment effect changes depending on the setting of the search range of
between-study variance, a rationale is needed for the prior setting of the
search range. The SJ method with M-estimator was the most conservative of the
comparable methods. We consider the reason to be to update the estimate of
between-study variance when estimating location parameters.

Figure 1: Coverage probabilities of the overall treatment effect by the DL,
SJ, Mi, and proposed mDL, mSJ, and mMi methods. The rows are within-study
variance settings, and the columns are between-study heterogeneity
$I^{2}=0.25,0.5,0.9$. Figure 2: Lengths of the overall treatment effect by the
DL, SJ, Mi, and proposed mDL, mSJ, and mMi methods. The rows are within-study
variance settings, and the columns are between-study heterogeneity
$I^{2}=0.25,0.5,0.9$.

###  4.3 Re-Analysis of Data from Small Meta-Analysis

We re-analyse the occurrence of a particular genetic variant, CCR5, in
patients with JIA compared to the general population (Hinks et al., 2010). The
investigation included three studies examining the association of JIA with a
biomarker via odds ratios. We apply the DL and SJ methods, and our M-estimator
approach to log odds ratios of the treatment effect from the JIA data. The
HKSJ approach with each method and Michael’s method are also applied to the
JIA data. The meta-analysis results are shown in Table 1. The between-study
variance by DL and SJ methods was estimated to be zero for these data.
Regarding a typical approach, the M-estimator approach had a wider CI than
typical DL and SJ methods because of the difference in handling between-study
variance. Although the HKSJ approach followed the same trend, the SJ method,
modified from the estimation of the between-study variance, was more affected
by the proposed M-estimator than the DL method. The Mi and mMi methods are the
conservative results as the CI includes zero. In this example, the M-estimator
approach allows analysis that considers the distribution of between-study
variance and avoids excessively conservative results.

Table 1: Means and $95\%$ CIs of the overall treatment effect for the JIA data. SE: standard error. CI: confidence interval. OR: odds ratio. HKSJ: Hartang-Knapp-Sidik-Jonkman approach, HKSJ-M: HKSJ approach with M-estimator. Mi: the method proposed by Michael et al. (2019). mMi: Mi method with M-estimator. JIA data |  |  |
---|---|---|---
Study | $\log\mathrm{OR}$ | SE |  $95\%$ CI
Prahalad et al. (2006) | -0.13 | 0.10 | (-0.33, 0.08)
Hinks et al. (2006) | -0.24 | 0.09 | (-0.41, -0.06)
Lindner et al. (2007) | -0.20 | 0.14 | (-0.47, 0.07)
Meta-analysis | DerSimonian-Laird | Sidik-Jonkman
Typical | -0.195 (-0.314, -0.075) | -0.195 (-0.317, -0.072)
M-estimator | -0.195 (-0.343, -0.046) | -0.195 (-0.347, -0.042)
HKSJ | -0.195 (-0.337, -0.053) | -0.195 (-0.320, -0.070)
HKSJ-M | -0.195 (-0.341, -0.046) | -0.195 (-0.342, -0.046)
Mi method | -0.195 (-0.446, 0.056) |
mMi method | -0.195 (-0.397, 0.007) |

##  5 Discussion

We proposed an estimation of the overall treatment effect that considers the
distribution of estimated between-study variance. We demonstrated the
distribution of between-study variance using general moment estimation and SJ
estimation. The distribution of between-study variance includes an unknown
location parameter. At the same time, we proposed the optimal estimated value
for the location parameter when the between-study variance estimated by each
method becomes the M-estimator. This approach enabled us to estimate the
overall treatment effect by considering the variability of the estimated
quantity even when unknown parameters were included in the distribution of the
estimated between-study variance.

The SJ method usually inserts another estimate for the unknown parameters
included in the estimate of the between-study variance. Still, it ignores the
variability of the inserted estimate. Considering the between-study variance
as an M-estimator, the variability is no longer present, and we can
concentrate on estimating its distribution.

The simulation study showed that our M-estimator approach can estimate more
conservatively than the corresponding typical method. In particular, the DL
and SJ methods with an M-estimator provide more accurate estimates than those
without the M-estimator. The application of the HKSJ approach gave similar
conservative results to typical methods with an M-estimator, but it is too
conservative in situations with one small study case. If the true distribution
of the test statistic is used, a more accurate estimation of the HKSJ approach
with the M-estimator may be possible. However, the actual distribution of the
HKSJ approach’s test statistic has yet to be discovered. The Mi method can
construct a conservative CI, but it cannot be computed in rare cases. The mMi
method extends the Mi method to calculate the CI at any study treatment
effects and within-study variances without prior setting the lower and upper
bounds of between-study variance. From the actual meta-analysis, for JIA data
comparison with a small number of studies, the typical and HKSJ approaches
with an M-estimator had longer CIs than those without an M-estimator despite a
between-study variance of zero.

This study has several limitations. While the distribution of between-study
variance was shown for general moment estimation, it was not provided for
restricted maximum likelihood estimation. This is because the latter
necessitates iterative computations, making it challenging to illustrate
explicitly the distribution of the between-study variance. Additionally, the
convex function for the M-estimator utilized only absolute value functions,
which are considered the most bias-robust estimators. It may be possible to
find more optimal optimization functions depending on the distribution of the
applied estimators.

In summary, we estimated the overall treatment effect considering the
distribution of the between-study variance by utilizing the optimal location
parameter. When existing methods based on asymptotic theory are not applicable
because of having only a few studies, employing our proposed approach for to
estimate the overall treatment effect enables a more accurate estimation. In
such circumstances, we recommend using the M-estimator approach.

## Acknowledgements

This work was supported by JSPS KAKENHI (Grant number JP24K14853).

## Supplementary Materials

Technical details of mathematical results and results from simulation studies
are available. Our M-estimator approach proposed in Section 2 and the random
sampling function of DL between-study variance described in Section 3 can be
installed from GitHub
(https://github.com/keisuke-hanada/metaMest/) as the metaMest package.

## References

  * Box (1954) Box, G. E. (1954).  Some theorems on quadratic forms applied in the study of analysis of variance problems, I. Effect of inequality of variance in the one-way classification.  The Annals of Mathematical Statistics pages 290–302.
  * DerSimonian and Laird (1986) DerSimonian, R. and Laird, N. (1986).  Meta-analysis in clinical trials.  Controlled Clinical Trials 7, 177–188.
  * Hanada and Sugimoto (2023) Hanada, K. and Sugimoto, T. (2023).  Inference using an exact distribution of test statistic for random-effects meta-analysis.  Annals of the Institute of Statistical Mathematics 75, 281–302.
  * Hartung (1999) Hartung, J. (1999).  An alternative method for meta-analysis.  Biometrical Journal: Journal of Mathematical Methods in Biosciences 41, 901–916.
  * Hartung and Knapp (2001a) Hartung, J. and Knapp, G. (2001a).  On tests of the overall treatment effect in meta-analysis with normally distributed responses.  Statistics in medicine 20, 1771–1782.
  * Hartung and Knapp (2001b) Hartung, J. and Knapp, G. (2001b).  A refined method for the meta-analysis of controlled clinical trials with binary outcome.  Statistics in medicine 20, 3875–3889.
  * Hinks et al. (2010) Hinks, A., Martin, P., Flynn, E., Eyre, S., Packham, J., Barton, A., Worthington, J., and Thomson, W. (2010).  Association of the ccr5 gene with juvenile idiopathic arthritis.  Genes & Immunity 11, 584–589.
  * Hinks et al. (2006) Hinks, A., Worthington, J., and Thomson, W. (2006).  The association of ptpn22 with rheumatoid arthritis and juvenile idiopathic arthritis.  Rheumatology 45, 365–368.
  * Huber (1992) Huber, P. J. (1992).  Robust estimation of a location parameter.  In Breakthroughs in statistics: Methodology and distribution, pages 492–518. Springer.
  * IntHout et al. (2014) IntHout, J., Ioannidis, J. P., and Borm, G. F. (2014).  The hartung-knapp-sidik-jonkman method for random effects meta-analysis is straightforward and considerably outperforms the standard dersimonian-laird method.  BMC medical research methodology 14, 1–12.
  * Li et al. (1994) Li, Y., Shi, L., and Daniel Roth, H. (1994).  The bias of the commonly-used estimate of variance in meta-analysis.  Communications in statistics-theory and methods 23, 1063–1085.
  * Lindner et al. (2007) Lindner, E., Nordang, G. B., Melum, E., Flatø, B., Selvaag, A. M., Thorsby, E., Kvien, T. K., Førre, Ø. T., and Lie, B. A. (2007).  Lack of association between the chemokine receptor 5 polymorphism ccr5delta32 in rheumatoid arthritis and juvenile idiopathic arthritis.  BMC medical genetics 8, 1–5.
  * Michael et al. (2019) Michael, H., Thornton, S., Xie, M., and Tian, L. (2019).  Exact inference on the random-effects model for meta-analyses with few studies.  Biometrics 75, 485–493.
  * Prahalad et al. (2006) Prahalad, S., Bohnsack, J. F., Jorde, L. B., Whiting, A., Clifford, B., Dunn, D., Weiss, R., Moroldo, M., Thompson, S. D., Glass, D. N., et al. (2006).  Association of two functional polymorphisms in the ccr5 gene with juvenile rheumatoid arthritis.  Genes & Immunity 7, 468–475.
  * Raudenbush (2009) Raudenbush, S. W. (2009).  Analyzing effect sizes: Random-effects models.  The handbook of Research Synthesis and Meta-analysis 2, 295–316.
  * Rousseeuw (1982) Rousseeuw, P. J. (1982).  Most robust m-estimators in the infinitesimal sense.  Zeitschrift für Wahrscheinlichkeitstheorie und verwandte Gebiete 61, 541–551.
  * Röver et al. (2015) Röver, C., Knapp, G., and Friede, T. (2015).  Hartung-knapp-sidik-jonkman approach and its modification for random-effects meta-analysis with few studies.  BMC medical research methodology 15, 1–7.
  * Sidik and Jonkman (2002) Sidik, K. and Jonkman, J. N. (2002).  A simple confidence interval for meta-analysis.  Statistics in medicine 21, 3153–3159.
  * Sidik and Jonkman (2005) Sidik, K. and Jonkman, J. N. (2005).  Simple heterogeneity variance estimation for meta-analysis.  Journal of the Royal Statistical Society Series C: Applied Statistics 54, 367–384.
  * Sidik and Jonkman (2007) Sidik, K. and Jonkman, J. N. (2007).  A comparison of heterogeneity variance estimators in combining results of studies.  Statistics in Medicine 26, 1964–1981.
  * Sidik and Jonkman (2022) Sidik, K. and Jonkman, J. N. (2022).  Quantifying uncertainty in method of moments estimates of the heterogeneity variance in random effects meta-analysis.  Biometrical Journal 64, 598–616.
  * Takeuchi (2020) Takeuchi, K. (2020).  Robust Estimation and Robust Parameter, pages 89–101.  Springer Japan, Tokyo.
  * Veroniki et al. (2016) Veroniki, A. A., Jackson, D., Viechtbauer, W., Bender, R., Bowden, J., Knapp, G., Kuss, O., Higgins, J. P., Langan, D., and Salanti, G. (2016).  Methods to estimate the between-study variance and its uncertainty in meta-analysis.  Research synthesis methods 7, 55–79.
  * Viechtbauer (2005) Viechtbauer, W. (2005).  Bias and efficiency of meta-analytic variance estimators in the random-effects model.  Journal of Educational and Behavioral Statistics 30, 261–293.
  * Viechtbauer (2010) Viechtbauer, W. (2010).  Conducting meta-analyses in r with the metafor package.  Journal of statistical software 36, 1–48.

##  Appendix A Supplemental Materials

###  A.1 Proof of Theorem 1

We have the following identity equation,

$\displaystyle
F_{\hat{\theta}}(q|\hat{\tau}_{M}^{2})-F_{\hat{\theta}|\hat{\tau}^{2}}(q|\hat{\tau}^{2})$
$\displaystyle=\int_{0}^{\infty}\left\\{F_{\hat{\theta}|\hat{\tau}^{2}}(q|x)-F_{\hat{\theta}|\hat{\tau}^{2}}(q|\hat{\tau}^{2})\right\\}f_{\hat{\tau}^{2}}(x|\hat{\tau}_{M}^{2})dx.$

Now, we consider Taylor expansion of $F_{\hat{\theta}|\hat{\tau}^{2}}(q|x)$
around $x=\hat{\tau}^{2}$. Then, there exists a constant $c$ such that

$\displaystyle
F_{\hat{\theta}}(q|\hat{\tau}_{M}^{2})-F_{\hat{\theta}|\hat{\tau}^{2}}(q|\hat{\tau}^{2})$
$\displaystyle=\int_{0}^{\infty}\left.\frac{d}{dx}F_{\hat{\theta}|\hat{\tau}^{2}}(q|x)\right|_{x=c}(x-\hat{\tau}^{2})f_{\hat{\tau}^{2}}(x|\hat{\tau}_{M}^{2})dx$
$\displaystyle=\left.\frac{d}{dx}F_{\hat{\theta}|\hat{\tau}^{2}}(q|x)\right|_{x=c}\int_{0}^{\infty}(x-\hat{\tau}^{2})f_{\hat{\tau}^{2}}(x|\hat{\tau}_{M}^{2})dx$
$\displaystyle=\left.\frac{d}{dx}F_{\hat{\theta}|\hat{\tau}^{2}}(q|x)\right|_{x=c}\left(E[\hat{\tau}^{2}|\hat{\tau}_{M}^{2}]-\hat{\tau}^{2}\right).$

Because we have $\hat{\theta}|\hat{\tau}^{2}\sim
N\left(\theta,\left(\sum_{k=1}^{K}w_{k}(x)\right)^{-1}\right)$, we can get

$\displaystyle\frac{d}{dx}F_{\hat{\theta}|\hat{\tau}^{2}}(q|x)$
$\displaystyle=\frac{d}{dx}\int_{-\infty}^{q}\frac{1}{\sqrt{2\pi(\sum_{k=1}^{K}w_{k}(x))^{-1}}}\exp\left\\{-\frac{(y-\theta)^{2}}{2(\sum_{k=1}^{K}w_{k}(x))^{-1}}\right\\}dy$
$\displaystyle=\int_{-\infty}^{q}\frac{d}{dx}\frac{1}{\sqrt{2\pi(\sum_{k=1}^{K}w_{k}(x))^{-1}}}\exp\left\\{-\frac{(y-\theta)^{2}}{2(\sum_{k=1}^{K}w_{k}(x))^{-1}}\right\\}dy$
$\displaystyle=\int_{-\infty}^{q}\left(1-(y-\theta)^{2}\sum_{k=1}^{K}w_{k}(x)\right)\frac{1}{\sqrt{2\pi(\sum_{k=1}^{K}w_{k}(x))^{-1}}}$
$\displaystyle\times\sum_{k=1}^{K}\frac{-w_{k}^{2}(x)}{2}\exp\left\\{-\frac{(y-\theta)^{2}}{2(\sum_{k=1}^{K}w_{k}(x))^{-1}}\right\\}dy$
$\displaystyle=-\frac{1}{2}\sum_{k=1}^{K}w_{k}^{2}(x)\int_{-\infty}^{(q-\theta)\sum_{k=1}^{K}w_{k}(x)}(1-z^{2})\frac{1}{\sqrt{2\pi}}e^{-\frac{z^{2}}{2}}dz.$
(9)

For any real number $c_{0}$, we have

$\displaystyle\int_{-\infty}^{c_{0}}(1-z^{2})\frac{1}{\sqrt{2\pi}}e^{-\frac{z^{2}}{2}}dz$
$\displaystyle=\int_{-\infty}^{\infty}(1-z^{2})\frac{1}{\sqrt{2\pi}}e^{-\frac{z^{2}}{2}}dz-\int_{c_{0}}^{\infty}(1-z^{2})\frac{1}{\sqrt{2\pi}}e^{-\frac{z^{2}}{2}}dz$
$\displaystyle=-\int_{c_{0}}^{\infty}(1-z^{2})\frac{1}{\sqrt{2\pi}}e^{-\frac{z^{2}}{2}}dz.$
(10)

Thus, we can obtain

$\displaystyle\int_{-\infty}^{0}(1-z^{2})\frac{1}{\sqrt{2\pi}}e^{-\frac{z^{2}}{2}}dz$
$\displaystyle=\frac{1}{2}\int_{-\infty}^{\infty}(1-z^{2})\frac{1}{\sqrt{2\pi}}e^{-\frac{z^{2}}{2}}dz=0.$

When $q>\theta$, we have

$\displaystyle\int_{-\infty}^{(q-\theta)\sum_{k=1}^{K}w_{k}(x)}(1-z^{2})\frac{1}{\sqrt{2\pi}}e^{-\frac{z^{2}}{2}}dz$
$\displaystyle=\int_{-\infty}^{0}(1-z^{2})\frac{1}{\sqrt{2\pi}}e^{-\frac{z^{2}}{2}}dz$
$\displaystyle+\int_{0}^{(q-\theta)\sum_{k=1}^{K}w_{k}(x)}(1-z^{2})\frac{1}{\sqrt{2\pi}}e^{-\frac{z^{2}}{2}}dz$
$\displaystyle=\int_{0}^{1}(1-z^{2})\frac{1}{\sqrt{2\pi}}e^{-\frac{z^{2}}{2}}dz$
$\displaystyle+\int_{1}^{(q-\theta)\sum_{k=1}^{K}w_{k}(x)}(1-z^{2})\frac{1}{\sqrt{2\pi}}e^{-\frac{z^{2}}{2}}dz$
$\displaystyle>\int_{0}^{1}(1-z^{2})\frac{1}{\sqrt{2\pi}}e^{-\frac{z^{2}}{2}}dz+\int_{1}^{\infty}(1-z^{2})\frac{1}{\sqrt{2\pi}}e^{-\frac{z^{2}}{2}}dz$
$\displaystyle=0.$

Hence, we can obtain as
$\frac{d}{dx}F_{\hat{\theta}|\hat{\tau}^{2}}(q|x)|_{x=c}<0~{}(q>\theta)$ from
equation (A.1). When $q<\theta$, we have
$\frac{d}{dx}F_{\hat{\theta}|\hat{\tau}^{2}}(q|x)|_{x=c}<0$ because of
equation (A.1). As a result, we obtain

$\displaystyle\left.\frac{d}{dx}F_{\hat{\theta}|\hat{\tau}^{2}}(q|x)\right|_{x=c}$
$\displaystyle<0\,\,\,\,\,(q>\theta)$
$\displaystyle\left.\frac{d}{dx}F_{\hat{\theta}|\hat{\tau}^{2}}(q|x)\right|_{x=c}$
$\displaystyle>0\,\,\,\,\,(q<\theta)$
$\displaystyle\left.\frac{d}{dx}F_{\hat{\theta}|\hat{\tau}^{2}}(q|x)\right|_{x=c}$
$\displaystyle=0\,\,\,\,\,(q=\theta).$

Thus, if $E[\hat{\tau}^{2}|\hat{\tau}_{M}^{2}]\leq\hat{\tau}^{2}$, the
$100(1-\alpha)\%$ CI of $\theta$ with $\hat{\tau}^{2}_{M}$ is more
conservative than the one with $\hat{\tau}^{2}$. From Proposition 2 in the
next Appendix A.2, if $\hat{\tau}^{2}\to_{p}\tau^{*2}$, then
$\hat{\tau}_{M}^{2}\to_{p}\tau^{*2}$. Thus, for any real number $q$, we have

$\displaystyle|F_{\hat{\theta}}(q|\hat{\tau}_{M}^{2})-F_{\hat{\theta}}(q|\tau^{*2})|$
$\displaystyle=\left|\left\\{\left.\frac{dF_{\hat{\theta}}}{d\tau^{2}}(q|\tau^{2})\right|_{\tau^{2}=\tau^{2}_{c}}\right\\}(\hat{\tau}_{M}^{2}-\tau^{*2})\right|$
$\displaystyle\leq\left|\left\\{\left.\frac{dF_{\hat{\theta}}}{d\tau^{2}}(q|\tau^{2})\right|_{\tau^{2}=\tau^{2}_{c}}\right\\}\right|\left|\hat{\tau}_{M}^{2}-\tau^{*2}\right|\to_{p}0\,\,\,(K\to\infty).$

Therefore, the coverage probability converges of the $100(1-\alpha)\%$ CI of
$\theta$ to its true value if the number of studies is sufficiently large.

$\square$

###  A.2 Some Properties of M-estimator

####  A.2.1 Proposition 1

Under the random-effect model (1) in the main manuscript, the variance of
M-estimator $\hat{\tau}^{2}_{M}$ is given as

$\displaystyle V[\hat{\tau}_{M}^{2}]$
$\displaystyle=\left(\frac{1}{r(\tau^{*2})}\right)^{2}V[\hat{\tau}^{2}|\tau^{2}=\tau^{*2}]+O_{p}(a_{k})\,\,\,\,\,(K\to\infty),$

where $r(x)=dR/d\tau^{2}(x)$, $a_{K}=r(c_{K})-r(\tau^{*2})$, $c_{K}$ is a
positive constant such that
$|c_{K}-\tau^{*2}|<|\hat{\tau}_{M}^{2}-\tau^{*2}|$, and $O_{p}(a_{k})$ denotes
the order less than or equal to that of $a_{K}$ in probability. Thus, the
variance estimator of $\hat{\tau}_{M}^{2}$ is

$\displaystyle\hat{V}[\hat{\tau}_{M}^{2}]$
$\displaystyle=\left(\frac{1}{r(\hat{\tau}_{M}^{2})}\right)^{2}V[\hat{\tau}^{2}|\tau^{2}=\hat{\tau}_{M}^{2}].$

Regarding the proof of Proposition 1, we first consider a Taylor expansion of
$\hat{\tau}^{2}=R(\hat{\tau}_{M}^{2})$ around $\hat{\tau}_{M}^{2}=\tau^{*2}$.
In this case, there exists a positive constant $c_{K}>0$ such that

$\displaystyle\hat{\tau}^{2}$
$\displaystyle=R(\hat{\tau}_{M}^{2})=R(\tau^{*2})+r(c_{K})(\hat{\tau}_{M}^{2}-\tau^{*2}),$

where $|c_{K}-\tau^{*2}|<|\hat{\tau}_{M}^{2}-\tau^{*2}|$ and
$r(c)=dR/d\tau^{2}(\tau^{2})|_{\tau^{2}=c}$. Thus, the variance of
$\hat{\tau}_{M}^{2}$ is written as

$\displaystyle V[\hat{\tau}_{M}^{2}]$
$\displaystyle=\left(\frac{1}{r(c_{K})}\right)^{2}V[\hat{\tau}^{2}]$
$\displaystyle=\left(\frac{1}{r(\tau^{*2})}\right)^{2}V[\hat{\tau}^{2}]+\left(\frac{r(c_{K})-r(\tau^{*2})}{r(\tau^{*2})}\right)^{2}V[\hat{\tau}^{2}]$
$\displaystyle=\left(\frac{1}{r(\tau^{*2})}\right)^{2}V[\hat{\tau}^{2}]+\left(\frac{\sqrt{V[\hat{\tau}^{2}]}}{r(\tau^{*2})}\right)^{2}a_{K}^{2}.$

Now, we have $r(x)>0$ because $R(x)$ is a monotonic increasing function,
$V[\hat{\tau}^{2}]<\infty$, and $c_{K}\to_{p}\tau^{*2}$, and then the second
term on the right-hand side can be calculated as follows

$\displaystyle\left(\frac{\sqrt{V[\hat{\tau}^{2}]}}{r(\tau^{*2})}\right)^{2}a_{K}^{2}$
$\displaystyle\leq Aa_{K}^{2}\to_{p}0\,\,\,(K\to 0),$

where $A(<\infty)$ is a positive constant. Hence, we can obtain

$\displaystyle V[\hat{\tau}_{M}^{2}]$
$\displaystyle=\left(\frac{1}{r(\tau^{*2})}\right)^{2}V[\hat{\tau}^{2}]+O_{p}(a_{K}).$

$\square$

####  A.2.2 Proposition 2

Under the random-effect model (1) in the main manuscript, if the estimator
$\hat{\tau}^{2}$ has consistency, then the M-estimator $\hat{\tau}_{M}^{2}$ is
also consistent. Thus,
$\hat{\tau}^{2}\to_{p}\tau^{*2}\,\,\,(K\to\infty)\Rightarrow\hat{\tau}^{2}_{M}\to_{p}\tau^{*2}\,\,\,(K\to\infty)$.
Regarding the proof of Proposition 2, if we assume
$\hat{\tau}^{2}\to_{p}\tau^{*2}$, then for any $\varepsilon>0$, we have

$\displaystyle\Pr\\{|\hat{\tau}^{2}-\tau^{*2}|>\varepsilon\\}\to
0\,\,\,(K\to\infty).$

Now, we evaluate the difference between $\hat{\tau}_{M}^{2}$ and $\tau^{*2}$
as

$\displaystyle|\hat{\tau}_{M}^{2}-\tau^{*2}|$
$\displaystyle\leq\left|R^{-1}(\hat{\tau}^{2})-R^{-1}(\tau^{*2})\right|+\left|R^{-1}(\tau^{*2})-\tau^{*2}\right|.$

The inverse of the estimation function $R^{-1}(\hat{\tau}^{2})$ can be
evaluated as

$\displaystyle|R^{-1}(\hat{\tau}^{2})-R^{-1}(\tau^{*2})|$
$\displaystyle=\left|\left\\{\left.\frac{dR^{-1}}{d\tau^{2}}(\tau^{2})\right|_{\tau^{2}=c}\right\\}(\hat{\tau}^{2}-\tau^{*2})\right|$
$\displaystyle\leq
M\left|\hat{\tau}^{2}-\tau^{*2}\right|\to_{p}0\,\,\,(K\to\infty),$

where $M(<\infty)$ is a positive constant. Thus, we have
$R^{-1}(\hat{\tau}^{2})\to_{p}R^{-1}(\tau^{*2})\,\,\,(K\to\infty)$. From the
consistency of $\hat{\tau}^{2}\to_{p}\tau^{*2}$ and the dominated convergence
theorem, we get

$\displaystyle\lim_{K\to\infty}R(\tau^{*2})$
$\displaystyle=\lim_{K\to\infty}\underset{t\geq
0}{\mathrm{argmin}}~{}\int_{0}^{\infty}\rho(x-t)f_{\hat{\tau}^{2}}(x|\tau^{*2})dx$
$\displaystyle=\underset{t\geq
0}{\mathrm{argmin}}~{}\rho(\tau^{*2}-t)=\tau^{*2}.$

The above equation implies that
$\lim_{K\to\infty}R^{-1}(\tau^{*2})=\tau^{*2}$. Thus, the first inequality is
convergence, $|\hat{\tau}_{M}^{2}-\tau^{*2}|\to_{p}0\,\,\,(K\to\infty)$, and
then $\hat{\tau}_{M}^{2}\to_{p}\tau^{*2}$.

$\square$

###  A.3 Proof of Theorem 2: Distribution of General Moment Estimator

From the theorem of Box (1954), $Q_{a}$ is distributed as
$\sum_{r=1}^{R}\lambda_{r}\chi_{r(1)}^{2}$, where
$\lambda_{1},\dots,\lambda_{R}$ are non-negative eigenvalues of matrix $VW$,
$V$ is the variance-covariance matrix of
$(\hat{\theta}_{1}-\hat{\theta}_{a},\dots,\hat{\theta}_{K}-\hat{\theta}_{a})$,
$W$ is a diagonal matrix with $a_{1},\dots,a_{K}$, and
$\chi_{1(1)}^{2},\dots,\chi_{R(1)}^{2}$ are independently and identically
distributed as chi-square distributions with 1 degree of freedom.

$\displaystyle E[(\hat{\theta}_{k}-\hat{\theta}_{a})^{2}]$
$\displaystyle=E\left[\hat{\theta}_{k}^{2}-2\hat{\theta}_{k}\hat{\theta}_{a}+\hat{\theta}_{a}^{2}\right]$
$\displaystyle=\left(\theta^{2}+(\sigma_{k}^{2}+\tau^{2})\right)-2\left(\theta^{2}+\frac{a_{k}(\sigma_{k}^{2}+\tau^{2})}{A_{1}^{2}}\right)+\left(\theta^{2}+\frac{1}{A_{1}^{2}}\sum_{i=1}^{K}a_{i}^{2}(\sigma_{i}^{2}+\tau^{2})\right)$
$\displaystyle=(\sigma_{k}^{2}+\tau^{2})-2\frac{a_{k}(\sigma_{k}^{2}+\tau^{2})}{A_{1}^{2}}+\frac{1}{A_{1}^{2}}\sum_{i=1}^{K}a_{i}^{2}(\sigma_{i}^{2}+\tau^{2})$
$\displaystyle=\left(\sigma_{k}^{2}+\frac{1}{A_{1}^{2}}\sum_{i=1}^{K}a_{i}^{2}\sigma_{i}^{2}-\frac{2a_{k}\sigma_{k}^{2}}{A_{1}}\right)+\left(1+\frac{A_{2}}{A_{1}^{2}}-\frac{2a_{k}}{A_{1}}\right)\tau^{2},$

where $A_{r}=\sum_{i=1}^{K}a_{i}^{r}$. We assume $k\neq k^{\prime}$,

$\displaystyle
E[(\hat{\theta}_{k}-\hat{\theta}_{a})(\hat{\theta}_{k^{\prime}}-\hat{\theta}_{a})]$
$\displaystyle=E\left[\hat{\theta}_{k}\hat{\theta}_{k^{\prime}}-\hat{\theta}_{k}\hat{\theta}_{a}-\hat{\theta}_{k^{\prime}}\hat{\theta}_{a}+\hat{\theta}_{a}^{2}\right]$
$\displaystyle=\theta^{2}-\left(\theta^{2}+\frac{a_{k}(\sigma_{k}^{2}+\tau^{2})}{A_{1}^{2}}\right)-\left(\theta^{2}+\frac{a_{k^{\prime}}(\sigma_{k^{\prime}}^{2}+\tau^{2})}{A_{1}^{2}}\right)$
$\displaystyle+\left(\theta^{2}+\frac{1}{A_{1}^{2}}\sum_{i=1}^{K}a_{i}^{2}(\sigma_{i}^{2}+\tau^{2})\right)$
$\displaystyle=-\frac{a_{k}(\sigma_{k}^{2}+\tau^{2})}{A_{1}^{2}}-\frac{a_{k^{\prime}}(\sigma_{k^{\prime}}^{2}+\tau^{2})}{A_{1}^{2}}+\frac{1}{A_{1}^{2}}\sum_{i=1}^{K}a_{i}^{2}(\sigma_{i}^{2}+\tau^{2})$
$\displaystyle=\left(-\frac{a_{k}\sigma_{k}^{2}+a_{k^{\prime}}\sigma_{k^{\prime}}^{2}}{A_{1}}+\frac{1}{A_{1}^{2}}\sum_{i=1}^{K}a_{i}^{2}\sigma_{i}^{2}\right)+\left(\frac{A_{2}}{A_{1}^{2}}-\frac{a_{k}+a_{k^{\prime}}}{A_{1}}\right)\tau^{2}.$

Thus, we can calculate the $(k,k^{\prime})$ element of $V$ as equation (6) in
the main manuscript. The distribution function of the general moment estimator
of between-study variance can be given as

$\displaystyle F_{\hat{\tau}_{G,u}^{2}}(x|\tau^{2})$
$\displaystyle=\Pr\left\\{\hat{\tau}^{2}\leq x\right\\}$
$\displaystyle=\Pr\left\\{\frac{Q_{a}-\left(\sum_{k=1}^{K}a_{k}\sigma_{k}^{2}-A_{1}^{-1}\sum_{k=1}^{K}a_{k}^{2}\sigma_{k}^{2}\right)}{A_{1}-A_{2}/A_{1}}\leq
x\right\\}$
$\displaystyle=\Pr\left\\{Q_{a}\leq(A_{1}-A_{2}/A_{1})x+\left(\sum_{k=1}^{K}a_{k}\sigma_{k}^{2}-A_{1}^{-1}\sum_{k=1}^{K}a_{k}^{2}\sigma_{k}^{2}\right)\right\\}$
$\displaystyle=F_{Q_{a}}\left((A_{1}-A_{2}/A_{1})x+\left(\sum_{k=1}^{K}a_{k}\sigma_{k}^{2}-A_{1}^{-1}\sum_{k=1}^{K}a_{k}^{2}\sigma_{k}^{2}\right)\right).$

Therefore, the density function of $\hat{\tau}_{G,u}^{2}$ can be given as

$\displaystyle f_{\hat{\tau}_{G,u}^{2}}(x)$
$\displaystyle=(A_{1}-A_{2}/A_{1})f_{Q_{a}}\left((A_{1}-A_{2}/A_{1})x+\left(\sum_{k=1}^{K}a_{k}\sigma_{k}^{2}-A_{1}^{-1}\sum_{k=1}^{K}a_{k}^{2}\sigma_{k}^{2}\right)\right),$

and $\hat{\tau}_{G,u}^{2}$ is distributed as

$\displaystyle\hat{\tau}_{G,u}\sim\frac{\sum_{r=1}^{R}\lambda_{r}\chi_{r(1)}^{2}-\left(\sum_{k=1}^{K}a_{k}\sigma_{k}^{2}-A_{1}^{-1}\sum_{k=1}^{K}a_{k}^{2}\sigma_{k}^{2}\right)}{A_{1}-A_{2}/A_{1}}.$

$\square$

###  A.4 Proof of Theorem 3: Distribution of SJ Estimator

We show that the SJ estimator $\hat{\tau}_{SJ}$ is distributed as
$\frac{\tau^{2}}{K-1}\chi_{K-1}^{2}$. From the theorem of Box (1954),
$\hat{\tau}_{SJ}^{2}$ is distributed as
$\sum_{r=1}^{R}\lambda_{r}\chi_{r(1)}^{2}$, where
$\lambda_{1},\dots,\lambda_{R}$ are non-negative eigenvalues of matrix $VW$,
$V$ is the variance-covariance matrix of
$(\hat{\theta}_{1}-\hat{\theta}_{v},\dots,\hat{\theta}_{K}-\hat{\theta}_{v})$,
$W$ is a diagonal matrix with
$(K-1)^{-1}v_{1}^{-1},\dots,(K-1)^{-1}v_{K}^{-1}$, and
$\chi_{1(1)}^{2},\dots,\chi_{R(1)}^{2}$ are independently and identically
distributed as chi-square distributions with 1 degree of freedom. Here, we can
calculate as below as well as the flow of proof of Theorem 2:

$\displaystyle E[(\hat{\theta}_{k}-\hat{\theta}_{v})^{2}]$
$\displaystyle=\sigma_{k}^{2}+\tau^{2}-\frac{1}{W_{1}(\tau^{2})}$
$\displaystyle
E[(\hat{\theta}_{k}-\hat{\theta}_{v})(\hat{\theta}_{k^{\prime}}-\hat{\theta}_{v})]$
$\displaystyle=-\frac{1}{W_{1}(\tau^{2})}.$

Thus,

$\displaystyle\det(\lambda-VW)$ $\displaystyle=\det\left(\lambda
I_{K}-\frac{\tau^{2}}{K-1}\left(I_{K}-\frac{1}{W_{1}(\tau^{2})}\left(\begin{array}[]{cccc}w_{1}(\tau^{2})&w_{2}(\tau^{2})&\cdots&w_{K}(\tau^{2})\\\
w_{1}(\tau^{2})&w_{2}(\tau^{2})&\cdots&w_{K}(\tau^{2})\\\
\vdots&\vdots&\vdots&\vdots\\\
w_{1}(\tau^{2})&w_{2}(\tau^{2})&\cdots&w_{K}(\tau^{2})\end{array}\right)\right)\right)$
$\displaystyle=\left|\begin{array}[]{cccc}\lambda-\frac{\tau^{2}}{K-1}\left(1-\frac{w_{1}(\tau^{2})}{W_{1}(\tau^{2})}\right)&-\frac{w_{2}(\tau^{2})}{W_{1}(\tau^{2})}&\cdots&-\frac{w_{K}(\tau^{2})}{W_{1}(\tau^{2})}\\\
-\frac{w_{1}(\tau^{2})}{W_{1}(\tau^{2})}&\lambda-\frac{\tau^{2}}{K-1}\left(1-\frac{w_{2}(\tau^{2})}{W_{1}(\tau^{2})}\right)&\cdots&-\frac{w_{K}(\tau^{2})}{W_{1}(\tau^{2})}\\\
\vdots&\vdots&\vdots&\vdots\\\
-\frac{w_{1}(\tau^{2})}{W_{1}(\tau^{2})}&-\frac{w_{2}(\tau^{2})}{W_{1}(\tau^{2})}&\cdots&\lambda-\frac{\tau^{2}}{K-1}\left(1-\frac{w_{K}(\tau^{2})}{W_{1}(\tau^{2})}\right)\end{array}\right|$
$\displaystyle=\lambda\left|\begin{array}[]{cccc}1&-\frac{w_{2}(\tau^{2})}{W_{1}(\tau^{2})}&\cdots&-\frac{w_{K}(\tau^{2})}{W_{1}(\tau^{2})}\\\
1&\lambda-\frac{\tau^{2}}{K-1}\left(1-\frac{w_{2}(\tau^{2})}{W_{1}(\tau^{2})}\right)&\cdots&-\frac{w_{K}(\tau^{2})}{W_{1}(\tau^{2})}\\\
\vdots&\vdots&\vdots&\vdots\\\
1&-\frac{w_{2}(\tau^{2})}{W_{1}(\tau^{2})}&\cdots&\lambda-\frac{\tau^{2}}{K-1}\left(1-\frac{w_{K}(\tau^{2})}{W_{1}(\tau^{2})}\right)\end{array}\right|$
$\displaystyle=\lambda\left|\begin{array}[]{cccc}1&-\frac{w_{2}(\tau^{2})}{W_{1}(\tau^{2})}&\cdots&-\frac{w_{K}(\tau^{2})}{W_{1}(\tau^{2})}\\\
0&\lambda-\frac{\tau^{2}}{K-1}&\cdots&0\\\ \vdots&\vdots&\vdots&\vdots\\\
0&0&\cdots&\lambda-\frac{\tau^{2}}{K-1}\end{array}\right|$
$\displaystyle=\lambda\left(\lambda-\frac{\tau^{2}}{K-1}\right)^{K-1}.$

Therefore, we can get $\lambda_{r}=\frac{\tau^{2}}{K-1}$ $(r=1,\dots,K-1)$,
and then the distribution of $\hat{\tau}_{SJ}^{2}$ can be written as

$\displaystyle\hat{\tau}_{SJ}^{2}$
$\displaystyle\sim\sum_{r=1}^{R}\lambda_{r}\chi_{r(1)}^{2}=\sum_{r=1}^{R}\frac{\tau^{2}}{K-1}\chi_{r(1)}^{2}=\frac{\tau^{2}}{K-1}\chi_{K-1}^{2}.$

$\square$

###  A.5 Simulation Result of HKSJ Approach

The coverage probabilities and lengths of the CI of $\theta$ by the HKSJ
approach with comparison methods and our proposed one are shown in Figure 3
and 4 in Appendix A.5. As with typical methods, our proposed method always had
a higher coverage probability than the HKSJ approach with the DL and SJ
estimators. However, the HKSJ approach with the proposed M-estimator is
sometimes too conservative in scenarios where the existing HKSJ approach
preserves nominal levels. When there was a small between-study heterogeneity,
the correction of the proposed M-estimator was quite effective. Still, when
there was a large between-study heterogeneity, the difference between the
existing method and the one with the proposed M-estimator became small.

Figure 3: Coverage probabilities of the overall treatment effect by the HKSJ
approach with the DL estimator (HKSJ-DL), the SJ estimator (HKSJ-SJ) and
proposed M-estimator method (HKSJ-mDL, HKSJ-mSJ). Rows are within-study
variance settings, and columns are between-study heterogeneity
$I^{2}=0.25,0.5,0.9$. Figure 4: Lengths of the overall treatment effect by the
HKSJ approach with the DL estimator (HKSJ-DL), the SJ estimator (HKSJ-SJ) and
proposed M-estimator method (HKSJ-mDL, HKSJ-mSJ). Rows are within-study
variance settings, and columns are between-study heterogeneity
$I^{2}=0.25,0.5,0.9$.

###  A.6 The CI by Michael’s Approach

The number of computational failures regarding the CI of overall treatment
effect by Michael et al. (2019) is illustrated in Figure 5 in Appendix A.6.
The CI of the overall treatment effect is feasible when the number of studies
is significantly large, but it becomes infeasible when the number of studies
is small. This occurs when the CI for the between-study variance cannot be
computed. The CI is used as the search range for the between-study variance.
Therefore, if the CI for the between-study variance cannot be calculated, the
CI for the overall treatment effect cannot be calculated either.

Figure 5: Number of computational failures regarding the CI of the overall
treatment effect by Michael et al. (2019) (Mi), and proposed M-estimator
method (mMi). Rows are within-study variance settings, and columns are
between-study heterogeneity $I^{2}=0.25,0.5,0.9$. The number of iterations is
$10,000$.