The current statistical procedures implemented in statistical software packages for pooling of diagnostic test accuracy data include hSROC (Rutter and Gatsonis 2001) regression and the bivariate random-effects meta-analysis model (BRMA) ((Reitsma et al. 2005), (Arends et al. 2008), (Chu and Cole 2006), (Riley, Abrams, Sutton, et al. 2007)). However, these models do not report the overall mean but rather the mean for a central study with random-effect equal to zero and have difficulties estimating the correlation between sensitivity and specificity when the number of studies in the meta-analysis is small and/or when the between-study variance is relatively large (Riley, Abrams, Lambert, et al. 2007).

This tutorial on advanced statistical methods for meta-analysis of diagnostic accuracy studies discusses and demonstrates Bayesian modeling using CopulaDTA (Nyaga 2023) package in R (R Core Team 2022) to fit different models to obtain the meta-analytic parameter estimates. The focus is on the joint modelling of sensitivity and specificity using copula based bivariate beta distribution. Essentially,

- we present the Bayesian approach which offers flexibility and ability to perform complex statistical modelling even with small data sets and
- include covariate information, and
- provide an easy to use code.

The statistical methods are illustrated by re-analysing data of two published meta-analyses.

Modelling sensitivity and specificity using the bivariate beta distribution provides marginal as well as study-specific parameter estimates as opposed to using bivariate normal distribution (e.g., in BRMA) which only yields study-specific parameter estimates. Moreover, copula based models offer greater flexibility in modelling different correlation structures in contrast to the normal distribution which allows for only one correlation structure.

In a systematic review of diagnostic test accuracy, the statistical analysis section aims at estimating the average (across studies) sensitivity and specificity of a test and the variability thereof, among other measures. There tends to be a negative correlation between sensitivity and specificity, which postulates the need for correlated data models. The analysis is statistically challenging because the user

- deals with two summary statistics,
- has to account for correlation between sensitivity and specificity,
- has to account for heterogeneity in sensitivity and specificity across the studies and
- should be allowed to incorporate covariates.

Currently, the HSROC (Rutter and Gatsonis 2001) regression or the bivariate random-effects meta-analysis model (BRMA) ((Reitsma et al. 2005), (Arends et al. 2008), (Chu and Cole 2006), (Riley, Abrams, Sutton, et al. 2007)) are recommended for pooling of diagnostic test accuracy data. These models fit a bivariate normal distribution which allows for only one correlation structure to the logit transformed sensitivity and specificity. The resulting distribution has no closed form and therefore the mean sensitivity and specificity is only estimated after numerical integration or other approximation methods, an extra step which is rarely taken.

Within the maximum likelihood estimation methods, the BRMA and HSROC models have difficulties in convergence and estimating the correlation parameter when the number of studies in the meta-analysis are small and/or when the between-study variances are relatively large (Takwoingi et al. 2015). When the correlation is close to the boundary of its parameter space, the between study variance estimates from the BRMA are upwardly biased as they are inflated to compensate for the range restriction on the correlation parameter (Riley, Abrams, Lambert, et al. 2007). This occurs because the maximum likelihood estimator truncates the between-study covariance matrix on the boundary of its parameter space, and this often occurs when the within-study variation is relatively large or the number of studies is small (Riley, Abrams, Sutton, et al. 2007).

The BRMA and HSROC assume that the transformed data is approximately normal with constant variance, however for sensitivity and specificity and proportions in generals, the mean and variance depend on the underlying probability. Therefore, any factor affecting the probability will change the mean and the variance. This implies that the in models where the predictors affect the mean but assume a constant variance will not be adequate.

Joint modelling of study specific sensitivity and specificity using existing or copula based bivariate beta distributions overcomes the above mentioned difficulties. Since both sensitivity and specificity take values in the interval space (0, 1), it is a more natural choice to use a beta distribution to describe their distribution across studies, without the need for any transformation. The beta distribution is conjugate to the binomial distribution and therefore it is easy to integrate out the random-effects analytically giving rise to the beta-binomial marginal distributions. Moreover no further integration is needed to obtain the meta-analytically pooled sensitivity and specificity. Previously, (Cong, Cox, and Cantor 2007) fitted separate beta-binomial models to the number of true positives and the number of false positives. While the model ignores correlation between sensitivity and specificity, (Cong, Cox, and Cantor 2007) reported that the model estimates are comparable to those from the SROC model (Moses, Shapiro, and Littenberg 1993), the predecessor of the HSROC model.

Ignoring the correlation would have negligible influence on the meta-analysis results when the within-study variability is large relative to the between-study variability (Riley 2009). By utilising the correlation, we allow borrowing strength across sensitivities and specificities resulting in smaller standard errors. The use of copula based mixed models within the frequentist framework for meta-analysis of diagnostic test accuracy was recently introduced by (Nikoloulopoulos 2015) who evaluated the joint density numerically.

This tutorial, presents and demonstrates hierarchical mixed models for meta-analysis of diagnostic accuracy studies. In the first level of the hierarchy, given sensitivity and specificity for each study, two binomial distributions are used to describe the variation in the number of true positives and true negatives among the diseased and healthy individuals, respectively. In the second level, we model the unobserved sensitivities and specificities using a bivariate distribution. While hierarchical models are used, the focus of meta-analysis is on the pooled average across studies and rarely on a given study estimate.

The methods are demonstrated using datasets from two previously published meta-analyses:

- on diagnostic accuracy of telomerase in urine as a tumour marker for the diagnosis of primary bladder cancer since it is a problematic dataset that has convergence issues caused by the correlation parameter being estimated to be -1 and has no covariate (Glas et al. 2003) and
- on the comparison of the sensitivity and specificity of human papillomavirus testing (using the HC2 assay) versus repeat cytology to triage women with minor cytological cervical lesions to detect underlying cervical precancer(Arbyn et al. 2013). The second dataset is used to demonstrate meta-regression with one covariate which can be naturally extended to include several covariates.

A bivariate copula function describes the dependence structure between two random variables. Two random variables \(X_1\) and \(X_2\) are joined by a copula function C if their joint cumulative distribution function can be written as \[F(x_1, x_2) = C(F_1 (x_1), F_2 (x_2 )), -\infty \le x_1, x_2 \le +\infty\]

According to the theorem of (Sklar 1959), there exists for every bivariate (multivariate in extension) distribution a copula representation C which is unique for continuous random variables. If the joint cumulative distribution function and the two marginals are known, then the copula function can be written as \[C(u, ~v) = F(F_1^{-1} (u), ~F_2^{-1} (v)),~ 0 \le~ u, ~v ~\le~ 1\]

A 2-dimensional copula is in fact simply a 2-dimensional cumulative function restricted to the unit square with standard uniform marginals. A comprehensive overview of copulas and their mathematical properties can be found in (Nelsen 2006). To obtain the joint probability density, the joint cumulative distribution \(F(x_1, x_2)\) should be differentiated to yield \[ f(x_1, ~x_2) = f_1(x_1) ~f_2(x_2 ) ~c(F_1(x_1), ~F_2 (x_2)) \] where \(f_1\) and \(f_2\) denote the marginal density functions and c the copula density function corresponding to the copula cumulative distribution function C. Therefore from equation above, a bivariate probability density can be expressed using the marginal and the copula density, given that the copula function is absolutely continuous and twice differentiable.

When the functional form of the marginal and the joint densities are known, the copula density can be derived as follows \[c(F_1(x_1), ~F_2(x_2)) = \frac{f(x_1, ~x_2)}{f_1 (x_1 )~ f_2 (x_2 )}\]

While our interest does not lie in finding the copula function, the equations above serve to show how one can move from the copula function to the bivariate density or vice-versa, given that the marginal densities are known. The decompositions allow for constructions of other and possible better models for the variables than would be possible if we limited ourselves to only existing standard bivariate distributions.

We finish this section by mentioning an important implication when Sklar’s theorem is extended to a meta-regression setting with covariates. According to (Patton 2006), it is important that the conditioning variable remains the same for both marginal distributions and the copula, as otherwise the joint distribution might not be properly defined. This implies that covariate information should be introduced in both the marginals and the association parameters of the model.

Since there are two sources of heterogeneity in the data, the within- and between-study variability, the parameters involved in a meta-analysis of diagnostic accuracy studies vary at two levels. For each study \(i\), \(i = 1, ..., n\), let \(Y_{i}~=~(Y_{i1},~ Y_{i2})\) denote the true positives and true negatives, \(N_{i}~=~( N_{i1},~ N_{i2})\) the diseased and healthy individuals respectively, and \(\pi_{i}~ =~ (\pi_{i1},~ \pi_{i2})\) represent the `unobserved’ sensitivity and specificity respectively.

Given study-specific sensitivity and specificity, two separate binomial distributions describe the distribution of true positives and true negatives among the diseased and the healthy individuals as follows \[Y_{ij}~ |~ \pi_{ij}, ~\textbf{x}_i~ \sim~ bin(\pi_{ij},~ N_{ij}), i~=~1, ~\dots ~n, ~j~=~1, ~2\]

where \(\textbf{x}_i\) generically
denotes one or more covariates, possibly affecting \(\pi_{ij}\). The equation above forms the
higher level of the hierarchy and models the within-study variability.
The second level of the hierarchy aims to model the between study
variability of sensitivity and specificity while accounting for the
inherent negative correlation thereof, with a bivariate distribution as
follows \[
\begin{pmatrix}
g(\pi_{i1})\\
g(\pi_{i2})
\end{pmatrix} \sim f(g(\pi_{i1}),~ g(\pi_{i2}))~ =~ f(g(\pi_{i1}))
~f(g(\pi_{i2})) ~c(F_1(g(\pi_{i1})),~ F_2(g(\pi_{i2}))),
\]

where \(g(.)\) denotes a transformation
that might be used to modify the (0, 1) range to the whole real line.
While it is critical to ensure that the studies included in the
meta-analysis satisfy the specified entry criterion, there are study
specific characteristics like different test thresholds and other
unobserved differences that give rise to the second source of
variability, the between-study variability. It is indeed the difference
in the test thresholds between the studies that gives rise to the
correlation between sensitivity and specificity. Including study level
covariates allows us to model part of the between-study variability. The
covariate information can and should be used to model the mean as well
as the correlation between sensitivity and specificity.

In the next section we give more details on different bivariate distributions \(f(g(\pi_{i1}),~g(\pi_{i2}))\) constructed using the logit or identity link function \(g(.)\), different marginal densities and/or different copula densities \(c\). We discuss their implications and demonstrate their application in meta-analysis of diagnostic accuracy studies. An overview of suitable parametric families of copula for mixed models for diagnostic test accuracy studies was recently given by Nikoloupolous (2015). Here, we consider five copula functions which can model negative correlation.

Given the density and the distribution function of the univariate and bivariate standard normal distribution with correlation parameter \(\rho \in (-1, 1)\), the bivariate Gaussian copula function and density is expressed (Meyer 2013) as \[C(u, ~v, ~\rho) = \Phi_2(\Phi^{-1}(u),~ \Phi^{-1}(v),~ \rho), \] \[c(u, ~v, ~\rho) =~ \frac{1}{\sqrt{1~-~\rho^2}} ~exp\bigg(\frac{2~\rho~\Phi^{-1}(u) ~\Phi^{-1}(v) - \rho^2~ (\Phi^{-1}(u)^2 + \Phi^{-1}(v)^2)}{2~(1 - \rho^2)}\bigg ) \]

The logit transformation is often used in binary logistic regression to relate the probability of “success” (coded as 1, failure as 0) of the binary response variable with the linear predictor model that theoretically can take values over the whole real line. In diagnostic test accuracy studies, the `unobserved’ sensitivities and specificities can range from 0 to 1 whereas their logits = \(log ( \frac{\pi_{ij}}{1~-~ \pi_{ij}})\) can take any real value allowing to use the normal distribution as follows \[ logit (\pi_{ij}) ~\sim~ N(\mu_j, ~\sigma_j) ~<=> ~logit (\pi_{ij}) ~=~ \mu_j ~+~ \varepsilon_{ij}, \] where, \(\mu_j\) is a vector of the mean sensitivity and specificity for a study with zero random effects, and \(\varepsilon_{i}\) is a vector of random effects associated with study \(i\). Now \(u\) is the normal distribution function of \(logit(\pi_{i1}\)) with parameters \(\mu_1\) and \(\sigma_1\), is the normal distribution function of \(logit(\pi_{i2})\) with parameters \(\mu_2\) and \(\sigma_2\), \(\Phi_2\) is the distribution function of a bivariate standard normal distribution with correlation parameter \(\rho \in (-1, ~1)\) and \(\Phi^{-1}\) is the quantile of the standard normal distribution. In terms of \(\rho\), Kendall’s tau is expressed as (\(\frac{2}{\pi}\))arcsin\((\rho)\).

With simple algebra the copula density with normal marginal distributions simplifies to \[ c(u, ~v, ~\rho) = \frac{1}{\sqrt{1 - \rho^2}} ~exp\bigg(\frac{1}{2 ~(1 ~-~ \rho^2)}~ \bigg( \frac{2~\rho~(x - \mu_1)~(y - \mu_2)}{\sigma_1 ~\sigma_2} ~-~ \rho^2 ~\bigg (\frac{{x ~-~ \mu_1}^2}{\sigma_1} ~+~ \frac{{y ~-~ \mu_2}^2}{\sigma_2}\bigg)\bigg ) \bigg ). \]

The product of the copula density above, the normal marginal of \(logit(\pi_{i1}\)) and \(logit(\pi_{i2}\)) form a bivariate normal
distribution which characterize the model by (Reitsma et al. 2005), (Arends et al. 2008), (Chu and Cole 2006), and (Riley, Abrams, Lambert, et al. 2007), the
so-called bivariate random-effects meta-analysis (BRMA) model,
recommended as the appropriate method for meta-analysis of diagnostic
accuracy studies. Study level covariate information explaining
heterogeneity is introduced through the parameters of the marginal and
the copula as follows \[\boldsymbol{\mu}_j =
\textbf{X}_j\textbf{B}_j^{\top}. \]

\(\boldsymbol{X}_j\) is a \(n \times p\) matrix containing the
covariates values for the mean sensitivity(\(j
= 1\)) and specificity(\(j =
2\)). For simplicity, assume that \(\boldsymbol{X}_1 = \boldsymbol{X}_2 =
\boldsymbol{X}\). \(\boldsymbol{B}_j^\top\) is a \(p \times\) 1$ vector of regression
parameters, and \(p\) is the number of
parameters. By inverting the logit functions, we obtain \[
\pi_{ij} = logit^{-1} (\mu_j + \varepsilon_{ij}).
\]

Therefore, the meta-analytic sensitivity and specificity obtained by
averaging over the random study effect, is given by, for \(j = 1, 2\) \[
E(\pi_{j} )~=~E(logit^{-1} (\mu_j ~+~ \varepsilon_{ij}))
~=~\int_{-\infty}^{\infty}logit^{-1}(\mu_j ~+~ \varepsilon_{ij})
f(\varepsilon_{ij},~\sigma_j)~d\varepsilon_{ij},
\]

assuming that \(\sigma_1^2 > 0\) and
\(\sigma_2^2 > 0\). The integration
above has no analytical expression and therefore needs to be numerically
approximated and the standard are not easily available. Using MCMC
simulation in the Bayesian framework the meta-analytic estimates can be
easily computed as well as a standard error estimate and a credible
intervals \(E(\pi_{j})\) with minimum
effort by generating predictions of the fitted bivariate normal
distribution.

In the frequentist framework, it is more convenient however to use
numerical averaging by sampling a large number \(M\) of random-effects \(\hat{\varepsilon}_{ij}\) from the fitted
distribution and to estimate the meta-analytic sensitivity and
specificity by, for \(j = 1, 2\) \[
\hat{E}(\pi_{j}) = \frac{1}{M} \sum_{i~ =~ 1}^{M}logit^{-1}(\hat{\mu}_j
+ \hat{\varepsilon}_{ij}).
\]

However, inference is not straightforward in the frequentist framework
since the standard errors are not available. When \(\varepsilon_{ij} = 0\), then \[
E(\pi_{j} ~|~ \varepsilon_{ij} = 0) ~=~ logit^{-1}(\mu_j).
\]

Inference for \(E(\pi_{j} ~|~\varepsilon_{ij}
~=~ 0)\) as expressed in the equation above can be done in both
Bayesian and frequentist framework. The equation represents the mean
sensitivity and specificity for a ``central” study with \(\varepsilon_{ij} ~=~ 0\). Researchers often
seem to confuse \(E(\pi_{j}
~|~\varepsilon_{ij} ~=~ 0)\) with \(E(\pi_{j})\) but due to the non-linear
logit transformations, they are clearly not the same parameter.

With the identity link function, no transformation on study-specific sensitivity and specificity is performed. A natural choice for and would be beta distribution functions with parameters (\(\alpha_1, ~\beta_1\)) and (\(\alpha_2, ~\beta_2\)) respectively. Since \(\pi_{ij} ~\sim ~ beta(\alpha_j, ~\beta_j)\), the meta-analytic sensitivity and specificity are analytically solved as follows \[ E(\pi_{j} ) = \frac{\alpha_j}{\alpha_j + \beta_j}, \]

After reparameterising the beta distributions using the mean (\(\mu_j ~=~ \frac{\alpha_j}{\alpha_j ~+~
\beta_j}\)) and certainty (\(\psi_j ~=~
\alpha_j ~+~ \beta_j\)) or dispersion (\(\varphi_j ~=~ \frac{1}{1~ +~
\alpha_j~+~\beta_j}\)) parameters different link functions
introduce covariate information to the mean, certainty/dispersion and
association (\(\rho\)) parameters. A
typical model parameterisation is \[
\boldsymbol{\mu}_j ~=~ logit^{-1}(\textbf{XB}_j^{\top}), \\
\boldsymbol{\psi}_j ~=~ g(\textbf{WC}_j^{\top}), \\
\boldsymbol{\alpha}_j ~=~ \boldsymbol{\mu}_{j} ~\circ~
\boldsymbol{\psi}_{j}, \\
\boldsymbol{\beta}_j ~=~ ( \textbf{1} ~-~ \boldsymbol{\mu}_{j}) ~\circ~
\boldsymbol{\psi}_{j}, \\
\boldsymbol{\rho} ~=~ tanh(\textbf{ZD}_j^{\top}) ~=~ \frac{exp(2 \times
\textbf{ZD}_j^{\top}) ~-~ 1}{exp(2\times \textbf{ZD}_j^{\top}) ~+~ 1}.
\]

\(\boldsymbol{X, W}\) and \(\boldsymbol{Z}\) are a \(n \times p\) matrices containing the
covariates values for the mean, dispersion and correlation which we will
assume has similar information and denoted by \(\boldsymbol{X}\) for simplicity purpose,
\(p\) is the number of parameters,
\(\textbf{B}_j^{\top}\) , \(\textbf{V}_j^{\top}\) and \(\textbf{D}_j^{\top}\) are a \(p \times 1\) vectors of regression
parameters relating covariates to the mean, variance and correlation
respectively. \(g(.)\) is the log link
to mapping \(\textbf{XC}_j^{\top}\) to
the positive real number line and \(\circ\) is the Hadamard product.

This flexible copula in the so-called family of Archimedean copulas by (Frank 1979). The functional form of the copula and the density is given by; \[ C(F(\pi_{i1} ), ~F(\pi_{i2} ),\theta) ~=~ - \frac{1}{\theta} ~log \bigg [ 1 + \frac{(e^{-\theta ~F(\pi_{i1})} - 1) (e^{-\theta~ F(\pi_{i2})} - 1)}{e^{-\theta} - 1} \bigg ], \\ c(F(\pi_{i1} ), ~F(\pi_{i2} ),\theta) ~=~ \frac {\theta~(1 - e^{-\theta})~e^{-\theta~(F(\pi_{i1}) ~ +~ F(\pi_{i2}))}}{[1 - e^{-\theta}-(1 - e^{-\theta ~F(\pi_{i1})})~(1 - e^{-\theta ~F(\pi_{i2})})]^2} . \]

Since \(\theta \in \mathbb{R}\),
both positive and negative correlation can be modelled, making this one
of the more comprehensive copulas. When \(\theta\) is 0, sensitivity and specificity
are independent. For \(\theta > 0\),
sensitivity and specificity exhibit positive quadrant dependence and
negative quadrant dependence when \(\theta
< 0\). The Spearman correlation \(\rho_s\) and Kendall’s tau \(\tau_k\) can be expressed in terms of \(\theta\) as \[
\rho_s = 1 - 12~\frac{D_2 (-\theta) ~-~ D_1(-\theta)}{\theta}, \\
\tau_k = 1 + 4~ \frac{D_1(\theta) ~-~ 1}{\theta},
\]

where \(D_j(\delta)\) is the Debye
function defined as \[
D_j(\delta) ~=~ \frac{j}{\delta^j}
\int_{\theta}^{\delta}\frac{t^j}{exp(t) ~-~ 1} ~dt, ~j~= 1, ~2.
\]

Covariate information is introduced in a similar manner as earlier. The identity link is used for the association parameter \(\theta\).

This popular copula by (Farlie 1960),
(Gumbel 1960) and (Morgenstern 1956) is defined as \[
C(F(\pi_{i1} ), ~F(\pi_{i2}),~\theta) ~= ~F(\pi_{i1})~F(\pi_{i2})[1 ~+~
\theta~(1 - F(\pi_{i1}))~(1 ~-~ F(\pi_{i2}))], \nonumber \\
c(F(\pi_{i1} ), ~F(\pi_{i2}),~\theta) ~=~ [1 ~+~ \theta~(2~F(\pi_{i1})
~-~ 1)~(2~F(\pi_{i2}) ~-~ 1)]. \]

Because \(\theta \in (-1, 1)\), the
Spearman correlation and Kendall’s tau are expressed in terms of \(\theta\) as \(\theta/3\) and \(2\theta/9\) respectively, making this
copula only appropriate for data with weak dependence since \(|\rho_s| ~\leq~ 1/3\). The logit link,
log/identity link and Fisher’s z transformation can be used to introduce
covariate information in modelling the mean, dispersion and association
parameter.

The Clayton copula function and density by (Clayton 1978) is defined as \[ C(F(\pi_{i1}), ~F(\pi_{i2}),~\theta) = [F(\pi_{i1})^{- \theta} ~+~ F(\pi_{i2})^{-\theta} ~-~ 1]^{\frac{- 1}{\theta}}, \nonumber \\ c(F(\pi_{i1}), ~F(\pi_{i2}),~\theta) = (1 ~+~ \theta)~F(\pi_{i1})^{- (1 ~+~ \theta)} ~F(\pi_{i2})^{- (1 + \theta)}~[F(\pi_{i1} )^{- \theta} + F(\pi_{i2})^{- \theta} ~-~ 1]^{\frac{- (2~\theta ~+~ 1)}{\theta}}. \]

Since \(\theta \in (0, ~\infty)\),
the Clayton copula typically models positive dependence; Kendall’s tau
equals \(\theta/(\theta ~+~ 2)\).
However, the copula function can be rotated by \(90^\circ\) or \(270^\circ\) to model negative dependence.
The distribution and density functions following such rotations are
given by \[
C_{90} (F(\pi_{i1}), ~F(\pi_{i2}), ~\theta) ~=~ F(\pi_{i2}) ~-~ C(1 ~-~
F(\pi_{i1}), ~F(\pi_{i2}), ~\theta), \nonumber \\
c_{90} (F(\pi_{i1}),~ F(\pi_{i2}),~\theta) ~=~ (1 ~+~ \theta)(1 ~-~
F(\pi_{i1}))^{- (1 ~+~ \theta)}~F(\pi_{i2})^{- (1 ~+~ \theta)} ~[(1 -
F(\pi_{i1}))^{-\theta} \nonumber \\
+~
F(\pi_{i2})^{- \theta} - 1]^{\frac{- (2~\theta ~+~ 1)}{\theta}},
\] and \[
C_{270} (F(\pi_{i1}), ~F(\pi_{i2}), ~\theta) ~= ~F(\pi_{i1})-
C(F(\pi_{i1}), ~1 ~-~ F(\pi_{i2}),\theta), \nonumber
\\
c_{270} (F(\pi_{i1}), ~F(\pi_{i2}), ~\theta) ~=~(1 ~+~ \theta)
~F(\pi_{i1} )^{- (1 ~+~ \theta)}~ (1 ~-~ F(\pi_{i2} ))^{- (1 ~+~
\theta)}~[F(\pi_{i1})^{- \theta} \nonumber \\
+
(1 ~-~ F(\pi_{i2}))^{- \theta} ~-~ 1]^{\frac{- (2~\theta ~+~
1)}{\theta}}. \]

The logit, log/identity and log/identity links can be used to introduce
covariate information in modelling the mean (\(\mu_j\)), certainty (\(\psi_j\))/dispersion (\(\varphi_j\)) and association (\(\theta\)) parameters respectively.

Of course other copula functions that allow for negative association can be chosen. It is also an option to use known bivariate beta distributions. However, it is not always straightforward and analytically attractive to derive the corresponding copula function for all bivariate distributions. The use of existing bivariate beta distributions in meta-analysis of diagnostic accuracy studies has been limited because these densities model positive association( e.g., (Libby and Novick 1982), (Olkin and Liu 2003)), or both positive and negative association but over a restricted range( e.g., (Sarmanov 1966)).

Within the Bayesian framework, the analyst updates a prior
opinion/information of a parameter based on the observed data whereas in
the frequentist framework, the analyst investigates the behaviour of the
parameter estimates in hypothetical repeated samples from a certain
population. Due to its flexibility and use of MCMC simulations, complex
modelling can often be implemented more easily within the Bayesian
framework. By manipulating the prior distributions, Bayesian inference
can circumvent identifiability problems whereas numerical approximation
algorithms in frequentist inference without prior distributions can
become stuck caused by identifiability problems. However, Bayesian
methods typically require statistical expertise and patience because the
MCMC simulations are computationally intensive. In contrast, most
frequentist methods have been wrapped up in standard “procedures” that
require less statistical knowledge and programming skills. Moreover
frequentist methods are optimized with maximum likelihood estimation
(MLE) that have much shorter run-times as opposed to MCMC simulations.
`CopulaREMADA`

(Nikoloulopoulos
2016) is such a MLE based `R`

package.

The `CopulaDTA`

package is an extension of
`rstan`

(Stan Development Team
2023), the `R`

interface to `Stan`

(Carpenter et al. 2017) for diagnostic test
accuracy data. Stan is a probabilistic programming language which has
implemented Hamilton Monte Carlo(MHC) and uses No-U-Turn sampler (NUTS)
(Hoffman and Gelman 2011). The package
facilitates easy application of complex models and their visualization
within the Bayesian framework with much shorter run-times.

JAGS (Martyn Plummer et al. 2003) is an
alternative extensible general purpose sampling engine to Stan.
Extending `JAGS`

requires knowledge of `C++`

to
assemble a dynamic link library(DLL) module. From experience,
configuring and building the module is a daunting and tedious task
especially in the Windows operation system. The above short-comings
coupled with the fact that `Stan`

tends to converge with
fewer iterations even from bad initial values than `JAGS`

made us prefer the `Stan`

MCMC sampling engine.

The `CopulaDTA`

package is available via the Comprehensive
`R`

Archive Network (CRAN) at https://CRAN.R-project.org/package=CopulaDTA. With a
working internet connection, the `CopulaDTA`

package is
installed and loaded in `R`

with the following commands

```
#install.packages("CopulaDTA", dependencies = TRUE)
library(CopulaDTA)
#> Loading required package: rstan
#> Loading required package: StanHeaders
#>
#> rstan version 2.26.15 (Stan version 2.26.1)
#> For execution on a local, multicore CPU with excess RAM we recommend calling
#> options(mc.cores = parallel::detectCores()).
#> To avoid recompilation of unchanged Stan programs, we recommend calling
#> rstan_options(auto_write = TRUE)
#> For within-chain threading using `reduce_sum()` or `map_rect()` Stan functions,
#> change `threads_per_chain` option:
#> rstan_options(threads_per_chain = 1)
#> Do not specify '-march=native' in 'LOCAL_CPPFLAGS' or a Makevars file
#>
#> Attaching package: 'CopulaDTA'
#> The following objects are masked from 'package:rstan':
#>
#> plot, traceplot
#> The following object is masked from 'package:graphics':
#>
#> plot
#> The following object is masked from 'package:base':
#>
#> plot
```

The `CopulaDTA`

package provide functions to fit bivariate
beta-binomial distributions constructed as a product of two beta
marginal distributions and copula densities discussed earlier. The
package also provides forest plots for a model with categorical
covariates or with intercept only. Given the chosen copula function, a
beta-binomial distribution is assembled up by the `cdtamodel`

function which returns a `cdtamodel`

object. The main
function `fit`

takes the `cdtamodel`

object and
fits the model to the given dataset and returns a `cdtafit`

object for which `print`

, `summary`

and
`plot`

methods are provided for.

To assess model convergence, mixing and stationarity of the chains, it is necessary to check the potential scale reduction factor \(\hat{R}\), effective sample size (ESS), MCMC error and trace plots of the parameters. When all the chains reach the target posterior distribution, the estimated posterior variance is expected to be close to the within chain variance such that the ratio of the two, \(\hat{R}\) is close to 1 indicating that the chains are stable, properly mixed and likely to have reached the target distribution. A large \(\hat{R}\) indicates poor mixing and that more iterations are needed. Effective sample size indicates how much information one actually has about a certain parameter. When the samples are auto correlated, less information from the posterior distribution of our parameters is expected than would be if the samples were independent. ESS close to the total post-warm-up iterations is an indication of less autocorrelation and good mixing of the chains. Simulations with higher ESS have lower standard errors and more stable estimates. Since the posterior distribution is simulated there is a chance that the approximation is off by some amount; the Monte Carlo (MCMC) error. MCMC error close to 0 indicates that one is likely to have reached the target distribution.

Watanabe-Alkaike Information Criterion (WAIC) (Watanabe 2010), a recent model comparison tool to measure the predictive accuracy of the fitted models in the Bayesian framework, will be used to compare the models. WAIC can be viewed as an improvement of Deviance Information Criterion(DIC) which, though popular, is known to be have some problems (M. Plummer 2008). WAIC is a fully Bayesian tool, closely approximates the Bayesian cross-validation, is invariant to reparameterisation and can be used for simple as well as hierarchical and mixture models.

(Glas et al. 2003) systematically
reviewed the sensitivity and specificity of cytology and other markers
including telomerase for primary diagnosis of bladder cancer. They
fitted a bivariate normal distribution to the logit transformed
sensitivity and specificity values across the studies allowing for
heterogeneity between the studies. From the included 10 studies, they
reported that telomerase had a sensitivity and specificity of
`0.75 [0.66, 0.74]`

and `0.86 [0.71, 0.94]`

respectively. They concluded that telomerase was not sensitive enough to
be recommended for daily use. This dataset is available within the
package and the following commands

```
data(telomerase)
telomerase
```

loads the data into the R enviroment and generates the following output

```
#> ID TP TN FN FP
#> 1 1 25 25 8 1
#> 2 2 17 11 4 3
#> 3 3 88 31 16 16
#> 4 4 16 80 10 3
#> 5 5 40 137 17 1
#> 6 6 38 24 9 6
#> 7 7 23 12 19 0
#> 8 8 27 18 6 2
#> 9 9 14 29 3 3
#> 10 10 37 7 7 22
```

`ID`

is the study identifier, `DIS`

is the
number of diseased, `TP`

is the number of true positives,
`NonDis`

is the number of healthy and `TN`

is the
number of true negatives.

(Arbyn et al. 2013) performed a
Cochrane review on the accuracy of human papillomavirus testing and
repeat cytology to triage of women with an equivocal Pap smear to
diagnose cervical precancer. They fitted the BRMA model in
`SAS`

using `METADAS`

on 10 studies where both
tests were used. They reported absolute sensitivity of
`0.909 [0.857, 0.944]`

and `0.715 [0.629, 0.788]`

for HC2 and repeat cytology respectively. The specificity was
`0.607 [0.539, 0.68]`

and `0.684 [0.599, 0.758]`

for HC2 and repeat cytology respectively. These data is used to
demonstrate how the intercept only model is extended in a
meta-regression setting. This dataset is also available within the
package and the following commands

```
data(ascus)
ascus
```

loads the data into the R environment and generates the following output

```
#> Test StudyID TP FP TN FN
#> 1 RepC Andersson 2005 6 14 28 4
#> 2 RepC Bergeron 2000 8 28 71 4
#> 3 RepC Del Mistro 2010 20 191 483 7
#> 4 RepC Kulasingam 2002 20 74 170 6
#> 5 RepC Lytwyn 2000 4 20 26 2
#> 6 RepC Manos 1999 48 324 570 15
#> 7 RepC Monsonego 2008 10 18 168 15
#> 8 RepC Morin 2001 14 126 214 5
#> 9 RepC Silverloo 2009 24 43 105 10
#> 10 RepC Solomon 2001 227 1132 914 40
#> 11 HC2 Andersson 2005 6 17 25 4
#> 12 HC2 Bergeron 2000 10 38 61 2
#> 13 HC2 Del Mistro 2010 27 154 566 2
#> 14 HC2 Kulasingam 2002 23 115 129 3
#> 15 HC2 Lytwyn 2000 4 19 33 1
#> 16 HC2 Manos 1999 58 326 582 7
#> 17 HC2 Monsonego 2008 22 110 72 2
#> 18 HC2 Morin 2001 17 88 253 2
#> 19 HC2 Silverloo 2009 34 65 81 2
#> 20 HC2 Solomon 2001 256 1050 984 11
```

`Test`

is an explanatory variable showing the type of
triage test, `StudyID`

is the study identifier,
`TP`

is the number of true positives, `FP`

is the
number of false positives, `TN`

is the number of true
negatives, `FN`

is the number of false negatives.

The `CopulaDTA`

package has five different correlation
structures that result to five different bivariate beta-binomial
distributions to fit to the data. The correlation structure is specified
by indicating `copula ="gauss"`

or `"fgm"`

or
`"c90"`

or `"270"`

or `"frank"`

in the
`cdtamodel`

function.

`.1 <- cdtamodel("gauss") gauss`

The code for the fitted model and the model options are displayed as follows

```
print(gauss.1)
str(gauss.1)
```

The Gaussian copula bivariate beta-binomial distribution is fitted to
the `telomerase`

data with the following code

```
.1 <- fit(
fitgauss.1,
gaussdata = telomerase,
SID = "ID",
iter = 28000,
warmup = 1000,
thin = 30,
seed = 3)
```

By default, `chains = 3`

and `cores = 3`

and
need not be specified unless otherwise. From the code above,
`28000`

samples are drawn from each of the `3`

chains, the first `1000`

samples are discarded and thereafter
every `30`

\(^{th}\) draw
kept such that each chain has 900 post-warm-up draws making a total of
2700 post-warm-up draws. The seed value, `seed = 3`

,
specifies a random number generator to allow reproducibility of the
results and `cores = 3`

allows for parallel-processing of the
chains by using `3`

cores, one core for each chain. They were
no initial values specified and in that case, the program randomly
generates random values satisfying the parameter constraints.

The trace plots below show satisfactory mixing of the chains and convergence.

`traceplot(fitgauss.1)`

Next, obtain the model summary estimates as follows

```
print(fitgauss.1, digits = 4)
#> Posterior marginal mean and median sensitivity and specificity
#> with 95% credible intervals
#>
#> Parameter Mean Lower Median Upper n_eff Rhat
#> MUse[1] Sensitivity 0.756762 6.904e-01 0.756036 0.81658 1196.9 1.001
#> MUsp[1] Specificity 0.798289 6.171e-01 0.813517 0.90640 704.4 1.004
#> ktau[1] Correlation -0.820176 -9.861e-01 -0.876343 -0.33334 269.2 1.015
#> Varse[1] Var(Sens) 0.006198 8.321e-06 0.005047 0.01947 165.7 1.007
#> Varsp[1] Var(Spec) 0.048111 1.357e-02 0.041060 0.12204 169.5 1.007
#>
#>
#> Model characteristics
#>
#> Copula function: gauss, sampling algorithm: NUTS(diag_e)
#>
#> Formula(1): MUse ~ 1
#> Formula(2): MUsp ~ 1
#> Formula(3): Omega ~ 1
#> 3 chain(s)each with iter=28000; warm-up=1000; thin=30.
#> post-warmup draws per chain=900;total post-warmup draws=2700.
#>
#> Predictive accuracy of the model
#>
#> Log point-wise predictive density (LPPD): -38.0607
#> Effective number of parameters: 7.5807
#> Watanabe-Akaike information Criterion (WAIC): 91.2828
```

From the output above, `n_eff`

and `Rhat`

both
confirm proper mixing of the chains with little autocorrelation. The
meta-analytic sensitivity `MUse[1]`

and specificity
`MUsp[1]`

is `0.7568 [0.6904, 0.8166]`

and
`0.7983 [0.6171, 0.9064]`

respectively. The between-study
variability in sensitivity and specificity is
`0.0062 [0, 0.0195]`

and
`0.0048 [0.0136, 0.1220]`

, respectively. The Kendall’s tau
correlation between sensitivity and specificity is estimated to be
`-0.8202 [-0.9861, -0.3333]`

.

The command below produces a series of forest plots.

```
plot(fitgauss.1)
#> $G1
```

```
#>
#> $G2
```

```
#>
#> $G3
```

`$G1`

is a plot of the study-specific sensitivity and
specificity (magenta points) and their corresponding 95 % exact
confidence intervals (black lines). `$G2`

is a plot of the
posterior study-specific (blue points) and marginal(blue diamond)
sensitivity and specificity and their corresponding 95 % credible
intervals (black lines).

`$G3`

is a plot of the posterior study-specific (blues
stars), and marginal(blue diamond) sensitivity and specificity and their
corresponding 95 % credible intervals (black lines). The study-specific
sensitivity and specificity (magenta points) and their corresponding 95
% exact confidence intervals (thick grey lines) are also presented.

As observed in plots above, the posterior study-specific sensitivity and specificity are less extreme and variable than the “observed” study-specific sensitivity and specificity. In other words, there is “shrinkage” towards the overall mean sensitivity and specificity as studies borrow strength from each other in the following manner: the posterior study-specific estimates depends on the global estimate and thus also on all other the studies.

The other four copula based bivariate beta distributions are fitted as follows

```
.1 <- cdtamodel(copula = "c90")
c90
.1 <- fit(c90.1,
fitc90data=telomerase,
SID="ID",
iter=28000,
warmup=1000,
thin=30,
seed=718117096)
.1 <- cdtamodel(copula = "c270")
c270
.1 <- fit(c270.1,
fitc270data=telomerase,
SID="ID",
iter=19000,
warmup=1000,
thin=20,
seed=3)
.1 <- cdtamodel(copula = "fgm")
fgm
.1 <- fit(fgm.1,
fitfgmdata=telomerase,
SID="ID",
iter=19000,
warmup=1000,
thin=20,
seed=3)
.1 <- cdtamodel(copula = "frank")
frank
.1 <- fit(frank.1,
fitfrankdata=telomerase,
SID="ID",
iter=19000,
warmup=1000,
thin=20,
seed=1959300748)
```

For comparison purpose, the current recommended model; the BRMA,
which uses normal marginals is also fitted to the data though it is not
part of the `CopulaDTA`

package. The model is first expressed
in `Stan`

modelling language in the code below and is stored
within `R`

environment as character string named
`BRMA1`

.

```
<- "
BRMA1 data{
int<lower = 0> N;
int<lower = 0> tp[N];
int<lower = 0> dis[N];
int<lower = 0> tn[N];
int<lower = 0> nondis[N];
}
parameters{
real etarho;
vector[2] mul;
vector<lower = 0>[2] sigma;
vector[2] logitp[N];
vector[2] logitphat[N];
}
transformed parameters{
vector[N] p[2];
vector[N] phat[2];
real MU[2];
vector[2] mu;
real rho;
real ktau;
matrix[2,2] Sigma;
rho = tanh(etarho);
ktau = (2/pi())*asin(rho);
for (a in 1:2){
for (b in 1:N){
p[a][b] = inv_logit(logitp[b][a]);
phat[a][b] = inv_logit(logitphat[b][a]);
}
mu[a] = inv_logit(mul[a]);
}
MU[1] = mean(phat[1]);
MU[2] = mean(phat[2]);
Sigma[1, 1] = sigma[1]^2;
Sigma[1, 2] = sigma[1]*sigma[2]*rho;
Sigma[2, 1] = sigma[1]*sigma[2]*rho;
Sigma[2, 2] = sigma[2]^2;
}
model{
etarho ~ normal(0, 10);
mul ~ normal(0, 10);
sigma ~ cauchy(0, 2.5);
for (i in 1:N){
logitp[i] ~ multi_normal(mul, Sigma);
logitphat[i] ~ multi_normal(mul, Sigma);
}
tp ~ binomial(dis,p[1]);
tn ~ binomial(nondis, p[2]);
}
generated quantities{
vector[N*2] loglik;
for (i in 1:N){
loglik[i] = binomial_lpmf(tp[i]| dis[i], p[1][i]);
}
for (i in (N+1):(2*N)){
loglik[i] = binomial_lpmf(tn[i-N]| nondis[i-N], p[2][i-N]);
}
}
"
```

Next, prepare the data by creating a list as follows

```
= list(
datalist tp = telomerase$TP,
dis = telomerase$TP + telomerase$FN,
tn = telomerase$TN,
nondis = telomerase$TN + telomerase$FP,
N = 10)
```

In the `data`

block the dimensions and names of variables
in the dataset are specified, here `Ns`

indicate the number
of studies in the dataset. The `parameters`

block introduces
the unknown parameters to be estimated. These are `etarho`

; a
scalar representing the Fisher’s transformed form of the association
parameter \(\rho\), `mul`

;a
\(2 \times 1\) vector representing the
mean of sensitivity and specificity on the logit scale for a central
study where the random-effect is zero, `sigma`

; a \(2 \times 1\) vector representing the
between study standard deviation of sensitivity and specificity on the
logit scale, `logitp`

; a \(Ns
\times 2\) array of study-specific sensitivity in the first
column and specificity in the second column on logit scale, and
`logitphat`

; a \(Ns \times
2\) array of predicted sensitivity in the first column and
predicted specificity in the second column on logit scale.

The parameters are further transformed in the
`transformed parameters`

block. Here, `p`

is a
\(2 \times Ns\) array of sensitivity in
the first column and specificity in the second column after inverse
logit transformation of `logitp`

, and `phat`

is a
\(2 \times Ns\) array of predicted
sensitivity in the first column and predicted specificity in the second
column after inverse logit transformation of `logitphat`

to
be used in computing the meta-analytic sensitivity and specificity.
`mu`

is a \(2 \times 1\)
vector representing the mean of sensitivity and specificity for a
certain study with a random effect equal to 0, `MU`

is a
\(2 \times 1\) vector containing the
meta-analytic sensitivity and specificity, `Sigma`

; a $2
\(\times\) 2$ matrix representing the
variance-covariance matrix of sensitivity and specificity on the logit
scale, `rho`

and `ktau`

are scalars representing
the Pearson’s and Kendall’s tau correlation respectively. The prior
distributions for the all parameters and data likelihood are defined in
the `model`

block. Finally, in the
`generated quantities`

block, `loglik`

is a \((2Ns) \times 1\) vector of the log
likelihood needed to compute the WAIC.

Next, call the function `stan`

from the `rstan`

package to translate the code into `C++`

, compile the code
and draw samples from the posterior distribution as follows

```
.1 <- stan(model_code = BRMA1,
brmadata = datalist,
chains = 3,
iter = 5000,
warmup = 1000,
thin = 10,
seed = 3,
cores = 3)
```

The parameter estimates are extracted and the chain convergence and autocorrelation examined further with the following code

```
print(brma.1, pars=c('MU', 'mu', 'rho', "Sigma"), digits=4, prob=c(0.025, 0.5, 0.975))
#> Inference for Stan model: 61572683b29d52354783115614fab729.
#> 3 chains, each with iter=5000; warmup=1000; thin=10;
#> post-warmup draws per chain=400, total post-warmup draws=1200.
#>
#> mean se_mean sd 2.5% 50% 97.5% n_eff Rhat
#> MU[1] 0.7525 0.0018 0.0517 0.6323 0.7562 0.8415 849 0.9999
#> MU[2] 0.7908 0.0034 0.1095 0.5273 0.8094 0.9539 1045 1.0008
#> mu[1] 0.7668 0.0015 0.0388 0.6869 0.7688 0.8369 701 0.9990
#> mu[2] 0.8937 0.0027 0.0753 0.6943 0.9115 0.9825 789 0.9992
#> rho -0.9311 0.0068 0.1353 -0.9996 -0.9813 -0.5626 401 1.0077
#> Sigma[1,1] 0.3376 0.0091 0.2918 0.0579 0.2554 0.9851 1026 1.0023
#> Sigma[1,2] -1.2291 0.0275 0.8765 -3.4195 -1.0031 -0.2724 1017 0.9991
#> Sigma[2,1] -1.2291 0.0275 0.8765 -3.4195 -1.0031 -0.2724 1017 0.9991
#> Sigma[2,2] 5.6827 0.1282 4.1931 1.4720 4.6330 16.9031 1070 1.0002
#>
#> Samples were drawn using NUTS(diag_e) at Mon Oct 09 09:19:55 2017.
#> For each parameter, n_eff is a crude measure of effective sample size,
#> and Rhat is the potential scale reduction factor on split chains (at
#> convergence, Rhat=1).
```

The meta-analytic sensitivity (`MU[1]`

) and specificity
(`MU[2]`

) and 95% credible intervals are
`0.7525[0.6323, 0.8415]`

and
`0.7908[0.5273, 0.9539]`

respectively. This differs from what
the authors published (0.75[0.66, 0.74] and 0.86[0.71, 0.94]) in two
ways. The authors fitted the standard bivariate normal distribution to
the logit transformed sensitivity and specificity values across the
studies allowing for heterogeneity between the studies and disregarded
the higher level of the hierarchical model. Because of this the authors
had to use a continuity correction of 0.5 since the seventh study had
“observed” specificity equal to 1, a problem not encountered in the
hierarchical model. Secondly the authors do not report the meta-analytic
values but rather report the mean sensitivity (`mu[1]`

) and
specificity (`mu[2]`

) for a particular, hypothetical study
with random-effect equal to zero, which in our case is
`0.7668[0.6869, 0.8369]`

and
`0.8937[0.6943, 0.9825]`

respectively and is comparable to
what the authors reported. This discrepancy between `MU`

and
`mu`

will indeed increase with increase in the between study
variability. Here, (`MU[1]`

) and (`mu[1]`

) are
similar because the between-study variability in sensitivity in the
logit scale (`Sigma[1,1]`

) is small
`0.333 [0.0579, 0.9851]`

. (`MU[2]`

) is
substantially smaller than (`mu[2]`

) as a result of the
substantial between-study heterogeneity in specificity in the logit
scale (`Sigma[2,2]`

) =
`5.6827 [1.4720, 16.9031]`

.

The figure below shows satisfactory chain mixing with little autocorrelation for most of the fitted models except the Clayton copula models.

```
.1 <- traceplot(fitgauss.1)
f1.2 <- traceplot(fitc90.1)
f1.3 <- traceplot(fitc270.1)
f1.4 <- traceplot(fitfgm.1)
f1.5 <- traceplot(fitfrank.1)
f1
#draws <- as.array(brma.1, pars=c('MU'))
.6 <- rstan::stan_trace(brma.1, pars=c('MU'))
f1
```

```
library(Rmisc)
#> Loading required package: lattice
#> Loading required package: plyr
multiplot(f1.1, f1.2, f1.3, f1.4, f1.5, f1.6, cols=2)
```

The mean sensitivity and specificity as estimated by all the fitted distributions are in the table below.

```
<- data.frame(Parameter=c("Sensitivity", "Specificity", "Correlation", "Var(lSens)", "Sigma[1,2]", "Sigma[2,1]","Var(lSpec)"),
brma.summary1 summary(brma.1, pars=c('MU', 'ktau', 'Sigma'))$summary[,c(1, 4, 6, 8:10)])
<- brma.summary1[-c(5,6),]
brma.summary1
names(brma.summary1)[2:5] <- c("Mean", "Lower", "Median", "Upper")
library(loo)
#> This is loo version 2.5.1
#> - Online documentation and vignettes at mc-stan.org/loo
#> - As of v2.0.0 loo defaults to 1 core but we recommend using as many as possible. Use the 'cores' argument or set options(mc.cores = NUM_CORES) for an entire session.
#> - Windows 10 users: loo may be very slow if 'mc.cores' is set in your .Rprofile file (see https://github.com/stan-dev/loo/issues/94).
#>
#> Attaching package: 'loo'
#> The following object is masked from 'package:rstan':
#>
#> loo
<- cbind(Model=rep(c("Gaussian", "C90", "C270", "FGM", "Frank", "BRMA"), each=5),
Table1 rbind(summary(fitgauss.1)$Summary,
summary(fitc90.1)$Summary,
summary(fitc270.1)$Summary,
summary(fitfgm.1)$Summary,
summary(fitfrank.1)$Summary,
brma.summary1),WAIC = t(data.frame(rep(c(summary(fitgauss.1)$WAIC[1],
summary(fitc90.1)$WAIC[1],
summary(fitc270.1)$WAIC[1],
summary(fitfgm.1)$WAIC[1],
summary(fitfrank.1)$WAIC[1],
::waic(extract_log_lik(brma.1, parameter_name="loglik"))$estimates[3,1]), each=5))))
loo#> Warning:
#> 6 (30.0%) p_waic estimates greater than 0.4. We recommend trying loo instead.
rownames(Table1) <- NULL
print(Table1, digits=4)
#> Model Parameter Mean Lower Median Upper n_eff
#> 1 Gaussian Sensitivity 0.756762 6.904e-01 0.756036 8.166e-01 1196.857
#> 2 Gaussian Specificity 0.798289 6.171e-01 0.813517 9.064e-01 704.379
#> 3 Gaussian Correlation -0.820176 -9.861e-01 -0.876343 -3.333e-01 269.179
#> 4 Gaussian Var(Sens) 0.006198 8.321e-06 0.005047 1.947e-02 165.705
#> 5 Gaussian Var(Spec) 0.048111 1.357e-02 0.041060 1.220e-01 169.508
#> 6 C90 Sensitivity 0.751379 6.913e-01 0.753546 8.098e-01 25.638
#> 7 C90 Specificity 0.807051 6.549e-01 0.821119 9.069e-01 119.897
#> 8 C90 Correlation -0.528340 -9.766e-01 -0.725178 -4.020e-18 4.111
#> 9 C90 Var(Sens) 0.004885 3.400e-04 0.003297 1.955e-02 11.615
#> 10 C90 Var(Spec) 0.045694 1.556e-02 0.038049 1.020e-01 137.149
#> 11 C270 Sensitivity 0.757528 6.877e-01 0.761163 8.210e-01 273.236
#> 12 C270 Specificity 0.803502 6.328e-01 0.811740 9.097e-01 1081.987
#> 13 C270 Correlation -0.697493 -9.827e-01 -0.808717 -3.332e-06 40.012
#> 14 C270 Var(Sens) 0.006662 2.667e-04 0.005293 2.027e-02 556.055
#> 15 C270 Var(Spec) 0.044767 1.268e-02 0.037922 1.112e-01 1098.815
#> 16 FGM Sensitivity 0.759407 6.891e-01 0.761931 8.174e-01 2475.208
#> 17 FGM Specificity 0.802588 6.453e-01 0.812498 9.045e-01 2293.332
#> 18 FGM Correlation -0.174538 -2.222e-01 -0.222221 2.222e-01 785.016
#> 19 FGM Var(Sens) 0.005390 7.425e-07 0.004181 1.813e-02 1019.633
#> 20 FGM Var(Spec) 0.041890 1.177e-02 0.036671 9.997e-02 2479.371
#> 21 Frank Sensitivity 0.756683 6.855e-01 0.758340 8.152e-01 2686.631
#> 22 Frank Specificity 0.808239 6.472e-01 0.818777 9.110e-01 1910.561
#> 23 Frank Correlation -0.706819 -8.550e-01 -0.692019 1.000e+00 2700.000
#> 24 Frank Var(Sens) 0.006678 5.896e-04 0.005280 2.140e-02 2699.766
#> 25 Frank Var(Spec) 0.042067 1.201e-02 0.035908 1.039e-01 1937.653
#> 26 BRMA Sensitivity 0.752531 6.323e-01 0.756181 8.415e-01 849.491
#> 27 BRMA Specificity 0.790796 5.273e-01 0.809420 9.539e-01 1044.902
#> 28 BRMA Correlation -0.822353 -9.824e-01 -0.876654 -3.804e-01 239.123
#> 29 BRMA Var(lSens) 0.337556 5.792e-02 0.255387 9.851e-01 1025.609
#> 30 BRMA Var(lSpec) 5.682692 1.472e+00 4.632967 1.690e+01 1070.481
#> Rhat WAIC
#> 1 1.0014 91.28
#> 2 1.0044 91.28
#> 3 1.0154 91.28
#> 4 1.0072 91.28
#> 5 1.0069 91.28
#> 6 1.1047 91.40
#> 7 1.0304 91.40
#> 8 1.3707 91.40
#> 9 1.1005 91.40
#> 10 1.0311 91.40
#> 11 1.0096 90.75
#> 12 1.0001 90.75
#> 13 1.0407 90.75
#> 14 1.0024 90.75
#> 15 0.9999 90.75
#> 16 0.9998 97.37
#> 17 0.9996 97.37
#> 18 1.0070 97.37
#> 19 1.0034 97.37
#> 20 0.9999 97.37
#> 21 0.9994 90.55
#> 22 0.9997 90.55
#> 23 NaN 90.55
#> 24 0.9990 90.55
#> 25 0.9992 90.55
#> 26 0.9999 86.76
#> 27 1.0008 86.76
#> 28 1.0218 86.76
#> 29 1.0023 86.76
#> 30 1.0002 86.76
```

The results are presented graphically as follows

```
<- ggplot2::ggplot(Table1[Table1$Parameter %in% c("Sensitivity", "Specificity"),],
g1 ::aes(x=Model,
ggplot2y= Mean)) +
::geom_point(size=3) +
ggplot2::theme_bw() +
ggplot2::coord_flip() +
ggplot2::facet_grid(~ Parameter, switch="x") +
ggplot2::geom_errorbar(ggplot2::aes(ymin=Lower,
ggplot2ymax=Upper),
linewidth=.75,
width=0.15) +
::theme(axis.text.x = ggplot2::element_text(size=13, colour='black'),
ggplot2axis.text.y = ggplot2::element_text(size=13, colour='black'),
axis.title.x = ggplot2::element_text(size=13, colour='black'),
strip.text = ggplot2::element_text(size = 13, colour='black'),
axis.title.y= ggplot2::element_text(size=13, angle=0, colour='black'),
strip.text.y = ggplot2::element_text(size = 13, colour='black'),
strip.text.x = ggplot2::element_text(size = 13, colour='black'),
plot.background = ggplot2::element_rect(fill = "white", colour='white'),
panel.grid.major = ggplot2::element_blank(),
panel.background = ggplot2::element_blank(),
strip.background = ggplot2::element_blank(),
axis.line.x = ggplot2::element_line(color = 'black'),
axis.line.y = ggplot2::element_line(color = 'black')) +
::scale_y_continuous(name="Mean [95% equal-tailed credible intervals]",
ggplot2limits=c(0.45,1.1),
breaks=c(0.5, 0.75, 1),
labels = c("0.5", "0.75", "1"))
g1
```

`Table1`

above showed that the correlation as estimated by
the BRMA model and the Gaussian copula bivariate beta are more extreme
but comparable to the estimates from the Frank, \(90^{\circ}\)- and \(270^{\circ}\)- Clayton copula. On the other
extreme is the estimate from the model FGM copula bivariate beta and
this is due to the constraints on the association parameter in the FGM
copula where values lie within |2/9|.

In the figure above `g1`

, the marginal mean sensitivity
and specificity from the five bivariate beta distributions are
comparable with subtle differences in the 95 percent credible intervals
despite differences in the correlation structure.

Glas et al. (2003) and Riley et al. (2007a) estimated the Pearson’s
correlation parameter in the BRMA model \(\rho\) as -1 within the frequentist
framework. Using maximum likelihood estimation, Riley et al. (2007b)
showed that the between-study correlation from the BRMA is often
estimated as +/-1. Without estimation difficulties, the table above
shows an estimated Pearson’s correlation of
`-0.8224[-0.9824, -0.3804]`

. This is because Bayesian methods
are not influenced by sample size and therefore able to handle cases of
small sample sizes with less issues.

Essentially, all the six models are equivalent in the first level of
hierarchy and differ in specifying the prior distributions for the
“study-specific” sensitivity and specificity. As thus, the models should
have the same number of parameters in which case it makes sense then to
compare the log predictive densities. Upon inspection, the log
predictive densities from the five copula-based models are practically
equivalent (`min=-38.77, max=-37.89`

) but the effective
number of parameters differed a bit (`min=7.25, max=9.92`

).
The BRMA had `6`

effective parameters but a lower log
predictive density of `-43.4`

. The last column of table above
indicates that the BRMA fits the data best based on the WAIC despite its
lower predictive density.

The `ascus`

dataset has `Test`

as a covariate.
The covariate is used as it is of interest to study its effect on the
joint distribution of sensitivity and specificity (including the
correlation). The following code fits the copula based bivariate
beta-binomial distribution to the data `ascus`

data.

```
.2 <- cdtamodel(copula = "fgm",
fgmmodelargs = list(formula.se = StudyID ~ Test + 0))
.2 <- fit(fgm.2,
fitfgmdata=ascus,
SID="StudyID",
iter=19000,
warmup=1000,
thin=20,
seed=3)
.2 <- cdtamodel(copula = "gauss",
gaussmodelargs = list(formula.se = StudyID ~ Test + 0))
.2 <- fit(gauss.2,
fitgaussdata=ascus,
SID="StudyID",
iter=19000,
warmup=1000,
thin=20,
seed=3)
.2 <- cdtamodel(copula = "c90",
c90modelargs = list(formula.se = StudyID ~ Test + 0))
.2 <- fit(c90.2,
fitc90data=ascus,
SID="StudyID",
iter=19000,
warmup=1000,
thin=20,
seed=3)
.2 <- cdtamodel(copula = "c270",
c270modelargs = list(formula.se = StudyID ~ Test + 0))
.2 <- fit(c270.2,
fitc270data=ascus,
SID="StudyID",
iter=19000,
warmup=1000,
thin=20,
seed=3)
.2 <- cdtamodel(copula = "frank",
frankmodelargs = list(formula.se = StudyID ~ Test + 0))
.2 <- fit(frank.2,
fitfrankdata=ascus,
SID="StudyID",
iter=19000,
warmup=1000,
thin=20,
seed=3)
```

The `BRMA`

is fitted by the code below

```
<- "
BRMA2 data{
int<lower=0> N;
int<lower=0> tp[N];
int<lower=0> dis[N];
int<lower=0> tn[N];
int<lower=0> nondis[N];
matrix<lower=0>[N,2] x;
}
parameters{
real etarho; // g(rho)
matrix[2,2] betamul;
vector<lower=0>[2] sigma;
vector[2] logitp[N];
vector[2] logitphat[N];
}
transformed parameters{
row_vector[2] mul_i[N];
vector[N] p[2];
vector[N] phat[2]; //expit transformation
matrix[2,2] mu;
real rho; //pearsons correlation
real ktau; //Kendall's tau
matrix[2,2] Sigma; //Variance-covariance matrix
matrix[2,2] MU;
row_vector[2] RR;
rho = tanh(etarho); //fisher z back transformation
ktau = (2/pi())*asin(rho);
for (j in 1:2){
for (k in 1:2){
mu[j, k] = inv_logit(betamul[j, k]);
}
}
for (a in 1:N){
mul_i[a] = row(x*betamul,a);
for (b in 1:2){
p[b][a] = inv_logit(logitp[a][b]);
phat[b][a] = inv_logit(logitphat[a][b]);
}
}
MU[1,1] = sum(col(x, 1).*phat[1])/sum(col(x, 1)); //sensitivity HC2
MU[1,2] = sum(col(x, 1).*phat[2])/sum(col(x, 1)); //specificity HC2
MU[2,1] = sum(col(x, 2).*phat[1])/sum(col(x, 2)); //sensitivity RepC
MU[2,2] = sum(col(x, 2).*phat[2])/sum(col(x, 2)); //specificity RepC
RR = row(MU, 2)./row(MU, 1); //RepC vs HC2
Sigma[1, 1] = sigma[1]^2;
Sigma[1, 2] = sigma[1]*sigma[2]*rho;
Sigma[2, 1] = sigma[1]*sigma[2]*rho;
Sigma[2, 2] = sigma[2]^2;
}
model{
#Priors
etarho ~ normal(0, 10);
sigma ~ cauchy(0, 2.5);
for (i in 1:2){
betamul[i] ~ normal(0, 10);
}
for (l in 1:N){
logitp[l] ~ multi_normal(mul_i[l], Sigma);
logitphat[l] ~ multi_normal(mul_i[l], Sigma);
}
//Likelihood
tp ~ binomial(dis,p[1]);
tn ~ binomial(nondis, p[2]);
}
generated quantities{
vector[N*2] loglik;
for (i in 1:N){
loglik[i] = binomial_lpmf(tp[i]| dis[i], p[1][i]);
}
for (i in (N+1):(2*N)){
loglik[i] = binomial_lpmf(tn[i-N]| nondis[i-N], p[2][i-N]);
}
}
"
<- list(
datalist2 tp = ascus$TP,
dis = ascus$TP + ascus$FN,
tn = ascus$TN,
nondis = ascus$TN + ascus$FP,
N = 20,
x = structure(.Data=c(2-as.numeric(as.factor(ascus$Test)),
rev(2-as.numeric(as.factor(ascus$Test)))),
.Dim=c(20, 2)))
.2 <- stan(model_code = BRMA2,
brmadata = datalist2,
chains = 3,
iter = 5000,
warmup = 1000,
thin = 10,
seed = 3,
cores=3)
```

The figure below shows the trace plots for all the six models fitted
to the `ascus`

data where all parameters, including the
correlation parameter(except the BRMA) are modeled as a function of the
covariate. There is proper chains mixing and convergence except for the
case of the Clayton copula based bivariate beta.

```
.1 <- traceplot(fitgauss.2)
f2.2 <- traceplot(fitc90.2)
f2.3 <- traceplot(fitc270.2)
f2.4 <- traceplot(fitfgm.2)
f2.5 <- traceplot(fitfrank.2)
f2
#draws <- as.array(brma.2, pars=c('MU'))
.6 <- rstan::stan_trace(brma.2, pars=c('MU')) f2
```

`multiplot(f2.1, f2.2, f2.3, f2.4, f2.5, f2.6, cols=2)`

The `n_eff`

in table below indicate substantial
autocorrelation in sampling the correlation parameters except in the
`Gaussian`

, `FGM`

and `Frank`

models.
From the copula based bivariate beta distributions, it is apparent that
the correlation between sensitivity and specificity in HC2 and repeat
cytology is different.

```
<- data.frame(Parameter=c(rep(c("Sensitivity", "Specificity"), times=2), "RSE", "RSP", "Var(lSens)",
brma.summary2 "Cov(lSens_Spec)", "Cov(lSens_Spec)", "Var(lSpec)", "Correlation"),
Test=c(rep(c("HC2", "Repc"), times=2), rep("Repc", 2), rep("Both", 5)),
summary(brma.2, pars=c('MU', "RR", 'Sigma', 'ktau'))$summary[,c(1, 4, 6, 8:10)])
names(brma.summary2)[3:6] <- c("Mean", "Lower", "Median", "Upper")
<- cbind(Model=rep(c("Gaussian", "C90", "C270", "FGM", "Frank"), each=14),
Table2 Test=rep(c("HC2", "Repc"), length.out=70),
rbind(summary(fitgauss.2)$Summary,
summary(fitc90.2)$Summary,
summary(fitc270.2)$Summary,
summary(fitfgm.2)$Summary,
summary(fitfrank.2)$Summary),
WAIC = t(data.frame(rep(c(summary(fitgauss.2)$WAIC[1],
summary(fitc90.2)$WAIC[1],
summary(fitc270.2)$WAIC[1],
summary(fitfgm.2)$WAIC[1],
summary(fitfrank.2)$WAIC[1]), each=14))))
$Parameter <- rep(rep(c("Sensitivity", "Specificity", "RSE", "RSP", "Correlation", "Var(Sens)", "Var(Spec)"), each=2), times=5)
Table2
<- rbind(Table2, cbind(Model=rep("BRMA", 11),
Table2
brma.summary2,WAIC=rep(loo::waic(extract_log_lik(brma.2, parameter_name="loglik"))$estimates[3,1], 11)))
#> Warning:
#> 19 (47.5%) p_waic estimates greater than 0.4. We recommend trying loo instead.
rownames(Table2) <- NULL
print(Table2[Table2$Parameter=="Correlation",], digits=4)
#> Model Test Parameter Mean Lower Median Upper n_eff
#> 9 Gaussian HC2 Correlation -0.43812 -0.9984 -6.959e-01 9.847e-01 154.238
#> 10 Gaussian Repc Correlation -0.91991 -0.9997 -9.643e-01 -6.103e-01 75.582
#> 23 C90 HC2 Correlation -0.06588 -0.7610 -1.039e-17 -7.624e-19 39.786
#> 24 C90 Repc Correlation -0.85157 -0.9804 -9.120e-01 -4.906e-01 11.437
#> 37 C270 HC2 Correlation -0.03038 -0.6452 -7.000e-18 -1.782e-18 104.680
#> 38 C270 Repc Correlation -0.77847 -0.9757 -7.058e-01 -5.394e-01 2.995
#> 51 FGM HC2 Correlation -0.07618 -0.2222 -2.215e-01 2.222e-01 2284.371
#> 52 FGM Repc Correlation -0.19819 -0.2222 -2.222e-01 1.894e-01 2550.145
#> 65 Frank HC2 Correlation -0.48806 -0.8140 -4.497e-01 1.000e+00 NaN
#> 66 Frank Repc Correlation -0.73784 -0.8627 -7.275e-01 1.000e+00 NaN
#> 81 BRMA Both Correlation -0.84808 -0.9839 -8.980e-01 -4.497e-01 118.687
#> Rhat WAIC
#> 9 1.0066 236.4
#> 10 1.0542 236.4
#> 23 1.0690 235.7
#> 24 1.1128 235.7
#> 37 1.0326 227.5
#> 38 1.4613 227.5
#> 51 1.0007 245.1
#> 52 0.9997 245.1
#> 65 NaN 238.3
#> 66 NaN 238.3
#> 81 1.0254 233.7
```

The `Clayton90`

model has the lowest WAIC even though
sampling from the posterior distribution was difficult as seen in their
trace plots in figure above and the `n_eff`

and
`Rhat`

in the table above. The difficulty in sampling from
the posterior could be signalling over-parameterisation of the
correlation structure. It would thus be interesting to re-fit the models
using only one correlation parameter and compare the models. WAIC is
known to fail in certain settings and this examples shows that it is
crucial to check the adequacy of the fit and plausibility of the model
and not blindly rely on an information criterion to select the best fit
to the data.

From the posterior relative sensitivity and specificity plotted below, all the models that converged generally agree that repeat cytology was less sensitive than HC2 without significant loss in specificity.

```
<- ggplot2::ggplot(Table2[Table2$Parameter %in% c("RSE", "RSP") & (Table2$Test == "Repc") ,],
g2 ::aes(x=Model, y= Mean, ymax=1.5)) +
ggplot2::geom_point(size=3) +
ggplot2::theme_bw() +
ggplot2::coord_flip() +
ggplot2::facet_grid(~ Parameter, switch="x") +
ggplot2::geom_errorbar(ggplot2::aes(ymin=Lower, ymax=Upper),linewidth=.75, width=0.15) +
ggplot2::geom_hline(ggplot2::aes(yintercept = 1),
ggplot2linetype = "longdash",
colour="blue") +
::theme(axis.text.x = ggplot2::element_text(size=13, colour='black'),
ggplot2axis.text.y = ggplot2::element_text(size=13, colour='black'),
axis.title.x = ggplot2::element_text(size=13, colour='black'),
strip.text = ggplot2::element_text(size = 13, colour='black'),
axis.title.y= ggplot2::element_text(size=13, angle=0, colour='black'),
strip.text.y = ggplot2::element_text(size = 13, colour='black'),
strip.text.x = ggplot2::element_text(size = 13, colour='black'),
panel.grid.major = ggplot2::element_blank(),
panel.background = ggplot2::element_blank(),
strip.background = ggplot2::element_blank(),
plot.background = ggplot2::element_rect(fill = "white", colour='white'),
axis.line.x = ggplot2::element_line(color = 'black'),
axis.line.y = ggplot2::element_line(color = 'black')) +
::scale_y_continuous(name="Mean [95% equal-tailed credible intervals]",
ggplot2limits=c(0.5, 2))
g2
```

Copula-based models offer great flexibility and ease but their use is not without caution. While the copulas used in this paper are attractive as they are mathematically tractable, (Mikosch 2006) and (Genest and Remillard 2006) noted that it might be difficult to estimate copulas from data. Furthermore, the concepts behind copula models is slightly more complex and therefore require statistical expertise to understand and program them as they are not yet available as standard procedure/programs in statistical software.

In this paper, several advanced statistical models for meta-analysis
of diagnostic accuracy studies were briefly discussed. The use of the
`R`

package `CopulaDTA`

within the flexible
`Stan`

interface was demonstrated and shows how complex
models can be implemented in a convenient way.

In most practical situations, the marginal mean structure is of primary interest and the correlation structure is treated a nuisance making the choice of copula less critical. Nonetheless, an appropriate correlation structure is critical in the interpretation of the random variation in the data as well as obtaining valid model-based inference for the mean structure.

When the model for the mean is correct but the true distribution is misspecified, the estimates of the model parameters will be consistent but the standard errors will be incorrect Agresti (2002). Nonetheless, the bivariate beta distribution has the advantage to allow direct joint modelling of sensitivity and specificity, without the need of any transformation, and consequently providing estimates with the appropriate meta-analytic interpretation but with the disadvantage of being more computationally intensive for some of the copula functions.

(Leeflang et al. 2013) showed that the sensitivity and specificity often vary with disease prevalence. The models presented above can easily be extended and implemented to jointly model prevalence, sensitivity and specificity using tri-variate copulas.

There were some differences between the models in estimating the meta-analytic sensitivity and specificity and the correlation. Therefore, further research is necessary to investigate the effect of certain parameters, such as the number of studies, sample sizes and misspecification of the joint distribution on the meta-analytic estimates.

The proposed Bayesian joint model using copulas to construct
bivariate beta distributions, provides estimates with both the
appropriate marginal as well as conditional interpretation, as opposed
to the typical BRMA model which estimates sensitivity and specificity
for specific studies with a particular value for the random-effects.
Furthermore, the models do not have estimation difficulties with small
sample sizes or large between-study variance because: **i**
the between-study variances are not constant but depends on the
underlying means and **ii** Bayesian methods are less
influenced by small samples sizes.

The fitted models generally agree that the mean specificity was slightly lower than what Glas et al. (2003) reported and based on this we conclude that telomerase was not sensitive and specific enough to be recommended for daily use.

In the ASCUS triage data, conclusion based on the fitted models is in line with what the authors conclude: that HC2 was considerably more sensitive but slightly and non-significantly less specific than repeat cytology to triage women with an equivocal Pap smear to diagnose cervical precancer.

While the BRMA had the lowest WAIC for both datasets, we still recommend modelling of sensitivity and specificity using bivariate beta distributions as they easily and directly provide meta-analytic estimates.

Arbyn, M, J Roelens, C Simoens, F Buntinx, E Paraskevaidis, P P L
Martin-Hirsch, and W Prendiville. 2013. “Human Papillomavirus
Testing Versus Repeat Cytology for Triage of Minor Cytological Cervical
Lesions.” *Cochrane Database of Systematic Reviews*,
31–201.

Arends, L R, T H Hamza, J C van Houwelingen, M H Heijenbrok-Kal, M G
Hunink, and T Stijnen. 2008. “Bivariate Random Effects
Meta-Analysis of ROC Curves.” *Medical Decision Making* 28
(5): 621–38.

Carpenter, Bob, Andrew Gelman, Matthew D. Hoffman, Daniel Lee, Ben
Goodrich, Michael Betancourt, Marcus Brubaker, Jiqiang Guo, Peter Li,
and Allen Riddell. 2017. “Stan : A Probabilistic Programming
Language.” *Journal of Statistical Software* 76 (1). https://doi.org/10.18637/jss.v076.i01.

Chu, H, and SR Cole. 2006. “Bivariate Meta-Analysis of Sensitivity
and Specificity with Sparse Data: A Generalized Linear Mixed Model
Approach.” *Journal of Clinical Epidemiology* 59 (12):
1331–32.

Clayton, David G. 1978. “A Model for Association in Bivariate Life
Tables and Its Application in Epidemiological Studies of Familial
Tendency in Chronic Disease Incidence.” *Biometrika* 65
(1): 141–51.

Cong, X, D D Cox, and S B Cantor. 2007. “Bayesian Meta-Analysis of
Papanicolaou Smear Accuracy.” *Gynecologic Oncology* 107
(1 Suppl 1): S133–37.

Farlie, D G J. 1960. “The Performance of Some Correlation
Coefficients for a General Bivariate Distribution.”
*Biometrika* 47: 307–23.

Frank, M J. 1979. “On the Simultaneous Associativity of f(x, y)
and x + y - f(x, y).” *Aequationes Mathematicae*, 194–226.

Genest, C, and B Remillard. 2006. “Comments on t. Mikosh’s Paper,
Copulas: Tales and Facts.” *Extremes* 9: 27–36.

Glas, A S, D Roos, M Deutekom, A H Zwinderman, P M M Bossuyt, and K H
Kurth. 2003. “Tumor Markers in the Diagnosis of Primary Bladder
Cancer. A Systematic Review.” *The Journal of Urology* 169
(6): 1975–82.

Gumbel, E J. 1960. “Bivariate Exponential Distributions.”
*Journal of the American Statistical Association* 55: 698–707.

Hoffman, M D, and A Gelman. 2011. “The No-u-Turn Sampler:
Adaptively Setting Path Lengths in Hamiltonian Monte Carlo.”
*ARXIV*, 11–30.

Leeflang, M M, AW Rutjes, JB Reitsma, L Hooft, and P M Bossuyt. 2013.
“Variation of a Test’s Sensitivity and Specificity with Disease
Drevalence.” *Canadian Medical Association Journal* 185
(11): E537–44.

Libby, D L, and R E Novick. 1982. “Multivariate Generalized Beta
Distributions with Applications to Utility Assessment.”
*Journal of Educational Statistics* 7 (4): 271–94.

Meyer, C. 2013. “The Bivariate Normal Copula.”
*Communications in Statistics - Theory and Methods* 42 (13):
2402–22.

Mikosch, T V. 2006. “Copulas: Tales and Facts. Discussion Paper
with a Rejoinder.” *Extremes* 9: 3–20, 55–62.

Morgenstern, D. 1956. “Einfache Beispiele Zweidimensionaler
Verteilungen.” *Mitteilungsblatt fürMathematische
Statistik* 8: 234–35.

Moses, Lincoln E, David Shapiro, and Benjamin Littenberg. 1993.
“Combining Independent Studies of a Diagnostic Test into a Summary
ROC Curve: Data-Analytic Approaches and Some Additional
Considerations.” *Statistics in Medicine* 12 (14):
1293–1316.

Nelsen, R B. 2006. *An Introduction to Copulas*. Springer-Verlag,
New York.

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.”
*ARXIV*, 1–23.

———. 2016. *: Copula Random Effects Model for Bivariate Meta-Analysis
of Diagnostic Test Accuracy Studies*. https://cran.r-project.org/package=CopulaREMADA.

Nyaga, Victoria N. 2023. *CopulaDTA: Copula Based Bivariate
Beta-Binomial Model for Diagnostic Test Accuracy Studies*. https://cran.r-project.org/package=CopulaDTA.

Olkin, I, and R Liu. 2003. “A Bivariate Beta Distribution.”
*Statistics and Probability Letters* 62: 407–12.

Patton, A J. 2006. “Modelling Asymmetric Exchange Rate
Dependence.” *International Economic Review* 47 (2):
527–56.

Plummer, M. 2008. “Penalized Loss Functions for Bayesian Model
Comparison.” *Biostatistics* 9: 523–39.

Plummer, Martyn et al. 2003. “JAGS: A Program for Analysis of
Bayesian Graphical Models Using Gibbs Sampling.” In
*Proceedings of the 3rd International Workshop on Distributed
Statistical Computing*, 124:125. Technische Universit at Wien Wien,
Austria.

R Core Team. 2022. * : A Language and Environment for Statistical
Computing*. Vienna, Austria: Foundation for Statistical Computing.
https://www.R-project.org/.

Reitsma, J B, A S Glas, A W Rutjes, R J Scholten, P M Bossuyt, and A H
Zwinderman. 2005. “Bivariate Analysis of Sensitivity and
Specificity Produces Informative Summary Measures in Diagnostic
Reviews.” *Journal of Clinical Epidemiology* 58 (10):
982–90.

Riley, R D. 2009. “Multivariate Meta-Analysis: The Effect of
Ignoring Within-Study Correlation.” *Journal of the Royal
Statistical Society* 172 (4).

Riley, R D, K R Abrams, P C Lambert, A J Sutton, and J R Thompson. 2007.
“An Evaluation of Bivariate Random-Effects Meta-Analysis for the
Joint Synthesis of Two Correlated Outcomes.” *Statistics in
Medicine* 26 (1): 78–97.

Riley, R D, K R Abrams, A J Sutton, P C Lambert, and J R Thompson. 2007.
“Bivariate Random-Effects Meta-Analysis and the Estimation of
Between-Study Correlation.” *BMC Medical Research
Methodology* 7 (3).

Rutter, C M, and C M Gatsonis. 2001. “A Hierarchical Regression
Approach to Meta-Analysis of Diagnostic Test Accuracy
Evaluations.” *Statistics in Medicine* 20: 2865–84.

Sarmanov, O. 1966. “Generalized Normal Correlation and
Two-Dimensional Fréchet Classes.” *Soviet Mathematics -
Doklady* 7: 596–99.

Sklar, A. 1959. “Fonctions de Répartition á n Dimensions Et Leurs
Marges.” *Publications de l’Institut de Statistique de
L’Université de Paris* 8: 229–31.

Stan Development Team. 2023. “RStan: The
R Interface to Stan.” https://mc-stan.org/.

Takwoingi, Yemisi, Boliang Guo, Richard D Riley, and Jonathan J Deeks.
2015. “Performance of Methods for Meta-Analysis of Diagnostic Test
Accuracy with Few Studies or Sparse Data.” *Statistical
Methods in Medical Research*, 0962280215592269.

Watanabe, S. 2010. “Asymptotic Equivalence of Bayes Cross
Validation and Widely Applicable Information Criterion in Singular
Learning Theory.” *Journal of Machine Learning Research*
11: 3571–94.