# Easy confidence interval formulas for network meta-analysis and adjustment
of confidence intervals for a small number of studies

Masahiro Kojima111Address: Biometrics Department, R&D Division, Kyowa Kirin
Co., Ltd. Otemachi Financial City Grand Cube, 1-9-2 Otemachi, Chiyoda-ku,
Tokyo, 100-004, Japan. Tel: +81-3-5205-7200 E-Mail:
<EMAIL_ADDRESS>Kyowa Kirin Co., Ltd  The Graduate
University for Advanced Studies

###### Abstract

We propose simple formulas of confidence intervals for the Wald statistic,
likelihood ratio statistic, and score statistic for a network meta-analysis.
In addition, we consider resolutions for concerns that network meta-analyses
with a small number of studies cannot hold a nominal confidence level. For a
bias adjustment in analyses with a small number of studies, a Bartlett-type
adjustment is a well-known method. Many Bartlett-type adjustment-type methods
are based on maximum likelihood estimators. However, the network meta-analysis
often uses the restricted maximum likelihood estimators that have not been
extensively discussed in Bartlett-type adjustment.
In this paper, we propose a Bartlett-type adjustment method for the Wald
statistic, likelihood ratio statistic, and score statistic when nuisance
parameters are estimated by not only the maximum likelihood method but also
the restricted maximum likelihood method. In addition, we propose a higher-
order adjustment by applying the bootstrap method to the Bartlett-type
adjusted statistics. Using a computer simulation, we confirmed that the
adjusted confidence intervals maintained a nominal confidence level. In
addition, we confirmed that the confidence interval of the likelihood ratio
test based on the restricted maximum likelihood method performs well without
further bootstrap adjustment. Finally, we demonstrated that confidence
intervals were adjusted for actual network meta-analysis.

Key words and phrases: Bartlett-type adjustments; Random effects meta-analysis
model; Parametric bootstrap

##  1 Introduction

Over the past few decades, many researchers have shown an interest in network
meta-analysis (NMA) which has been actively used in a wide range of fields
including ecology, economics, medical science, and medicine. As Gurevitch et
al.[1] indicated, studies using the meta-analysis began in the 1970s; in the
recent years, the number of studies has rapidly increased and currently
exceeds 200,000 papers. The number of papers using network meta-analyses has
also grown[2]. With the number of clinical trials increasing each year[3],
studies using the network meta-analysis seem likely to be more numerous in the
future.

NMA synthesizes treatment effects obtained from high-quality studies via
systematic reviews[4]. There may be discrepancies in the results of the
studies to be synthesized, and such discrepancies are called heterogeneity. To
address the problem of the heterogeneity, a random effect is used to adjust
heterogeneities between studies. In the NMA, the random effects meta-analysis
model[5] is often used[6, 7, 8].

Here, we consider the NMA with a small number of studies. The problem in such
cases is that the confidence intervals cannot hold a nominal confidence level
($100-{\alpha}$)$\%$. Consequently, the confidence intervals become wider or
narrower than the true confidence intervals, potentially resulting in
misleading decision-making. To resolve this problem with small numbers of
studies, the Bartlett-type adjustment method is often used. Because the
Bartlett-type adjustment can adjust the second-order bias in type I error[9,
10], the accuracy of confidence intervals is improved. Many Bartlett-type
adjustment methods are based on maximum likelihood estimators (MLE). However,
the NMA often uses the restricted maximum likelihood estimators (REMLE), which
have not been extensively discussed in this context. Noma et al.[11] proposed
the Bartlett-type adjustment using $m$ bootstrap samples for the NMA with a
small number of studies. However, there is a problem with the bootstrap
samples, i.e., the calculation load is very high because the MLE or REMLE is
derived by numerically solving a likelihood equation that is not closed-form.
In addition, because the likelihood ratio (LR) statistic and score statistic
on the REMLE are not closed-form, those confidence intervals also need to be
derived numerically. When using the bootstrap method, we do not know which
bias terms cause the bias. Kojima and Kubokawa[12] analytically derived
Bartlett-type adjustment terms in general linear models via an analytical
approach. However, they did not propose a Bartlett-type adjustment method
based on the restricted likelihood function. van Aert and Jackson[13] proposed
using t-quantiles instead of normal quantiles to improve the actual coverage
probability to be closer to the nominal level. However, their proposed method
cannot be applied the NMA straightforwardly.

In this paper, we propose novel Bartlett-type adjustment terms for confidence
intervals based on the Wald, LR, and score using the MLE and REMLE. Because
the Bartlett adjustment terms can be calculated analytically, our proposed
method can calculate the adjusted confidence intervals at 1/$m$ of the speed
of the bootstrap method. In addition, we propose the closed-form confidence
intervals of LR and score statistic based on a novel constrained likelihood
function under a null hypothesis. Thus, our method does not need to use
numerical analysis except for the derivation of the MLE. We focus on the
likelihood-based inference methods to derive the equations for adjustment.

We give further details of the Bartlett adjustment and Bartlett-type
adjustment. Here, ”Bartlett adjustment” refers to the adjustment only for the
likelihood ratio statistics. In general, the Bartlett adjustment terms can be
derived from the expectation of the LR statistic. By contrast, the ”Bartlett-
type adjustment” is an extension of the Bartlett adjustment that also supports
the adjustment of the Wald statistic and score statistic. In other words, the
Bartlett-type adjustment method is more generalized method. The Bartlett-type
adjustment terms can be derived from the asymptotic expansion of the
statistic.

In Section 2, we introduce the random effects meta-analysis model and the
three statistics. In addition, we provide the formulas for confidence
intervals. In Section 3, we obtain the Bartlett-type adjustment terms via the
analytical approach and the adjusted confidence intervals. In Section 4, we
propose the higher-order adjustment using the parametric bootstrap method. In
Section 5, we use a computer simulation to confirm that the adjusted
confidence intervals maintained a nominal confidence level. In Section 6, we
demonstrate that confidence intervals were adjusted for actual network meta-
analysis. Section 7 is the discussion. All technical issues are given in the
Appendix in the supplemental material.

##  2 Methods

##  3 Three statistics in the random effects meta-analysis model

###  3.1 Random effects meta-analysis model

We assume that the number of studies is $n$ and the number of treatment
effects of interest is $p$. The treatment effects $p$ are differences of
outcomes between treatment groups (e.g., any pairwise comparison) and a
reference treatment (e.g., a placebo). The effect measures are a mean
difference, log odds ratio difference, and log hazard ratio[14, 15, 16, 17].
In the NMA, because the every study may not include the all treatment groups,
we define that the number of outcomes of the $i$-th study is $p_{i}$
independent contrast estimates and $p_{i}\leq p$ satisfies. We assume $p$ is a
finite, $N=\sum^{n}_{i=1}p_{i}$, and $N>p$.

We consider the random effects meta-analysis model consisting of the within-
study model and between-studies model[4, 11],

$\displaystyle\mbox{(within-study model) }{\text{\boldmath$y$}}_{i}=$
$\displaystyle{\text{\boldmath$X$}}_{i}{\text{\boldmath$\beta$}}_{i}+{\text{\boldmath$e$}}_{i},{}$
$\displaystyle\mbox{(between-studies model) }{\text{\boldmath$\beta$}}_{i}=$
$\displaystyle{\text{\boldmath$\mu$}}+{\text{\boldmath$v$}}_{i},{}$

where ${\text{\boldmath$y$}}_{i}$ is a $p_{i}\times 1$ vector of a known
estimator of a treatment effect within the $i$-th study,
${\text{\boldmath$X$}}_{i}$ is a known $p_{i}\times p$ design matrix in which
the element of a treatment corresponding to $y_{ij}$ is 1 and the others are
0, ${\text{\boldmath$\beta$}}_{i}$ is a $p\times 1$ vector of a treatment
effect within the $i$-th study, ${\text{\boldmath$e$}}_{i}$ is a $p_{i}\times
1$ vector of an error within the $i$-th
study$,{\text{\boldmath$e$}}_{i}\sim{\cal
N}({\bf\text{\boldmath$0$}},{\text{\boldmath$R$}}_{i})$,
${\text{\boldmath$R$}}_{i}$ is a known $p_{i}\times p_{i}$ covariance matrix
of the $i$-th study, $\mu$ is a $p\times 1$ parameter vector of the average
treatment effect across all studies, ${\text{\boldmath$v$}}_{i}$ is a $p\times
1$ vector of an error between studies, ${\text{\boldmath$v$}}_{i}\sim{\cal
N}({\bf\text{\boldmath$0$}},{\text{\boldmath$V$}}({\text{\boldmath$\theta$}}))$,
${\text{\boldmath$V$}}({\text{\boldmath$\theta$}})$ is a $p\times p$
covariance matrix, and $\theta$ is a $q\times 1$ nuisance unknown parameter
vector. The covariance structure of
${\text{\boldmath$V$}}({\text{\boldmath$\theta$}})$ is known. The treatment
effects ${\text{\boldmath$\beta$}}_{i}$ within the $i$-th study are extracted
$p_{i}$ parameters by the design matrix ${\text{\boldmath$X$}}_{i}$.

To easily calculate the Bartlett-type adjustment, we consider the marginal
model of the random effects meta-analysis model as

${\text{\boldmath$y$}}={\text{\boldmath$X$}}{\text{\boldmath$\mu$}}+{\text{\boldmath$\epsilon$}},\hskip
10.0pt{\text{\boldmath$\epsilon$}}\sim{\cal
N}({\bf\text{\boldmath$0$}},{\text{\boldmath${\Sigma}$}}({\text{\boldmath$\theta$}})),\hskip
10.0pt{\text{\boldmath${\Sigma}$}}({\text{\boldmath$\theta$}})={\text{\boldmath$R$}}+{\text{\boldmath$G$}}({\text{\boldmath$\theta$}}),$
(3.1)

where
${\text{\boldmath$y$}}=({\text{\boldmath$y$}}_{1}^{T},\ldots,{\text{\boldmath$y$}}_{n}^{T})^{T}$
is an $N\times 1$ vector,
${\text{\boldmath$X$}}=({\text{\boldmath$X$}}_{1}^{T},\ldots,{\text{\boldmath$X$}}_{n}^{T})^{T}$
is an $N\times p$ matrix,
$R$=diag$({\text{\boldmath$R$}}_{1},\ldots,{\text{\boldmath$R$}}_{n})$ is a
$N\times N$ matrix,
${\text{\boldmath$G$}}({\text{\boldmath$\theta$}})={\text{\boldmath$Z$}}$diag$({\text{\boldmath$V$}}({\text{\boldmath$\theta$}}),\ldots,{\text{\boldmath$V$}}({\text{\boldmath$\theta$}})){\text{\boldmath$Z$}}^{T}$
is an $N\times N$ matrix,
$Z$=diag$({\text{\boldmath$X$}}_{1},\ldots,{\text{\boldmath$X$}}_{n})$ is an
$N\times np$ matrix. diag denotes the block diagonal matrix. We assume
rank$({\text{\boldmath$X$}})=p$.

We provide a simple example of the model. There are three studies in the NMA.
The first study has three treatments (drug A, drug B, and placebo). The second
study has two treatments (drug A and placebo). The third study has two
treatments (drug B and placebo). The difference of A and P is $y_{A-P}$, the
difference of B and P is $y_{B-P}$, the standard error of the difference of A
and P is $s_{A-P}$, the standard error of the difference of A and P is
$s_{B-P}$. For simplicity, this example does not consider the correlation of
$A-P$ and $B-P$ and the difference between A and B. In the first study,
${\text{\boldmath$y$}}_{1}=(y_{A-P},y_{B-P})^{T}$,
${\text{\boldmath$X$}}_{1}=((1,0),(0,1))$,
${\text{\boldmath$R$}}_{1}=((s^{2}_{A-P},s_{A-P}s_{B-P}/2),(s_{A-P}s_{B-P}/2,s^{2}_{B-P}))$,
${\text{\boldmath$V$}}({\text{\boldmath$\theta$}})=(({\theta}_{1},{\theta}_{3}),({\theta}_{3},{\theta}_{2}))$.
In the second study, $y_{2}=y_{A-P}$, ${\text{\boldmath$X$}}_{2}=(1,0)$,
${\text{\boldmath$R$}}_{2}=s^{2}_{A-P}$. In the third study, $y_{3}=y_{B-P}$,
${\text{\boldmath$X$}}_{3}=(0,1)$, ${\text{\boldmath$R$}}_{3}=s^{2}_{B-P}$.
The marginal model of the random effects meta-analysis model is obtained by
combining the information of each study.

In general, the NMA makes evaluations using the confidence intervals of
$\mu_{i}$. The confidence intervals based on LR statistic and Score statistic
need the null hypothesis[11, 18]. The null hypothesis is defined,

$H_{0}\ :\
{\text{\boldmath$r$}}^{T}{\text{\boldmath$\mu$}}={\text{\boldmath$r$}}^{T}{\text{\boldmath$\mu$}}_{0},$
(3.2)

where $r$ and ${\text{\boldmath$\mu$}}_{0}$ are known $p\times 1$ constant
vectors, the $i$-th element of $r$ is $1$ and the others are $0$. The element
$1$ of $r$ means to indicate the confidence interval for the parameter of
interest. If we are interested in differences between treatments, we set
${\text{\boldmath$\mu$}}={\bf\text{\boldmath$0$}}$. If you are interested in
whether the confidence interval is included for some clinically meaningful
threshold, put the clinically meaningful threshold in
${\text{\boldmath$\mu$}}_{0}$. If we are interested in whether or not a
confidence interval includes for some clinically meaningful value, we put a
clinically meaningful threshold in ${\text{\boldmath$\mu$}}_{0}$.

###  3.2 Wald, LR, and Score statistic

We introduce the Wald, LR, and Score statistics. The three statistics consist
of two types of the log likelihood function, the general log likelihood
function and the constrained log likelihood function under the null
hypothesis.

First of all, the log likelihood function is

$l^{L}({\widehat{\text{\boldmath$\mu$}}}({\text{\boldmath$\theta$}}),{\text{\boldmath$\theta$}})=C_{1}-\frac{1}{2}\mbox{ln}|{\Sigma}({\text{\boldmath$\theta$}})|-\frac{1}{2}\left({\text{\boldmath$y$}}-{\text{\boldmath$X$}}{\widehat{\text{\boldmath$\mu$}}}({\text{\boldmath$\theta$}})\right)^{T}{\Sigma}^{-1}({\text{\boldmath$\theta$}})\left({\text{\boldmath$y$}}-{\text{\boldmath$X$}}{\widehat{\text{\boldmath$\mu$}}}({\text{\boldmath$\theta$}})\right),$
(3.3)

where $C_{1}$ is a constant value and
${\widehat{\text{\boldmath$\mu$}}}({\text{\boldmath$\theta$}})=\left({\text{\boldmath$X$}}^{\prime}{\Sigma}^{-1}({\text{\boldmath$\theta$}}){\text{\boldmath$X$}}\right)^{-1}{\text{\boldmath$X$}}^{\prime}{\Sigma}^{-1}({\text{\boldmath$\theta$}}){\text{\boldmath$y$}}$
is the maximum likelihood estimator derived by maximizing
$l^{L}({\text{\boldmath$\mu$}},{\text{\boldmath$\theta$}})$ about $\mu$. The
MLE of $\theta$ is derived by maximizing (3.3). Next, we consider the
constrained log likelihood function under the null hypothesis, which can be
derived using the method of the Lagrange multipliers about
$l^{L}({\text{\boldmath$\mu$}},{\text{\boldmath$\theta$}})$. The constrained
log likelihood function under the null hypothesis has
$l^{L}({\widetilde{\text{\boldmath$\mu$}}}({\text{\boldmath$\theta$}}),{\text{\boldmath$\theta$}})$,
where

${\widetilde{\text{\boldmath$\mu$}}}({\text{\boldmath$\theta$}})={\widehat{\text{\boldmath$\mu$}}}({\text{\boldmath$\theta$}})-\frac{\left({\text{\boldmath$r$}}^{T}{\widehat{\text{\boldmath$\mu$}}}({\text{\boldmath$\theta$}})-{\text{\boldmath$r$}}^{T}{\text{\boldmath$\mu$}}_{0}\right)}{{\text{\boldmath$r$}}^{T}({\text{\boldmath$X$}}^{T}{\Sigma}^{-1}({\text{\boldmath$\theta$}}){\text{\boldmath$X$}})^{-1}{\text{\boldmath$r$}}}({\text{\boldmath$X$}}^{T}{\Sigma}^{-1}({\text{\boldmath$\theta$}}){\text{\boldmath$X$}})^{-1}{\text{\boldmath$r$}}.$
(3.4)

The constrained log likelihood function under the null hypothesis can be
expressed by replacing
${\widehat{\text{\boldmath$\mu$}}}({\text{\boldmath$\theta$}})$ of (3.3) with
${\widetilde{\text{\boldmath$\mu$}}}({\text{\boldmath$\theta$}})$, and this
relationship is important to derive the Bartlett-type adjustment terms. The
MLE of $\theta$ under the null hypothesis is derived by maximizing
$l^{L}({\widetilde{\text{\boldmath$\mu$}}}({\text{\boldmath$\theta$}}),{\text{\boldmath$\theta$}})$.

In the NMA, the REMLE is often used. We assume that $K$ is a function of the
$X$ satisfying
${\text{\boldmath$K$}}{\text{\boldmath$X$}}={\text{\boldmath$O$}}$ and
rank($K$)= $N-$rank$({\text{\boldmath$X$}})$. The restricted likelihood of
data $y$ is defined as the likelihood of data $K$$y$, the restricted log
likelihood function is shown as
$l^{RL}({\widehat{\text{\boldmath$\mu$}}}({\text{\boldmath$\theta$}}),{\text{\boldmath$\theta$}})=C_{2}-\frac{1}{2}\mbox{ln}|X^{T}{\Sigma}^{-1}({\text{\boldmath$\theta$}})X|+l^{L}({\widehat{\text{\boldmath$\mu$}}}({\text{\boldmath$\theta$}}),{\text{\boldmath$\theta$}})$[19,
20], where $C_{2}$ is a constant value. The REMLE of $\theta$ is derived by
maximizing
$l^{RL}({\widehat{\text{\boldmath$\mu$}}}({\text{\boldmath$\theta$}}),{\text{\boldmath$\theta$}})$.
We consider the constrained log likelihood function under the null hypothesis.
Since the restricted log likelihood function does not include $\mu$, we cannot
define the constrained log likelihood function under the null hypothesis in
the same way as the log likelihood function
$l^{L}({\widetilde{\text{\boldmath$\mu$}}}({\text{\boldmath$\theta$}}),{\text{\boldmath$\theta$}})$.
Welham and Thompson[19] show the two types of constrained log likelihood
function under the null hypothesis. The first is the sub-model and The second
is dropped terms relating to the null hypothesis. However, the both restricted
log likelihood functions cannot lead to the closed-form confidence intervals
for the LR and Score statistics. In this paper, we refer to Welham and
Thompson[19]’s function dropping terms, which assumes a part of $\mu$ being
$0$. We define a new extension of Welham and Thompson[19]’s function for
leading the closed-form confidence intervals. Now, we substitute $\mu$ into
${\widehat{\text{\boldmath$\mu$}}}({\text{\boldmath$\theta$}})$ of
$l^{RL}({\widehat{\text{\boldmath$\mu$}}}({\text{\boldmath$\theta$}}),{\text{\boldmath$\theta$}})$
and we obtain the constrained log likelihood function under the null
hypothesis in the same manner as the log likelihood function. Then, the new
constrained log likelihood function under the null hypothesis can be expressed
as
$l^{RL}({\widetilde{\text{\boldmath$\mu$}}}({\text{\boldmath$\theta$}}),{\text{\boldmath$\theta$}})$,
and the REMLE of $\theta$ under the null hypothesis can be obtained. The REMLE
of $\theta$ is confirmed to have consistency in Appendix A.10.

From here, when it is not necessary to distinguish $l^{L}$ from $l^{RL}$, we
write just $l$. In addition, we define that
${\widehat{\text{\boldmath$\theta$}}}$ is the MLE or REMLE, and that
${\widetilde{\text{\boldmath$\theta$}}}$ is the MLE or REMLE under the null
hypothesis.

We introduce the three statistics using the difference between two log
likelihood functions.

[1] Wald statistic

$\displaystyle
W({\widehat{\text{\boldmath$\theta$}}})=-2\left(l({\widetilde{\text{\boldmath$\mu$}}}({\widehat{\text{\boldmath$\theta$}}}),{\widehat{\text{\boldmath$\theta$}}})-l({\widehat{\text{\boldmath$\mu$}}}({\widehat{\text{\boldmath$\theta$}}}),{\widehat{\text{\boldmath$\theta$}}})\right)=\frac{\left({\text{\boldmath$r$}}^{T}{\widehat{\text{\boldmath$\mu$}}}({\widehat{\text{\boldmath$\theta$}}})-{\text{\boldmath$r$}}^{T}{\text{\boldmath$\mu$}}_{0}\right)^{2}}{h({\widehat{\text{\boldmath$\theta$}}})^{2}},{}$

where
$h({\text{\boldmath$\theta$}})=\sqrt{{\text{\boldmath$r$}}^{T}({\text{\boldmath$X$}}^{T}{\text{\boldmath${\Sigma}$}}({\text{\boldmath$\theta$}})^{-1}{\text{\boldmath$X$}})^{-1}{\text{\boldmath$r$}}}$.
The derivation of the last term is shown in Appendix A.3. From the last term,
$W({\widehat{\text{\boldmath$\theta$}}})$ is the closed-form. Thus, we can
easily derive the confidence interval.

[2] LR statistic

$\displaystyle
LR({\widehat{\text{\boldmath$\theta$}}},{\widetilde{\text{\boldmath$\theta$}}})=$
$\displaystyle-2\left(l({\widetilde{\text{\boldmath$\mu$}}}({\widetilde{\text{\boldmath$\theta$}}}),{\widetilde{\text{\boldmath$\theta$}}})-l({\widehat{\text{\boldmath$\mu$}}}({\widehat{\text{\boldmath$\theta$}}}),{\widehat{\text{\boldmath$\theta$}}})\right)=-2{\delta}({\widehat{\text{\boldmath$\theta$}}},{\widetilde{\text{\boldmath$\theta$}}})+W({\widetilde{\text{\boldmath$\theta$}}}),{}$

because
$l({\widetilde{\text{\boldmath$\mu$}}}({\widetilde{\text{\boldmath$\theta$}}}),{\widetilde{\text{\boldmath$\theta$}}})-l({\widehat{\text{\boldmath$\mu$}}}({\widehat{\text{\boldmath$\theta$}}}),{\widehat{\text{\boldmath$\theta$}}})={\delta}({\widehat{\text{\boldmath$\theta$}}},{\widetilde{\text{\boldmath$\theta$}}})+\left(l({\widetilde{\text{\boldmath$\mu$}}}({\widetilde{\text{\boldmath$\theta$}}}),{\widetilde{\text{\boldmath$\theta$}}})-l({\widehat{\text{\boldmath$\mu$}}}({\widetilde{\text{\boldmath$\theta$}}}),{\widetilde{\text{\boldmath$\theta$}}})\right)$
and Appendix A.1.
${\delta}({\widehat{\text{\boldmath$\theta$}}},{\widetilde{\text{\boldmath$\theta$}}})=l({\widehat{\text{\boldmath$\mu$}}}({\widetilde{\text{\boldmath$\theta$}}}),{\widetilde{\text{\boldmath$\theta$}}})-l({\widehat{\text{\boldmath$\mu$}}}({\widehat{\text{\boldmath$\theta$}}}),{\widehat{\text{\boldmath$\theta$}}})=O_{p}(N^{-1})$
since it can be proved in Appendix A.7 and [21]. The
${\text{\boldmath$\mu$}}_{0}$ is only included in
$W({\widetilde{\text{\boldmath$\theta$}}})$,
$W({\widetilde{\text{\boldmath$\theta$}}})$ is the closed-form about
${\text{\boldmath$\mu$}}_{0}$. Thus, we can derive the closed-form confidence
interval of
$LR({\widehat{\text{\boldmath$\theta$}}},{\widetilde{\text{\boldmath$\theta$}}})$.
Hardy and Thompson[22] propose a confidence interval of likelihood ratio
statistic using a likelihood function. However, The confidence interval is
calculated to satisfy
$l({\text{\boldmath$\mu$}},{\widehat{\text{\boldmath$\theta$}}})>l({\widehat{\text{\boldmath$\mu$}}}({\widehat{\text{\boldmath$\theta$}}}),{\widehat{\text{\boldmath$\theta$}}})-3.84/2$,
which increases the computational load. We consider the simple confidence
interval of likelihood ratio statistic based on the null hypothesis
${\text{\boldmath$r$}}^{T}{\text{\boldmath$\mu$}}={\text{\boldmath$r$}}^{T}{\text{\boldmath$\mu$}}_{0}$.

[3] Score statistic

$S({\widetilde{\text{\boldmath$\theta$}}})=-2\left(l({\widetilde{\text{\boldmath$\mu$}}}({\widetilde{\text{\boldmath$\theta$}}}),{\widetilde{\text{\boldmath$\theta$}}})-l({\widehat{\text{\boldmath$\mu$}}}({\widetilde{\text{\boldmath$\theta$}}}),{\widetilde{\text{\boldmath$\theta$}}})\right)=\frac{({\text{\boldmath$r$}}^{T}{\widehat{\text{\boldmath$\mu$}}}({\widetilde{\text{\boldmath$\theta$}}})-{\text{\boldmath$r$}}^{T}{\text{\boldmath$\mu$}}_{0})^{2}}{h({\widetilde{\text{\boldmath$\theta$}}})^{2}}=W({\widetilde{\text{\boldmath$\theta$}}}).{}$

The Score statistic can be defined by replacing
${\widehat{\text{\boldmath$\theta$}}}$ of the Wald statistic with
${\widetilde{\text{\boldmath$\theta$}}}$. Therefore, because
$S({\widetilde{\text{\boldmath$\theta$}}})$ can be expressed as closed-formed
$W(\cdot)$, we can formulate the three confidence intervals.

We formulate the confidence intervals. These three statistics are
asymptotically distributed as $\chi^{2}_{1}$. The coverage probability of
confidence interval of the Wald statistic is shown as

$\displaystyle P\left(W({\text{\boldmath$\theta$}})\leq
z_{{\alpha}\%,\chi^{2}}\right)$
$\displaystyle=P\left(-\sqrt{z_{{\alpha}\%,\chi^{2}}}\leq\frac{{\text{\boldmath$r$}}^{T}{\widehat{\text{\boldmath$\mu$}}}({\widehat{\text{\boldmath$\theta$}}})-{\text{\boldmath$r$}}^{T}{\text{\boldmath$\mu$}}}{h({\widehat{\text{\boldmath$\theta$}}})}\leq\sqrt{z_{{\alpha},\chi^{2}}}\right){}$
$\displaystyle=P\left({\text{\boldmath$r$}}^{T}{\widehat{\text{\boldmath$\mu$}}}({\widehat{\text{\boldmath$\theta$}}})-z_{{\alpha}\%/2}h({\widehat{\text{\boldmath$\theta$}}})\leq{\text{\boldmath$r$}}^{T}{\text{\boldmath$\mu$}}\leq{\text{\boldmath$r$}}^{T}{\widehat{\text{\boldmath$\mu$}}}({\widehat{\text{\boldmath$\theta$}}})+z_{{\alpha}/2}h({\widehat{\text{\boldmath$\theta$}}})\right),{}$

where $z_{{\alpha}\%,\chi^{2}}$ is the upper ${\alpha}\%$ of $\chi^{2}_{1}$,
$z_{{\alpha}\%/2}$ is the upper ${\alpha}\%/2$ of the standard normal
distribution and $z_{{\alpha}\%,\chi^{2}}=z_{{\alpha}\%/2}^{2}$. Since the
confidence intervals are not based on the null hypothesis,
${\text{\boldmath$\mu$}}_{0}$ in each statistic must be replaced with $\mu$.
Then, based on $W(\cdot)$, the confidence intervals of the other statistics
are also given below.

Confidence interval based on the Wald statistic:

$\displaystyle{\text{\boldmath$r$}}^{T}{\widehat{\text{\boldmath$\mu$}}}({\widehat{\text{\boldmath$\theta$}}})-z_{{\alpha}/2}h({\widehat{\text{\boldmath$\theta$}}})\leq{\text{\boldmath$r$}}^{T}{\text{\boldmath$\mu$}}\leq{\text{\boldmath$r$}}^{T}{\widehat{\text{\boldmath$\mu$}}}({\widehat{\text{\boldmath$\theta$}}})+z_{{\alpha}/2}h({\widehat{\text{\boldmath$\theta$}}}).$
(3.5)

Confidence interval based on the LR statistic:

$\displaystyle{\text{\boldmath$r$}}^{T}{\widehat{\text{\boldmath$\mu$}}}({\widetilde{\text{\boldmath$\theta$}}})-\sqrt{z_{{\alpha}/2}^{2}+2{\delta}({\widehat{\text{\boldmath$\theta$}}},{\widetilde{\text{\boldmath$\theta$}}})}h({\widetilde{\text{\boldmath$\theta$}}})\leq{\text{\boldmath$r$}}^{T}{\text{\boldmath$\mu$}}\leq{\text{\boldmath$r$}}^{T}{\widehat{\text{\boldmath$\mu$}}}({\widetilde{\text{\boldmath$\theta$}}})+\sqrt{z_{{\alpha}/2}^{2}+2{\delta}({\widehat{\text{\boldmath$\theta$}}},{\widetilde{\text{\boldmath$\theta$}}})}h({\widetilde{\text{\boldmath$\theta$}}}).$
(3.6)

Confidence interval based on the Score statistic:

$\displaystyle{\text{\boldmath$r$}}^{T}{\widehat{\text{\boldmath$\mu$}}}({\widetilde{\text{\boldmath$\theta$}}})-z_{{\alpha}/2}h({\widetilde{\text{\boldmath$\theta$}}})\leq{\text{\boldmath$r$}}^{T}{\text{\boldmath$\mu$}}\leq{\text{\boldmath$r$}}^{T}{\widehat{\text{\boldmath$\mu$}}}({\widetilde{\text{\boldmath$\theta$}}})+z_{{\alpha}/2}h({\widetilde{\text{\boldmath$\theta$}}}).$
(3.7)

We show numerical examples in the supplementary materials.

Noma et al.[11] do not show the closed-form confidence interval of the LR and
Score statistics. However, we can derive the closed-form confidence intervals
of all statistics based on $W(\cdot)$. On the other hand, the three confidence
intervals cannot hold the nominal confidence level when the number of studies
is small because the approximation accuracy of the statistics to
$\chi^{2}_{1}$ depends on the number of studies. To be strict, distribution
functions of the statistics have a bias term of the second order. In other
words, in the case of Wald, the distribution function is
$P\left(W({\widehat{\text{\boldmath$\theta$}}})\leq
z_{{\alpha},\chi^{2}}\right)=1-{\alpha}+O\left(N^{-1}\right)$, which has the
bias term($O\left(N^{-1}\right)$). In the next chapter, we derive the
asymptotic expanded distribution function and the Bartlett adjustment-type
terms of the second-order bias.

##  4 Bartlett-type adjustment

The Bartlett-type adjustment is a method that removes the bias of the second-
order. If we apply the Bartlett-type adjustment to the statistic,
$P\left(W_{ad}({\widehat{\text{\boldmath$\theta$}}})\leq
z_{{\alpha},\chi^{2}}\right)=1-{\alpha}+O\left(N^{-3/2}\right)$ where
$W_{ad}({\widehat{\text{\boldmath$\theta$}}})$ is the adjusted statistic. We
show the technical notation and assumption of
${\widehat{\text{\boldmath$\theta$}}}$ in the supplemental material.

Firstly, we consider to expanding asymptotically the distribution function of
the Wald statistic. Kojima and Kubokawa[12] have already derived the
asymptotic expansion of the distribution function of the Wald statistic.
However, it was expressed as a matrix form and is complicated to understand.
In this study, we re-calculate it in a vector form to simplify some equations
and make it easy to understand.

The original Wald statistic $W({\widehat{\text{\boldmath$\theta$}}})$ is
difficult to asymptotic, thus we divide it into three parts,

$W({\widehat{\text{\boldmath$\theta$}}})=\frac{(a+b)^{2}}{1-c},{}$

where

$\displaystyle
a=\frac{{\text{\boldmath$r$}}^{T}{\widehat{\text{\boldmath$\mu$}}}({\widehat{\text{\boldmath$\theta$}}})-{\text{\boldmath$r$}}^{T}{\widehat{\text{\boldmath$\mu$}}}({\text{\boldmath$\theta$}})}{h({\text{\boldmath$\theta$}})},\hskip
2.0ptb=\frac{{\text{\boldmath$r$}}^{T}{\widehat{\text{\boldmath$\mu$}}}({\text{\boldmath$\theta$}})-{\text{\boldmath$r$}}^{T}{\text{\boldmath$\mu$}}_{0}}{h({\text{\boldmath$\theta$}})},\hskip
2.0ptc=-\frac{h({\widehat{\text{\boldmath$\theta$}}})^{2}-h({\text{\boldmath$\theta$}})^{2}}{h({\text{\boldmath$\theta$}})^{2}}.{}$

From the assumption and Taylor series, we have $a=a_{1}+O_{p}(N^{-1})$, where
$a_{1}=O_{p}(N^{-1/2})$, and $c=c_{1}+c_{2}+O_{p}(N^{-3/2})$, where
$c_{1}=O_{p}(N^{-1/2})$, $c_{2}=O_{p}(N^{-1})$. We show the technical
description to derive the Bartlett-type adjustment in the supplementary
materials. The Bartlett-type adjustment $w_{ad}(x)$ satisfies
$P\left(W({\widehat{\text{\boldmath$\theta$}}})\leq
w_{ad}(x)x\right)=F_{\chi^{2}_{1}}(x)+O(N^{-3/2})$ where

$\displaystyle
w_{ad}(x)=1+\widehat{E[a_{1}^{2}]}+\widehat{E[c_{2}]}+\frac{1}{4}(1+x)\widehat{E[c_{1}^{2}]}.$
(4.1)

The $\widehat{E[a_{1}^{2}]}$, $\widehat{E[c_{2}]}$, and
$\widehat{E[c_{1}^{2}]}$ are plugged ${\widehat{\text{\boldmath$\theta$}}}$ to
$E[a_{1}^{2}]$, $E[c_{2}]$, and $E[c_{1}^{2}]$, respectively.
$\widehat{E[a_{1}^{2}]}$ means the bias of
${\widehat{\text{\boldmath$\beta$}}}({\widehat{\text{\boldmath$\theta$}}})$.
$\widehat{E[c_{2}]}$ and $\widehat{E[c_{1}^{2}]}$ mean the biases of
$h({\widehat{\text{\boldmath$\theta$}}})$. To understand the bias terms, we
consider the simplest case which is
${\text{\boldmath${\Sigma}$}}({\text{\boldmath$\theta$}})={\theta}_{1}{\text{\boldmath$I$}}_{N}$.
In this case, $\widehat{E[a_{1}^{2}]}$ is zero because
${\widehat{\text{\boldmath$\mu$}}}({\text{\boldmath$\theta$}})$ does not
include $\theta$, $\widehat{E[c_{2}]}$ is zero because
$h({\text{\boldmath$\theta$}})^{2}$ is a linear function about $\theta$, and
$\widehat{E[c_{1}^{2}]}$ is $\frac{1}{2N}$. Therefore, from (4.1), the Wald
statistic has a bias of $\frac{1}{8N}(1+x)$. However, in this study,
${\text{\boldmath${\Sigma}$}}({\text{\boldmath$\theta$}})$ becomes more
complicated, and each bias term may be larger depending on $\theta$.

Next, we derive the Bartlett-type adjustment terms of the LR statistic. By
using the Taylor expansion of
$l({\widetilde{\text{\boldmath$\mu$}}}({\widehat{\text{\boldmath$\theta$}}}),{\widehat{\text{\boldmath$\theta$}}})$
of the LR statistic around ${\widetilde{\text{\boldmath$\theta$}}}$, the
expanded LR statistic is obtained (details in Appendix A.9). The
characteristic function of the LR statistic is given using Appendix A.10. The
expanded distribution function of the LR statistic is

$\displaystyle
P(LR({\widehat{\text{\boldmath$\theta$}}},{\widetilde{\text{\boldmath$\theta$}}})\leq
x)=$ $\displaystyle
F_{\chi^{2}_{1}}(x)-x\left(E[c_{2}]+\frac{1}{4}E[c_{1}^{2}]\right)f_{\chi^{2}_{1}}(x)+O(N^{-3/2}).{}$

There is no bias term of
${\widehat{\text{\boldmath$\mu$}}}({\widehat{\text{\boldmath$\theta$}}})$ in
the second-order. In addition, the weight of the bias $E[c_{1}^{2}]$ is
reduced to $\frac{x}{4}$ compared with the Wald statistic. Therefore, the LR
statistic is expected to be less biased than the Wald statistic. Besides,
because ${\widehat{\text{\boldmath$\theta$}}}^{REML}$ has been shown to be
less the bias asymptotically than ${\widehat{\text{\boldmath$\theta$}}}^{ML}$
by Kubokawa[21], the LR statistic with the REMLE is considered to be less
biased.

We obtain the Bartlett-type adjustment terms,

$\displaystyle
lr_{ad}(x)=1+\widehat{E[c_{2}]}+\frac{1}{4}\widehat{E[c_{1}^{2}]}.$ (4.2)

Finally, we derive the Bartlett-type adjustment terms of the Score statistic.
The Score statistic is asymptotically expressed by the Taylor expansion
$l({\widetilde{\text{\boldmath$\mu$}}}({\widehat{\text{\boldmath$\theta$}}}),{\widehat{\text{\boldmath$\theta$}}})$
of the LR statistic around ${\widetilde{\text{\boldmath$\theta$}}}$, yielding
$S({\widetilde{\text{\boldmath$\theta$}}})=W({\widetilde{\text{\boldmath$\theta$}}})=W({\widehat{\text{\boldmath$\theta$}}})-\frac{1}{2}\sum_{i,j=1}^{q}{\hat{\theta}}_{1i}{\hat{\theta}}_{1j}(2a_{1i}b-b^{2}c_{1i})(2a_{1j}b-b^{2}c_{1j})+O_{p}(N^{-3/2})$.
The characteristic function is obtained using Appendix A.10. Thus, calculating
the distribution function using the inversion formula, we obtain the Bartlett-
type adjustment terms where

$\displaystyle
s_{ad}(x)=1-\widehat{E[a^{2}_{1}]}+\widehat{E[c_{2}]}+\frac{1}{4}(1-x)\widehat{E[c_{1}^{2}]}.$
(4.3)

The Score statistic also includes the bias $\widehat{E[a_{1}^{2}]}$ of
${\widehat{\text{\boldmath$\mu$}}}({\widehat{\text{\boldmath$\theta$}}})$, and
the extra bias of the variance is similar to that of the Wald statistic. We
consider why Wald and Score statistics include the extra biases relative to
the LR statistic. The Wald and Score statistics have the mismatched estimator
${\widetilde{\text{\boldmath$\mu$}}}({\widehat{\text{\boldmath$\theta$}}})$ or
${\widehat{\text{\boldmath$\mu$}}}({\widetilde{\text{\boldmath$\theta$}}})$,
which have mixed estimators calculated under the non-null and null hypotheses.
Therefore, the extra biases occur. On the other hand, the LR statistic is that
there is no problem in the combination of estimators, and the bias is small.
We show the formulas of each adjustment terms in Appendix 14.

Theorem 1.
Following the assumption in appendix A.2, the Bartlett-type adjusted three
statistics are (4.1), (4.2), and (4.3), respectively. Then, we obtain the
adjusted closed-form confidence intervals below. Confidence interval of the
Bartlett-type adjusted Wald statistic.
$\displaystyle{\text{\boldmath$r$}}^{T}{\widehat{\text{\boldmath$\mu$}}}({\widehat{\text{\boldmath$\theta$}}})-z_{{\alpha}/2}h({\widehat{\text{\boldmath$\theta$}}})\sqrt{w_{ad}(z_{{\alpha}/2}^{2})}\leq{\text{\boldmath$r$}}^{T}{\text{\boldmath$\mu$}}\leq{\text{\boldmath$r$}}^{T}{\widehat{\text{\boldmath$\mu$}}}({\widehat{\text{\boldmath$\theta$}}})+z_{{\alpha}\%/2}h({\widehat{\text{\boldmath$\theta$}}})\sqrt{w_{ad}(z_{{\alpha}/2}^{2})}$
(4.4) Confidence interval of the Bartlett-type adjusted LR statistic.
$\displaystyle{\text{\boldmath$r$}}^{T}{\widehat{\text{\boldmath$\mu$}}}({\widetilde{\text{\boldmath$\theta$}}})-\sqrt{LR_{ad}}h({\widetilde{\text{\boldmath$\theta$}}})\leq{\text{\boldmath$r$}}^{T}{\text{\boldmath$\mu$}}\leq{\text{\boldmath$r$}}^{T}{\widehat{\text{\boldmath$\mu$}}}({\widetilde{\text{\boldmath$\theta$}}})+\sqrt{LR_{ad}}h({\widetilde{\text{\boldmath$\theta$}}}),$
(4.5) where
$LR_{ad}=z_{{\alpha}/2}^{2}lr_{ad}(z_{{\alpha}/2}^{2})+2{\delta}({\widehat{\text{\boldmath$\theta$}}},{\widetilde{\text{\boldmath$\theta$}}})$.
Confidence interval of the Bartlett-type adjusted Score statistic.
$\displaystyle{\text{\boldmath$r$}}^{T}{\widehat{\text{\boldmath$\mu$}}}({\widetilde{\text{\boldmath$\theta$}}})-z_{{\alpha}/2}h({\widetilde{\text{\boldmath$\theta$}}})\sqrt{s_{ad}(z_{{\alpha}/2}^{2})}\leq{\text{\boldmath$r$}}^{T}{\text{\boldmath$\mu$}}\leq{\text{\boldmath$r$}}^{T}{\widehat{\text{\boldmath$\mu$}}}({\widetilde{\text{\boldmath$\theta$}}})+z_{{\alpha}/2}h({\widetilde{\text{\boldmath$\theta$}}})\sqrt{s_{ad}(z_{{\alpha}/2}^{2})}$
(4.6) The coverage probabilities in the above adjusted confidence intervals
are $1-{\alpha}+O(N^{-3/2})$, thus the second-order bias is removed.

We show an example of the Bartlett-type adjustment terms by using the
assumptions used in the the confidence interval example in the section 3.
$\widehat{E[a^{2}_{1}]}=0$,
$\widehat{E[c_{1}^{2}]}=\frac{n}{(1+{\hat{\theta}})^{4}}$,
$\widehat{E[c_{2}]}=\frac{n^{2}}{2(1+{\hat{\theta}})^{4}}$ with the MLE, and
$\widehat{E[c_{2}]}=0$ with the REMLE.

##  5 Higher order adjustment via parametric bootstrap

We consider whether further adjustments can be made to the adjusted
statistics. Kojima and Kubokawa[12] and Noma et al.[11] propose a bias
adjustment method using a parametric bootstrap. The parametric bootstrap
generates the bootstrap samples of $m$ sets from a parametric distribution.
The statistics using each bootstrap sample are calculated, we obtain the upper
${\alpha}\%$ of the sorted $m$ statistics as the rejection point. By using
this rejection point, the accuracy is impoved.

In this paper, we consider to use the parametric bootstrap method to
$P\left(W_{ad}({\widehat{\text{\boldmath$\theta$}}})\leq
z_{{\alpha}\%,\chi^{2}}\right)=1-{\alpha}+O\left(N^{-3/2}\right)$, thus we
have $P\left(W_{ad}({\widehat{\text{\boldmath$\theta$}}})\leq
z_{{\alpha}\%,w}^{*}\right)=1-{\alpha}+o\left(N^{-3/2}\right)$ where $*$
refers to bootstrap samples and $z_{{\alpha}\%,w}^{*}$ is the upper
${\alpha}\%$ of $W^{*}_{ad}({\widehat{\text{\boldmath$\theta$}}}^{*})$. We
propose the procedure for the bootstrap adjustment method below.

[Step 1.]
Generate the parametric bootstrap samples ${\text{\boldmath$y$}}^{*}$ of $m$
sets from
${\text{\boldmath$X$}}{\widetilde{\text{\boldmath$\mu$}}}({\widehat{\text{\boldmath$\theta$}}})+{\text{\boldmath$\epsilon$}}^{*}$
where [12] and [11] propose to generate using
${\widetilde{\text{\boldmath$\mu$}}}$ constrained under the null hypothesis.

[Step 2.]
Calculate the $m$ Bartlett-type adjusted statistics.

[Step 3.]
Select the upper ${\alpha}\%$ of m of sorted statistics, and find the square
root.

Theorem 2.
Following the assumption in appendix A.2, the Bartlett-type adjusted three
statistics are (4.1), (4.2), and (4.3), respectively. The upper ${\alpha}\%$
of the three statistics via the parametric bootstrap method are
$z_{{\alpha}\%,w}^{*}$, $z_{{\alpha},LR}^{*}$, and $z_{{\alpha},S}^{*}$. Then,
we obtain the adjusted closed-dorm confidence interval using the bootstrap
adjustment method, Bootstrap Confidence interval of the Bartlett-type adjusted
Wald statistic
$\displaystyle{\text{\boldmath$r$}}^{T}{\widehat{\text{\boldmath$\mu$}}}({\widehat{\text{\boldmath$\theta$}}})-h({\widehat{\text{\boldmath$\theta$}}})\sqrt{z_{{\alpha}\%,w}^{*}w_{ad}(z_{{\alpha}/2}^{2})}\leq{\text{\boldmath$r$}}^{T}{\text{\boldmath$\mu$}}\leq{\text{\boldmath$r$}}^{T}{\widehat{\text{\boldmath$\mu$}}}({\widehat{\text{\boldmath$\theta$}}})+z_{{\alpha}/2}h({\widehat{\text{\boldmath$\theta$}}})\sqrt{z_{{\alpha}\%,w}^{*}w_{ad}(z_{{\alpha}/2}^{2})}{}$
Bootstrap Confidence interval of the Bartlett-type adjusted LR statisticR
$\displaystyle{\text{\boldmath$r$}}^{T}{\widehat{\text{\boldmath$\mu$}}}({\widetilde{\text{\boldmath$\theta$}}})-\sqrt{LR^{\ast}_{ad}}h({\widetilde{\text{\boldmath$\theta$}}})\leq{\text{\boldmath$r$}}^{T}{\text{\boldmath$\mu$}}\leq{\text{\boldmath$r$}}^{T}{\widehat{\text{\boldmath$\mu$}}}({\widetilde{\text{\boldmath$\theta$}}})+\sqrt{LR^{\ast}_{ad}}h({\widetilde{\text{\boldmath$\theta$}}}),{}$
where
$LR^{\ast}_{ad}=z_{{\alpha},lr}^{*}lr_{ad}(z_{{\alpha}/2}^{2})+2{\delta}({\widehat{\text{\boldmath$\theta$}}},{\widetilde{\text{\boldmath$\theta$}}})$.
Bootstrap Confidence interval of the Bartlett-type adjusted Score statistic
$\displaystyle{\text{\boldmath$r$}}^{T}{\widehat{\text{\boldmath$\mu$}}}({\widetilde{\text{\boldmath$\theta$}}})-h({\widetilde{\text{\boldmath$\theta$}}})\sqrt{z_{{\alpha},s}^{*}s_{ad}(z_{{\alpha}/2}^{2})}\leq{\text{\boldmath$r$}}^{T}{\text{\boldmath$\mu$}}\leq{\text{\boldmath$r$}}^{T}{\widehat{\text{\boldmath$\mu$}}}({\widetilde{\text{\boldmath$\theta$}}})+h({\widetilde{\text{\boldmath$\theta$}}})\sqrt{z_{{\alpha},s}^{*}s_{ad}(z_{{\alpha}/2}^{2})}{}$
The coverage probabilities of (4.4), (4.5), and (4.6) are
$1-{\alpha}+O(N^{-3/2})$ The coverage probabilities in the bootstrap
confidence intervals are $1-{\alpha}+o(N^{-3/2})$, thus the third-order bias
is removed.

##  6 Simulation study

We perform a computer simulation to confirm the accuracy of the Bartlett-type
adjusted statistics through two simulation studies. We generate the datasets
using the same procedure as Sidik and Jonkman[23], Noma et al. [11], and
Nagashima et al. [24].

###  6.1 Simulation 1

There are four treatments (A-drug, B-drug, C-drug, and placebo [P]); the
network diagram and the number of studies are shown in Figure 1 and Table 1.
We are interested in a comparison between the outcome of P and each of the
three drugs.

Figure 1: Diagram of simulation study 1

Table 1: Setting of simulation study 1 Design | # of Studies | # of Studies
---|---|---
P vs A | 3 | 9
P vs B | 2 | 6
P vs C | 2 | 6
A vs B | 2 | 6
A vs C | 1 | 3
Total | 10 | 30

The first step is to generate the number of subjects in the studies
($n=n_{A}=n_{B}=n_{C}=n_{P}$) from the uniform distribution ${\cal U}(30,300)$
and the probability $p_{P}$ of occurrence of events in the placebo group from
the uniform distribution ${\cal U}(0.05,0.65)$. We set
${\text{\boldmath$\mu$}}=(\mu_{1},\mu_{2},\mu_{3})=(\mu_{A-P},\mu_{B-P},\mu_{C-P})=(0.398,0.702,0.866)$
and ${\theta}={\theta}_{1}=\ldots={\theta}_{6}=0.3$. All off-diagonal elements
of ${\text{\boldmath$V$}}({\theta})$ are ${\theta}/2$[25]. The next step is to
generate the treatment effect within the study from ${\cal
N}({\text{\boldmath$\mu$}},{\text{\boldmath$V$}}({\theta}))$. The
probabilities of groups A, B, and C group can be derived from $p_{P}$ and the
treatment effect. The final step is to generate the number of events from the
binomial distribution ${\cal B}(p_{i},n)$, where $i\in\\{P,A,B,C\\}$. The null
hypotheses are $H_{0}:\mu_{1}=0.398$ for $A-P$, $H_{0}:\mu_{2}=0.702$ for
$B-P$, and $H_{0}:\mu_{3}=0.866$ for $C-P$. The type 1 error is 5%, i.e., we
will confirm to be held the confidence interval to 95%. We count the number of
the confidence intervals that include the true value. The number of
simulations is 10,000 replications. When we use the bootstrap method, it is
computed on 1,000 replications and 1001 bootstrap samples[26], as the
numerical calculations take a great deal of time.

Table 2: Result of a simulation study (three comparisons)

N=10
---
estimator |  | W | WBC | Wboot | W${}^{BC}_{boot}$ | LR | LRBC | LRboot | LR${}^{BC}_{boot}$ | S | SBC | Sboot | S${}^{BC}_{boot}$
MLE | $\mu_{1}$ | 86.1 | 93.6 | 93.1 | 92.9 | 88.9 | 94.4 | 94.0 | 96.8 | 93.3 | 95.5 | 94.8 | 95.1
MLE | $\mu_{2}$ | 86.8 | 93.8 | 92.9 | 92.9 | 89.2 | 94.8 | 93.5 | 97.6 | 93.7 | 96.1 | 95.1 | 95.1
MLE | $\mu_{3}$ | 86.7 | 93.9 | 91.3 | 91.1 | 89.5 | 94.8 | 92.0 | 96.8 | 93.5 | 95.9 | 92.9 | 93.1
REMLE | $\mu_{1}$ | 91.3 | 94.6 | 96.2 | 96.2 | 94.8 | 95.5 | 95.8 | 95.8 | 98.2 | 96.3 | 95.8 | 95.8
REMLE | $\mu_{2}$ | 91.7 | 94.6 | 95.9 | 96.1 | 95.0 | 95.9 | 96.1 | 96.1 | 98.5 | 96.7 | 95.7 | 95.5
REMLE | $\mu_{3}$ | 91.9 | 94.6 | 94.8 | 94.8 | 95.1 | 95.9 | 95.2 | 95.0 | 98.5 | 96.7 | 94.8 | 95.0
N=30
estimator |  | W | WBC | Wboot | W${}^{BC}_{boot}$ | LR | LRBC | LRboot | LR${}^{BC}_{boot}$ | S | SBC | Sboot | S${}^{BC}_{boot}$
MLE | $\mu_{1}$ | 92.3 | 94.6 | 95.9 | 95.9 | 93.2 | 94.8 | 95.7 | 96.8 | 94.3 | 94.9 | 95.5 | 94.5
MLE | $\mu_{2}$ | 92.3 | 94.6 | 94.6 | 94.7 | 93.1 | 94.8 | 94.6 | 96.3 | 94.4 | 95.0 | 94.5 | 94.5
MLE | $\mu_{3}$ | 92.6 | 95.0 | 95.9 | 95.8 | 93.6 | 95.1 | 95.9 | 96.8 | 94.7 | 95.2 | 96.1 | 96.1
REMLE | $\mu_{1}$ | 93.8 | 94.9 | 95.7 | 95.7 | 94.7 | 94.9 | 95.4 | 95.3 | 95.7 | 95.0 | 95.3 | 95.3
REMLE | $\mu_{2}$ | 93.9 | 94.9 | 94.5 | 95.5 | 94.8 | 95.0 | 94.5 | 94.5 | 95.6 | 95.1 | 94.5 | 94.5
REMLE | $\mu_{3}$ | 94.2 | 95.1 | 96.0 | 95.9 | 95.0 | 95.3 | 95.8 | 95.8 | 96.1 | 95.3 | 96.0 | 96.0

  * •

W: Wald; WBC: Bartlett-type adjusted Wald; Wboot: Wald using bootstrap method;
W${}^{BC}_{boot}$: Bartlett-type adjusted Wald using bootstrap method; LR:
Likelihood Ratio; LRBC: Bartlett-type adjusted LR; LRboot: LR using bootstrap
method; LR${}^{BC}_{boot}$: Bartlett-type adjusted LR bootstrap method; S:
Score; SBC: Bartlett-type adjusted Score; Sboot: Score using bootstrap method;
S${}^{BC}_{boot}$: Bartlett-type adjusted Score using bootstrap method

The simulation results are in Table2. The all adjusted confidence intervals
have been improved and closed to 95%. The coverage probability of the Wald
statistic is under 95%, thus the confidence interval of the Wald statistic is
narrower (details in Appendix A.12 Average confidence intervals in computer
simulations). On the other hand, the confidence interval of the Score
statistic is wider, since the coverage probability exceed 95%. The confidence
interval of the LR statistic is relatively stable because the bias terms are
fewer than the Wald and Score statistics. In the case of MLE under $N=10$, the
adjusted confidence interval based on the Wald statistic have been
significantly improved. The confidence interval of the LR statistic using REML
estimator seems less biased because the percent change in the adjusted
statistic based on the unadjusted statistic in Appendix 12 in the supplemental
material is very small as 3%. As for the average simulation time of method,
when $N=10$, our method took $42.3$ minutes (0.3 seconds per one simulation),
whereas the bootstrap method took over about $4,000$ minutes (4 minutes per
one simulation). Moreover, when under $N=30$, our method took $49.8$ minutes
(0.3 seconds per one simulation), whereas the bootstrap method took over
$5,000$ minutes (5 minutes per one simulation).

###  6.2 Simulation 2

We assume that there seven treatments (A-drug, B-drug, C-drug, D-drug, E-drug,
F-drug, and placebo [P]); the network diagram and the number of studies are
shown in Figure 2 and Table 3. We are interested in a comparison between an
outcome of P and all six drugs. We assume the almost same assumption as in
Simulation 1, the probabilities of events are
$p_{A}=0.6,p_{B}=0.65,p_{C}=0.7,p_{D}=0.75,p_{E}=0.8,p_{F}=0.85,p_{P}=0.6$,
and bootstrap samples are generated 501 times. The null hypotheses are
$H_{0}:\mu_{1}=0$ for $A-P$, $H_{0}:\mu_{2}=0.093$ for $B-P$,
$H_{0}:\mu_{3}=0.192$ for $C-P$, $H_{0}:\mu_{4}=0.301$ for $D-P$,
$H_{0}:\mu_{5}=0.426$ for $E-P$, and $H_{0}:\mu_{6}=0.577$ for $F-P$. The
network diagram and the number of studies are shown in Figure 2 and Table 2.
The total number of studies to $N=20$ instead of $N=10$.

Figure 2: Diagram of simulation study 2

Table 3: Setting of simulation study 2 Comparison | # of Studies | # of Studies
---|---|---
P vs A | 6 | 4
P vs B | 3 | 2
P vs C | 3 | 2
A vs B | 6 | 4
A vs C | 3 | 2
A vs D | 3 | 2
B vs E | 3 | 2
C vs F | 3 | 2
Total | 30 | 20

Table 4: Covarage probability (%) in a simulation study (six comparisons)

N=20
---
estimator |  | W | WBC | Wboot | W${}^{BC}_{boot}$ | LR | LRBC | LRboot | LR${}^{BC}_{boot}$ | S | SBC | Sboot | S${}^{BC}_{boot}$
MLE | $\mu_{1}$ | 87.3 | 93.3 | 95.2 | 95.5 | 88.8 | 93.8 | 95.2 | 98.3 | 91.4 | 94.6 | 94.5 | 94.5
MLE | $\mu_{2}$ | 87.6 | 93.6 | 92.3 | 92.3 | 89.3 | 94.0 | 92.3 | 97.2 | 91.5 | 94.7 | 92.7 | 92.7
MLE | $\mu_{3}$ | 87.4 | 93.4 | 94.1 | 93.4 | 89.2 | 94.2 | 93.8 | 95.8 | 91.7 | 94.9 | 93.8 | 93.8
MLE | $\mu_{4}$ | 87.4 | 93.4 | 93.1 | 93.1 | 89.0 | 94.1 | 93.4 | 96.9 | 91.5 | 95.0 | 93.8 | 93.8
MLE | $\mu_{5}$ | 87.7 | 93.0 | 92.0 | 92.0 | 89.3 | 93.8 | 92.7 | 96.9 | 91.6 | 94.6 | 92.7 | 92.7
MLE | $\mu_{6}$ | 87.4 | 93.1 | 95.5 | 95.5 | 89.0 | 93.7 | 96.5 | 98.6 | 91.4 | 94.6 | 97.9 | 97.6
REMLE | $\mu_{1}$ | 93.2 | 94.7 | 95.8 | 95.5 | 95.0 | 95.3 | 95.8 | 95.8 | 96.9 | 95.9 | 95.8 | 95.8
REMLE | $\mu_{2}$ | 93.4 | 94.9 | 93.7 | 93.7 | 95.1 | 95.5 | 94.1 | 94.1 | 96.8 | 95.9 | 94.4 | 94.8
REMLE | $\mu_{3}$ | 93.6 | 95.1 | 95.5 | 95.5 | 95.4 | 95.7 | 95.5 | 95.5 | 97.2 | 96.2 | 94.4 | 94.4
REMLE | $\mu_{4}$ | 93.6 | 95.1 | 93.8 | 93.8 | 95.4 | 95.7 | 94.8 | 94.8 | 97.0 | 96.3 | 94.8 | 94.8
REMLE | $\mu_{5}$ | 93.2 | 94.6 | 94.4 | 94.4 | 95.0 | 95.4 | 94.8 | 94.8 | 96.7 | 95.9 | 95.5 | 95.5
REMLE | $\mu_{6}$ | 93.2 | 94.6 | 97.6 | 97.6 | 94.9 | 95.4 | 97.6 | 97.6 | 96.9 | 95.9 | 97.9 | 97.9
N=30
estimator |  | W | WBC | Wboot | W${}^{BC}_{boot}$ | LR | LRBC | LRboot | LR${}^{BC}_{boot}$ | S | SBC | Sboot | S${}^{BC}_{boot}$
MLE | $\mu_{1}$ | 90.2 | 93.9 | 94.8 | 94.8 | 91.3 | 94.3 | 95.1 | 96.7 | 92.7 | 94.8 | 95.3 | 95.1
MLE | $\mu_{2}$ | 90.2 | 94.0 | 93.6 | 93.7 | 91.2 | 94.3 | 93.8 | 95.9 | 92.6 | 94.6 | 93.8 | 93.8
MLE | $\mu_{3}$ | 90.0 | 94.0 | 93.7 | 93.7 | 91.1 | 94.3 | 93.7 | 96.3 | 92.6 | 94.7 | 94.2 | 94.1
MLE | $\mu_{4}$ | 89.8 | 93.7 | 93.7 | 93.7 | 90.8 | 94.1 | 93.8 | 96.7 | 92.3 | 94.3 | 94.6 | 94.4
MLE | $\mu_{5}$ | 90.3 | 94.0 | 93.1 | 93.1 | 91.4 | 94.4 | 93.5 | 95.9 | 92.7 | 94.8 | 93.2 | 93.1
MLE | $\mu_{6}$ | 90.1 | 93.8 | 93.7 | 93.8 | 91.0 | 94.1 | 94.1 | 96.6 | 92.4 | 94.5 | 94.6 | 94.6
REMLE | $\mu_{1}$ | 93.8 | 94.9 | 94.9 | 95.0 | 94.8 | 95.1 | 95.1 | 95.1 | 96.1 | 95.3 | 95.3 | 95.3
REMLE | $\mu_{2}$ | 93.6 | 94.8 | 93.8 | 93.7 | 94.8 | 95.0 | 93.8 | 93.8 | 95.9 | 95.2 | 93.8 | 93.8
REMLE | $\mu_{3}$ | 93.7 | 94.8 | 94.2 | 94.1 | 94.8 | 95.1 | 94.5 | 94.5 | 96.1 | 95.2 | 94.2 | 94.4
REMLE | $\mu_{4}$ | 93.5 | 94.5 | 94.9 | 94.9 | 94.5 | 94.9 | 95.3 | 95.0 | 95.9 | 95.2 | 95.0 | 95.0
REMLE | $\mu_{5}$ | 93.8 | 94.8 | 93.9 | 93.8 | 94.8 | 95.2 | 93.7 | 93.7 | 96.1 | 95.4 | 93.9 | 93.9
REMLE | $\mu_{6}$ | 93.5 | 94.7 | 94.6 | 94.6 | 94.7 | 95.0 | 94.2 | 94.2 | 96.0 | 95.3 | 94.2 | 94.2

The simulation results are shown in Table 4. Since the number of parameters
has increased, the coverage probabilities of the confidence intervals have
more bias. However, the adjusted confidence interval can be close to 95%.

As for the average simulation time of each method, when $N=20$, our method
took $230.1$ minutes (1.38 seconds per one simulation), whereas the bootstrap
method took over about $23,000$ minutes (23 minutes per one simulation).
Moreover, when $N=30$, our method took $283.7$ minutes (1.70 seconds per one
simulation), whereas the bootstrap method took over $28,000$ minutes (28
minutes per one simulation).

We overview the entire results of simulation. Since the LR statistic is less
biased compared to the Wald and score statistics and the adjusted LR statistic
is the most stable, we recommend to use the adjusted LR statistic. The REMLE
is more stable, We think it is better to use the REMLE than MLE.

##  7 Applications to actual network meta-analysis

###  7.1 Dual Antiplatelet therapy (DAPT)

Palmerini et al.[27] examined the mortality rate of three different extended
durations of the dual Antiplatelet therapy (DAPT) for the definite or probable
stent thrombosis rate via network meta-analysis of eight studies. The network
diagram is shown in Figure 2. The study showed that short duration of DAPT
decreases the mortality rate. We reanalyze the rate of the definite or
probable stent thrombosis using the network meta-analysis. Let the reference
group be 6 months, the variance-covariance matrix of random effect be the
compound symmetry, and the sets of bootstraps samples be 10,001. The null
hypotheses are set $H_{0}:\mu_{1}=0$ for 6 months vs 1 year, $H_{0}:\mu_{2}=0$
for 6 months vs 1 year over, and $H_{0}:\mu_{3}=0$ for 1 year vs 1 year over,
because we are interested whether there are differences.

We show the results of re-analysis in Figure 4. From the simulation results,
the naive confidence intervals are too narrow or too wide due to bias effects.
We confirmed the adjusted confidence intervals of the Wald statistic are
wider. The adjusted confidence intervals of the LR statistic are also wider.
The adjusted confidence intervals of the score statistic are narrower.

Figure 3: Palmerini2015 DAPT Figure 4: Palmerini2015 DAPT, definite or
probable stent thrombosis

##  8 Discussion

The problem with NMA using a small number of studies is that the confidence
intervals cannot hold a nominal confidence level. To resolve this problem, we
applied the Bartlett-type adjustment. Many Bartlett-type adjustment methods
are based on the MLE. However, the NMA often uses the REMLE, which has not
been extensively discussed in this context. We derived the novel Bartlett-type
adjustment terms (4.1) - (4.3) using not only MLE but also REMLE.

Noma et al.[11] proposed the Bartlett-type adjustment using $m$ bootstrap
samples for NMA with a small number of studies. However, there is a problem
when using bootstrap samples: the calculation load is very high, as the MLE or
REMLE is derived by numerically solving a likelihood equation that is not
closed-form. In addition, since the LR statistic and score statistic on the
REMLE are not closed-form, and those confidence intervals also need to be
derived numerically. When using the bootstrap method, we do not know which
bias terms cause the bias. We proposed the closed-form confidence intervals of
LR and score statistics and clarified the bias terms by deriving the Bartlett
adjustment term analytically. In addition, we proposed a higher-order
adjustment via the parametric bootstrap method.

In the simulation results, because the LR statistic is less biased than the
Wald and score statistics, and the adjusted LR statistic is the most stable,
we recommend using the adjusted LR statistic. Because the REMLE is more
stable, we think that it is preferable to use the REMLE rather than the MLE.
In particular, the LR test statistic with the REMLE performed very well
because it had less bias than the other statistics. We recommend using the LR
test statistic with the REMLE, but just in case, we recommend applying the
Bartlett-type adjustment to completely remove the bias. From demonstrations of
actual studies, we confirmed that the adjusted confidence intervals were
improved compared to the naive confidence intervals.

In the future, we believe that NMA will become more popular[2]. Consequently,
the Bartlett-type adjusted method will be used more frequently for the NMA
with a small number of in which the number of studies is small.

Acknowledgments. The author is grateful to the editor, the associate editor
and the referee for their valuable comments and helpful suggestions. The
author would like to thank Associate Professor Hisashi Noma for his
encouragement and helpful suggestions. This study was supported by a Grants-
in-Aid for Scientific Research from the Japan Society for the Promotion of
Science (Grant numbers: JP19H04074).

Disclosures.

None

## References

  * [1] Gurevitch J., Koricheva J., Nakagawa S., et al . Meta-analysis and the science of research synthesis. Nature. 2018;555:175-182.
  * [2] Li L, Catala-Lopez F, Alonso-Arroyo A, et al . The Global Research Collaboration of Network Meta-Analysis: A Social Network Analysis. PLOS ONE. 2016;11(9):1-19.
  * [3] Viergever RF, Li K. Trends in global clinical trial registration: an analysis of numbers of registered clinical trials in different parts of the world from 2004 to 2013. BMJ Open. 2015;5:1-15.
  * [4] Jackson D, Riley R, White IR. Multivariate meta-analysis: Potential and promise. Statistics in Medicine. 2011;30:2481–2498.
  * [5] White IR, Barrett JK, Jackson D, et al . Consistency and inconsistency in network meta-analysis: model estimation using multivariate meta-regression. Res. Syn. Meth.. 2012;3:111-125.
  * [6] Tricco CA, Thomas MS, Veroniki AA, et al . Comparisons of Interventions for Preventing Falls in Older Adults. JAMA. 2017;318(17):1687-1699.
  * [7] Brunoni RA, Chaimani A, Moffa HA, et al . Repetitive Transcranial Magnetic Stimulation for the Acute Treatment of Major Depressive Episodes. JAMA Psychiatry. 2017;74(2):143-153.
  * [8] Mitra S, Florez DI, Tamayo EM, et al . Association of Placebo, Indomethacin, Ibuprofen, and Acetaminophen With Closure of Hemodynamically Significant Patent Ductus Arteriosus in Preterm Infants. JAMA. 2018;319(12):1221-1238.
  * [9] Rothenberg TJ. Hypothesis Testing in Linear Models when the Error Covariance Matrix is Nonscalar. Econometrica. 1984;52(4):827-842.
  * [10] Cordeiro GM, Cribari-Neto F. An Introduction to Bartlett Correction and Bias Reduction.  Springer-Verlag Berlin Heidelberg; 2014.
  * [11] Noma H, Nagashima K, Maruo K, et al . Bartlett-type corrections and bootstrap adjustments of likelihood‐based inference methods for network meta-analysis. Statistics in Medicine. 2018;37:1178-1190.
  * [12] Kojima M, Kubokawa T. Bartlett-type adjustments for hypothesis testing in linear models with general error covariance matrices. Journal of Multivariate Analysis. 2013;122:162-174.
  * [13] Aert R.C.M, Jackson D.. A new justification of the Hartung-Knapp method for random-effects meta-analysis based on weighted least squares regression. Res. Syn. Meth.. 2019;10:515-527.
  * [14] Bruce PM, Thomas T, Curt FD, et al . Health Outcomes Associated With Various Antihypertensive Therapies Used as First-Line Agents. JAMA. 2003;289(19):2534-2544.
  * [15] Stettler C, Wandel S, Allemann S, Kastrati A, et al . Outcomes associated with drug-eluting and bare-metal stents: a collaborative network meta-analysis. Lancet. 2007;370(937-948).
  * [16] Eliott WJ, Meyer PM. Incident diabetes in clinical trials of antihypertensive drugs: a network meta-analysis. Lancet. 2007;369:201-207.
  * [17] Lee SWH, Chan CKY, Chua SS, Chaiyakunapruk N. Comparative effectiveness of telemedicine strategies on type 2 diabetes management: A systematic review and network meta-analysis. Scientific Reports. 2017;7(12680):1-11.
  * [18] Sugasawa S, Noma H. A unified method for improved inference in random effects meta-analysis. Biostatistics. 2021;22(1):114-130.
  * [19] Welham SJ, Thompson R. Likelihood Ratio Tests for Fixed Model Terms using Residual Maximum Likelihood. J. R. Statist. Soc. B. 1997;59(3):701-714.
  * [20] Jennrich IR, Schluchter DM. Unbalanced Repeated-Measures Models with Structured Covariance Matrices. Biometrics. 1986;42(4):805-820.
  * [21] Kubokawa T. On Measuring Uncertainty Of Small Area Estimators With Higher Order Accuracy. J. Japan Statist. Soc.. 2011;41(2):93-119.
  * [22] Hardy RJ, SG Thompson. A likelihood approach to meta‐analysis with random effects. Statistics in Medicine. 1996;15(6):619-629.
  * [23] Sidik K, Jonkman NJ. A comparison of heterogeneity variance estimators in combining results of studies. Statistics in Medicine. 2007;26:1964–1981.
  * [24] Nagashima K, Noma H, Furukawa T. Prediction intervals for random-effects meta-analysis: A confidence distribution approach. Statistical Methods in Medical Research. 2019;28(6):1689-1702.
  * [25] Higgins JPT. Borrowing strength from external trials in a meta-analysis. Statistics in Medicine. 1996;15:2733-2749.
  * [26] Davidson R, MacKinnon GJ. Bootstrap Tests: How Many Bootstraps?. Econometric Reviews. 2000;19(1):55-68.
  * [27] Palmerini T, Benedetto U, Bacchi-Reggiani L, et al . Mortality in patients treated with extended duration dual antiplatelet therapy after drug-eluting stent implantation: a pairwise and Bayesian network meta-analysis of randomised trials. Lancet. 2015;385:2371-82.