---
abstract: |
    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 diagnostic studies which compare the accuracy of more than one 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, 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 flexible reflection asymmetric tail and non-linear dependence. 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.

    **Key Words:** Diagnostic tests; factor copulas; multivariate meta-analysis; mixed models; sensitivity; specificity.
author:
- 'Aristidis K. Nikoloulopoulos [^1]'
title: |
    An one-factor copula mixed model for joint meta-analysis\
    of multiple diagnostic tests
---

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 or more accurate than the current 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 considerable 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 (e.g., @noma-etal-2020-biometrics). Most of the existing literature has mainly focused on a single test as acknowledged in a recent review by [@Ma-etal-2013].

However, as the understanding of a particular disease increases, along with technological advances, the development of multiple tests is apparent. Therefore, many studies that evaluate the accuracy of several tests instead of the accuracy of a single test have been conducted (e.g., @Nishimura-etal2007 [@Bennett-etal2007; @Kodama-etal2013]). 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 diagnostic studies which compared the accuracy of two tests .

[@trikalinos-etal-2014-rsm] were the first who developed a model for the joint meta-analysis of studies comparing two diagnostic tests. They proposed a multinomial generalized linear mixed model 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-etal2016] extended the bivariate model of [@Reitsma-etal-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.

As the information on the agreement between the two tests is usually not available from all the primary studies, proposed a model that is solely based on the information from the two (one per test) $2\times2$ 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 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-2018-smmr] 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-etal-2013). The extension of the aforementioned models 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, 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. We assume that the same individuals receive all tests, hence the test results are correlated. As a valid statistical model should account for this dependence, for the random effects distribution, we employ an one-factor copula . The one-factor copula can provide flexible reflection asymmetric tail and non-linear dependence as it is a truncated canonical vine copula rooted at the latent variable/factor. [@joeetal10] have shown that by choosing bivariate linking copulas appropriately, vine copulas can have a flexible range of lower/upper tail dependence (dependence among extreme values), 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 $\binom{T}{2}$, but more importantly its derivation requires bi-dimensional instead of $2T$-dimensional integration.

The remainder of the paper proceeds as follows. Section \[copula-mixed-model-sec\] introduces the one-factor copula mixed model for the comparison of multiple diagnostic tests and discusses its relationship with existing models. Section \[miss-section\] contains small-sample efficiency calculations to investigate the effect of misspecifying the random effects distribution on parameter estimators and standard errors. Section \[app-sec\] illustrates our methodology. We conclude with some discussion in Section \[discussion\].

\[copula-mixed-model-sec\]The one-factor copula mixed model 
============================================================

We first introduce the notation used in this paper. The data are $y_{ijkt},\, i = 1, . . . ,N,\, j=0,1,\,k=0,1,\,t=1,\ldots,T$, where $i$ is an index for the individual studies, $j$ is an index for the test outcome (0:negative; 1:positive), $k$ is an index for the disease outcome (0: non-diseased; 1: diseased) and $t$ is an index for the diagnostic test. The “classic" $2\times 2$ table is extended to a $2\times 2T$ table (Table \[2times2T\]), that is $T$ $2\times 2$ tables. Each cell in Table \[2times2T\] provides the cell frequency corresponding to a combination of index test and disease outcomes in study $i$ for test $t$. We assume that the gold standard is the same for the $T$ tests, i.e. $y_{i+01}=\cdots=y_{i+0T}$ and $y_{i+11}=\cdots=y_{i+1T}$.

  ------------ -- ------------ ------------ ---------- -------------- -- ------------ ------------ ---------- -------------- -- ------------ ------------ -- -- -- --
                                                                                                                                                                   
  [Test 1]{}          $-$          $+$       $\ldots$   [Test $t$]{}     $-$              $+$       $\ldots$   [Test $T$]{}     $-$          $+$                   
  $-$              $y_{i001}$   $y_{i011}$   $\ldots$       $-$          $y_{i00t}$    $y_{i01t}$   $\ldots$       $-$          $y_{i00T}$   $y_{i01T}$            
  $+$              $y_{i101}$   $y_{i111}$   $\ldots$       $+$          $y_{i10t}$    $y_{i11t}$   $\ldots$       $+$          $y_{i10T}$   $y_{i11T}$            
  Total            $y_{i+01}$   $y_{i+11}$   $\ldots$      Total         $y_{i+0t}$    $y_{i+1t}$   $\ldots$      Total         $y_{i+0T}$   $y_{i+1T}$            
  ------------ -- ------------ ------------ ---------- -------------- -- ------------ ------------ ---------- -------------- -- ------------ ------------ -- -- -- --

  : \[2times2T\] Data from an individual study in a $2\times 2T$ table.

The within-study model assumes that the number of true positives $Y_{i11t}$ and true negatives $Y_{i00t}$ for $t=1,\ldots,T$ are conditionally independent and binomially distributed given $(\X_{1t},\X_{0t})=(\x_{1t},\x_{0t})$, where $\X_t=(X_{1t},X_{0t})$ denotes the bivariate latent pair of (transformed) sensitivity and specificity for the test $t$. That is $$\begin{aligned}
\label{withinBinom}
Y_{i11t}|X_{1t}=x_{1t}&\sim& \mbox{Binomial}\Bigl(y_{i+1t},l^{-1}(x_{1t})\Bigr);\nonumber\\
Y_{i00t}|X_{0t}=x_{0t}&\sim& \mbox{Binomial}\Bigl(y_{i+0t},l^{-1}(x_{0t})\Bigr),\end{aligned}$$ for $t=1,\ldots, T$, where $l(\cdot)$ is a link function.

The stochastic representation of the between studies model takes the form

$$\begin{aligned}
\label{copula-between-norm}
&&\Bigl(F\bigl(X_{11};l(\pi_{11}),\de_{11}\bigr),\ldots,
F\bigl(X_{1t};l(\pi_{1t}),\de_{1t}\bigr),\ldots,F\bigl(X_{1T};l(\pi_{1T}),\de_{1T}\bigr)\nonumber\\
&&F\bigl(X_{01};l(\pi_{01}),\de_{01}\bigr),\ldots,
F\bigl(X_{0t};l(\pi_{0t}),\de_{0t}\bigr),\ldots,F\bigl(X_{0T};l(\pi_{0T}),\de_{0T}\bigr)\Bigr)\sim C(\cdot;\thbf),\end{aligned}$$

where $F\bigl(\cdot;l(\pi),\de\bigr)$ is the cdf of the univariate distribution of the random effect and $C(\cdot;\thbf)$ is an one-factor copula with dependence parameter vector $\thbf=(\th_{11},\ldots,\th_{1t},\ldots,\th_{1T},\th_{11},\ldots,\th_{1t},\ldots,\th_{1T})$. That is, we assume that $U_{11},\ldots,U_{1T},U_{01},\ldots,U_{0T}$, where $U_{kt}=F\bigl(X_{kt};l(\pi_{kt}),\de_{kt}\bigr),\,k=0,1,\,t=1,\ldots,T$, are conditionally independent given another standard uniform latent variable/factor $U$. Then we can easily get the joint density of the one-factor copula using 1-dimensional numerical integration, viz., $$c(u_{11},\ldots,u_{1t},\ldots,u_{1T},u_{01},\ldots,u_{0y},\ldots,u_{0T};\thbf)=\int_0^1\biggr\{\prod_{t=1}^T\biggr[c(u_{1t},u;\th_{1t})c(u_{0t},u;\th_{0t})\biggl]\biggl\}du.$$ It is seen that the $2T$-variate density decomposes in a product of $2T$ bivariate copula densities that link the random effects $U_{kt},\,k=0,1,\,t=1,\ldots,T$ with $U$. The parameters $(\pi_{11},\ldots,\pi_{1t},\ldots,\pi_{1T})$ $:=\bfpi_1$ and $(\pi_{01},\ldots,\pi_{0t},\ldots,\pi_{0T}):=\bfpi_0$ are those of main interest denoting the meta-analytic parameters for the sensitivities and specificities, respectively, while the univariate parameters $(\de_{11},$ $\ldots,\de_{1t},\ldots,\de_{1T}):=\debf_1$ and $(\de_{01},\ldots,\de_{0t},$ $\ldots,\de_{0T}):=\debf_0$ are of secondary interest denoting the between-study variabilities for the sensitivities and specificities, respectively. The copula parameter vector $\thbf$ has parameters of the random effects model and they are separated from the univariate parameters $\bfpi_1,\bfpi_0,\debf_1,\debf_0$.

The models in (\[withinBinom\]) and (\[copula-between-norm\]) together specify an one-factor copula mixed model with joint likelihood

$$\begin{aligned}
\label{mixed-cop-likelihood}
&L(\bfpi_1,\bfpi_0,\debf_1,\debf_0,\thbf)\nonumber\\
&
=\prod_{i=1}^N\int_{[0,1]^{2T}}\biggl\{\prod_{t=1}^T\biggl[ g\biggl(y_{i11t};y_{i+1t},l^{-1}\Bigl(F^{-1}\bigl(u_{1t};l(\pi_{1t}),\de_{1t}\bigr)\Bigr)\biggr) g\biggl(y_{i00t};y_{i+0t},
l^{-1}\Bigl(F^{-1}\bigl(u_{0t};\nonumber\\&l(\pi_{0t}),\de_{0t}\bigr)\Bigr)\biggr)\biggr]
\int_0^1\biggl\{\prod_{t=1}^T\biggl[c(u_{1t},u;\th_{1t})c(u_{0t},u;\th_{0t})\biggr]\biggr\}du\biggr\}du_{11},\ldots,du_{1T}\,du_{01},\ldots,du_{0T}\nonumber\\
&=
\prod_{i=1}^N\int_0^1\biggl\{\prod_{t=1}^T\biggl[ \int_0^1 \biggl\{g\biggl(y_{i11t};y_{i+1t},l^{-1}\Bigl(F^{-1}\bigl(u_{1t};l(\pi_{1t}),\de_{1t}\bigl)\Bigr)\biggr)c(u_{1t},u;\th_{1t})\biggr\}du_{1t}\nonumber\\
&\int_0^1\biggl\{g\biggl(y_{i00t};y_{i+0t},l^{-1}\Bigl(F^{-1}\bigl(u_{0t};l(\pi_{0t}),\de_{0t}\bigl)\Bigr)\biggr)c(u_{0t},u;\th_{0t})\biggr\}du_{0t}\biggr]\biggr\}du\end{aligned}$$

where $g\bigl(y;n,\pi\bigr)=\binom{n}{y}\pi^y(1-\pi)^{n-y},\quad y=0,1,\ldots,n,\quad 0<\pi<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 [@Nikoloulopoulos2015b; @Nikoloulopoulos2015c; @Nikoloulopoulos-2016-SMMR; @Nikoloulopoulos2018-AStA; @nikoloulopoulos-2018-smmr; @Nikoloulopoulos-2018-3dmeta-NE; @Nikoloulopoulos-2018-4dmeta-NE] we use

-   bivariate parametric copulas with different tail dependence behaviour, namely the bivariate normal (BVN) with intermediate tail dependence, Frank with tail independence, and Clayton with positive lower tail dependence. For the latter we also use its rotated versions to provide negative upper-lower tail dependence (Clayton rotated by 90$^\circ$), positive upper tail dependence (Clayton rotated by 180$^\circ$) and negative lower-upper tail dependence (Clayton rotated by 270$^\circ$);

-   the choices of $F\bigl(\cdot;l(\pi),\de\bigr)$ and $l$ that are given in Table \[choices\]. With a beta distribution we work on the original scale of the latent sensitivities and specificities.

   $F\bigl(\cdot;l(\pi),\de\bigr)$     $l$          $\pi$        $\de$
  --------------------------------- ---------- --------------- ----------
             $N(\mu,\s)$              logit     $l^{-1}(\mu)$     $\s$
         Beta$(\pi,\gamma)$          identity       $\pi$       $\gamma$

  : \[choices\]The choices of the $F\bigl(\cdot;l(\pi),\de\bigr)$ and $l$ in the one-factor copula mixed model.

Relationship with the $2T$-variate GLMM
---------------------------------------

We show what happens when all the bivariate copulas are BVN and the univariate distribution of the random effects is the $N(\mu,\s)$ distribution. One can easily deduce that the within-study model in (\[withinBinom\]) is the same as in the $2T$-variate GLMM.

Furthermore, when all the bivariate copulas are BVN copulas with correlation parameters $\rhobf=(\rho_{11},\ldots,$ $\rho_{1t},\ldots,\rho_{1T},\rho_{01},\ldots,\rho_{0t},\ldots,\rho_{0T})$, the resulting random effects distribution is the $2T$-variate normal distribution with mean vector $\bmu=\bigl(l(\bfpi_1),l(\bfpi_0)\bigr)$ and variance-covariance matrix $$\Sbf=\begin{pmatrix}
\s_{11}^2&\cdots&\rho_{11}\rho_{1T}\s_{11}\s_{1T}&\rho_{11}\rho_{01}\s_{11}\s_{01}&\cdots&\rho_{11}\rho_{0T}\s_{11}\s_{0T}\\
\vdots &\ddots&\vdots& \vdots&\vdots&\vdots\\
\rho_{1T}\rho_{11}\s_{1T}\s_{11}&\cdots&\s_{1T}^2&\rho_{1T}\rho_{01}\s_{1T}\s_{01}&\cdots&\rho_{1T}\rho_{0T}\s_{1T}\s_{0T}\\
\rho_{01}\rho_{11}\s_{01}\s_{11}&\cdots&\rho_{01}\rho_{1T}\s_{01}\s_{1T}&\s_{01}^2&\cdots&\rho_{1T}\rho_{0T}\s_{1T}\s_{0T}\\
\vdots &\vdots&\vdots& \vdots&\ddots&\vdots\\
\rho_{0T}\rho_{11}\s_{0T}\s_{11}&\cdots&\rho_{0T}\rho_{1T}\s_{0T}\s_{1T}&\rho_{0T}\rho_{01}\s_{0T}\s_{01}&\cdots&\s_{0T}^2\\
\end{pmatrix}$$ Hence, the proposed model has as special case the $2T$-variate GLMM.

\[computation\]Maximum likelihood estimation and computational details
----------------------------------------------------------------------

Estimation of the model parameters $(\bfpi_1,\bfpi_0,\debf_1,\debf_0,\thbf_1,\thbf_0)$ can be approached by the standard maximum likelihood (ML) method, by maximizing the logarithm of the joint likelihood in (\[mixed-cop-likelihood\]). The estimated parameters can be obtained by using a quasi-Newton [@nash90] 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 (\[mixed-cop-likelihood\]), numerical evaluation of the joint pmf can be achieved with the following steps:

1.  Calculate Gauss-Legendre quadrature points $\{u_q: q=1,\ldots,n_q\}$ and weights $\{w_q: q=1,\ldots,n_q\}$ in terms of standard uniform.

2.  Numerically evaluate the joint pmf $$\begin{aligned}
    &\int_0^1\biggl\{\prod_{t=1}^T\biggl[ \int_0^1 \biggl\{g\biggl(y_{i11t};y_{i+1t},l^{-1}\Bigl(F^{-1}\bigl(u_{1t};l(\pi_{1t}),\de_{1t}\bigl)\Bigr)\biggr)c(u_{1t},v;\th_{1t})\biggr\}du_{1t}\\
    &\int_0^1\biggl\{g\biggl(y_{i00t};y_{i+0t},l^{-1}\Bigl(F^{-1}\bigl(u_{0t};l(\pi_{0t}),\de_{0t}\bigl)\Bigr)\biggr)c(u_{0t},v;\th_{0t})\biggr\}du_{0t}\biggr]\biggr\}dv \end{aligned}$$ in a double sum:

    $$\begin{aligned}
    &\sum_{q_1=1}^{n_q}\biggl\{w_{q_1}\prod_{t=1}^T\biggl[\sum_{q_2=1}^{n_q} \biggl\{w_{q_2}g\biggl(y_{i11t};y_{i+1t},l^{-1}\Bigl(F^{-1}\bigl(C^{-1}(u_{q_2}|u_{q_1};\th_{1t});l(\pi_{1t}),\de_{1t}\bigl)\Bigr)\biggr)\biggr\}\\
    &\sum_{q_2=1}^{n_q}\biggl\{w_{q_2}g\biggl(y_{i00t};y_{i+0t},l^{-1}\Bigl(F^{-1}\bigl(C^{-1}(u_{q_2}|u_{q_1};\th_{0t});l(\pi_{0t}),\de_{0t}\bigl)\Bigr)\biggr)\biggr\}\biggr]\biggr\},\end{aligned}$$

    where $C^{-1}(v|u;\th)$ is the inverse conditional bivariate copula cdf. Note that the independent quadrature points $\{u_{q_1}: q_1=1,\ldots,n_q\}$ and $\{u_{q_2}: q_2=1,\ldots,n_q\}$ have converted to dependent quadrature points that have an one-factor copula distribution $C(\cdot;\thbf)$.

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 [@nash90]. Our one-factor 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.

\[miss-section\]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 \[true\], we gauge the small-sample efficiency of the ML method in Section \[computation\] and investigate the misspecification of either the parametric margin or bivariate copula of the random effects distribution. In Section \[miss\], we investigate the mixed model misspecification by using the D-vine copula mixed model proposed by [@nikoloulopoulos-2018-smmr] as the true model.

\[true\]Misspecification of the parametric margin or bivariate pair-copulas
---------------------------------------------------------------------------

We randomly generated $B = 10^4$ 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.

The simulation process is as below:

1.  Simulate $(u_{11},\ldots,u_{1T},u_{01},\ldots,u_{0T})$ from an factor copula $C(;\taubf)$; $\taubf$ is converted to the copula parameter vector $\thbf$ of BVN, Frank or (rotated) Clayton copulas via the relations $$\label{tauBVN}
    \tau=\frac{2}{\pi}\arcsin(\th),$$ $$\label{tauFrank}
    \tau=\left\{\begin{array}{ccc}
    1-4\theta^{-1}-4\theta^{-2}\int_\theta^0\frac{t}{e^t-1}dt &,& \th<0\\
    1-4\theta^{-1}+4\theta^{-2}\int^\theta_0\frac{t}{e^t-1}dt &,& \th>0\\
    \end{array}\right.,$$ and $$\label{tauCln}
    \tau=\left\{\begin{array}{rcl}
    \th/(\th+2)&,& \mbox{by 0$^\circ$ or 180$^\circ$ }\\
    -\th/(\th+2)&,& \mbox{by 90$^\circ$ or 270$^\circ$}\\
    \end{array}\right.,$$ in [@HultLindskog02], [@genest87], and , respectively.

2.  For each test $t$ in $1,\ldots,T$ convert to proportions via $$\begin{aligned}
    x_{1t}&=&l^{-1}\biggl(F^{-1}\Bigl(u_{1t};l(\pi_{1t}),\de_{1t}\Bigr)\biggr) \\
    x_{0t}&=&l^{-1}\biggl(F^{-1}\Bigl(u_{0t};l(\pi_{0t}),\de_{0t}\Bigr)\biggr)\end{aligned}$$

3.  Simulate the study size $n$ from a shifted gamma distribution, i.e., $n\sim \mbox{sGamma}(\a=1.2,\b=0.01,\mbox{lag}=30)$ and round off to the nearest integer.

4.  Draw the number of diseased $n_{1}$ from a $B(n,0.4)$ distribution and set $n_0=n-n_1$ .

5.  For each test $t$ in $1,\ldots,T$ generate $y_{11t}$ and $y_{00t}$ from a $B(n_1,x_{1t})$ and $B(n_0,x_{0t})$ distribution, respectively, and set $y_{01t}=n_1-y_{11t},y_{10t}=n_0-y_{00t}.$

Representative summaries of findings on the performance of the ML method in Section \[computation\] are given in Table \[sim1\] for 6-dimensional ($T = 3$) one-factor copula models with normal and beta margins. The true (simulated) copula distributions are the Clayton and the Clayton copula rotated by $270^\circ$ to handle the positive and negative dependencies, respectively. True sensitivity $\bfpi_1$ and specificity $\bfpi_2$ vectors are $(0.8,0.7,0.8)$ and $(0.7,0.8,0.7)$, the variability parameter vectors are $\debf_1=\debf_2=(1,1,1)$ or $\debf_1=\debf_2=(0.1,0.1,0.1)$ for normal or beta margin, respectively, and the Kendall’s $\taubf=(0.6,0.7,0.5,-0.3,-0.4,-0.2)$. Under each margin, 10,000 meta-analysis data sets are simulated with $N=50$ studies in each data set.

The table contains the resultant biases, root mean square errors (RMSEs) and standard deviations (SDs), along with average theoretical variances, scaled by 100, for the MLEs of the main parameters of interest under different copula choices and margins. The theoretical variances of the MLEs are obtained via the gradients and the Hessian that were computed numerically during the maximization process. The results for the variability and dependence parameters are presented in the Supplementary material.

[llllccccccccccccc]{} & & &\

& & & & & &\

& margin & copula & & $\pi_{11}$ & $\pi_{12}$ & $\pi_{13}$ & $\pi_{01}$ & $\pi_{02}$ & $\pi_{03}$ & & $\pi_{11}$ & $\pi_{12}$ & $\pi_{13}$ & $\pi_{01}$ & $\pi_{02}$ & $\pi_{03}$\
Bias & normal & BVN & & -0.35 & -0.43 & -0.32 & -0.03 & 0.04 & -0.04 & & 2.93 & 1.97 & 2.95 & 2.24 & 3.25 & 2.24\
& & Cln{$0^\circ,270^\circ$} & & -0.31 & -0.40 & -0.31 & 0.05 & 0.08 & 0.02 & & 3.06 & 2.14 & 3.08 & 2.26 & 3.21 & 2.26\
& & Cln{$180^\circ,90^\circ$} & & -1.38 & -1.98 & -1.12 & 0.30 & 0.50 & 0.08 & & 2.35 & 1.15 & 2.53 & 2.44 & 3.57 & 2.31\
& & Frank & & -3.71 & -4.92 & -3.42 & 2.40 & 2.28 & 1.64 & & 0.59 & -0.94 & 0.80 & 3.80 & 4.70 & 3.33\
& beta & BVN & & -3.99 & -3.25 & -4.01 & -3.38 & -4.16 & -3.33 & & 0.02 & -0.04 & 0.01 & -0.05 & -0.05 & -0.02\
& & Cln{$0^\circ,270^\circ$} & & -4.28 & -3.57 & -4.31 & -3.21 & -3.92 & -3.24 & & 0.02 & -0.04 & 0.00 & -0.02 & 0.01 & 0.00\
& & Cln{$180^\circ,90^\circ$} & & -5.42 & -5.02 & -5.21 & -3.28 & -4.14 & -3.30 & & -0.92 & -1.13 & -0.76 & 0.06 & 0.02 & 0.02\
& & Frank & & -7.03 & -7.04 & -6.82 & -1.21 & -1.86 & -1.88 & & -2.27 & -2.72 & -2.10 & 1.45 & 1.52 & 1.00\
SD & normal & BVN & & 2.51 & 3.25 & 2.50 & 3.13 & 2.45 & 3.14 & & 2.06 & 2.47 & 2.05 & 2.40 & 2.00 & 2.40\
& & Cln{$0^\circ,270^\circ$} & & 2.55 & 3.31 & 2.55 & 3.15 & 2.47 & 3.16 & & 2.07 & 2.48 & 2.07 & 2.41 & 2.00 & 2.42\
& & Cln{$180^\circ,90^\circ$} & & 2.90 & 3.77 & 2.84 & 3.29 & 2.59 & 3.23 & & 2.26 & 2.71 & 2.21 & 2.49 & 2.11 & 2.46\
& & Frank & & 3.58 & 4.49 & 3.52 & 3.65 & 2.81 & 3.53 & & 2.75 & 3.23 & 2.72 & 2.81 & 2.32 & 2.72\
& beta & BVN & & 2.46 & 2.86 & 2.46 & 2.77 & 2.44 & 2.79 & & 2.00 & 2.27 & 1.98 & 2.20 & 1.92 & 2.21\
& & Cln{$0^\circ,270^\circ$} & & 2.54 & 2.95 & 2.56 & 2.81 & 2.48 & 2.82 & & 2.00 & 2.28 & 1.99 & 2.20 & 1.94 & 2.22\
& & Cln{$180^\circ,90^\circ$} & & 2.92 & 3.31 & 2.87 & 2.81 & 2.47 & 2.82 & & 2.34 & 2.59 & 2.28 & 2.23 & 1.96 & 2.24\
& & Frank & & 3.47 & 3.88 & 3.43 & 3.10 & 2.63 & 3.06 & & 2.79 & 3.03 & 2.76 & 2.49 & 2.13 & 2.44\
$\sqrt{\bar V}$ & normal & BVN & & 2.40 & 3.09 & 2.41 & 3.08 & 2.40 & 3.07 & & 1.90 & 2.35 & 1.91 & 2.33 & 1.91 & 2.32\
& & Cln{$0^\circ,270^\circ$} & & 2.32 & 2.95 & 2.34 & 2.95 & 2.25 & 2.99 & & 1.86 & 2.29 & 1.87 & 2.29 & 1.86 & 2.30\
& & Cln{$180^\circ,90^\circ$} & & 2.66 & 3.45 & 2.64 & 3.18 & 2.48 & 3.14 & & 2.08 & 2.58 & 2.06 & 2.40 & 1.97 & 2.38\
& & Frank & & 2.66 & 3.27 & 2.68 & 3.00 & 2.22 & 3.06 & & 2.10 & 2.46 & 2.11 & 2.29 & 1.79 & 2.32\
& beta & BVN & & 2.29 & 2.70 & 2.30 & 2.69 & 2.29 & 2.67 & & 1.89 & 2.19 & 1.90 & 2.16 & 1.90 & 2.15\
& & Cln{$0^\circ,270^\circ$} & & 2.17 & 2.54 & 2.19 & 2.55 & 2.13 & 2.59 & & 1.82 & 2.10 & 1.83 & 2.11 & 1.83 & 2.12\
& & Cln{$180^\circ,90^\circ$} & & 2.56 & 2.99 & 2.55 & 2.74 & 2.36 & 2.72 & & 2.12 & 2.42 & 2.10 & 2.22 & 1.96 & 2.20\
& & Frank & & 2.45 & 2.77 & 2.48 & 2.60 & 2.10 & 2.65 & & 2.05 & 2.26 & 2.07 & 2.10 & 1.76 & 2.14\
RMSE & normal & BVN & & 2.54 & 3.28 & 2.52 & 3.13 & 2.45 & 3.14 & & 3.58 & 3.16 & 3.59 & 3.28 & 3.81 & 3.28\
& & Cln{$0^\circ,270^\circ$} & & 2.57 & 3.33 & 2.57 & 3.15 & 2.47 & 3.16 & & 3.69 & 3.28 & 3.71 & 3.30 & 3.78 & 3.31\
& & Cln{$180^\circ,90^\circ$} & & 3.21 & 4.26 & 3.05 & 3.31 & 2.64 & 3.23 & & 3.27 & 2.94 & 3.36 & 3.49 & 4.15 & 3.38\
& & Frank & & 5.15 & 6.66 & 4.91 & 4.37 & 3.62 & 3.89 & & 2.81 & 3.36 & 2.83 & 4.72 & 5.24 & 4.30\
& beta & BVN & & 4.69 & 4.33 & 4.70 & 4.37 & 4.83 & 4.34 & & 2.00 & 2.27 & 1.98 & 2.20 & 1.92 & 2.21\
& & Cln{$0^\circ,270^\circ$} & & 4.97 & 4.63 & 5.01 & 4.27 & 4.64 & 4.29 & & 2.00 & 2.28 & 1.99 & 2.20 & 1.94 & 2.22\
& & Cln{$180^\circ,90^\circ$} & & 6.15 & 6.01 & 5.95 & 4.32 & 4.83 & 4.34 & & 2.51 & 2.83 & 2.40 & 2.23 & 1.96 & 2.24\
& & Frank & & 7.84 & 8.03 & 7.64 & 3.33 & 3.23 & 3.59 & & 3.60 & 4.07 & 3.47 & 2.88 & 2.61 & 2.64\

Cln{$\omega_1^\circ,\omega_2^\circ$}: The bivariate copulas are the Clayton rotated by $\omega_1^\circ$ and $\omega_2^\circ$ to handle the positive and negative dependencies, respectively.

Conclusions from the values in the table and the Tables 1 and 2 in the Supplementary material 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 $\bfpi_1,\bfpi_2$ are not robust to margin misspecification, e.g., in Table \[sim1\] where the true univariate margins are normal (beta) the scaled biases for the MLEs of $\pi_{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 SDs of $\bfpi_1,\bfpi_2$ are not robust to margin misspecification, e.g., in Table \[sim1\] where the true univariate margins are normal the scaled SDs for the MLEs of $\pi_{01}$ for the various one-factor copula mixed models with normal margins range from $3.13$ to $3.61$, while for the various one-factor copula mixed models with beta margins range from $2.77$ to $3.10$.

-   The MLEs of $\bfpi_1,\bfpi_2$ 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 \[sim1\] the scaled biases for the MLEs of $\pi_{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 SDs of $\bfpi_1,\bfpi_2$ are rather robust to bivariate copula misspecification, but increase when the assumed bivariate copulas have different tail dependence behaviour, e.g., in Table \[sim1\] the scaled SDs for the MLEs of $\pi_{11}$ for the various one-factor copula mixed models with normal (beta) margins increase to $2.90$ (2.34) and $3.58$ (2.79) when rotated Clayton copulas with opposite direction tail dependence and Frank copulas with tail independence, respectively, are called.

\[miss\]Misspecification of the copula-mixed model
--------------------------------------------------

We randomly generate 10,000 samples from the D-vine copula mixed model with both normal and beta margins using the algorithm in [@nikoloulopoulos-2018-smmr]. We set the sample size $N$, the study size $n$, the true univariate and Kendall’s $\tau$ parameters, and the disease prevalence to mimic the rheumatoid arthritis data in [@Nishimura-etal2007]. The true (simulated) D-vine copula mixed model uses a Clayton rotated by 270$^\circ$ and Clayton copula to join the (latent) sensitivity and specificity of the first and second test, respectively.

Table \[sim2\] contains the resultant biases, root mean square errors (RMSEs) and standard deviations (SDs), along with average theoretical variances, scaled by 100, for the MLEs of the main parameters of interest under different copula choices and margins. The theoretical variances of the MLEs are obtained via the gradients and the Hessian that were computed numerically during the maximization process. The results for the variability parameters are presented in Table 3 and Table 4 of the Supplementary material.

  ----------------- -------- --------------------------- -- ------------ ------------ ------------ ------------ -- ------------ ------------ ------------ ------------
                                                                                                                                                          
                                                                                                                                                          
                    margin   copula                          $\pi_{11}$   $\pi_{12}$   $\pi_{01}$   $\pi_{02}$      $\pi_{11}$   $\pi_{12}$   $\pi_{01}$   $\pi_{02}$
  Bias              normal   BVN                               -0.07        -0.29        -0.21        -0.05            1.95         1.74         5.08         0.78
                             Cln{$0^\circ,270^\circ$}          -0.10        -0.29        -0.32        -0.06            2.00         1.79         4.99         0.78
                             Cln{$180^\circ,90^\circ$}         -0.68        -0.87        -0.10        -0.03            1.43         1.23         5.17         0.80
                             Frank                             -1.33        -1.47         0.06         0.03            0.53         0.34         5.29         0.84
                    beta     BVN                               -1.58        -1.71        -4.31        -1.22            0.11        -0.05         0.03        -0.02
                             Cln{$0^\circ,270^\circ$}          -1.64        -1.75        -4.46        -1.24            0.13        -0.03        -0.06        -0.02
                             Cln{$180^\circ,90^\circ$}         -2.40        -2.49        -4.26        -1.23           -0.68        -0.83         0.08        -0.01
                             Frank                             -2.64        -2.73        -4.14        -1.15           -1.18        -1.33         0.25         0.07
  SD                normal   BVN                                3.61         3.48         3.39         0.84            3.69         3.64         3.18         0.69
                             Cln{$0^\circ,270^\circ$}           3.62         3.50         3.42         0.84            3.70         3.65         3.23         0.69
                             Cln{$180^\circ,90^\circ$}          3.89         3.76         3.49         0.86            3.96         3.91         3.27         0.71
                             Frank                              4.06         3.93         3.58         0.88            4.17         4.15         3.31         0.72
                    beta     BVN                                3.35         3.25         3.58         1.06            3.38         3.35         3.27         0.74
                             Cln{$0^\circ,270^\circ$}           3.37         3.28         3.65         1.10            3.37         3.33         3.32         0.75
                             Cln{$180^\circ,90^\circ$}          3.72         3.62         3.61         1.09            3.71         3.69         3.31         0.75
                             Frank                              3.74         3.67         3.83         1.13            3.83         3.80         3.42         0.76
  $\sqrt{\bar V}$   normal   BVN                                3.38         3.31         3.22         0.80            3.48         3.46         2.88         0.64
                             Cln{$0^\circ,270^\circ$}           3.22         3.15         3.04         0.77            3.28         3.25         2.67         0.62
                             Cln{$180^\circ,90^\circ$}          3.26         3.19         3.15         0.81            3.39         3.37         2.79         0.64
                             Frank                              3.34         3.25         3.11         0.78            3.45         3.41         2.77         0.62
                    beta     BVN                                3.14         3.08         3.17         0.94            3.21         3.19         3.07         0.72
                             Cln{$0^\circ,270^\circ$}           2.97         2.91         2.89         0.88            3.00         2.98         2.74         0.70
                             Cln{$180^\circ,90^\circ$}          3.04         2.98         3.06         0.95            3.13         3.12         2.99         0.73
                             Frank                              3.07         3.01         3.01         0.90            3.15         3.13         2.86         0.69
  RMSE              normal   BVN                                3.61         3.49         3.40         0.84            4.17         4.03         5.99         1.04
                             Cln{$0^\circ,270^\circ$}           3.62         3.51         3.43         0.85            4.21         4.07         5.94         1.04
                             Cln{$180^\circ,90^\circ$}          3.95         3.86         3.49         0.86            4.21         4.10         6.12         1.07
                             Frank                              4.27         4.19         3.58         0.88            4.20         4.16         6.24         1.11
                    beta     BVN                                3.71         3.67         5.60         1.62            3.38         3.35         3.27         0.74
                             Cln{$0^\circ,270^\circ$}           3.74         3.71         5.76         1.65            3.37         3.33         3.32         0.75
                             Cln{$180^\circ,90^\circ$}          4.43         4.40         5.59         1.64            3.77         3.78         3.32         0.75
                             Frank                              4.58         4.57         5.64         1.61            4.01         4.03         3.43         0.76
  ----------------- -------- --------------------------- -- ------------ ------------ ------------ ------------ -- ------------ ------------ ------------ ------------

  : \[sim2\]Small sample of sizes $N = 22$ simulations ($10^4$ replications, $n_q=25$) from the D-vine copula mixed model and biases, root mean square errors (RMSEs) and standard deviations (SDs), along with the square root of the average theoretical variances ($\sqrt{\bar V}$), scaled by 100, for the MLEs of the main parameters of interest of the one-factor copula mixed model under different copula choices and margins. The true (simulated) D-vine copula mixed model uses a Clayton rotated by 270$^\circ$ and Clayton copula to join the (latent) sensitivity and specificity of the first and second test, respectively.

Cln{$\omega_1^\circ,\omega_2^\circ$}: The bivariate copulas are the Clayton rotated by $\omega_1^\circ$ and $\omega_2^\circ$ to handle the positive and negative dependencies, respectively.

From the table it is seen that the one-factor copula mixed model with normal (beta) margins led to unbiased and efficient estimates when the true model is the quadrivariate vine copula mixed model with normal (margins) when the bivariate copulas are a combination of Clayton and rotated Clayton by $270^\circ$ to model the 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. This is due the fact that the one-factor copula can be explained as a truncated Canonical vine rooted at the latent variable . That is the 1-factor copula has bivariate copulas in the 1st level of the vine and independence copulas in all the remaining levels of the vine (truncated after the 1st level). [@joeetal10] show 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,\ldots,2T-1$ to 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. 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 \[times\] summarizes the computing times (averaged over 20 replications) in seconds. Clearly the D-vine copula mixed approach requires a much higher computing time. Note in passing that for the estimation of the D-vine copula mixed model we even used parallel programming and 22 (one core for each observation) instead of 1 core that used in the estimation of the one-factor copula mixed model. 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.

   $n_q$   factor    vine
  ------- -------- ---------
    15      35.4     799.9
    30      65.7    7355.4
    50     126.2    42997.6

  : \[times\]Small sample of size $N=22$ simulations (20 replications) from the quadrivariate D-vine copula mixed model and computing times (averaged over 20 replications) in seconds of the one-factor and quadrivariate D-vine copula mixed approaches.

\[app-sec\] Application 
========================

We apply the one-factor copula mixed model in order to determine whether anti-cyclic citrullinated peptide (anti-CCP) antibody identifies more accurately patients with rheumatoid arthritis than rheumatoid factor (RF) does. We include $N=22$ studies that assessed both RF and anti-CCP2 antibody for diagnosing rheumatoid arthritis [@Nishimura-etal2007; @dimou-etal2016; @nikoloulopoulos-2018-smmr].

We fit the one-factor copula mixed model for all choices of parametric families of copulas and margins. Because the number of parameters is the same between the models, we use the log-likelihood at estimates as a rough diagnostic measure for goodness of fit between the models. To make it easier to compare strengths of dependence, we convert from the BVN, Frank and (rotated) Clayton $\theta$’s to $\tau$’s via the relations in (\[tauBVN\]), (\[tauFrank\]), and (\[tauCln\]).

[lccccccccccccccc]{} & & & & & & & & & &\
& & Est. & SE & & Est. & SE & & Est. & SE & & Est. & SE & & Est. & SE\
\
$\pi_{11}$ & & 0.681 & 0.034 & & 0.660 & 0.033 & & 0.678 & 0.033 & & 0.681 & 0.036 & & 0.676 & 0.034\
$\pi_{12}$ & & 0.684 & 0.034 & & 0.655 & 0.031 & & 0.673 & 0.032 & & 0.675 & 0.034 & & 0.674 & 0.035\
$\pi_{01}$ & & 0.825 & 0.033 & & 0.834 & 0.032 & & 0.827 & 0.033 & & 0.826 & 0.033 & & 0.827 & 0.033\
$\pi_{02}$ & & 0.960 & 0.008 & & 0.962 & 0.000 & & 0.960 & 0.008 & & 0.960 & 0.008 & & 0.960 & 0.008\
$\s_{11}$ & & 0.685 & 0.128 & & 0.698 & 0.134 & & 0.691 & 0.122 & & 0.722 & 0.133 & & 0.687 & 0.129\
$\s_{12}$ & & 0.697 & 0.124 & & 0.675 & 0.123 & & 0.657 & 0.112 & & 0.687 & 0.121 & & 0.722 & 0.134\
$\s_{01}$ & & 1.028 & 0.181 & & 1.028 & 0.177 & & 1.037 & 0.183 & & 1.029 & 0.181 & & 1.027 & 0.178\
$\s_{02}$ & & 0.790 & 0.175 & & 0.795 & 0.164 & & 0.794 & 0.184 & & 0.792 & 0.170 & & 0.797 & 0.171\
$\tau_{11}$ & & 0.644 & 0.168 & & 0.680 & 0.119 & & 0.719 & 0.137 & & 0.716 & 0.124 & & 0.818 & 0.223\
$\tau_{12}$ & & 0.802 & 0.395 & & 0.839 & 0.152 & & 0.750 & 0.149 & & 0.826 & 0.144 & & 0.466 & 0.136\
$\tau_{01}$ & & -0.125 & 0.168 & & -0.218 & 0.160 & & -0.149 & 0.161 & & -0.213 & 0.148 & & -0.227 & 0.162\
$\tau_{02}$ & & -0.201 & 0.182 & & -0.289 & 0.183 & & -0.228 & 0.333 & & -0.272 & 0.203 & & -0.278 & 0.221\
$-\log(L)$ & & & & & & & & & &\
Beta margins & & & & & & & & & & & & & & &\
$\pi_{11}$ & & 0.667 & 0.031 & & 0.648 & 0.032 & & 0.664 & 0.033 & & 0.665 & 0.031 & & 0.661 & 0.032\
$\pi_{12}$ & & 0.670 & 0.032 & & 0.646 & 0.032 & & 0.661 & 0.032 & & 0.661 & 0.030 & & 0.658 & 0.033\
$\pi_{01}$ & & 0.782 & 0.034 & & 0.789 & 0.033 & & 0.783 & 0.034 & & 0.784 & 0.033 & & 0.785 & 0.033\
$\pi_{02}$ & & 0.949 & 0.009 & & 0.950 & 0.009 & & 0.949 & 0.009 & & 0.949 & 0.009 & & 0.949 & 0.009\
$\s_{11}$ & & 0.087 & 0.028 & & 0.092 & 0.030 & & 0.089 & 0.030 & & 0.097 & 0.029 & & 0.089 & 0.029\
$\s_{12}$ & & 0.091 & 0.028 & & 0.092 & 0.027 & & 0.083 & 0.028 & & 0.091 & 0.026 & & 0.098 & 0.032\
$\s_{01}$ & & 0.132 & 0.039 & & 0.132 & 0.039 & & 0.133 & 0.039 & & 0.132 & 0.039 & & 0.130 & 0.039\
$\s_{02}$ & & 0.025 & 0.012 & & 0.026 & 0.013 & & 0.025 & 0.012 & & 0.026 & 0.013 & & 0.027 & 0.013\
$\tau_{11}$ & & 0.635 & 0.226 & & 0.937 & 0.004 & & 0.723 & 0.140 & & 0.731 & 0.128 & & 0.815 & 0.231\
$\tau_{12}$ & & 0.849 & 0.644 & & 0.651 & 0.103 & & 0.764 & 0.168 & & 0.811 & 0.126 & & 0.497 & 0.134\
$\tau_{01}$ & & -0.111 & 0.169 & & -0.175 & 0.167 & & -0.120 & 0.164 & & -0.217 & 0.144 & & -0.234 & 0.173\
$\tau_{02}$ & & -0.203 & 0.179 & & -0.195 & 0.187 & & -0.212 & 0.290 & & -0.248 & 0.192 & & -0.278 & 0.221\
$-\log(L)$ & & & & & & & & & &\

Cln{$\omega_1^\circ,\omega_2^\circ$}: The bivariate copulas are the Clayton rotated by $\omega_1^\circ$ and $\omega_2^\circ$ to handle the positive and negative dependencies, respectively.

The log-likelihoods showed that an one-factor copula mixed model with Clayton and Clayton rotated by 270$^\circ$ degrees copulas with normal margins to join the factor with each of the sensitivities/specificities provides the best fit (Table \[arthritis-res\]). 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-2018-smmr], 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.

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$^\circ$ reveals that there is tail dependence among the latent sensitivities and specificities.

\[discussion\]Discussion
========================

We have proposed an one-factor copula mixed model for joint meta-analysis and comparison of multiple diagnostic tests. This is the most general meta-analytic model, with univariate parameters separated from dependence parameters. Our general model includes the $2T$-variate GLMM as a special case and can provide an improvement over the latter as the random effects distribution is expressed via an one-factor copula that allows for flexible dependence modelling, different from assuming simple linear correlation structures and normality.

We have proposed a fast and efficient ML estimation technique based on dependent Gauss-Legendre quadrature points that have an one-factor copula distribution. 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. The proposed model can form the vehicle for conducting meta-analysis of comparative accuracy studies with three or more tests.

In practice, comparative accuracy studies with paired designs where each test is applied to the same patients should report the data as separate $2\times 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\times 2$ tables of test results as in Table \[2times2T\]. Comparative accuracy studies should rightly be designed so that patients receive each test in order to reduce biases and ensure the clinical relevance of the resulting inferences [@trikalinos-etal-2014-rsm].

As the number of tests $T$ increases, it becomes less likely that each patient will receive each test and have a recorded result. This will create missing data. For a single test reporting continuous results, primary studies usually provide test performance at multiple but often different thresholds . Hence, similarly with the meta-analysis of studies evaluating multiple tests, meta-analysis of studies evaluating multiple thresholds has multiple pairs of sensitivities and specificities with the different thresholds also creating missing data. [@Riley-etal-2015] and [@ensor-etal-2018-rsm] proposed methods that impute missing threshold results by randomly selecting from the set of all possible discrete combinations which lie between the results for 2 known bounding thresholds. Future research will focus on developing methods for handling the missing data for the case of different multiple diagnostic tests.

Software {#software .unnumbered}
========

[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-2018-CopulaREMADA].

Acknowledgements {#acknowledgements .unnumbered}
================

The simulations presented in this paper were carried out on the High Performance Computing Cluster supported by the Research and Specialist Computing Support service at the University of East Anglia.

Bennett, C., Guo, M., and Dharmage, S. (2007). Hba1c as a screening tool for detection of type 2 diabetes: A systematic review. , 24(4):333–343.

Chu, H. and Cole, S. R. (2006). Bivariate meta-analysis of sensitivity and specificity with sparse data: a generalized linear mixed model approach. , 59(12):1331–1332.

Dimou, N. L., Adam, M., and Bagos, P. G. (2016). A multivariate method for meta-analysis and comparison of diagnostic tests. , 35(20):3509–3523.

Dukic, V. and Gatsonis, C. ([2003]{}). . , [59]{}([4]{}):[936–946]{}.

Ensor, J., Deeks, J., Martin, E., and Riley, R. (2018). Meta-analysis of test accuracy studies using imputation for partial reporting of multiple thresholds. , 9(1):100–115.

Genest, C. (1987). Frank’s family of bivariate distributions. , 74(3):549–555.

Genest, C. and MacKay, J. (1986). The joy of copulas: bivariate distributions with uniform marginals. , 40(4):280–283.

Hoyer, A. and Kuss, O. (2018). Meta-analysis for the comparison of two diagnostic tests to a common gold standard: A generalized linear mixed model approach. , 27(5):1410–1421.

Hult, H. and Lindskog, F. ([2002]{}). . , [34]{}:[587–608]{}.

Joe, H., Li, H., and Nikoloulopoulos, A. K. (2010). Tail dependence functions and vine copulas. , 101:252–270.

Kodama, S., Horikawa, C., Fujihara, K., Hirasawa, R., Yachi, Y., Yoshizawa, S., Tanaka, S., Sone, Y., Shimano, H., Iida, K., Saito, K., and Sone, H. (2013). Use of high-normal levels of haemoglobin a1c and fasting plasma glucose for diabetes screening and for prediction: A meta-analysis. , 29(8):680–692.

Krupskii, P. and Joe, H. (2013). Factor copula models for multivariate data. , 120:85–101.

Ma, X., Nie, L., Cole, S. R., and Chu, H. (2016). Statistical methods for multivariate meta-analysis of diagnostic tests: An overview and tutorial. , 25:1596–1619.

Nash, J. (1990). . Hilger, New York. 2nd edition.

Nikoloulopoulos, A. K. (2015). A mixed effect model for bivariate meta-analysis of diagnostic test accuracy studies using a copula representation of the random effects distribution. , 34:3842–3865.

Nikoloulopoulos, A. K. (2017). A vine copula mixed effect model for trivariate meta-analysis of diagnostic test accuracy studies accounting for disease prevalence. , 26(5):2270–2286.

Nikoloulopoulos, A. K. (2018a). Hybrid copula mixed models for combining case-control and cohort studies in meta-analysis of diagnostic tests. , 27(8):2540–2553.

Nikoloulopoulos, A. K. (2018b). On composite likelihood in bivariate meta-analysis of diagnostic test accuracy studies. , 102:211–227.

Nikoloulopoulos, A. K. (2019a). . package version 1.3. URL: <http://CRAN.R-project.org/package=CopulaREMADA>.

Nikoloulopoulos, A. K. (2019b). A [D]{}-vine copula mixed model for joint meta-analysis and comparison of diagnostic tests. , 28(10-11):3286–3300.

Nikoloulopoulos, A. K. (2020a). An extended trivariate vine copula mixed model for meta-analysis of diagnostic studies in the presence of non-evaluable outcomes. . [DOI:10.1515/ijb-2019-0107](http://dx.doi.org/10.1515/ijb-2019-0107).

Nikoloulopoulos, A. K. (2020b). A multinomial quadrivariate [D]{}-vine copula mixed model for meta-analysis of diagnostic studies in the presence of non-evaluable subjects. . [DOI:10.1177/0962280220913898](http://dx.doi.org/10.1177/0962280220913898).

Nikoloulopoulos, A. K. and Joe, H. (2015). Factor copula models for item response data. , 80:126–150.

Nishimura, K., Sugiyama, D., Kogata, Y., Tsuji, G., Nakazawa, T., Kawano, S., Saigo, K., Morinobu, A., Koshiba, M., Kuntz, K., Kamae, I., and Kumagai, S. (2007). Meta-analysis: Diagnostic accuracy of anti-cyclic citrullinated peptide antibody and rheumatoid factor for rheumatoid arthritis. , 146(11):797–808.

Noma, H., Nagashima, K., and Furukawa, T. A. (2020). Permutation inference methods for multivariate meta-analysis. , 76(1):337–347.

Reitsma, J. B., Glas, A. S., Rutjes, A. W. S., Scholten, R. J. P. M., Bossuyt, P. M., and Zwinderman, A. H. (2005). Bivariate analysis of sensitivity and specificity produces informative summary measures in diagnostic reviews. , 58(10):982–990.

Riley, R. D., Ahmed, I., Ensor, J., Takwoingi, Y., Kirkham, A., Morris, R. K., Noordzij, J. P., and Deeks, J. J. ([2015]{}). . , [4]{}:[12]{}.

Stroud, A. H. and Secrest, D. (1966). . Prentice-Hall, Englewood Cliffs, NJ.

Takwoingi, Y., Leeflang, M., and Deeks, J. (2013). Empirical evidence of the importance of comparative studies of diagnostic test accuracy. , 158(7):544–554.

Trikalinos, T. A., Hoaglin, D. C., Small, K. M., Terrin, N., and Schmid, C. H. (2014). Methods for the joint meta-analysis of multiple tests. , 5(4):294–312.

[^1]: <a.nikoloulopoulos@uea.ac.uk>, School of Computing Sciences, University of East Anglia, Norwich NR4 7TJ, U.K.
