An one-factor copula mixed model for joint meta-analysis of multiple diagnostic tests

As the meta-analysis of more than one diagnostic tests can impact clinical decision making and patient health, there is an increasing body of research in models and methods for meta-analysis of studies comparing multiple diagnostic tests. The application of the existing models to compare the accuracy of three or more tests suffers from the curse of multi-dimensionality, i.e., either the number of model parameters increase rapidly or high dimensional integration is required. To overcome these issues in joint meta-analysis of studies comparing $T>2$ diagnostic tests in a multiple tests design with a gold standard, we propose a model that assumes the true positives and true negatives for each test are conditionally independent and binomially distributed given the $2T$-variate latent vector of sensitivities and specificities. For the random effects distribution, we employ an one-factor copula that provides tail dependence or tail asymmetry. Maximum likelihood estimation of the model is straightforward as the derivation of the likelihood requires bi-dimensional instead of $2T$-dimensional integration. Our methodology is demonstrated with an extensive simulation study and an application example that determines which is the best test for the diagnosis of rheumatoid arthritis.


Introduction
The identification of the most accurate diagnostic test for a particular disease contributes to the prevention of unnecessary risks to patients and healthcare costs. Diagnostic test accuracy studies aim to identify a new diagnostic test that is as accurate as the current perfect reference standard, also known as gold standard, yet less expensive or invasive.
Clinical and policy decisions are usually made on the basis of the results from many diagnostic test accuracy studies on the same research question. The considerably large number of diagnostic test accuracy studies has led to the use of meta-analysis. The purpose of a meta-analysis of diagnostic test accuracy studies is to combine information over different studies, and provide an integrated analysis that will have more statistical power to detect an accurate diagnostic test than an analysis based on a single study. As the accuracy of a diagnostic test is commonly measured by a pair of indices such as sensitivity and specificity, synthesis of diagnostic test accuracy studies is the most common medical application of multivariate meta-analysis (Jackson et al., 2011).
Most of the existing meta-analysis models and methods, when a perfect reference standard is available, have mainly focused on a single test (e.g., Rutter and Gatsonis 2001;Reitsma et al. 2005;Chu and Cole 2006).
However, as the understanding of a particular disease increases, along with technological advances, the comparative test accuracy of more than one diagnostic tests is apparent. As summarized by Takwoingi et al. (2013), diagnostic test accuracy studies can be comparative when they assess two or more tests or noncomparative when they assess one diagnostic test. Estimates of comparative test accuracy can be obtained from either category of studies, but the ones from the latter category are confounded by study setting. The robust comparative studies of diagnostic test accuracy use either a multiple test (also called paired or crossover) design, in which all patients undergo all tests together with the perfect reference standard, or more rarely, a randomised (also called parallel) design, in which all patients undergo the perfect reference standard test but are randomly allocated to have only one of the other tests. A multiple test design is statistically much more efficient, in that one needs much smaller sample sizes to detect a given difference in test accuracy, compared with a randomized design.
As the meta-analysis of more than one diagnostic tests can impact clinical decision making and patient health there is an increasing body of research that focus on the development of meta-analysis models and methods for the synthesis of studies comparing multiple diagnostic tests. Trikalinos et al. (2014) were the first who developed a model for the joint meta-analysis of studies comparing two diagnostic tests in a multiple tests design with a gold standard. They proposed a multinomial generalized linear mixed model (GLMM) which assumes independent multinomial distributions for the counts of each combination of test results in diseased patients, and, the counts of each combination of test results in non-diseased patients, conditional on the transformed latent true positive rate (TPR) and false positive rate (FPR) for each test, and latent joint TPR and FPR, which capture information on the agreement between the two tests in each study. Dimou et al. (2016) extended the bivariate model of Reitsma et al. (2005), which jointly meta-analyses the study-estimates of sensitivity and specificity for the case of a single test, to the case of two tests. They modelled the transformed study-estimates of TPR and FPR of the two tests using a quadrivariate normal distribution, with the information on the agreement between the two tests incorporated in the calculation of the within-study covariance matrix which is assumed fixed. Nikoloulopoulos (2020c) proposed a multinomial truncated D-vine copula mixed model for the joint meta-analysis of studies comparing two diagnostic tests, which assumes independent multinomial distributions for the counts of each combination of test results in diseased and non-diseased patients, conditional on the latent vector of probabilities of each combination of test results in diseased and non-diseased patients.
Their proposed model includes the multinomial GLMM (Trikalinos et al., 2014) as a special case, but can also operate on the original scale of the latent proportions.
As the information on the agreement between the two tests is usually not available from all the primary studies, Hoyer and Kuss (2018) proposed a model that is solely based on the information from the two (one per test) 2 × 2 tables with the number of true positives, true negatives, false negatives and false positives per study. They extended the bivariate generalized mixed model (GLMM) proposed by Chu and Cole (2006) to the quadrivariate case. The proposed quadrivariate GLMM assumes that the true positives and true negatives from the two tests are conditionally independent and binomially distributed given the bivariate latent pairs of transformed sensitivity and specificity, which are quadrivariate normally distributed. Nikoloulopoulos (2019b) generalised the quadrivariate GLMM by proposing a model that instead links the four random effects using a quadrivariate D-vine copula rather than the quadrivariate normal distribution.
However, for a particular disease there may be three (or more) diagnostic tests developed, where each of the tests is subject to several studies (e.g., Takwoingi et al. 2013). The extension of the aforementioned models (Trikalinos et al., 2014;Dimou et al., 2016;Hoyer and Kuss, 2018;Nikoloulopoulos, 2019bNikoloulopoulos, , 2020c to compare the accuracy of more than two tests suffers from the curse of multi-dimensionality, i.e., either the number of model parameters increase rapidly or high dimensional integration is required. In this paper to overcome the drawbacks in existing models for the joint meta-analysis of studies comparing T > 2 diagnostic tests in a multiple test design with a gold standard, we propose a model that assumes the true positives and true negatives for each test are conditionally independent and binomially distributed given the 2T -variate latent (random) vector of (transformed) sensitivities and specificities. For the random effects distribution, we employ an one-factor copula (Krupskii and Joe, 2013;Nikoloulopoulos and Joe, 2015;Kadhem and Nikoloulopoulos, 2021). The one-factor copula can provide, with appropriately chosen linking copulas, asymmetric dependence structure as well as tail dependence (dependence among extreme values) as it is an 1-truncated C-vine copula (Brechmann et al., 2012) rooted at the latent variable/factor. Joe et al. (2010) have shown that by choosing bivariate linking copulas appropriately, vine copulas can have a flexible range of lower/upper tail dependence, and different lower/upper tail dependence parameters for each bivariate margin.
With an one-factor copula, dimension reduction is achieved as the dependence among the latent sensitivities and specificities is explained by one other latent variable/factor. Hence, the proposed model has 2T dependence parameters instead of T (2T − 1), but more importantly its derivation requires bi-dimensional instead of 2T -dimensional integration.
The remainder of the paper proceeds as follows. Section 2 introduces the one-factor copula mixed model for the comparison of multiple diagnostic tests in a multiple tests design with a gold standard and Section 3 discusses its relationship with the 2T -variate GLMM. Section 4 deduces summary receiver operating characteristic (SROC) curves from the proposed model through quantile regression techniques. Section 5 provides a fast and efficient maximum likelihood (ML) estimation technique based on dependent Gauss-Legendre quadrature points that have an one-factor copula distribution and Section 6 contains small-sample efficiency calculations to investigate the effect of misspecifying the random effects distribution on parameter estimators and standard errors. Section 7 applies our methodology to data from a meta-analysis of diagnostic tests for rheumatoid arthritis. We conclude with some discussion in Section 8, followed by a brief section with software details.
2 The one-factor copula mixed model We first introduce the notation used in this paper. Let i be an index for the individual studies, j an index for the test outcome (0:negative; 1:positive), k an index for the disease outcome (0: non-diseased; 1: diseased) and t an index for the diagnostic test. The frequency data y ijkt , i = 1, ..., N, j = 0, 1, k = 0, 1, t = 1, . . . , T , corresponding to a combination of index test and disease outcomes in study i for test t, form a 2 × 2T table (Table 1), that is T "classic" 2 × 2 tables. We assume that the gold standard is the same for the T tests, i.e. y i+01 = · · · = y i+0T and y i+11 = · · · = y i+1T .
The one-factor copula can be explained as an 1-truncated C-vine rooted at the latent variable V (Krupskii and Joe, 2013;Nikoloulopoulos and Joe, 2015;Kadhem and Nikoloulopoulos, 2021). 2T -dimensional C-vine copulas can cover flexible dependence structures through the specification of 2T bivariate marginal copulas at level 1 and T (2T − 1) bivariate conditional copulas at higher levels (Nikoloulopoulos et al., 2012). For the 2Tdimensional one-factor copula, the pairs at level 1 are U, U kt , for k = 0, 1, t = 1, . . . , T , and for higher levels the (conditional) copula pairs are set to independence. That is the 1-factor copula has 2T bivariate copulas C kt,V (·; θ kt ) that link U kt , k = 0, 1, t = 1, . . . , T with V in the 1st level of the vine and independence copulas in all the remaining levels of the vine (truncated after the 1st level). Figure 1 depicts the graphical representation of the 1-factor copula model. Joe et al. (2010) have shown that in order for a vine copula to have (tail) dependence for all bivariate margins, it is only necessary for the bivariate copulas in level 1 to have (tail) dependence and it is not necessary for the conditional bivariate copulas in levels 2, . . . , 2T to have (tail) dependence.
The models in (1) and (4) together specify an one-factor copula mixed model with joint likelihood where g y; n, π = n y π y (1 − π) n−y , y = 0, 1, . . . , n, 0 < π < 1, is the binomial probability mass function (pmf). It is shown that the joint likelihood is represented as an one-dimensional integral of a function which in turn is a product of 2T one-dimensional integrals. As a result, 2T -dimensional numerical integration has been avoided.
Our general statistical model allows for selection of bivariate copulas and univariate margins independently, i.e., there are no constraints in the choices of parametric bivariate copulas and univariate margins. In line with our previous contributions in copula mixed models (Nikoloulopoulos, 2015(Nikoloulopoulos, , 2017(Nikoloulopoulos, , 2018a(Nikoloulopoulos, ,b, 2019b(Nikoloulopoulos, , 2020a we use • the choices of F ·; l(π), δ and l that are given in Table 2. With a beta distribution we work on the original scale of the latent sensitivities and specificities. The choices of the F ·; l(π), δ and l in the one-factor copula mixed model.
We show what happens when all the bivariate copulas C kt,V (; θ kt ) are BVN and the univariate distribution of the random effects is the N (µ, σ) distribution.
One can easily deduce that the within-study model in (1) is the same as in the 2T -variate GLMM. Furthermore, when C kt,V (; θ kt ) are all BVN copulas, then (2) becomes the copula of the multivariate normal distribution with an one-factor correlation structure. Let C kt,V (; θ kt ) be the BVN copula with correlation parameter θ kt . Let Φ and φ denote the standard normal cdf and density function, and let Φ 2 (·; ρ) be the BVN cdf , to get a 2T -variate distribution with N (0, 1) margins. Then This model is the same as the 2T -variate normal model with an one-factor correlation structure This occurs because the multivariate cdf in (6) comes from the representation where W, ǫ kt are i.i.d. N (0, 1) random variables (Krupskii and Joe, 2013;Nikoloulopoulos and Joe, 2015).
The resulting random effects distribution for (X 11 , . . . , X 1t , . . . , X 1T , X 01 , . . . , X 0t , . . . , X 0T ) is the 2Tvariate normal distribution with mean vector µ = l(π 1 ), l(π 0 ) and variance-covariance matrix Hence, the proposed model has as special case the 2T -variate GLMM with an one-factor correlation structure that has a latent additive structure as seen in (8). Nevertheless, if other bivariate copulas are used, then the one-factor copula mixed model has a latent structure that is non-additive.

Summary receiver operating characteristic curves
Though typically the focus of meta-analysis has been to derive the summary-effect estimates, there is increasing interest in alternative summary outputs, such as summary receiver operating characteristic (SROC) curves (e.g., Arends et al. 2008;Rücker and Schumacher 2009). In this section we derive the SROC curves from the onefactor copula mixed model.
For the one-factor copula mixed model, the model parameters (including dependence parameters), the choice of the copula, and the choice of the margin will affect the shape of the SROC curve. Let the joint cdf of (U 1t , U 0t ) be given by the copula C 1t,0t (·; θ 1t,0t ). The copula parameters θ 1t,0t , t = 1, . . . , T can be derived using the following steps: 1. Convert the copula parameters θ 1t and θ 0t of BVN, Frank or (rotated) Clayton copulas to Kendall's τ 1t and τ 0t via the relations in Hult and Lindskog (2002), Genest (1987), or Genest and MacKay (1986), respectively.
Then, the SROC curves for the latent pair (X 1t , X 0t ) can be deduced through the quantile regression techniques proposed by Nikoloulopoulos (2015): 1. Set C 1t|0t (u 1t |u 0t ; θ 1t,0t ) = q; 2. Solve for the quantile regression curve u 1t := u 1t (u 0t , q; θ 1t,0t ) = C −1 1t|0t (q|u 0t ; θ 1t,0t ); 3. Replace u kt by F x kt ; l(π kt ), δ kt ; As there is no priori reason to regress X 1t on X 0t instead of the other way around (Arends et al., 2008), quantile regression curves of X 0t on X 1t are also derived in a similar manner. We use the median regression curves (q = 0.5), along with the quantile regression curves with a focus on high (q = 0.99) and low quantiles (q = 0.01), which are strongly associated with the upper and lower tail dependence, respectively, imposed from each parametric family of bivariate copulas. These can be seen as confidence regions, as per the terminology in Rücker and Schumacher (2009), of the median regression curves. Finally, in order to reserve the nature of a bivariate response instead of a univariate response along with a covariate, we plot the corresponding contour graph of the bivariate copula density. The contour plot can be seen as the predictive region (analogously to Reitsma et al. 2005) of the estimated pair (π 1t , π 0t ) of the meta-analytic parameters of sensitivity and specificity at test t.

Maximum likelihood estimation and computational details
Estimation of the model parameters (π 1 , π 0 , δ 1 , δ 0 , θ) can be approached by the standard maximum likelihood (ML) method, by maximizing the logarithm of the joint likelihood in (5). The estimated parameters can be obtained by using a quasi-Newton (Nash, 1990) method applied to the logarithm of the joint likelihood. This numerical method requires only the objective function, i.e., the logarithm of the joint likelihood, while the gradients are computed numerically and the Hessian matrix of the second order derivatives is updated in each iteration. The standard errors (SE) of the ML estimates can be also obtained via the gradients and the Hessian computed numerically during the maximization process.
For one-factor copula mixed models of the form with joint likelihood as in (5), numerical evaluation of the joint pmf can be achieved with the following steps: 1. Calculate Gauss-Legendre (Stroud and Secrest, 1966) quadrature points {u q : q = 1, . . . , N q } and weights {w q : q = 1, . . . , N q } in terms of standard uniform.
With Gauss-Legendre quadrature, the same nodes and weights are used for different functions; this helps in yielding smooth numerical derivatives for numerical optimization via quasi-Newton (Nash, 1990). Our onefactor copula mixed model for meta-analysis of multiple diagnostic tests is straightforward computationally as it requires the calculation of a double summation over the quadrature points.
6 Small-sample efficiency -Misspecification of the random effects distribution In this section, we study the small-sample efficiency and robustness of the ML estimation of the one-factor copula mixed model. In Section 6.1, we gauge the small-sample efficiency of the ML method in Section 5 and investigate the misspecification of either the parametric margin or bivariate copula of the random effects distribution. In Section 6.2, we investigate the mixed model misspecification by using the D-vine copula mixed model proposed by Nikoloulopoulos (2019b) as the true model. That is we include a sensitivity analysis to the conditional independence assumption.

Misspecification of the parametric margin or bivariate pair-copulas
We randomly generated 10,000 samples of size n = 20, 50, 100 from an one-factor copula mixed model with both normal and beta margins that jointly meta-analyses T = {2, 3, 4} diagnostic tests.
4. Draw the number of diseased n 1 from a B(n, 0.4) distribution and set n 0 = n − n 1 .
Representative summaries of findings on the performance of the ML method in Section 5 are given in Tables   3 and 4 for 6-dimensional (T = 3) one-factor copula models with normal and beta margins, respectively. The true (simulated) bivariate copulas are the Clayton and the Clayton copula rotated by 270 • to handle the positive and negative dependencies, respectively. True sensitivity π 1 and specificity π 0 vectors are (0.8, 0.7, 0.8) and Conclusions from the values in the Tables 3 and 4 are the following: • ML with the true one-factor copula mixed model is highly efficient according to the simulated biases, SDs and RMSEs.
• The MLEs of π 1 , π 0 are not robust to margin misspecification, e.g., in Table 3 (Table 4) where the true univariate margins are normal (beta) the scaled biases for the MLEs of π 02 for the various one-factor copula mixed models with beta (normal) margins range from −4.16 (3.21) to −1.86 (4.70).  • The MLEs of π 1 , π 0 are rather robust to bivariate copula misspecification, but their biases increase when the assumed bivariate copulas have different tail dependence behaviour. For example, in Table 3 (Table   4) the scaled biases for the MLEs of π 11 for the various one-factor copula mixed models with normal (beta) margins increase to −1.38 (-0.92) and −3.71 (-2.27) when rotated Clayton copulas with opposite direction tail dependence and Frank copulas with tail independence, respectively, are called.
• The MLEs of δ 1 , δ 0 are rather robust to bivariate copula misspecification, but their biases increase when the assumed bivariate copula has tail dependence of opposite direction from the true bivariate copula. For example, in Table 3 (Table 4) the scaled biases for the MLEs of σ 02 (γ 02 ) for the various one-factor copula mixed models with normal (beta) margins range from −0.83 ( −0.25 ) to 1.75 (0.05), but the scaled bias increases to 6.85 (1.01) when rotated Clayton copulas with opposite direction tail dependence are called.
• The ML estimates of τ 's are robust to margin misspecification, as the copula remains invariant under any series of strictly increasing transformations of the components of the random vector, e.g., in Table 3 the scaled bias ofτ 13 is 1.81 for the true one-factor copula mixed model and 1.57 for an one-factor copula mixed model with the true bivariate copulas but beta margins.

Misspecification of the copula-mixed model -Sensitivity analysis to the conditional independence
We show a sensitivity analysis to the conditional independence assumption. We randomly generate 10,000 samples from the D-vine copula mixed model with both normal (Table 5) and beta (Table 6) margins using the algorithm in Nikoloulopoulos (2019b). We set the sample size N , the study size n, the true univariate and Kendall's τ parameters, and the disease prevalence to mimic the rheumatoid arthritis data in Nishimura et al. (2007). The D-vine copula mixed model assumes full dependence among the tests as the D-vine copula is not truncated, i.e., there are bivariate copulas not only at level 1 of the D-vine. Figure   U 01 U 02 |U 12 Level 3

Figure 2: Graphical representation of the 4-dimensional D-vine copula model with 3 levels.
From Table 5 (Table 6) it is seen that the one-factor copula mixed model with normal (beta) margins led to unbiased and efficient estimates when the bivariate copulas are a combination of Clayton and rotated Clayton by 270 • to model the positive and negative dependencies, respectively. These are the same with the true (simulated) copulas of the D-vine copula mixed model which imply that the sensitivity and specificity of each test have tail dependence. Hence, the tail dependence between the factor and each of the latent sensitivities/specificities is inherited to the tail dependence between the latent sensitivities and specificities, and thus, the conditional independence assumption has no impact on the estimation of the meta-analytic parameters of sensitivity and specificity of each test when this assumption is violated. This is due the fact that the one-factor copula can be explained as an 1-truncated C-vine rooted at the factor (Krupskii and Joe, 2013;Nikoloulopoulos and Joe, 2015;Kadhem and Nikoloulopoulos, 2021). Note also that in line with the results in the preceding subsection, the biases of the estimates increase when the assumed bivariate copulas have tail dependence of opposite direction from the true copulas or tail independence. When the BVN copulas with intermediate tail dependence are used to link the factor with the latent sensitivities/specificities, the estimates are robust to misspecification of the copula mixed model as long as the univariate margins are correctly specified.
Finally in order to study the relative performance of the one-factor copula mixed model over the quadrivariate vine copula mixed model as the number of quadrature points increase we randomly generated B = 20 samples of size N = 22 from the D-vine copula mixed model. The model parameters are set as before. The simulations were carried out on a Broadwell E5-2680 v4@2.40GHz. Table 7 summarizes the computing times (averaged over 20 replications) in seconds. Clearly the D-vine copula mixed approach requires a much higher computing time. Hence it is demonstrated that even for the case of T = 2 tests, the computational improvement of the one-factor copula mixed model is substantial, as one has to calculate numerically bivariate integrals instead of much more difficult quadrivariate integrals.  margin copula π 11 π 12 π 01 π 02 γ 11 γ 12 γ 01 γ 02  7 Application Nishimura et al. (2007) contacted a systematic review and summarized data of rheumatoid factor (RF) and anticyclic citrullinated peptide (anti-CCP) antibodies for diagnosing rheumatoid arthritis. They included N = 22 studies that assessed both RF and anti-CCP2 antibody for diagnosing rheumatoid arthritis and used the 1987 revised American College of Rheumatology (ACR) criteria as the perfect reference standard of rheumatoid arthritis (Arnett et al., 1988). These data have been frequently used as an example for methodological papers on joint meta-analysis of diagnostic accuracy studies in a multiple tests design with a gold standard (e.g., Dimou et al. 2016;Nikoloulopoulos 2019b). Liu et al. (2015) in one of their examples deal with the same data, but as they propose models for the meta-analysis of the accuracy of a diagnostic test under evaluation and an imperfect reference test, they use only the the RF test as the index test for detection of rheumatoid arthritis and assume that the ACR 1987 revised criteria are an imperfect reference test for classification. Their analysis confirmed that the ACR 1987 revised criteria are a prefect reference test as the estimates of sensitivity and specificity of the ACR 1987 criteria (reference test) were 1, suggesting that such reference test is in fact a gold standard.
We use the one-factor copula mixed model in order to determine whether anti-CCP antibody identifies more accurately patients with rheumatoid arthritis than RF does. We fit the one-factor copula mixed model for all choices of parametric families of copulas and margins. To make it easier to compare strengths of dependence, we convert from the BVN, Frank and (rotated) Claytonθ's toτ 's via the relations in (9), (10), and (11).
Because the number of parameters is the same between the models, we use the log-likelihood at the maximum likelihood estimates as a rough diagnostic measure for model selection between the models. For vine copulas (one-factor copula is an 1-truncated C-vine copula), Dissmann et al. (2013) found that pair-copula selection based on likelihood seems to be better than even using bivariate goodness-of-fit tests. The goodness-of-fit procedures involve a global distance measure between the model-based and empirical distribution, hence they might not be sensitive to tail behaviours and are not diagnostic in the sense of suggesting improved parametric models in the case of small p-values (Joe, 2014, page 254). A larger likelihood value indicates a model that better approximates both the dependence structure of the data and the strength of dependence in the tails.
The log-likelihoods showed that an one-factor copula mixed model with Clayton and Clayton rotated by 270 • degrees copulas with normal margins to join the factor with each of the sensitivities/specificities provides the best fit (Table 8). For this particular example it is revealed that an one-factor copula mixed model with the sensitivities and specificities on the transformed scale provides better fit than an one-factor copula mixed model with beta margins, which models the sensitivity and specificity on the original scale.
The resultant sensitivities and specificities show that the anti-CCP2 antibody is better compared with RF.
Both tests have fairly similar sensitivity but the anti-CCP2 is much more specific. On the one hand, the estimated univariate parameters and standard errors are in line with the ones in Nikoloulopoulos (2019b), but the implementation of the proposed model is much faster, since a numerically time-consuming four-dimensional integral calculation is replaced with a numerically fast two-dimensional integral calculation on the other.

Rheumatoid Factor
Anti-CCP2 antibody  Figure 3: Contour plots (predictive region) and quantile regression curves from the best fitted one-factor copula mixed model for the rheumatoid arthritis data. Red and green lines represent the quantile regression curves x1t := x1t(x0t, q) and x0t := x0t(x1t, q), respectively; for q = 0.5 solid lines and for q ∈ {0.01, 0.99} dotted lines (confidence region). The axes are in logit scale since we also plot the estimated contour plot of the random effects distribution as predictive region; this has been estimated for the logit pair of (Sensitivity, Specificity) for each test.
x 1t := x 1t (x 0t , q = 0.5) x 0t := x 0t (x 1t , q = 0.5) From the Kendall's tau estimates and standard errors there is strong evidence of dependence between the two diagnostic tests. The fact that the best-fitting bivariate copulas are Clayton and Clayton rotated by 270 • reveals that there is tail dependence among the latent sensitivities and specificities. This can be further seen trough the predictive region of the SROC curves. Figure 3 depicts the SROC curves and summary operating points (a pair of average sensitivity and specificity) with a confidence region and a predictive region for each test from the best fitted one-factor copula mixed model. Sharper corners in the predictive region indicate tail dependence. Figure 4 provides a direct and visual comparison between the two competing diagnostic tests and reveals that the anti-CCP2 antibody is better compared with RF.

Discussion
We have proposed an one-factor copula mixed model for joint meta-analysis and comparison of multiple diagnostic tests in a multiple tests design with a gold standard. This is a parsimonious meta-analytic model that (a) has the 2T -variate GLMM with an additive latent structure as a special case when the BVN copulas are used, (b) can have a latent structure that is not additive if other than BVN copulas are called, (c) can model the latent sensitivities and specificities on the original scale rather than a transformed scale as in the 2T -variate GLMM (d) enables the meta-analytic parameters of interest to be separated from the copula (dependence) parameters which are interpretable as dependence of the latent sensitivity/specificity with another latent variable, (e) avoids the curse of multi-dimensionality and (f) models adequately the dependence among the latent sensitivities and specificities as it can be explained as an 1-truncated C-vine copula.
Our model can provide an improvement over the 2T -variate GLMM with an additive latent structure as the random effects distribution is expressed via an one-factor copula that provides a wide range of dependence with 2T dependence parameters and allow for different types of tail behaviour, different from assuming simple linear correlation structures, normality and tail independence. This strength of multivariate meta-analysis approaches that use copulas has been pointed out by Jackson and White (2018) and Jackson et al. (2020) and it has also been exploited in network meta-analysis (Phillippo et al., 2020).
The 2T -variate D-vine copula mixed model, which it has as special case the 2T -variate GLMM with an unstructured correlation structure, provides full dependence, but it is intractable as the number of competing tests increases. The 2T -variate one-factor copula mixed model solves this problem since the joint likelihood reduces to an one-dimensional integral of a function which in turn is a product of 2T one-dimensional integrals, hence the method avoids 2T -dimensional integration which is time consuming even for T = 2 tests. Its parsimony is not a distributional concern about the dependence between the tests due to the main result in Joe et al. (2010): all the bivariate margins of the vine copula have (tail) dependence if the bivariate copulas at level 1 have (tail) dependence. This is satisfied by the one-factor copula as it is an 1-truncated C-vine. Hence, the proposed model can form the vehicle for conducting meta-analysis of comparative accuracy studies with three or more tests.
When the focus is on estimates of the meta-analytic univariate parameters of interest, the outgrowth of joint analysis is modest, in that the differences in the summary estimates and standard errors from separate meta-analyses for each test are not that distinct. The most striking differences between separate and joint meta-analyses arise when one deduces comparative diagnostic accuracy, i.e., an SROC curve. An SROC curve makes much more sense and will help decision makers to assess the actual diagnostic accuracy of the competing diagnostic tests. In an era of evidence-based medicine, decision makers need high-quality procedures such as the SROC curves to support decisions about whether or not to use a diagnostic test in a specific clinical situation and, if so, which test. We have deduced SROC curves from the one-factor copula mixed model. The model parameters (including dependence parameters), the choice of the copula, and the choice of the margin affect the shape of the SROC curves. A series of independence models cannot be used to produce the SROC curves, since the dependence parameters affect the shape of the SROC curve and these are set to independence.
Comparative accuracy studies with paired designs where each test is applied to the same patients should report the data as separate 2 × 2 tables. Authors of primary studies of diagnostic accuracy that assess three or more tests in the same patients should be encouraged to report sufficient data to extract separate 2 × 2 tables of test results as in Table 1. Comparative accuracy studies should rightly use a multiple tests designs so that patients receive each test in order to reduce biases and ensure the clinical relevance of the resulting inferences (Trikalinos et al., 2014).
Nevertheless, in practice there exist comparative studies in a randomized design or even non-comparative studies (Takwoingi et al., 2013) and for some of them the reference test might be imperfect. Future research will focus on extending the one-factor copula mixed model to incorporate randomised designs and non-comparative studies with or without a gold standard. Ma et al. (2018) and Lian et al. (2019) proposed methods for comparing multiple diagnostic tests that can incorporate studies with different designs and studies with our without gold standard. As their methods assume that the between-studies model is the multivariate normal distribution that suffers for the curse of multidimensionality when the numbers of tests increases, we will exploit the use of the one-factor copula distribution. The one-factor copula distribution will provide computational and distributional improvements when adopted to the setting of Ma et al. (2018) and Lian et al. (2019).
Software R functions to implement the one-factor copula mixed model for meta-analysis of multiple diagnostic tests will be part of the next major release of the R package CopulaREMADA (Nikoloulopoulos, 2019a).