Improved Confidence Regions in Meta-analysis of Diagnostic Test Accuracy

TL;DR

This paper introduces a new, simple, and accurate confidence region method for meta-analysis of diagnostic test accuracy, addressing the limitations of standard methods that underestimate errors.

Contribution

The paper proposes an asymptotic expansion-based confidence region that improves inference accuracy without complex computations like bootstrap or Monte Carlo methods.

Findings

The new method provides more accurate coverage probabilities.

Simulation studies confirm its effectiveness over standard methods.

Applied to alcohol screening test data, it yields more reliable confidence regions.

Abstract

Meta-analyses of diagnostic test accuracy (DTA) studies have been gathering attention in research in clinical epidemiology and health technology development, and bivariate random-effects model is becoming a standard tool. However, standard inference methods usually underestimate statistical errors and possibly provide highly overconfident results under realistic situations since they ignore the variability in the estimation of variance parameters. To overcome the difficulty, a new improved inference method, namely, an accurate confidence region for the meta-analysis of DTA, by asymptotically expanding the coverage probability of the standard confidence region. The advantage of the proposed confidence region is that it holds a relatively simple expression and does not require any repeated calculations such as Bootstrap or Monte Carlo methods to compute the region, thereby the proposed method can be easily carried out in practical applications. The effectiveness of the proposed method is demonstrated through simulation studies and an application to meta-analysis of screening test accuracy for alcohol problems.

Tables

Table 1. Table 1: Averaged values of the heterogeneity measure I 2 ( % ) I^{2}(\%) based on 1000 replications.

		${\tau }^2$
$n$		0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8
8		27.2	63.4	75.1	81.4	85.0	87.2	89.1	90.8
16		49.0	74.9	83.1	87.4	89.7	91.6	92.7	93.7
24		56.2	78.0	85.4	88.9	91.2	92.7	93.7	94.6

Equations

$$y_i|{\mu }_i\sim N_2({\mu }_i,S_i),{\mu }_i\sim N_2(\beta ,\Sigma ),i=1,\dots ,n,$$

$$\left\{\beta :(\hat{\beta }(\hat{\Sigma })-\beta {)}^{t}V(\hat{\Sigma }{)}^{-1}(\hat{\beta }(\hat{\Sigma })-\beta )\le {\chi }_2^2(\alpha )\right\},$$

$${\hat{\Sigma }}_0=\frac{1}{n}\sum _{i=1}^n\left\{(y_i-X_i{\hat{\beta }}^{\mathrm{OLS}})(y_i-X_i{\hat{\beta }}^{\mathrm{OLS}}{)}^t-S_i\right\},$$

$$\hat{\beta }(\Sigma )=\{X^t(I_n\otimes \Sigma +S{)}^{-1}X\}{X}^t(I_n\otimes \Sigma +S{)}^{-1}y.$$

$$\bigg{\{}\beta:({\widehat{\beta}}({\widehat{\Sigma}})-{\beta})^{t}V({\widehat{\Sigma}})^{-1}({\widehat{\beta}}({\widehat{\Sigma}})-{\beta})\leq x(1+h({\widehat{\Sigma}}))\bigg{\}},$$

$$\displaystyle P\big{\{}({\widehat{\beta}}({\widehat{\Sigma}})-{\beta})^{t}V({\widehat{\Sigma}})^{-1}({\widehat{\beta}}({\widehat{\Sigma}})-{\beta})\leq x(1+h)\big{\}}$$

$$=F_k(x)+hx{f}_k(x)+\left(\frac{B_1}{4}-\frac{B_2}{2}+2B_3\right)f_{k+2}(x)$$

$$-\left(\frac{B_1}{4}+\frac{B_2}{2}\right)f_{k+4}(x)+O(n^{-3/2}),$$

$$B_{1}={\rm E}\Big{[}{\rm tr}(K({\widehat{\Sigma}},\Sigma))^{2}\Big{]},\ \ \ \ B_{2}={\rm tr}\Big{(}{\rm E}[K({\widehat{\Sigma}},\Sigma)^{2}]\Big{)},\ \ \ \ B_{3}={\rm tr}\Big{(}{\rm E}[K({\widehat{\Sigma}},\Sigma)]\Big{)},$$

$$hx{f}_k(x)+\left(\frac{B_1}{4}-\frac{B_2}{2}+2B_3\right)f_{k+2}(x)-\left(\frac{B_1}{4}+\frac{B_2}{2}\right)f_{k+4}(x)=0.$$

$$h(\Sigma )=\frac{1}{k}\left(\frac{B_1}{4}-\frac{B_2}{2}+2B_3\right)-\frac{x}{k(k+2)}\left(\frac{B_1}{4}+\frac{B_2}{2}\right).$$

$$\begin{split}B_{1}^{\ast}=&\frac{2}{n^{2}}\sum_{i=1}^{n}\sum_{j=1}^{n}\sum_{k=1}^{n}{\rm tr}\Big{(}VU_{jik}VU_{kij}\Big{)},\\ B_{2}^{\ast}=&\frac{1}{n^{2}}\sum_{i=1}^{n}{\rm tr}\Big{(}V\sum_{j=1}^{n}U_{jij}^{2}\Big{)}+\frac{1}{n^{2}}\sum_{i=1}^{n}\sum_{j=1}^{n}\sum_{k=1}^{n}{\rm tr}^{2}\Big{(}VU_{ijk}\Big{)},\\ B_{3}^{\ast}=&B_{2}-\frac{1}{n^{2}}\sum_{i=1}^{n}\sum_{j=1}^{n}{\rm tr}\Big{(}VU_{iji}D_{j}D_{i}^{-1}\Big{)}\\ &-\frac{1}{n^{2}}\sum_{i=1}^{n}\sum_{j=1}^{n}{\rm tr}\Big{(}D_{i}^{-1}D_{j}\Big{)}{\rm tr}\Big{(}VU_{iji}\Big{)},\end{split}$$

$$\mathrm{E}[Py(\tilde{\beta }-\beta {)}^t]V^{-1}$$

$$\hat{\beta }(\hat{\Sigma }(y+XT),y+XT)-\tilde{\beta }(\Sigma ,y+XT)=\hat{\beta }(\hat{\Sigma }(y),y)-\tilde{\beta }(\Sigma ,y),$$

$$\mathrm{E}[\hat{V}]-V=-\mathrm{E}[(\hat{\beta }-\tilde{\beta })(\hat{\beta }-\tilde{\beta }{)}^t]+O(n^{-5/2}).$$

$$\hat{V}-V=$$

$$=$$

$$+(\hat{V}-V)X^t({\hat{\Sigma }}_n+S{)}^{-1}({\hat{\Sigma }}_n-{\Sigma }_n)({\Sigma }_n+S{)}^{-1}XV$$

$$\equiv$$

$$I_1$$

$$-V{X}^t({\Sigma }_n+S{)}^{-1}({\hat{\Sigma }}_n-{\Sigma }_n)({\Sigma }_n+S{)}^{-1}({\hat{\Sigma }}_n-{\Sigma }_n)({\Sigma }_n+S{)}^{-1}XV$$

$$+O_p(n^{-5/2}),$$

$$I_2=$$

$$\times ({\Sigma }_n+S{)}^{-1}({\hat{\Sigma }}_n-{\Sigma }_n)({\Sigma }_n+S{)}^{-1}XV+O_p(n^{-5/2}).$$

$$\mathrm{E}[$$

$$=$$

$$-VR({\hat{\Sigma }}_n-{\Sigma }_n)({\Sigma }_n+S{)}^{-1}({\hat{\Sigma }}_n-{\Sigma }_n)R^{t}V]+O(n^{-5/2})$$

$$\begin{array}{rl}=& \mathrm{E}[VR({\hat{\Sigma }}_n-{\Sigma }_n)R^{t}VR({\hat{\Sigma }}_n-{\Sigma }_n)R^{t}V\\ & -VR({\hat{\Sigma }}_n-{\Sigma }_n)({\Sigma }_n+S{)}^{-1}({\hat{\Sigma }}_n-{\Sigma }_n)R^{t}V]+O(n^{-5/2}),\end{array}$$

$$\hat{\beta }-\tilde{\beta }=$$

$$+V{X}^t\{({\hat{\Sigma }}_n+S{)}^{-1}-({\Sigma }_n+S{)}^{-1}\}(y-X\beta )$$

$$=$$

$$J_1=$$

$$J_2=$$

$$\mathrm{E}[J_1{J}_1^t]$$

$$\mathrm{E}[J_2{J}_2^t]$$

$$\mathrm{E}[J_1{J}_2^t]$$

$$-\mathrm{E}[$$

$$=$$

Full text

Improved Confidence Regions in Meta-analysis of Diagnostic Test Accuracy

Tsubasa Ito1 and Shonosuke Sugasawa2

1M&D Data Science Center, Tokyo Medical and Dental University

2Center for Spatial Information Science, The University of Tokyo

Abstract

Meta-analyses of diagnostic test accuracy (DTA) studies have been gathering attention in research in clinical epidemiology and health technology development, and bivariate random-effects model is becoming a standard tool. However, standard inference methods usually underestimate statistical errors and possibly provide highly overconfident results under realistic situations since they ignore the variability in the estimation of variance parameters. To overcome the difficulty, a new improved inference method, namely, an accurate confidence region for the meta-analysis of DTA, by asymptotically expanding the coverage probability of the standard confidence region. The advantage of the proposed confidence region is that it holds a relatively simple expression and does not require any repeated calculations such as Bootstrap or Monte Carlo methods to compute the region, thereby the proposed method can be easily carried out in practical applications. The effectiveness of the proposed method is demonstrated through simulation studies and an application to meta-analysis of screening test accuracy for alcohol problems.

Key words: Asymptotic expansion; Bias correction; Confidence region; random-effects model

Introduction

Evidence synthesis methods have been gathering attention in diagnostic test accuracy (DTA) studies in clinical epidemiology and health technology development (Leeflang et al., 2008). In this meta-analysis, the summary statistics in each study are two primary correlated outcomes of diagnostic, sensitivity and false positive rate ($1-{\rm specificity}$), and we are typically interested in summary receiver operating characteristic curve. Moreover, DTA from different sources for studies are generally heterogeneous due to various factors, which should be adequately addressed to avoid underestimation of statistical errors and misleading conclusions (Higgins and Green, 2011). Due to the potential correlations between two summary measures and potential heterogeneity, the bivariate random-effects models is adopted as the standard method for the meta-analysis (Reitsma et al., 2005; Harbord et al., 2007).

In the bivariate random-effects meta-analyses, standard inference methods depend on large sample approximations for the number of studies synthesized, for example the extended DerSimonian-Laird methods (Chen et al., 2012; Jackson et al., 2010, 2013) and restricted maximum likelihood (REML) estimation (Reitsma et al., 2005; Jackson et al., 2011), but the numbers of trials are often moderate or small in practice. In this situation, validity of the inference methods can be violated, which may lead over-confidence results, that is, coverage probabilities of the confidence regions or intervals cannot retain their nominal confidence levels and also the type-I error probabilities of the corresponding tests can be inflated. Such problem with random-effects models was well recognized in the context of both univariate and multivariate meta-analysis, even when the models are completely specified (Veroniki et al., 2019). Recently, several refined methods have been proposed to improve confidence intervals in multivariate meta-analysis. For example, Noma et al. (2018) developed improved confidence intervals in network meta-analysis using Bartlett-type corrections, and Noma et al. (2020) and Sugasawa and Noma (2020) developed a unified method for computing accurate confidence intervals and regions in general random-effects meta-analysis. However, these methods require computationally very intensive methods based on Monte Carlo or Bootstrap methods. Also these methods considered confidence intervals or regions by inverting statistical hypothesis tests, thereby feasible ways to construct confidence regions are not necessarily obvious. On the other hand, there are a few analytical approaches to improve the standard approaches. Noma (2011) and Guolo (2012) considered higher order likelihood inference in the univariate meta-analysis, which cannot be directly applicable to more complicated multivariate meta-analysis. As more general approaches, Zucker et al. (2000) proposed an improved likelihood test in general linear mixed models through asymptotic expansions of the (restricted) maximum likelihood estimators, but the results include tedious algebraic expressions and are not useful in practice.

In this paper, we propose an improved confidence region for the bivariate random-effects meta-analysis for DTA, which does not require any repeated calculation methods and has relatively simple analytical expressions, thereby the proposed method could be easily employed in practical applications. The key mathematical tool is the distributional properties between the ordinary least squares estimator and residuals, and define a class of estimators of variance parameters in random-effects models. Then, we find a relatively simple formula for asymptotic approximation of the coverage probability of the crude Wald-type confidence intervals and regions, and construct a second order accurate confidence region. We carry out extensive simulation studies to compare the performance of the proposed confidence region with that of the standard REML method, and demonstrate that the proposed method shows quite reasonable empirical coverage than REML while the computational cost in both methods are almost identical. We also demonstrate the proposed method through an application to meta-analysis of screening test accuracy for alcohol problems.

This paper is set out as follows. In Section 2, we describe the proposed confidence region under bivariate random-effects models. In Section 3, we numerically demonstrate the proposed confidence region together with existing methods through extensive simulation studies and an application with real dataset. We conclude with a short discussion in Section 4. R code implementing the proposed method is available at GitHub repository (https://github.com/sshonosuke/CCR-BMA).

Improved Confidence Regions in Meta-analysis for Diagnostic Test Accuracy

Bivariate random-effects models and confidence region

There has been increasing interest in systematic reviews and meta-analyses of data from diagnostic accuracy studies. For this purpose, a bivariate random-effect model (Reitsma et al., 2005; Harbord et al., 2007) is widely used. Following Reitsma et al. (2005), we define $\mu_{Ai}$ and $\mu_{Bi}$ as the logit-transformed true sensitivity and specificity, respectively, in the $i$th study. Let $y_{Ai}$ and $y_{Bi}$ be the observed logit-transformed sensitivity and specificity, and $s_{Ai}$ and $s_{Bi}$ are associated standard errors. The bivariate model assumes that $\mu_{i}=(\mu_{Ai},\mu_{Bi})^{t}$ and $y_{i}=(y_{Ai},y_{Bi})^{t}$ follow bivariate normal distributions:

[TABLE]

where $\beta=(\beta_{A},\beta_{B})^{t}$ is a vector of the average logit-transformed sensitivity and specificity, and $S_{i}={\rm diag}(s_{Ai},s_{Bi})$. Note that there is no correlation between $y_{Ai}$ and $y_{Bi}$ given $\mu_{i}$ since sensitivity and specificity are calculated based on individuals identified as positive and negative, respectively. Here ${\Sigma}$ is unstructured, so that it allows correlation between $\mu_{Ai}$ and $\mu_{Bi}$. Let $y=(y_{1}^{t},\ldots,y_{n}^{t})^{t}\in{\mathbb{R}}^{2n}$, $X=(X_{1},\ldots,X_{n})^{t}$ with $X_{i}=I_{2}$ and $S={\rm diag}(S_{1},\ldots,S_{n})$. Then, the model (1) is equivalent to $y\sim N(X\beta,I_{n}\otimes\Sigma+S)$.

Our primary interest is a confidence region of $\beta$. Hence, the variance-covariance matrix $\Sigma$ is a nuisance parameter. These parameters are typically estimated via (restricted) maximum likelihood methods based on the model assumption (1). For summarizing the results of the meta-analysis, we typically employ confidence region of $\beta$ rather than separate confidence intervals since sensitivity and specificity could be highly correlated. Reitsma et al. (2005) suggested the $100(1-\alpha)\%$ confidence region for $\beta$ as the interior points of the ellipse defined as

[TABLE]

where ${\widehat{\beta}}(\Sigma)$ is the generalized least squares estimator of $\beta$ and ${\widehat{\Sigma}}$ is the restricted maximum likelihood estimator of $\Sigma$, $V({\Sigma})=\{X^{t}(I_{n}\otimes\Sigma+S)X\}^{-1}$ is the variance-covariance matrix of ${\widehat{\beta}}$, ${\widehat{\Sigma}}$ is the restricted maximum likelihood estimator and $\chi_{2}^{2}(\alpha)$ is the upper $100\alpha\%$ point of the $\chi^{2}$ distribution with $2$ degrees of freedom. The joint confidence region (2) is approximately valid, that is, the coverage error converges to $1-\alpha$ as the number of studies $n$ goes to infinity. However, when $n$ is not sufficiently large, the coverage error is not negligible, and the region (2) would under-cover the true $\beta$.

Improved confidence region

In this work, we derive an improved confidence region whose coverage error is $o(n^{-1})$, which has higher order accuracy than the standard confidence region (2). The main idea is to derive an approximation formula of the coverage probabilities of the confidence region of the form (2) with a certain class of estimators for $\Sigma$, and derive an improved confidence region in an analytical form.

We consider a class of estimators ${\widehat{\Sigma}}(y)$ satisfying the following conditions:

• (C1)

${\widehat{\Sigma}}$ is an even function of $y$ and translation invariant, that is, ${\widehat{\Sigma}}(y)={\widehat{\Sigma}}(-y)$, and ${\widehat{\Sigma}}(y+c)={\widehat{\Sigma}}(y)$ for any $c\in{\mathbb{R}}^{2n}$.

• (C2)

${\widehat{\Sigma}}$ is $\sqrt{n}$-consistent and ${\widehat{\Sigma}}$ is second-order unbiased, namely ${\widehat{\Sigma}}-\Sigma=O(n^{-1/2})$ and ${\rm E}[{\widehat{\Sigma}}]=\Sigma+o(n^{-1})$.

• (C3)

${\widehat{\Sigma}}$ is a function of $Py$ with $P=I_{2n}-X(X^{t}X)^{-1}X^{t}$.

The first condition (C1) is typically satisfied by typical estimators including (restricted) maximum likelihood estimator and moment-based estimators. The $\sqrt{n}$-consistency in (C2) is also a standard condition, but second order unbiasedness of ${\widehat{\psi}}$ is not always satisfied. For example, the maximum likelihood (ML) estimator does not necessarily hold the property. The condition (C3) requires that the estimator should be function of residuals based on ordinary least squares estimator of $\beta$, which is a key assumption in constructing the proposed confidence region. The condition (C3) enables us to get a relatively simple form of the corrected confidence region. Note that the typical estimators (e.g. REML) does not satisfy the condition (C3). As a specific estimator satisfying all the above conditions, we employ the following moment-based estimator:

[TABLE]

where ${\widehat{\beta}}^{{\rm OLS}}=(X^{t}X)^{-1}X^{t}y$ is the ordinary least squares estimator. Since this estimator is not second-order unbiased, let ${\widehat{\Sigma}}$ be a bias corrected version, that is, ${\widehat{\Sigma}}={\widehat{\Sigma}}_{0}-{\rm Bias}_{{\widehat{\Sigma}}_{0}}({\widehat{\Sigma}})$ with ${\rm Bias}_{{\widehat{\Sigma}}_{0}}({\Sigma})=-n^{-2}\sum_{i=1}^{n}(\Sigma+S_{i})$, which satisfies all the conditions (C1)$\sim$(C3). We also note that given the estimator of $\Sigma$, the parameter $\beta$ can be estimated via the generalized least squares estimator given by

[TABLE]

In order to improve the coverage accuracy of the confidence region (2), we consider a class of confidence regions of the form

[TABLE]

where $h(\cdot)$ is a function with order $O(n^{-1})$. When $h(\Sigma)=0$ and $x=\chi_{2}^{2}(\alpha)$, the confidence region (3) reduces to (2), thereby the function $h$ can be regarded as an adjustment function to achieve reasonable coverage properties. If ${\widehat{\Sigma}}$ satisfies the conditions (C1)$\sim$(C3), the approximation formula of coverage probability of the confidence region (3) can be obtained in a relatively simple form, as summarized in the following theorem.

Theorem 1**.**

Suppose that ${\widehat{\Sigma}}$ satisfies the conditions (C1)$\sim$(C3), and $h\equiv h(\Sigma)$ is a function with order $O(n^{-1})$. Then, it follows that

[TABLE]

where $F_{k}(\cdot)$ and $f_{k}(\cdot)$ are the cumulative distribution and density function of the chi-squared distribution with degrees of freedom $k$, respectively, and $B_{1},B_{2}$ and $B_{3}$ are $O(n^{-1})$ quantities given by

[TABLE]

with K({\widehat{\Sigma}},\Sigma)=\big{\{}V({\widehat{\Sigma}})-V(\Sigma)\big{\}}V(\Sigma)^{-1}.

From Theorem 1, It turned out that the coverage probability of the confidence region (3) is a simple functional of the djustment function $h(\cdot)$. Hence, to achieve higher accuracy of the confidence region, it suffices to choose $h(\cdot)$ such that

[TABLE]

Since $f_{k+2}(x)/f_{k}(x)=x/k$ and $f_{k+4}(x)/f_{k}(x)=x^{2}/k(k+2)$, the solution with respect to $h$ is given by

[TABLE]

We also note that $h({\widehat{\Sigma}})=h(\Sigma)+o_{p}(n^{-1})$ since $h(\cdot)=O(n^{-1})$. Then, the confidence region given in (3) with $h(\cdot)$ given in (5) holds the second-order accuracy as shown in the following theorem.

Theorem 2**.**

Let ${\rm CCI}_{\alpha}$ be the confided region of the form (3) with $h(\cdot)$ given in (5). Then, it follows that $P(\beta\in{\rm CCI}_{\alpha})=1-\alpha+o(n^{-1})$.

It is notable that the derived confidence region has analytical expressions, so that it does not require any computationally intensive methods such as bootstrap and Monte Carlo integration as used in Sugasawa and Noma (2020) and Noma et al. (2020). For practical implementation, we need to obtain the expressions of $B_{1}$, $B_{2}$ and $B_{3}$ given in (4). We here provide approximation formulas. we can obtain $B_{\ell}^{\ast},\ \ell=1,2,3$ which satisfy $B_{\ell}^{\ast}=B_{\ell}+o(n^{-1})$, where

[TABLE]

for $U_{ijk}=D_{i}^{-1}D_{j}D_{k}^{-1}$ and $D_{i}=\Sigma+S_{i}$. The detailed derivation is given in the Supplementary Material. Note that using $B_{\ell}^{\ast}$ instead of $B_{\ell}$ in the derived confidence region does not change the coverage accuracy shown in Theorem 2 since the difference between $B_{\ell}^{\ast}$ and $B_{\ell}$ is only $o(n^{-1})$.

Numerical Studies

Simulation study

We carried out extensive simulation studies to assess the finite sample performance of the proposed confidence region (3) together with the approximate confidence region (2) by Reitsma et al. (2005). In this study, we do not consider possible competitors by Noma et al. (2020); Sugasawa and Noma (2020) due to two reasons; the coverage performance has been already confirmed in their papers, and calculation of sizes of their confidence regions are so intensive that it is not feasible to repeatedly calculate them in our simulation study. Hence, the following simulation study is supposed to compare the performance of the proposed and standard methods, both of which have almost the same computational time.

In the model (1), we set $\beta=(0,0)$ and $\Sigma_{11}=\Sigma_{22}=\tau^{2}$ and $\Sigma_{12}=\Sigma_{21}=\tau^{2}\rho$. We considered 8 scenarios of the between study variances $\tau^{2}\in\{0.1,0.2,\ldots,0.8\}$ and 5 scenarios of the between study correlations $\rho\in\{0,0.2,\ldots,0.8\}$. Following, Jackson and Riley (2014), for each simulation, two within-study variances $s_{Ai}$ and $s_{Bi}$ were simulated from a scaled chi-squared distribution with 1 degree of freedom, multiplied by 0.25, and truncated to lie within the interval $[0.009,0.6]$, so the expected values of the variance is $0.20$. We changed the number of studies $n$ over 8,16 and 24, and set the nominal level $\alpha$ to $0.05$. Based on 1000 replications, we evaluated empirical coverage probabilities of 95% confidence regions of the true parameters vector $\beta$ obtained from the proposed corrected (CCR) method as well as the standard naive (NCR) method. For simplicity, we evaluated coverage rates assessing rejection rates of the test of null hypothesis for the true parameters. Since areas of the corrected confidence region is approximately $1+h({\widehat{\Sigma}})$ times larger than those of naive ones, we also computed median values of $h({\widehat{\Sigma}})$ among 1000 replications. To see the degree of heterogeneity depending on $n$ and $\tau^{2}$, we computed heterogeneity measure given by $I^{2}=\tau^{2}/(Q+\tau^{2})$ with $Q=(n-1)\sum_{i=1}^{n}w_{i}/\{(\sum_{i=1}^{n}w_{i})^{2}-\sum_{i=1}^{n}w_{i}^{2}\}$ and $w_{i}=s_{Ai}^{-1}$ in each iteration, which were averaged over 1000 replications. Note that $I^{2}\in(0,1)$ and lager value of $I^{2}$ indicates more significant heterogeneity in the data.

The averaged values of $I^{2}$ are reported in Table 1, which indicates that our simulation scenarios contain a wide range of heterogeneity. The obtained coverage probabilities and the median values of $h({\widehat{\Sigma}})$ are shown in Figures 1 and 2, respectively. From Figure 1, it is observed that the simulated coverage probabilities of the standard NCR seriously smaller than the nominal level ($95\%$), especially in the case with the small number of studies ($n=8$), possibly because of the naive approximation in (2). On the other hand, the proposed CCR provides considerably better performance than NCR as the coverage probabilities are relatively close to the nominal level. Although the coverage probability of CCR tend to be larger than the nominal level when $\tau$ is small and/or $\rho$ is large, such a conservative property would be much more desirable than the over-confident property that NCR shows. From Figure 2, we can see that the area of CCR is much larger than that of NCR since CCR takes account of additional variability due to the estimation of the variance-covariance matrix, so it is quite reasonable that $h({\widehat{\Sigma}})$ decreases as $n$ increases. Moreover, we can also observe that the areas of CCR decreases as $\tau^{2}$ increases and increases as $\rho$ increases, which are consistent to the results of the overage probabilities shown in Figure 1.

Example: screening test accuracy for alcohol problems

Here we provide a re-analysis of the dataset given in Kriston et al. (2008), including $n=14$ studies regarding a short screening test for alcohol problems. Following Reitsma et al. (2005), we used logit-transformed values of sensitivity and specificity, denoted by $y_{Ai}$ and $y_{Bi}$, respectively, and associated standard errors $s_{Ai}$ and $s_{Bi}$. For the bivariate summary data, we first fitted the bivariate models (1) using the restricted maximum likelihood method and found that $\widehat{\rho}=0.854$, and the heterogeneity measure $I^{2}$ for sensitivity and specificity are respectively given by $94.8\%$ and $98.9\%$, so there seems notable heterogeneity in the data. We then computed $95\%$ CRs of $\beta$ based on NCR (2) given in Reitsma et al. (2005) and the proposed CCR. Following Reitsma et al. (2005), the obtained two CRs of $\beta$ were transformed to the scale $({\rm logit}(\beta_{A}),1-{\rm logit}(\beta_{B}))$, where ${\rm logit}(\beta_{A})$ and $1-{\rm logit}(\beta_{B})$ are the sensitivity and false positive rate, respectively. The obtained two CRs are presented in Figure 3 with a plot of the observed data, summary points ${\widehat{\beta}}$, and the summary receiver operating curve. The approximate CR is smaller than the proposed CR, which may indicate that the approximation method underestimates the variability of estimating nuisance variance parameters. In Figure 3, we also reported the confidence region based on Sugasawa and Noma (2020) using Monte Carlo simulation to compute accurate $p$-values of likelihood ratio statistics. The two regions based on the proposed method and Sugasawa and Noma (2020) are slightly different but both are clearly wider than the naive confidence region. On the other hand, the computation time of the proposed method was less than 1 second while the inference method by Sugasawa and Noma (2020) took more than 12 hours, where the program was run on a PC with a 3 GHz 8-Core Intel Xeon E5 8 Core Processor with approximately 16GB RAM.

Discussion

In this paper, we presented an improved confidence region for random effects meta-analysis for diagnostic test accuracy without using repeated calculations such as Monte Carlo or Bootstrap methods. The proposed confidence region has relatively simple form and they are shown to have second order accurate coverage probability while the standard inference methods (e.g. REML) have significant coverage errors. In simulation studies, we demonstrated that possible under-coverage properties of the standard methods under the small number of studies to be synthesized while the proposed method provides reasonable coverage properties.

A possible limitation of the proposed method might be that the coverage accuracy still depends on the number of studies. On the other hand, inference methods that does not rely on large sample approximation have been recently proposed (e.g. Noma et al., 2020; Sugasawa and Noma, 2020), which are computationally intensive, so they would not be necessarily practical. Then, the proposed method would be regarded as a reasonable compromise between methods with exact empirical coverage and computational efficiency.

Acknowledgments

This research was supported by Japan Society of Promotion of Science KAKENHI (grant number: 18K12757).

Key lemmas

We first introduce lemmas which play important roles in the proof of Theorems 1. The first lemma is used for deriving the conditional distribution of ${\widehat{\beta}}$.

Lemma S1**.**

Under the conditions (C1)-(C3)* given in the main document, ${\widetilde{\beta}}$ is independent of $Py$ for $P=I_{n}-X(X^{t}X)^{-1}X^{t}$. Also, ${\widehat{\beta}}-{\widetilde{\beta}}$ is a function of $Py$, and independent of ${\widetilde{\beta}}$.*

Proof.

Let ${\widetilde{\varepsilon}}=y-X{\beta}$, which is distributed as $N(0,\Sigma_{n}+S)$. Since ${\widetilde{\beta}}-{\beta}=VX^{t}(\Sigma_{n}+S)^{-1}{\widetilde{\varepsilon}}$, it holds that

[TABLE]

Since $V$ is a full-rank matrix, we have ${\rm E}[Py({\widetilde{\beta}}-{\beta})^{t}]=0$, that is the covariance of $Py$ and ${\widetilde{\beta}}$ is [math], which implies that ${\widetilde{\beta}}$ is independent of $Py$ from the normality assumption. Now, we write ${\widetilde{\beta}}$ as ${\widetilde{\beta}}({\Sigma},y)$ and ${\widehat{\beta}}$ as ${\widehat{\beta}}({\widehat{\Sigma}}(y),y)$. Since ${\widetilde{\beta}}({\Sigma},y+XT)={\widetilde{\beta}}({\Sigma},y)+T$ and ${\widehat{\beta}}({\widehat{\Sigma}}(y+XT),y+XT)={\widehat{\beta}}({\widehat{\Sigma}}(y),y)+T$ from (C3), we have

[TABLE]

which implies that ${\widehat{\beta}}-{\widetilde{\beta}}$ is invariance with respect to the translation $y\to y+XT$. Moreover, $Py$ is maximal invariant with respect to the translation $y\to y+XT$ since $P(y+XT)=Py$ and $Py_{1}=Py_{2}$ implies that $y_{1}=y_{2}+XT^{\prime}$ for $T^{\prime}=(X^{t}X)^{-1}X^{t}(y_{1}-y_{2})$. Then, ${\widehat{\beta}}-{\widetilde{\beta}}$ is a function of $Py$ from Theorem 2 in Berger (1985), p.403. ∎

In the next lemma, we show the first order bias of the plug-in estimator $\widehat{V}$ is approximately the same as the negative covariance of ${\widehat{\beta}}-{\widetilde{\beta}}$.

Lemma S2**.**

Under the conditions (C1)-(C3), it holds that

[TABLE]

Proof.

We will show the Lemma by directly comparing both sides of the equation in the Lemma. Noting that $V=\{X^{t}({\Sigma}_{n}+S)X\}^{-1}$ and $A^{-1}-B^{-1}=-A^{-1}(A-B)B^{-1}$ for some non-singular matrices $A$ and $B$, we have

[TABLE]

Since $V=O(n^{-1})$ and $\widehat{V}-V=O_{p}(n^{-1/2})$ from the condition (C2), we have

[TABLE]

and

[TABLE]

Then, for $R=X^{t}({\Sigma}_{n}+S)^{-1}$ we have

[TABLE]

where the last equality holds since ${\widehat{\Sigma}}$ is a second-order unbiased estimator of ${\Sigma}$.

Next, we evaluate the first term of the right side of the equation in the Lemma. We can write ${\widehat{\beta}}-{\widetilde{\beta}}$ as

[TABLE]

In order to approximate the covariance of ${\widehat{\beta}}-{\widetilde{\beta}}$ up to the order $O(n^{-5/2})$, we expand $J_{1}$ and $J_{2}$ as

[TABLE]

The straightforward calculation shows that

[TABLE]

thereby we have

[TABLE]

which has the same expression as (S1). ∎

Proof of Theorem 1

From Lemma S1, the conditional distribution of ${\widehat{\beta}}-{\beta}$ given $Py$ is $N_{2}({\widehat{\beta}}-{\widetilde{\beta}},V)$. Let $w=V^{-1/2}\{({\widehat{\beta}}-{\beta})-({\widehat{\beta}}-{\widetilde{\beta}})\}$. It is noted that $\widehat{V}-V=O_{p}(n^{-3/2})$. Then, the conditional distribution of $w$ given $Py$ is $w\sim N_{k}(0,I_{k})$, and the Mahalanobis’ distance is approximated via Taylor series expansion as

[TABLE]

where

[TABLE]

From (S2), the characteristic function $\varphi(t)={\rm E}[\exp\{it({\widehat{\beta}}-{\beta})^{t}\widehat{V}^{-1}({\widehat{\beta}}-{\beta})\}]$ is approximated as

[TABLE]

because $G_{1}=O_{p}(n^{-1/2})$, $g_{2}=O_{p}(n^{-1/2})$ and $g_{3}=O_{p}(n^{-1})$. From the law of iterated expectations and the conditional normality of $w$, the above equation reduces to

[TABLE]

For some deterministic matrix $A$ and $w\sim N_{k}(0,I_{k})$, it holds that

[TABLE]

Using these equalities, from the law of iterated expectations, we have

[TABLE]

For notational simplicity, let $J={\rm E}[{\rm tr}^{2}(G_{1})]+2{\rm tr}({\rm E}[G_{1}^{2}])$. Let $s=(1-2it)^{-1}$, or $it=(s-1)/(2s)$. Then, $(1-2it)^{k/2}\varphi(t)-1$ can be written as

[TABLE]

We shall evaluate the moments in (S3). First, $G_{1}$ can be expressed as

[TABLE]

thereby it holds that

[TABLE]

for $K=V^{-1/2}(\widehat{V}-V)V^{-1/2}$. Noting that the first term in (S4) is $O_{p}(n^{-1/2})$ and the second term is $O(n^{-1})$, we can expand $G_{1}^{2}$ and ${\rm tr}^{2}(G_{1})$ as

[TABLE]

which lead to ${\rm E}[G_{1}^{2}]={\rm E}[K^{2}]+O(n^{-3/2})$ and ${\rm E}[{\rm tr}^{2}(G_{1})]={\rm E}[{\rm tr}^{2}(K)]+O(n^{-3/2})$. Thus,

[TABLE]

It can be also observed that

[TABLE]

Then, from Lemma S2 we have

[TABLE]

Combining (S5), (S6) and (S7), we can see that the characteristic function of $({\widehat{\beta}}-{\widetilde{\beta}})^{t}\widehat{V}^{-1}({\widehat{\beta}}-{\widetilde{\beta}})$ can be written as

[TABLE]

for $B_{1},B_{2}$ and $B_{3}$ are defined in the main document. From the fact that the characteristic function of the chi-squared distribution with degrees of freedom $k+2h$ is given by $(1-2it)^{-k/2-h}=(1-2it)^{-k/2}s^{h}$, it follows that the asymptotic expansion of the cumulative distribution function of $({\widehat{\beta}}-{\beta})^{t}\widehat{V}^{-1}({\widehat{\beta}}-{\beta})$ is

[TABLE]

where $F_{k}(x)$ is the cumulative distribution function of the chi-squared distribution with degrees of freedom $k$. Note that $F_{k+r-2}(x)-F_{k+r}(x)=2f_{k+r}(x)$, where $f_{k}(x)$ is the density function of the chi-squared distribution with degrees of freedom $k$. Then, it holds that

[TABLE]

thereby, for a function $h=h({\Sigma})$ with order $O(n^{-1})$, we have

[TABLE]

which completes the proof.

Derivation of the equation (6)

We write functions given in Section 2 as functions of ${\Sigma}$ since the unknown parameter is ${\Sigma}$ in this example. For $V=(\sum_{i=1}^{n}D_{i}^{-1})^{-1}$ and $D_{i}={\Sigma}+S_{i}$, $\widehat{V}-V$ can be expanded as

[TABLE]

Since the first term on the right side of the above equation is of order $O_{p}(n^{-3/2})$, we only need to consider this term to derive the expressions given in (6).

At first, we evaluate $B_{1}^{\ast}$. It is noted that we have

[TABLE]

where $u_{j}$ are independently distributed as the standard normal distribution. Then, we have

[TABLE]

Next, we evaluate $B_{2}^{\ast}$. It is noted that we have

[TABLE]

and that ${\widehat{\Sigma}}-{\Sigma}$ can be written as

[TABLE]

Then, for $u_{i}$ for $i=1,\ldots,n$ which are independently distributed as the multivariate standard normal distribution, it holds that for $\ell,m=1,\ldots,n$,

[TABLE]

for $L_{i\ell m}=D_{i}^{1/2}D_{\ell}^{-1}VD_{m}^{-1}D_{i}^{1/2}$. Then, we have

[TABLE]

Finally, we evaluate $B_{3}^{\ast}$. From the equation (S1), for $V=(\sum_{i=1}^{n}D_{i}^{-1})^{-1}$ we have

[TABLE]

The trace of the first term in the above equation is exactly the same with $B_{2}^{\ast}$ and is given in (S9). To evaluate the second term, it is noted that

[TABLE]

Then, the trace of the second term in the above equation is given by

[TABLE]

Equation (S8), (S9) and (S10) lead to the expression given in (4) in the main document.