A mixed effect model for bivariate meta-analysis of diagnostic test accuracy studies using a copula representation of the random effects distribution

A mixed effect model for bivariate meta-analysis of diagnostic test accuracy studies using a copula representation of the random effects distribution

Aristidis K. Nikoloulopoulos1
11A.Nikoloulopoulos@uea.ac.uk, School of Computing Sciences, University of East Anglia, Norwich NR4 7TJ, UK
Abstract

Diagnostic test accuracy studies typically report the number of true positives, false positives, true negatives and false negatives. There usually exists a negative association between the number of true positives and true negatives, because studies that adopt less stringent criterion for declaring a test positive invoke higher sensitivities and lower specificities. A generalized linear mixed model (GLMM) is currently recommended to synthesize diagnostic test accuracy studies. We propose a copula mixed model for bivariate meta-analysis of diagnostic test accuracy studies. Our general model includes the GLMM as a special case and can also operate on the original scale of sensitivity and specificity. Summary receiver operating characteristic curves are deduced for the proposed model through quantile regression techniques and different characterizations of the bivariate random effects distribution. Our general methodology is demonstrated with an extensive simulation study and illustrated by re-analysing the data of two published meta-analyses. Our study suggests that there can be an improvement on GLMM in fit to data and makes the argument for moving to copula random effects models. Our modelling framework is implemented in the package CopulaREMADA within the open source statistical environment R.

Keywords: copula models; diagnostic tests; multivariate meta-analysis; random effects models; SROC, sensitivity/specificity.

1 Introduction

Synthesis of diagnostic test accuracy studies is the most common medical application of multivariate meta-analysis [21, 30]. Meta-analysis is broadly defined as the quantitative review of the results of related but independent studies [41]. 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. Accurate diagnosis plays an important role in the disease control and prevention [29].

Diagnostic test accuracy studies observe the result of a gold standard procedure which defines the presence or absence of a decease and the result of a diagnostic test. They typically report the number of true positives (diseased people correctly diagnosed), false positives (non-diseased people incorrectly diagnosed as diseased), true negatives and false negatives. As the sensitivity (proportion of those with the disease) and specificity (proportion of those without the disease) are estimated from different samples in each study (diseased and non-diseased patients), they can be assumed to be independent so that the within-study correlations are set to zero [30]. However, there may be a negative between-studies association which should be accounted for. A negative association between these quantities across studies is likely because studies that adopt less stringent criterion for declaring a test positive invoke higher sensitivities and lower specificities [21].

In situations where studies compare a diagnostic test with its gold standard, heterogeneity arises between studies due to the differences in disease prevalence, study design as well as laboratory and other characteristics [7]. Because of this heterogeneity, a generalized linear mixed model (GLMM) has been recommended in the biostatistics literature [4, 1, 14, 29] to synthesize information. Note in passing that it is equivalent with the hierarchical summary receiver operating characteristic model in Rutter and Gatsonis [46] for the case without covariates [15, 5]. The GLMM assumes independent binomial distributions for the true positives and true negatives, conditional on the latent pair of transformed (via a link function) sensitivity and specificity in each study. The random effects (latent pair of transformed sensitivity and specificity) are jointly analysed with a bivariate normal (BVN) distribution.

Chu et al. [7] propose an alternative mixed model which operates on the original scale of sensitivity and specificity. The random effects follow the bivariate Sarmanov’s [47] family of distributions with beta margins [28]. However, this random effects distribution has a limited range of dependence and is inappropriate for general modelling unless the responses are weakly dependent. Hence, this model is too restrictive in the context of diagnostic accuracy studies where strong (negative) dependence is likely.

We propose a copula mixed model as an extension of the GLMM and mixed model in Chu et al. [7] by rather using a copula representation of the random effects distribution with normal and beta margins, respectively. Copulas are a useful way to model multivariate data as they account for the dependence structure and provide a flexible representation of the multivariate distribution. The theory and application of copulas have become important in finance, insurance and other areas, in order to deal with dependence in the joint tails. Here, we indicate that this can also be important in meta-analysis of diagnostic test accuracy studies. Diagnostic test accuracy studies is a prime area of application for copula models, as the traditional assumption of multivariate normality is invalid in this context.

A copula approach for meta-analysis of diagnostic accuracy studies was recently proposed by Kuss et al. [27] who explored the use of a copula model for observed discrete variables (number of true positives and true negatives) which have beta-binomial margins. This model is actually an approximation of a copula mixed model with beta margins for the latent pair of sensitivity and specificity. Although, this approximation can only be used under the unrealistic case that the number of observations in the respective study group of healthy and diseased probands is the same for each study. In real data applications, the number of true positives and negatives do not have a common support over different studies, hence, one cannot conclude that there is a copula. The natural replicability is in the random effects probability for sensitivity and specificity.

The remainder of the paper proceeds as follows. Section 2 summarizes the standard GLMM for synthesis of diagnostic test accuracy studies. Section 3 has a brief overview of relevant copula theory and then introduces the copula mixed model for diagnostic test accuracy studies and discusses its relationship with existing mixed models. Section 4 discusses suitable parametric families of copulas for the copula mixed model, deduces summary receiver operating characteristic curves for the proposed model through quantile regression techniques and different characterizations of the bivariate random effects distribution, and demonstrates that they can show the effect of different model assumptions. Section 5 contains small-sample efficiency calculations to investigate the effect of misspecifying the random effects distribution on parameter estimators and standard errors and compare the proposed methodology to existing methods. Section 6 summarizes the assessment of the proposed models using the Vuong’s statistic [53], which is based on sample difference in Kullback-Leibler divergence between two models and can be used to differentiate two parametric models which could be non-nested. Section 7 presents applications of our methodology to four data frames with diagnostic accuracy data from binary test outcomes. We conclude with some discussion in Section 8, followed by a section with the software details and a technical Appendix.

2 The standard GLMM

We first introduce the notation used in this paper. The focus is on two-level (within-study and between-studies) cluster data. The data are are , where is an index for the within study measurements and is an index for the individual studies. The data, for study , can be summarized in a table with the number of true positives (), true negatives (), false negatives (), and false positives (); see Table 1.

The standard two-level model of meta-analysing diagnostic test accuracy studies [4, 15, 1, 14, 29] lies in the framework of mixed models [8]. The within-study model assumes that the number of true positives and true negatives are conditionally independent and binomially distributed given , where denotes the bivariate latent (random) pair of transformed sensitivity and specificity. That is

 Yi1|X1=x1 ∼ Binomial(ni1,l−1(x1)); Yi2|X2=x2 ∼ Binomial(ni2,l−1(x2)), (1)

where is a link function such as the commonly used logit. The between studies model assumes that is BVN distributed with mean vector and variance covariance matrix . That is

 X∼BVN(μ,Σ). (2)

The models in (2) and (2) together specify a GLMM with joint likelihood

 L(π1,π2,σ1,σ2,ρ)=N∏i=1∫∫2∏j=1g(yij;nij,l−1(xj))ϕ12(x1,x2;μ,Σ)dx1dx2,

where

 g(y;n,π)=(ny)πy(1−π)n−y,y=0,1,…,n,0<π<1,

is the binomial probability mass function (pmf) and is the BVN density with mean vector and variance covariance matrix . The parameters and are those of actual interest denoting the meta-analytic parameters for the sensitivity and specificity, respectively, while the univariate parameters and are of secondary interest denoting the variability between studies.

3 The copula mixed model for diagnostic test accuracy studies

In this section, we introduce the copula mixed model for diagnostic test accuracy studies and discuss its relationship with existing mixed models. Before that, the first subsection has some background on copula models. In Subsection 3.2 and Subsection 3.3 a copula representation of the random effects distribution with normal and beta margins respectively is presented. We complete this section with details on maximum likelihood estimation.

3.1 Overview and relevant background for copulas

A copula is a multivariate cumulative distribution function (cdf) with uniform margins [22, 33, 25]. If is a bivariate cdf with univariate margins , then Sklar’s [51] theorem implies that there is a copula such that

 F12(x1,x2)=C(F1(x1),F2(x2)). (3)

The copula is unique if are continuous, but not if some of the have discrete components. If is continuous and , then the unique copula is the distribution of leading to

 C(u1,u2)=F12(F−11(u1),F−12(u2)),0≤uj≤1,j=1,2,

where are inverse cdfs. In particular, if is the BVN cdf with correlation and standard normal margins, and is the univariate standard normal cdf, then the BVN copula is

 C(u1,u2)=Φ12(Φ−1(u1),Φ−1(u2);ρ).

The power of copulas for dependence modelling is due to the dependence structure being considered separate from the univariate margins; see e.g., [22, Section 1.6]. If is a parametric family of copulas and is a parametric model for the th univariate margin, then

 C(F1(x1;η1),F2(x2;η2);θ)

is a bivariate parametric model with univariate margins . For copula models, the variables can be continuous or discrete [38].

3.2 The copula mixed model for the latent pair of transformed sensitivity and specificity

Here we generalize the GLMM by proposing a model that links the two random effects using a copula function rather than the BVN distribution.

The within-study model is the same as in the standard GLMM; see (2). The stochastic representation of the between studies model takes the form

 (Φ(X1;l(π1),σ21),Φ(X2;l(π2),σ22))∼C(⋅;θ), (4)

where is a parametric family of copulas with dependence parameter and is the cdf of the N() distribution. The joint density of the transformed latent proportions can be derived as a double partial derivative of the cdf in (3)

 f12(x1,x2;π1,π2,σ1,σ2,θ)=∂C(Φ(x1;l(π1),σ21),Φ(x2;l(π2),σ22);θ)∂x1∂x2 (5) =c(Φ(x1;l(π1),σ21),Φ(x2;l(π2),σ22);θ)ϕ(x1;l(π1),σ21)ϕ(x2;l(π2),σ22),

where and is the copula and N() density, respectively. The models in (2) and (4) together specify a copula mixed model with joint likelihood

 L(π1,π2,σ1,σ2,θ) = = N∏i=1∫10∫102∏j=1g(yij;nij,l−1(Φ−1(uj;l(πj),σ2j)))c(u1,u2;θ)du1du2.

It is important to note that the copula parameter is a parameter of the random effects model and it is separated from the univariate parameters. The univariate parameters and are those of actual interest denoting the meta-analytic parameters for the sensitivity and specificity, while the univariate parameters and are of secondary interest expressing the variability between studies.

3.2.1 Relationship with the GLMM

In this subsection, we show what happens when the bivariate copula is the BVN copula. The resulting model is the same as the GLMM.

The BVN copula density is

 c(u1,u2;ρ)=1√1−ρ2exp(z21+z22−2ρz1z22√1−ρ2)exp(z21+z222),

where . Then for we have . Hence, the joint density in (5) becomes

 f12(x1,x2;π1,π2,σ1,σ2,ρ) = 12πσ1σ2√1−ρ2exp[12√1−ρ2{(x1−l(π1))22σ21+ (x2−l(π2))22σ22−2ρ(x1−l(π1))(x2−l(π2))σ1σ2}],

which apparently is the BVN density .

3.3 The copula mixed model for the latent pair of sensitivity and specificity

The within-study model also assumes that the number of true positives and true negatives are conditionally independent and binomially distributed given , where denotes the bivariate latent random pair of sensitivity and specificity. That is

 Yi1|X1=x1 ∼ Binomial(ni1,x1); Yi2|X2=x2 ∼ Binomial(ni2,x2). (7)

So one does not have to transform the latent sensitivity and specificity and can work on the original scale. The Beta() distribution can be used for the marginal modeling of the latent proportions and its density is

 f(x;α,β)=xα−1(1−x)β−1B(α,β),00.

In the sequel we will use the Beta() parametrization, where (mean parameter) and (dispersion parameter).

The stochastic representation of the between studies model is

 (F(X1;π1,γ1),F(X2;π2,γ2))∼C(⋅;θ), (8)

where is a parametric family of copulas with dependence parameter and is the cdf of the the Beta() distribution. The models in (3.3) and (8) together specify a copula mixed model with joint likelihood

 L(π1,π2,γ1,γ2,θ) = N∏i=1∫10∫102∏j=1g(yij;nij,xj)c(F(x1;π1,γ1),F(x2;π2,γ2);θ) ×2∏j=1f(xj;πj,γj)dx1dx2 = N∏i=1∫10∫102∏j=1g(yij;nij,F−1(uj;πj,γj))c(u1,u2;θ)du1du2.

As before, the copula parameter is a parameter of the random effects model and it is separated from the univariate parameters, the univariate parameters and are the meta-analytic parameters for the sensitivity and specificity, but, now and express the variability between studies.

3.3.1 Relationship with existing models

Chu et al. [7], instead of using a copula for the random effects distribution or a copula density for in (3.3), use the Sarmanov’s [47] family of bivariate densities

 f12(x1,x2)=f1(x1)f2(x2)(1+θψ1(x1)ψ2(x2)),

where is the marginal density of , is a bounded non-constant function such as , and for all . For the Sarmanov’s densities if one uses , then the Farlie–Gumbel–Morgenstern copula (density) is obtained. However in [7], “kernels” of the type are considered as in [28]. The advantage of this choice is that the corresponding likelihood function has a closed form, since the product of integrals can be evaluated analytically. The joint likelihood takes the form

 L(π1,π2,γ1,γ2,θ) = N∏i=1∫10∫102∏j=1g(yij;nij,xj)f(xj;πj,γj)(1+θ2∏j=1(xj−πj))dx1dx2 = N∏i=12∏j=1h(yij;nij,πj,γj)(1+θ2∏j=1yij−nijπjγ−1j+nij−1),

where

is the pmf of a Beta-Binomial() distribution with mean and variance . The disadvantage of this mixed model is that the Sarmanov’s density with beta margins in [28] has a limited range of dependence and is inappropriate for general modeling unless the responses are weakly dependent.

Kuss et al. [27] proposed a copula model with beta-binomial margins in this context. This model is actually an approximation of the copula mixed model with beta margins for the latent pair of sensitivity and specificity in (3.3) and (8). They attempt to approximate the likelihood in (3.3) with the likelihood of a copula model for observed discrete variables which have beta-binomial margins.

The approximation that they suggest is

 L(π1,π2,γ1,γ2,θ)≈N∏i=1c(H(yi1;ni1,π1,γ1),H(yi2;ni2,π2,γ2);θ)2∏j=1h(yij;nij,πj,γj),

where is the cdf of the the Beta-Binomial() distribution. In their approximation the authors also treat the observed variables which have beta-binomial distributions as being continuous, and model them under the theory for copula models with continuous margins. Kuss et al. [27], referring to Genest and Nešlehová [9], claim that there are problems on applying copula to discrete data especially in extreme cases with very small numbers of support points for the discrete marginal distributions. Genest and Nešlehová [9] only warn against estimation for discrete-margined copula models using rank-based methods, instead recommending maximum likelihood estimation. Essentially, Genest et al. [10] apply copula models to multivariate binary data (the extreme case of discreteness) and call on composite likelihood techniques for estimation. Multivariate copulas for discrete response data have been in use for a considerable length of time, e.g., in Joe [22], and earlier for some simple copula models. Several examples of copula models for multivariate discrete data can be found in the literature; see e.g., [40] for an application in biostatistics and [34] for a survey of copula models and methods for multivariate discrete response data.

However, the main problem in [27] is that the approximation (even if treating the observed variables which have beta-binomial distributions as being discrete) can only be used under the unrealistic case that the number of observations in the respective study group of healthy and diseased probands is the same for each study . In real data applications, the discrete do not have a common support over different studies or , hence, one cannot conclude that there is a copula for that applies when the vary with different studies . The natural replicability is in the random effects probability for sensitivity and specificity.

3.4 Maximum likelihood estimation and computational details

Estimation of the model parameters and can be approached by the standard maximum likelihood (ML) method, by maximizing the logarithm of the joint likelihood in (3.2) and (3.3), respectively. The estimated parameters can be obtained by using a quasi-Newton [32] 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. Assuming that the usual regularity conditions [49] for asymptotic maximum likelihood theory hold for the bivariate model as well as for its margins we have that ML estimates are asymptotically normal. Therefore one can build Wald tests to statistically judge any effect.

For mixed models of the form with joint likelihood as in (3.2) and (3.3), numerical evaluation of the joint pmf is easily done with the following steps:

1. Calculate Gauss-Legendre quadrature points and weights in terms of standard uniform; see e.g., [52].

2. Convert from independent uniform random variables and to dependent uniform random variables and that have distribution . The inverse of the conditional distribution corresponding to the copula is used to achieve this.

3. Numerically evaluate the joint pmf, e.g.,

 ∫10∫102∏j=1g(yj;nj,F−1(uj;πj,γj))c(u1,u2;θ)du1du2

in a double sum:

 nq∑q1=1nq∑q2=1wq1wq2g(y1;n,F−1(uq1;πj,γj))g(y2;n,F−1(C−1(uq2|uq1;θ);πj,γj)).

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 [32]. Our comparisons show that is adequate with good precision to at least at four decimal places; hence it also provides the advantage of fast implementation.

To sum up, our mixed effect model for meta-analysis of diagnostic test accuracy studies using a copula representation of the random effects distribution with a double integral is straightforward computationally. Note in passing that the linear mixed model in [43] can also provide handy computations, but it has limitations due to the use of continuity correction and normal approximation [4, 29].

4 Choices of parametric families of copulas

In our candidate set, families that have different strengths of tail behaviour (see e.g., [17]) are included. In the descriptions below, a bivariate copula is reflection symmetric if its density satisfies for all . Otherwise, it is reflection asymmetric often with more probability in the joint upper tail or joint lower tail. Upper tail dependence means that as and lower tail dependence means that as . If for a bivariate copula , then , where is the survival (or rotated by 180 degrees) copula of ; this “reflection” of each uniform random variable about changes the direction of tail asymmetry.

• Reflection symmetric copulas with tail independence satisfying and as , such as the Frank copula with inverse conditional cdf

 C−1(v|u;θ)=−1θlog[1−(1−e−θ)(v−1−1)e−θu+1],θ∈(−∞,∞)∖{0}.
• Reflection symmetric copulas with intermediate tail dependence [18] such as the BVN copula, which satisfies as with inverse conditional cdf

 C−1(v|u;θ)=Φ(√1−ρ2Φ−1(v)+ρΦ−1(u)),θ∈[−1,1].
• Reflection asymmetric copulas with lower tail dependence only such as the Clayton copula with inverse conditional cdf

 C−1(v|u;θ)={(v−θ/(1+θ)−1)u−θ+1}−1/θ,θ∈(0,∞).
• Reflection asymmetric copulas with upper tail dependence only such as the rotated by 180 degrees Clayton copula with inverse conditional cdf

 C−1(v|u;θ)=1−[{(1−v)−θ/(1+θ)−1}(1−u)−θ+1]−1/θ,θ∈(0,∞).

The Frank and BVN copulas interpolate from the Fréchet lower (perfect negative dependence) to the Fréchet upper (perfect positive dependence) bound, and, thus they are sufficient from bivariate studies on diagnostic accuracy where negative dependence between the number of true positives and true negatives is expected. The Clayton copula belongs in the Archimedean class of copulas. Archimedean copulas, see e.g. [22], have the form,

 C(u1,u2;θ)=ϕ(ϕ−1(u1;θ)+ϕ−1(u2;θ);θ), (10)

where the generator is the Laplace transform (LT) of a univariate family of distributions of positive random variables indexed by the parameter , such that and its inverse have closed forms. The Clayton copula interpolates from the independence () to the Fréchet upper (comonotonic copula) bound (). For extension of the Laplace transform for , the Clayton family extends to countermonotonicity (). However this extension is not generally useful for applications because the support of (10) is not all of [22, page 109]. Negative dependence in Clayton copulas can be introduced by applying decreasing transformations to the “oppositely” ordered variables. If where is a copula with positive dependence, one could always get some negative dependence, by supposing is the copula of (rotation by 90 degrees) or the copula of (rotation by 270 degrees). So it is worthwhile to rotate the Clayton copula by 90 and 270 degrees to model negative dependence. These rotated copulas interpolate from the Fréchet lower (perfect negative dependence) () to the independence (). Negative upper-lower tail dependence means that as and negative lower-upper tail dependence means that as [24]. So in order to model negative (tail) dependence the choices are:

• Reflection asymmetric copula family with negative upper-lower tail dependence, such as the rotated by 90 degrees Clayton copula with inverse conditional cdf

 C−1(v|u;θ)={(vθ/(1−θ)−1)(1−u)θ+1}1/θ,θ∈(0,∞).
• Reflection asymmetric copula family with negative lower-upper tail dependence, such as the as the rotated by 270 degrees Clayton copula with inverse conditional cdf

 C−1(v|u;θ)=1−[{(1−v)θ/(1−θ)−1}uθ+1]1/θ,θ∈(0,∞).

For this paper, the above copula families are sufficient for the applications in Section 7, since tail dependence is a property to consider when choosing amongst different families of copulas and the concept of upper/lower tail dependence is one way to differentiate families. Nikoloulopoulos and Karlis [39] have shown that it is hard to choose a copula with similar properties from real data, since copulas with similar (tail) dependence properties provide similar fit. Kuss et al. [27] used, in addition to these copulas, the Placket copula. Plackett copula is a reflection symmetric copula [22, pages 221-22] with tail independence [33, page 215] (not reflection asymmetric copula as stated in [27]) and is not used here since we have included another choice of copulas with similar properties i.e., the Frank copula.

4.1 Summary receiver operating characteristic curves

Rutter and Gatsonis [46] proposed a hierarchical summary receiver operating characteristic (SROC) curve which for some cases is the same with the corresponding GLMM SROC curve [5]. For the GLMM model, the model parameters control the shape of the SROC curve. The GLMM SROC curve can be obtained through a characterization of the estimated bivariate normal distribution by a line [5, 6, 7]. Based on the bivariate normality of the random effects, the expected sensitivity for a chosen specificity in the transformed scale is given in a closed form:

 E[X1|X2=x2]=[l(π1)−ρl(π2)σ1/σ2]+ρl(x2)σ1/σ2. (11)

In general, however, is not in closed form and thus does not have simple expressions in terms of distribution functions and copulas.

An alternative to the mean for specifying “typical” values of for each value of is the median, which leads to the notion of median regression of on . For in range of , let denote a solution to the equation . Then the scatter plot of and is the median regression curve of on .

For copula models, median regression curves [33, pages 217–218] can be easily calculated, since

 Pr(X1≤x1|X2=x2)=Pr(U1≤F1(x1)|U2=F2(x2))=C(u1|u2)∥∥∥u1=F1(x1)u2=F2(x2),

but their shape also depends on the choice of bivariate copulas. Furthermore, as emphasized in [45], since there is no unique definition of a SROC curve, it is preferable and will make more sense to deduce confidence regions as well. To this end in addition of using just median regression curves we will also exploit the use of quantile regression curves with a focus on high () and low quantiles () which are strongly associated with the upper and lower tail dependence imposed from each parametric family of copulas. These can be also seen as confidence regions of the median regression SROC curve. Note that Kendall’s tau only accounts for the dependence dominated by the middle of the data, and it is expected to be similar amongst different families of copulas. However, the tail dependence varies, as explained in Section 4, and is a property to differentiate amongst different families of copulas.

To find the quantile regression curves:

1. Set .

2. Solve for the quantile regression curve .

3. For replace by for beta margins or for normal margins.

4. Plot versus .

Of course, the quantile regression curve of on is defined similarly. However, there is no priori reason to regress on instead of the other way around [1]. In fact, if one wants to reserve the nature of a bivariate response instead of a univariate response along with a covariate, then a contour graph can be easily plotted. The contour plot can be seen as the predictive region (analogously to [43]) of the estimated pair of sensitivity and specificity. However, the resulted shape of the prediction region is not depended on the assumption of bivariate normality for the random effects.

To depict the different shapes of the SROC curves, in Figures 1 and 2 we plot them from the copula representation of the random effects distribution with normal and beta margins, respectively, and BVN, Frank Clayton by 90 and 270 copulas with the same model parameters and , respectively. We convert from to the BVN, Frank and rotated Clayton copula parameter via the relations

 τ=2πarcsin(θ), (12)
 τ=⎧⎪⎨⎪⎩1−4θ−1−4θ−2∫0θtet−1dt,θ<01−4θ−1+4θ−2∫θ0tet−1dt,θ>0, (13)
 τ={θ/(θ+2),by 0 or 180 degrees−θ/(θ+2),by 90 or 270 degrees (14)

in [19], [11], and [12] respectively.

5 Small-sample efficiency – Misspecification of the random effects distribution

An extensive simulation study is conducted (a) to gauge the small-sample efficiency of the ML and approximation in Kuss et al. [27]’s (hereafter KHS) methods, and (b) to investigate in detail the misspecification of the parametric margin or family of copulas of the random effects distribution.

To simulate the data we have used the generation process in [42] to get heterogeneous study sizes; the simulation steps follow:

1. Simulate the study size from a shifted gamma distribution, i.e., and round off to the nearest integer.

2. Simulate from a parametric family of copulas ; is converted to BVN, Frank and Clayton rotated by 90/270 dependence parameter via the relations in (12), (13), and (14).

3. Convert to beta realizations via or normal realizations via for ; for the latter convert to proportions via .

4. Draw the number of diseased from a distribution.

5. Set , and then round for .

We randomly generated samples of size from the Clayton rotated by 270 degrees copula mixed model with beta margins. Table 2 contains the resultant biases, root mean square errors (RMSE), and standard deviations (SD), along with average theoretical variances scaled by , for the MLEs under different copula choices and margins. The theoretical variances of the MLEs are obtained via the gradients and the Hessian computed numerically during the maximization process. We also provide biases, RMSEs and SDs for the KHS estimates under the ‘true’ model, i.e., the Clayton rotated by 270 degrees copula mixed model with beta margins.

Conclusions from the values in the table are the following:

• ML with the the ‘true’ copula mixed model is highly efficient according to the simulated biases and standard deviations.

• The MLEs of the meta-analytic parameters are slightly underestimated under copula misspecification. That is, there is some downward bias for these parameters, especially if the “working” model is not close to Kullback-Liebler distance with the “true” model, i.e., it is misspecified. For example in the table there is more bias for the Clayton rotated by 90 degrees and Frank copulas since they have different tail dependence from the ‘true’ model, i.e., the rotated Clayton by 270 degrees. An interesting result is that the BVN copula performed rather well under misspecification.

• The SDs are rather robust to the copula misspecification.

• The meta-analytic MLEs and SDs are not robust to the margin misspecification, while the MLE of and its SD is.

• The KHS approximation method yields to biased univariate parameter estimates.

• The efficiency of the KHS approximation method is low for the dependence parameter . The parameter is substantially underestimated.

The simulation results indicate that the KHS approximation method in [27] is an inefficient; hence flawed method. This was expected, since theoretically there are serious problems on modelling assumptions under the case of heterogeneous study sizes. If the number of true positives and negatives do not have a common support over different studies, then one cannot conclude that there is a copula. To make our study complete, we perform theoretical calculations, similarly to [23, 35, 37], to investigate the accuracy of the approximate copula likelihood method in [27] for the special case of a constant size of groups of diseased and healthy people in the single studies and show whether or not this leads to consistent estimate of the parameters of the bivariate random effects distribution. As shown in the Appendix, the KHS method leads to asymptotic bias for both the univariate and copula parameters and hence there is no consistency. Also given the resultant substantial asymptotic downward bias for the dependence parameter, the approximation deteriorates, and, hence cannot be used e.g., for prediction purposes via SROC curves. To this end, the KHS approximation method is not used in the sequel, since its inefficiency has been shown, and, it should be avoided for bivariate meta-analysis of diagnostic test accuracy studies.

The effect of misspecifying the copula choice can be seen as minimal for both the univariate parameters and Kendall’s tau. However, note that (a) the meta-analytic parameters are a univariate inference, and hence it is the univariate marginal distribution that matters and not the type of the copula, and, (b) as previously emphasized Kendall’s tau only accounts for the dependence dominated by the middle of the data (sensitivities and specificities), and it is expected to be similar amongst different families of copulas. However, the tail dependence varies, as explained in Section 4, and is a property to consider when choosing amongst different families of copulas, and, hence affects the shape of SROC curves, i.e., prediction. SROC will essentially show the effect of different model (random effect distribution) assumptions, since it is an inference that depends on the joint distribution.

6 Vuong’s test for model comparison

In this section we provide a methodology for the comparison of non-nested models. It would be used as a tool to show if the copula mixed model provides better fit than the standard GLMM. We will call a test proposed by Vuong [53]. The Vuong’s test is the sample version of the difference in Kullback-Leibler divergence between two models and can be used to differentiate two parametric models which could be non-nested.This test has been used extensively in the copula literature to compare copula models, see e.g., [2, 3, 25]

Assume that we have Models 1 and 2 with parametric densities and respectively. We can compare

 Δ1f✠=N−1[∑i{Ef✠[logf✠(Y1,Y2)]−Ef✠[logf(1)(Y1,Y2;θ(1))]}],

and

 Δ2f✠=N−1[∑i{Ef✠[logf✠(Y1,Y2)]−Ef✠[logf(2)(Y1,Y2;θ(2))]}],

where are the parameters in Models 1 and 2 respectively that lead to the closest Kullback-Leibler divergence to the true ; equivalently they are the limits in probability of the MLEs based on models 1 and 2 respectively. Model 1 is closer to the true , i.e., is the better fitting model if , and Model 2 is the better fitting model if . The sample version of with MLEs is

 ¯D=N∑i=1Di/N,

where . Vuong [53] has shown that asymptotically

 √N¯D/s∼N(0,1),

where .

7 Illustrations

We illustrate the use of the copula mixed model by re-analysing the data of two published meta-analysis [48, 13]. These data have been frequently used as an example for methodological papers on meta-analysis of diagnostic accuracy studies [43, 4, 46, 44, 16, 27].

We fit the copula mixed model for all different choices of parametric families of copulas and margins. To make it easier to compare strengths of dependence, we convert the copula parameters to Kendall’s ’s via the relations in (12), (13), and (14) for BVN, Frank and rotated Clayton copulas, respectively. Since 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. We further compute the Vuong’s tests with Model 1 being the BVN copula mixed model with normal margins, i.e., the standard GLMM, to reveal if any other copula mixed model provides better fit than the standard GLMM.

Finally, we demonstrate SROC curves and summary operating points (a pair of average sensitivity and specificity) with a confidence region and a predictive region as deduced in Section 4.1.

7.1 The telomerase and computed tomography data

In Glas et al. [13] the telomerase marker for the diagnosis of bladder cancer is evaluated using 10 studies. The size in each study ranges from 35 to 195. The interest was to define if this non-invasive and cheap marker could replace the standard of cystoscopy or histopathology. Riley et al. [44] applied the GLMM with different starting values and all produced a between-study correlation estimate of but with different meta-analytic parameter point estimates and standard errors. Clearly at this example, it is not possible to estimate the correlation between the logit sensitivity and specificity, and the maximum likelihood estimator should truncate the correlation to the left boundary of its parameter space, i.e., . In [27] it is acknowledged that the copula model for observed discrete variables which have beta-binomial margins yields sensible estimates for the dependence parameter and its standard error. This result is in error and due to the fact that the KHS method underestimates the dependence parameter as emphasized in Section 5 and shown in the Appendix for .

Fitting the copula mixed model for all different choices of parametric families of copulas and margins, the resultant estimate of the dependence parameter was close to the boundary of its parameter space. If the dependence parameter estimate is so large (on absolute value), the copula should be set to countermonotonic (Fréchet lower bound), and, then optimize over the remaining (univariate) parameters. With other words there exists negative perfect dependence, and thus there is only one copula: the countermonotonic copula. This is a limiting case for all the parametric families of copulas, listed in Section 4, when the dependence parameter is fixed to the left boundary of its parameter space.

This was also the case for the analysis of the data on 17 studies of computed tomography (CT) for the diagnosis of lymph node metastasis in women with cervical cancer, one of three imaging techniques in the meta-analysis in [48]. The size in each study ranges from 20 to 253. Diagnosis of metastatic disease by CT relies on nodal enlargement.