Robust bivariate random-effects model for accommodating outlying and influential studies in meta-analysis of diagnostic test accuracy studies

Abstract

Due to the inevitable inter-study correlation between test sensitivity (Se) and test specificity (Sp), mostly because of threshold variability, hierarchical or bivariate random-effects models are widely used to perform a meta-analysis of diagnostic test accuracy studies. Conventionally, these models assume that the random-effects follow the bivariate normal distribution. However, the inference made using the well-established bivariate random-effects models, when outlying and influential studies are present, may lead to misleading conclusions, since outlying or influential studies can extremely influence parameter estimates due to their disproportional weight. Therefore, we developed a new robust bivariate random-effects model that accommodates outlying and influential observations and gives robust statistical inference by down-weighting the effect of outlying and influential studies. The marginal model and the Monte Carlo expectation-maximization algorithm for our proposed model have been derived. A simulation study has been carried out to validate the proposed method and compare it against the standard methods. Regardless of the parameters varied in our simulations, the proposed model produced robust point estimates of Se and Sp compared to the standard models. Moreover, our proposed model resulted in precise estimates as it yielded the narrowest confidence intervals. The proposed model also generated a similar point and interval estimates of Se and Sp as the standard models when there are no outlying and influential studies. Two published meta-analyses have also been used to illustrate the methods.

Keywords

Influential studies meta-analysis outlying studies robust random-effects model sensitivity specificity

1 Introduction

A meta-analysis of diagnostic test accuracy (DTA) allows for the synthesis of diagnostic test characteristics, such as test Se and Sp, across several independent studies. Due to the inevitable correlation between Se and Sp, hierarchical random-effects model¹ or bivariate random-effects models^2,3 are widely used to perform meta-analysis of DTA studies. Conventionally, these models assume that the random-effects, which describe the variation among the true study-specific sensitivities and specificities, follow the bivariate normal distribution. However, in the presence of a few outlying and/or influential studies, this assumption may not adequately describe the variation in true Se and Sp across studies. Moreover, due to their weights, outlying and/or influential studies have the potential to extremely influence the estimated overall Se and Sp. Therefore, simply fitting the current standard bivariate random-effects models when outlying and/or influential studies are present may lead to misleading inferences.⁴

In the conventional meta-analysis of intervention studies, a substantial number of research articles^5–10 which advocate for the use of robust statistical methods in the presence of outlying or influential studies have been published. Methods of identifying outlying and/or influential studies had been proposed independently by Viechtbauer and Cheung⁷ and Gumedze and Jackson.⁸ Lee and Thompson⁶ proposed several heavy-tailed distributions including the t, skew-normal and skew-t for the random-effects in the presence of outlying studies. On the other hand, long-tailed distributions like the t and the Beta have been proposed to accommodate outlying studies by Baker and Jackson.⁵ A model that uses a mixture of normal distributions for the random-effects, identifies and down-weights the effect of outlying studies, was developed by Beath.⁹ As an improvement on their previous work, Baker and Jackson¹⁰ proposed new and mathematically tractable skewed marginal distributions which are also simpler to fit than the models proposed by Baker and Jackson.⁵ In their recent review article, Jackson and White¹¹ recommended the use of robust meta-analytic models that relaxes the assumption of normally distributed random-effects. Most recently, Negeri and Beyene¹² proposed and evaluated a skew-normal random-effects meta-analytic model that accounts for possible skewness in DTA data.

Despite the attention it has gained in the conventional meta-analysis of interventions, tackling the effect of outlying and/or influential studies was not studied in the context of a meta-analysis of DTA studies. To the best of our knowledge, the proposal of advanced statistical methods for detecting the presence of outlying or influential studies in a meta-analysis of DTA studies was only recently studied by Negeri and Beyene.⁴ In their article, it has been demonstrated that the well-established bivariate random-effects models are significantly affected by outlying or influential studies and may lead to misleading inferences. Although having a sound statistical method for identifying outlying and/or influential studies is a promising and necessary start, there is no clear methodology that guides how to deal with outlying and influential studies in a meta-analysis of DTA studies. There has been a general understanding among researchers that outlying or influential observations should not be discarded from the analyses due to their nature. However, a more robust statistical method that can accommodate those outlying and/or influential studies by keeping the effect of such studies to a minimal is desired.

Since the current hierarchical and bivariate random-effects models for meta-analysis of DTA studies do not account for outlying and/or influential studies and may lead to unreliable conclusions as pointed out by,⁴ we aim to fill this gap by developing a new alternative bivariate random-effects model which can accommodate outlying and/or influential studies by down-weighting their effect and gives a robust statistical inference. More specifically, we propose the bivariate Laplace distribution¹³ for the random-effects to model the variation in Se and Sp across studies with the aim of accounting for potentially outlying and/or influential studies. For two reasons we chose the bivariate Laplace distribution over the frequently used distributions like the t. First, the Laplace distribution has not been implemented in the conventional meta-analysis of treatments or in a meta-analysis of DTA studies, although there have been other proposed robust models that accommodate outlying studies in the former literature as discussed above. Second, since the bivariate Laplace distribution has exactly the same number of parameters as the bivariate normal, it leads to a meaningful model comparison using information criteria approaches. By assuming the bivariate normal distribution to model the within-study variation in Se and Sp, we derived our marginal model to be the bivariate normal-Laplace (NL) distribution. That led us to implement the Monte Carlo expectation-maximization (MCEM) algorithm to obtain the maximum likelihood estimates of our model parameters. The standard errors of the parameter estimators and information criteria for model selection are obtained at the last iteration of the MCEM algorithm once the algorithm has converged.

The rest of this article is organized as follows. Section 2 motivates the proposed model using two published datasets. Section 3 gives a comprehensive description of the standard as well as the proposed methods. Section 4 presents the simulation study results. In Section 5, we illustrate the methods using two real DTA meta-analyses, and, we close the article by a discussion and conclusion in Section 6.

2 Motivating examples

This section discusses two published meta-analyses of DTA studies to motivate the proposed model. These data are selected because they contain one or more outlying and/or influential studies according to, among others, the methods proposed by Negeri and Beyene.⁴ The performance of the different methods when fitted to these datasets will be discussed in Section 5.

2.1 Naked eye examination data

The first example, which is based on the naked eye examination (NEE) test for the diagnosis of primary melanoma compared to the dermoscopy test, is obtained from Vestergaard et al.¹⁴ The forest plot for the NEE data is given in Figure 1. The authors¹⁴ have identified and discussed the impact of two outlying studies on the diagnostic odds ratio of dermoscopy relative to the NEE test. Both studies 2 and 6 were detected as outlying and influential in specificity, and study 3 is identified to be outlying and influential in sensitivity.

Figure 1.

Forest plot for sensitivity (left) and specificity (right) of the NEE data.

2.2 Mini-mental state examination data

The second meta-analysis is from the review of Arevalo et al.¹⁵ about the mini-mental state examination (MMSE) test for detecting Alzheimer’s disease and dementia in people with mild cognitive impairment. Figure 2 displays the forest plot of the MMSE data. Studies 1 and 2 are detected as outlying and influential in specificity.

Figure 2.

Forest plot for sensitivity (left) and specificity (right) of the MMSE data.

3 Methods

In this section, we present the three bivariate random-effects models used to synthesize a pair of Se and Sp obtained from a set of k independent DTA studies. Before discussing the models, we define the following notations. Let $y_{i}$ = ${(y_{1 i}, y_{2 i})}^{T}$ = {logit $({\hat{S e}}_{i})$ , logit( ${\hat{S p}}_{i}$ ) $}'$ be a vector of the observed logit sensitivities and logit specificities for each study, i = 1,…,k. Let $μ$ and $Σ$ denote the location and scale parameters of a bivariate distribution, respectively.

3.1 The bivariate normal–normal model

The standard bivariate random-effects model of Reitsma et al.² assumes that the within-study variation in sensitivity and specificity is modelled according to

Y_{i} | b_{i} \sim N_{2} (μ + b_{i}, Ψ_{i}), i = 1, \dots, k

(1)

where

μ

is the fixed-effect and represents the pooled or overall logit-transformed Se and Sp,

b_{i}

denotes the study-specific true underlying effect (random-effects), and

Ψ_{i} = (\begin{matrix} ψ_{1 i}^{2} & 0 \\ 0 & ψ_{2 i}^{2} \end{matrix}), i = 1, \dots, k

represents the assumed to be known within-study covariance matrix.

On the second level, the bivariate normal–normal (BNN) model assumes that the true effect, $b_{i}$ , varies across studies, and models this variability as

b_{i} \sim N_{2} (0, Σ), i = 1, \dots, k .

(2)

where

Σ = (\begin{matrix} σ_{1}^{2} & σ_{12} \\ σ_{12} & σ_{2}^{2} \end{matrix})

denotes the between-study covariance matrix.

As in any hierarchical model, maximum likelihood estimation of parameters of the BNN model is also based on the marginal model. The marginal model obtained from equations (1) and (2) by integrating out the random-effects is given as

Y_{i} \sim N_{2} (μ, Σ_{i}), i = 1, \dots, k; where Σ_{i} = Σ + Ψ_{i}

(3)

The likelihood and log-likelihood function of the BNN model which follows from equation (3) is given as

L (μ, Σ | y) = \prod_{i = 1}^{k} \frac{1}{2 π | Σ_{i} |^{\frac{1}{2}}} \exp {- \frac{1}{2} {(y_{i} - μ)}^{T} Σ_{i}^{- 1} (y_{i} - μ)} and l (μ, Σ | y) = - k log (2 π) - \frac{1}{2} \sum_{i = 1}^{k} log | Σ_{i} | - \frac{1}{2} \sum_{i = 1}^{k} {(y_{i} - μ)}^{T} Σ_{i}^{- 1} (y_{i} - μ),

(4)

respectively.

From equation (4), it follows that the maximum likelihood estimates (MLE) of $μ$ and $Σ$ do not have closed form solutions and should be obtained iteratively using numerical methods. R packages such as the mada¹⁶ and the mvmeta¹⁷ can be used to fit the BNN model.

Once the parameters are estimated, inference regarding the parameter of interest, $μ$ , is made by first deriving its variance (or standard error). From equation (4), the solution to the score function with respect to $μ$ is given as

μ = {(\sum_{i = 1}^{k} W_{i})}^{- 1} \sum_{i = 1}^{k} W_{i} Y_{i} where W_{i} = {(Σ + Ψ_{i})}^{- 1} = Σ_{i}^{- 1}

(5)

Therefore, the estimated variance-covariance matrix of the pooled logit(Se) and logit(Sp) estimator, $μ$ , can easily be obtained from equation (5) as

\hat{var (μ)} = {(\sum_{i = 1}^{k} {\hat{W}}_{i})}^{- 1} = {(\sum_{i = 1}^{k} {\hat{Σ}}_{i}^{- 1})}^{- 1}

(6)

where

{\hat{Σ}}_{i} = \hat{Σ} + Ψ_{i}

, and

\hat{Σ}

is the maximum likelihood estimate of

Σ

obtained by maximizing equation (4). Finally, the desired

(1 - α) 100 %

Wald-type confidence interval for the unknown pooled logit(Se),

μ_{1}

, and pooled logit(Sp),

μ_{2}

, can be obtained using equations (5) and (6).

3.2 The bivariate binomial-normal model

The bivariate generalized random-effects model of Chu and Cole and Chu et al.^3,18 is defined by assuming the binomial distribution for modelling the within-study variability and the bivariate normal distribution for modelling the between-study variability in sensitivity and specificity. That is, the bivariate binomial-normal (BBN) assumes

\begin{array}{l} T P_{i} | b_{1 i} \sim Binomial (n_{1 i}, S e_{i}); y_{1 i} = μ_{1} + b_{1 i} \\ T N_{i} | b_{2 i} \sim Binomial (n_{2 i}, S p_{i}); y_{2 i} = μ_{2} + b_{2 i} \\ b_{i} \sim N_{2} (0, Σ) \end{array}

where TP_i and TN_i are the study-specific number of true positives and true negatives,

n_{1 i}

and

n_{2 i}

are the study-specific total number of diseased and non-diseased subjects, respectively;

y_{1 i}, y_{2 i}, b_{1 i}, b_{2 i}, μ_{1}, μ_{2}

are such that

y_{i} = {(y_{1 i}, y_{2 i})}^{T}, b_{i} = {(b_{1 i}, b_{2 i})}^{T}

, and

μ = {(μ_{1}, μ_{2})}^{T}

The marginal likelihood function of the BBN model has no closed-form expression since the integral in equation (7) does not have a closed form solution. However, the R package Metatron¹⁹ has recently implemented the BBN model by using the adaptive Gaussian quadrature algorithm^3,20 to approximate the likelihood numerically.

L (μ, Σ | y) = \int_{ℝ^{2}} \prod_{i = 1}^{k} f_{y_{i} | b_{i}} (y_{i} | b_{i}, μ) f_{b_{i}} (b_{i} | Σ) d b_{i}

(7)

3.3 The proposed model

This section is devoted to our proposed bivariate random-effects model. Before that, we review the symmetric multivariate Laplace (or double-exponential) distribution.

3.3.1 Review of the symmetric multivariate Laplace distribution

Kotz et al.¹³ defined a p-variate random vector X to have a symmetric multivariate Laplace distribution with a location parameter $0$ and a p × p positive-definite dispersion matrix $Σ$ , say $X \sim L_{p} (0, Σ)$ if it has a characteristic function

φ_{X} (t) = \frac{1}{1 + \frac{1}{2} t' Σ t}

and density function

f (x | 0, Σ) = \frac{2}{{(2 π)}^{\frac{p}{2}} | Σ |^{\frac{1}{2}}} {(\frac{x^{'} Σ^{- 1} x}{2})}^{ν / 2} K_{ν} (\sqrt{2 x' Σ^{- 1} x})

where

ν = (2 - p) / 2

and

K_{ν} (.)

is the modified Bessel function of the third kind (see Kotz et al.¹³ for detail).

As discussed in Eltoft et al.²¹ and Wang et al.,²² the multivariate Laplace distribution can be approximated by

f (x | 0, Σ) \sim \frac{\exp (- \sqrt{2 x' Σ^{- 1} x})}{2^{\frac{2 p - 1}{4}} π^{\frac{p - 1}{2}} | Σ |^{\frac{1}{2}} {(x^{'} Σ^{- 1} x)}^{\frac{1}{4}}}

(8)

upon using the asymptotic formula for the Bessel function:

K_{ν} (y) \sim \sqrt{\frac{π}{2 y}} \exp (- y) as y \to \infty

The mean and covariance matrix of X found from its characteristics function is given, respectively, to be

E (X_{i}) = 0 and var (X_{i}) = Σ

3.3.2 The bivariate normal-Laplace model

Our proposed bivariate normal-Laplace (BNL) model generalizes the BNN model by proposing a heavy-tailed distribution instead of the commonly used normal distribution for the random-effects. More specifically, we model the within-study variability by the bivariate normal distribution and the between-study variability by the symmetric bivariate Laplace distribution. Assuming that $b_{i}$ and $e_{i}$ are independent for each study, $i = 1, \dots, k$ , the BNL model can be defined as

\begin{array}{l} Y_{i} = μ + b_{i} + e_{i}, i = 1, \dots, k \\ b_{i} \sim L_{2} (0, Σ), i = 1, \dots, k, is the random-effects \\ e_{i} \sim N_{2} (0, Ψ_{i}), i = 1, \dots, k, is the random-error term \end{array}

(9)

where

μ, Σ

and

Ψ_{i}

are as defined earlier.

Deriving the marginal distribution of the response variable in equation (9) would let one to estimate parameters of the BNL model by the maximum likelihood method. To do so, we define a new random vector $η_{i}$ = $b_{i} + e_{i}$ and rewrite equation (9) as

Y_{i} = μ + η_{i}

Then using the assumption that $b_{i}$ and $e_{i}$ are independent for each study and by employing properties of characteristic functions, we obtain the characteristic function of $η_{i}$ to be

\begin{array}{l} φ_{η_{i}} (t) = φ_{b_{i} + e_{i}} (t) = φ_{b_{i}} (t) φ_{e_{i}} (t) \\ = (\frac{1}{1 + \frac{1}{2} t' Σ t}) exp (- \frac{1}{2} t' Ψ_{i} t) \\ = \frac{exp (- \frac{1}{2} t' Ψ_{i} t)}{1 + \frac{1}{2} t' Σ t} \end{array}

(10)

Using the uniqueness of characteristic functions, equation (10) is recognized as the characteristic function of the symmetric bivariate normal-Laplace distribution, $N L_{2} (0, Σ, Ψ_{i})$ . That is, $η_{i} \sim N L_{2} (0, Σ, Ψ_{i})$ . Therefore, using Property 1 of Jose and Thomas,²³ $Y_{i} \sim N L_{2} (μ, Σ, Ψ_{i})$ with characteristic function

φ_{y_{i}} (t) = \frac{exp (i t' μ - \frac{1}{2} t' Ψ_{i} t)}{1 + \frac{1}{2} t' Σ t}

(11)

As shown in Appendix A (available on-line as a supplementary material), the mean vector and covariance matrix of the marginal model of our BNL random-effects model can easily be obtained from equation (11) to be

E (Y_{i}) = μ, and var (Y_{i}) = Σ + Ψ_{i}

(12)

We want to note that our marginal model has no explicit density function since there is no closed-form density function for the multivariate NL distribution in general. Therefore, we cannot directly maximize the likelihood function to get estimates of our model parameters. That leads us to the next section—methods of estimating the parameters of our BNL model.

3.3.3 Maximum likelihood estimation of the parameters of the BNL model

Our BNL model can also be defined in a hierarchical way as

\begin{array}{l} Y_{i} | b_{i} \sim N_{2} (μ + b_{i}, Ψ_{i}), i = 1, \dots, k \\ b_{i} \sim L_{2} (0, Σ), i = 1, \dots, k \end{array}

leading to the following likelihood function for the marginal model

L (μ, Σ; y_{i}) = \prod_{i = 1}^{k} \int_{ℝ^{2}} f_{y_{i} | b_{i}} (y_{i} | b_{i}; μ, Ψ_{i}) f_{b_{i}} (b_{i}; Σ) d b_{i}

(13)

Due to the above-mentioned reasons, the likelihood in equation (13) does not have a closed-form expression. Therefore, by treating the random-effects ( $b_{i}$ ) as missing data, we propose the Monte Carlo (MC) version of the expectation-maximization (EM) algorithm known as the MCEM algorithm, to estimate the five parameters ( $μ_{1}, μ_{2}, σ_{1}^{2}, σ_{12}, σ_{2}^{2}$ ) of our BNL model.

3.3.4 The MCEM algorithm

Suppose that $y = {(y_{1}, y_{2}, \dots, y_{k})}^{T}$ and $b = {(b_{1}, b_{2}, \dots, b_{k})}^{T}$ denote the observed data and latent or unobserved data, respectively, and let $y_{c} = (y, b)$ represent the complete data. Let $θ = (μ, Σ)$ be the set of parameters of the BNL model. The complete-data likelihood and log-likelihood function are respectively given as

\begin{array}{l} L_{c} (θ; y_{c}) = \prod_{i = 1}^{k} f_{y | b} (y_{i} | b_{i}; μ, Ψ_{i}) f_{b} (b_{i}; Σ), and \\ l_{c} (θ; y_{c}) = \sum_{i = 1}^{k} [log {f_{y | b} (y_{i} | b_{i}; μ, Ψ_{i})} + log {f_{b} (b_{i}; Σ)}] \end{array}

The EM algorithm of Dempster et al.²⁴ is an iterative process that involves computing the Q-function—the expected value of the complete-data log-likelihood given the observed data and the current estimate of $θ$ at the E-step; and iteratively maximizing the resultant Q-function until convergence at the M-step of the algorithm.

Thus, the expected complete-data log-likelihood function given the observed data and the current estimate of $θ$ , denoted by $Q (θ; \hat{θ})$ is given by

\begin{array}{l} Q (θ; \hat{θ}) = E {l_{c} (θ; y_{c}) | y, \hat{θ}} \\ = E [\sum_{i = 1}^{k} log {f_{y | b} (y_{i} | b_{i}; μ, Ψ_{i})} | y, \hat{θ}] + E [\sum_{i = 1}^{k} log {f_{b} (b_{i}; Σ)} | y, \hat{θ}] \\ = \sum_{i = 1}^{k} E {log (f_{y | b} (y_{i} | b_{i}; μ, Ψ_{i})) | y, \hat{θ}} + \sum_{i = 1}^{k} E {log (f_{b} (b_{i}; Σ)) | y, \hat{θ}} \\ = \sum_{i = 1}^{k} \int_{ℝ^{2}} log {f_{y | b} (y_{i} | b_{i}; μ, Ψ_{i})} f_{b | y} (b_{i} | y_{i}; \hat{θ}) d b_{i} + \\ \sum_{i = 1}^{k} \int_{ℝ^{2}} log {f_{b} (b_{i}; Σ)} f_{b | y} (b_{i} | y_{i}; \hat{θ}) d b_{i} \end{array}

(14)

However, as shown in equation (14), we cannot directly compute the Q-function since the expectation (or integral) requires the knowledge of the marginal likelihood function of the response variable y—which does not exist in closed-form. Therefore, we propose the following MCEM algorithm where we approximate the E-step of the deterministic EM algorithm using MC methods.

Specifically, for a given Monte Carlo sample size R, we generate samples $b_{i, 1}, b_{i, 2}, \dots, b_{i, R}, i = 1, . ., k$ , using Metropolis–Hastings algorithm from the target density $f_{b | y} (b_{i} | y_{i}; \hat{θ})$ by specifying some candidate density say $f (b_{i}; .)$ . The acceptance probability of the Metropolis–Hastings algorithm has a nice form when the marginal density of the random-effects $f_{b_{i}} (b_{i}; \hat{Σ})$ is used as a candidate density as shown in Appendix B available on-line as a supplementary material.

Once the MC samples $b_{i, 1}, b_{i, 2}, \dots, b_{i, R}, i = 1, . ., k$ are obtained, the Q-function in the E-step of our MCEM algorithm can be approximated by

Q_{R} (θ, \hat{θ}) = \frac{1}{R} \sum_{i = 1}^{k} \sum_{r = 1}^{R} log {f_{y | b} (y_{i} | b_{i, r}; μ, Ψ_{i})} + \frac{1}{R} \sum_{i = 1}^{k} \sum_{r = 1}^{R} log {f_{b} (b_{i, r}; Σ)} .

(15)

By the law of large numbers, the estimator in equation (15) converges to the theoretical expectation given in equation (14).

In the maximization step of the algorithm, we maximize $Q_{R} (θ, \hat{θ})$ with respect to $μ$ and $Σ$ . Maximizing $Q_{R} (θ, \hat{θ})$ with respect to $μ$ yields a closed form solution

{\hat{μ}}^{(j)} = {(\sum_{i = 1}^{k} Ψ_{i}^{- 1})}^{- 1} \frac{1}{R} \sum_{i = 1}^{k} \sum_{r = 1}^{R} Ψ_{i}^{- 1} (y_{i} - b_{i, r}^{(j)})

(16)

at the j-th iteration of the algorithm. Maximizing

Q_{R} (θ, \hat{θ})

with respect to

Σ

requires application of numerical methods. In this sequel, we used the optim() function in the R-programming language to obtain the MLE,

{\hat{Σ}}^{(j)}

, of

Σ

at the j-th iteration. Point estimates of the five parameters of the BNN model, obtained after fitting the BNN model to DTA data by excluding the outlying and influential studies, are used to initiate the MCEM algorithm. If the between-study covariance matrix from the BNN model is not positive-definite, we use the 2 × 2 identity matrix as a starting value for

Σ

instead.

3.3.5 Monte Carlo sample size and convergence criteria

When implementing the MCEM algorithm, the MC sample size, R, has been given special consideration in the literature. Tanner²⁵ recommends to start the MCEM algorithm with a small R, say R = 10, when the approximations ${\hat{θ}}^{(j)}$ are not close to the MLE $\hat{θ}$ and to increase R when the iterates become closer to the MLE. Different researchers used different approaches when increasing R from iteration to iteration. McCulloch²⁶ employed a predetermined value of R for different categories of iterations. He suggested to use R = 50, 200 and 500 for iterations 1–19, 20–39 and 40–50, respectively, and to terminate the algorithm after the 50-th iteration. Booth and Hobert,²⁷ Levine and Casella²⁸ and Levine and Fan²⁹ all developed an automated approach that would increase the MC sample size only if the previous estimate ${\hat{θ}}^{(j - 1)}$ lies in an approximate confidence interval about the current estimate ${\hat{θ}}^{(j)}$ .

Though increasing the MC sample size, R, at each iteration of the MCEM algorithm is needed as pointed out by Neath et al.,³⁰ it would be practical to come up with an inexpensive approach, both in terms of programming and computational time. Robert and Casella,³¹ Section 5.4.4, demonstrated the MCEM algorithm by increasing the current MC sample size linearly by doubling the value of R in the previous iteration (i.e. $R^{(j)} = 2 \times R^{(j - 1)}$ ). Although this approach is straightforward to apply and easy to use in practice, the MC sample size will grow exponentially and leads to wastage of the MC sample sizes in the early stages of the MCEM algorithm where the current iterates ${\hat{θ}}^{(j)}$ are far from the true MLE value $\hat{θ}$ .²⁸

In this paper, we suggest increasing the MC sample size at every iteration of the algorithm in a linear way, say $R^{(j)} = R^{(j - 1)} + ⌊ R^{(j - 1)} / c ⌋$ , where $⌊ . ⌋$ means the integer part and c is a positive constant. This approach has been previously used by Booth and Hobert²⁷ and Levine and Casella²⁸ using c = 3, 4, 5 but only when the approximate confidence interval of the current estimate ${\hat{θ}}^{(j)}$ includes the previous estimate value ${\hat{θ}}^{(j - 1)}$ . We have successfully implemented our approach using a starting MC sample size of R = 10 and c = 5.

Another point that needs to be discussed is when to stop the algorithm and declare the current parameter estimates as the MLEs. Although several methods have been suggested in the literature, the criteria based on relative error suggested by Booth and Hobert²⁷ and Neath et al.³⁰ is appealing to be used in this sequel. For some user specified positive constants δ and ϵ, the ordinary EM algorithm will report convergence at the j-th iteration if the following relative error criterion is satisfied

{\max_{1 \leq i \leq 5}}_{} (\frac{| θ_{i}^{(j)} - θ_{i}^{(j - 1)} |}{| θ_{i}^{(j)} | + δ}) < ϵ

(17)

For the MCEM algorithm, due to the introduction of MC error in the E-step of the algorithm, both Booth and Hobert²⁷ and Neath et al.³⁰ recommend to declare convergence only if criteria (17) is met for three successive iterations. In this paper, we used criteria (17) with $δ = 0.001$ and $ϵ = 0.005$ following Booth and Hobert²⁷; and we would declare convergence if criteria (17) is met for three consecutive iterations, covariance matrix of the estimates (see section 3.3.6) is positive-definite and maximum number of iterations is not reached. Appendix C, which is available as a supplementary material, summarizes the MCEM algorithm used in this article.

3.3.6 Estimation of standard errors and model selection criteria

Since the EM algorithm does not provide standard errors as a byproduct, several researchers have developed alternative methods that approximate standard errors of parameters estimated by the EM algorithm. Among these methods, the Oakes’ approach³² is relatively simple to apply. It is well understood that, asymptotically, the inverse of the observed Fisher’s information matrix (FIM), denoted by I( $θ$ ), approximates the variance of the MLE estimator. That is

var (\hat{θ}) = I {(\hat{θ})}^{- 1} = {[- \frac{\partial^{2} log {L (θ | y)}}{\partial θ^{2}}]}^{- 1} | θ = \hat{θ}

Oakes³² provided an approximation to the matrix of second-order derivatives by approximating the Q-function using the following MC sum

\begin{array}{l} \frac{\partial^{2} log {L (θ | y)}}{\partial θ^{2}} = \frac{1}{R} \sum_{r = 1}^{R} \frac{\partial^{2} log {L (θ' | y; b_{r})}}{\partial θ^{2}} + \\ \frac{1}{R} \sum_{r = 1}^{R} {[\frac{\partial log {L (θ^{'} | y; b_{r})}}{\partial θ} - \frac{1}{R} \sum_{r = 1}^{R} \frac{\partial log {L (θ^{'} | y; b_{r})}}{\partial θ}]}^{2} \end{array}

(18)

using the random-effects

b_{i, 1}, b_{i, 2}

,…,

b_{i, R}, i = 1, \dots, k,

which have been already simulated at the last iteration of the MCEM algorithm.³³ In Appendix D (available as an on-line supplementary material), we have provided the first- and second-order partial derivatives of the complete-data log-likelihood required to compute equation (18).

We may compare goodness-of-fit of the proposed model against the well-established ones via AIC and BIC by approximating the marginal log-likelihood $L (\hat{θ}; y)$ through Monte Carlo. That is, since

\begin{array}{l} log L (\hat{θ}; y) = log \prod_{i = 1}^{k} f (y_{i}; \hat{θ}) = log \prod_{i = 1}^{k} \int_{ℝ^{2}} f (y_{i} | b_{i}; \hat{θ}) f (b_{i}; \hat{θ}) d b_{i} \\ = log \prod_{i = 1}^{k} E {f (y_{i} | b_{i}; \hat{θ})} \end{array}

at the last iteration of the MCEM algorithm, following Chen et al.,³⁴ we can generate large MC samples

b_{i, 1}, b_{i, 2}, \dots, b_{i, R}, i = 1, \dots, k

from

f (b_{i}; \hat{θ})

and approximate the marginal likelihood by the MC average

log L (\hat{θ}; y) = \sum_{i = 1}^{k} log {\frac{1}{R} \sum_{r = 1}^{R} f (y_{i} | b_{i, r}; \hat{θ})} .

Using the value of R and the MC samples $b_{i, 1}, b_{i, 2}, \dots, b_{i, R}, i = 1, \dots, k$ obtained at the last iteration of the MCEM algorithm, we compute the AIC and BIC as

\begin{array}{l} AIC = - 2 log L (\hat{θ}; y) - 2 m, and \\ BIC = - 2 log L (\hat{θ}; y) - log (2 k) m \end{array}

where m, k and

log L (\hat{θ}; y)

are the number of free parameters to be estimated, number of studies in the meta-analysis, and the maximized log-likelihood value, respectively. The smaller the AIC or BIC, the better the model.

4 Simulation study

4.1 Simulation design

A simulation study is performed to demonstrate the performance of the methods discussed in this sequel. A total of 56 settings has been considered by varying the random-effects distribution: bivariate normal or bivariate Laplace, model parameters: $μ$ and $Σ$ , and data characteristics: n and k. To inform our simulation, we reanalyzed 165 empirical meta-analyses obtained from the Cochrane database for diagnostic test reviews from 2011 to 2015. Accordingly, the following set of parameters is considered in our simulation study:

The overall Se and Sp for simulating DTA data with small mean vectors: (Se=0.5, Sp=0.6).

The overall Se and Sp for generating DTA data with large mean vectors: (Se=0.90, Sp=0.95).

The between-study variances and covariance parameter for simulating DTA data with no outlying studies: ( $σ_{1}^{2}$ =0.5, $σ_{2}^{2}$ =0.6, σ₁₂= −0.2).

The between-study variances and covariance parameter for generating DTA data with outlying studies: ( $σ_{1}^{2}$ =(0.2,1.5), $σ_{2}^{2}$ =(0.25,1.6), σ₁₂= −0.2).

The number of subjects in the diseased (n₁) and non-diseased (n₂) category: (n₁= 100, n₂= 200).

The number of studies in the meta-analysis: (k=10, 20, 30, and 40).

Seven different scenarios have been considered to systematically introduce outlying and influential studies by varying $μ$ and $Σ$ . Our definition of ‘small’ and ‘large’ between-study variances follows those of Hamza et al.³⁵ and Negeri et al.³⁶ We define outlying or influential studies as those studies in the meta-analysis which have large Se or Sp, and small or large between-study variances than the rest of the studies. In the first and second scenario, studies that are outlying and influential both in Se and Sp with large ( $σ_{1}^{2} = 1.5, σ_{2}^{2} = 1.6$ ) and small ( $σ_{1}^{2} = 0.2, σ_{2}^{2} = 0.25$ ) between-study variances, respectively, are considered. In the third and fourth scenario, we introduced studies which are outlying and influential only in Se with large and small between-study variances, respectively. In the fifth and sixth scenario, studies that are outlying and influential only in Sp with large and small between-study variances, respectively, are considered. Finally, studies that are outlying and influential neither in Se nor in Sp are generated. These seven scenarios are considered to illustrate how robust the proposed method is regardless of the presence or absence of outlying studies. We introduced 2, 3, 4, and 5 outlying and influential studies into the meta-analyses with 10, 20, 30, and 40 number of studies, respectively. These seven scenarios for incorporating outlying studies, along with the four k values, and the two distributional assumptions for the random-effects constitute the 56 scenarios we considered in our simulation study.

4.2 Simulation results

When outlying and influential studies are introduced both in Se and Sp, the between-studies variances are large, and random-effects are generated from the bivariate normal distribution (Figure 3), the proposed BNL model produced the most robust estimates, whereas the standard models yielded an inflated overall sensitivity and specificity estimates. For these scenarios, the absolute (relative) change in sensitivity between the proposed model and the BNN model ranges between 7% (13%) and 10% (21%), and the absolute (relative) difference in sensitivity between the proposed model and the BBN varies from 0% (0%) to 8% (16%). The absolute (relative) difference in specificity varies from 3% (4%) to 10% (16%), between the proposed BNL model and BNN, and from 1% (2%) to 5% (8%) between the BNL model and BBN model. For the same scenario but when the random-effects are generated from the bivariate Laplace distribution, the absolute (relative) change in Se between the proposed BNL and BNN model, and between the BNL and BBN model varies from 6% (12%) to 19% (38%), and from 2% (3%) to 12% (26%), respectively. On the other hand, the absolute (relative) change in Sp between the BNL and BNN model, and between the BNL and BBN model, respectively, ranges between 4% (5%) and 10% (16%), and between 1% (1%) and 6% (9%). Similar results have been observed when the small between-study variances have been employed.

Figure 3.

Point estimates of sensitivity (left) and specificity (right) of the three methods when data are generated from the bivariate normal distribution, $σ_{12} = - 0.2, n_{1} = 100, n_{2} = 200$ , outlying and influential studies are generated using large between-study variances ( $σ_{1}^{2} = 1.5, σ_{2}^{2} = 1.6$ ) and are introduced both in sensitivity and specificity.

Figure 4 presents the results when outlying and influential studies are introduced only in Se, the between-study variances are large, and random-effects are generated from the bivariate normal distribution. The proposed BNL model yielded robust estimates of sensitivity and specificity compared to the standard models’ estimates. The absolute (relative) change in sensitivity of the BNL model and the BNN model varies from 1% (1%) to 14% (23%), and that of the BNL model and the BBN model ranges between 1% (2%) and 6% (12%). When the bivariate Laplace distribution is used to generate the random-effects, the absolute (relative) change in Se between the BNL and BNN model, and the BNL and BBN model, respectively, ranges from 1% (1%) to 8% (16%), and from 1% (2%) to 6% (13%). Moreover, the robustness of our proposed BNL model is justified as it yielded similar estimates of specificity to the standard models. We observed similar results when the small between-study variances are used.

Figure 4. Point estimates of sensitivity (left) and specificity (right) of the three methods when data are generated from the bivariate normal distribution, $σ_{12} = - 0.2, n_{1} = 100, n_{2} = 200$ , outlying and influential studies are generated using large between-study variances ( $σ_{1}^{2} = 1.5, σ_{2}^{2} = 1.6$ ) and are introduced only in sensitivity.

The results when outlying and influential studies are incorporated only in Sp, the between-study variances are large, and random-effects are simulated from the bivariate normal distribution are displayed in Figure 5. The absolute (relative) change in specificity of the BNL model and the BNN model varies from 4% (6%) to 6% (10%), and that of between the BNL model and the BBN model ranges between 1% (1%) and 3% (4%). On the other hand, the absolute (relative) change in Sp between the BNL and BNN model, and the BNL and BBN model was between 2% (2%) and 9% (15%), and between 0% (0%) and 10% (13%), respectively, when the random-effects are simulated from the bivariate Laplace distribution. Once again our proposed model produced similar estimates in sensitivity to the standard models suggesting that the BNL model is robust to outlying and influential studies. Analogous results had been noticed when the small between-study variances are assumed.

Figure 5. Point estimates of sensitivity (left) and specificity (right) of the three methods when data are generated from the bivariate normal distribution, $σ_{12} = - 0.2, n_{1} = 100, n_{2} = 200$ , outlying and influential studies are generated using large between-study variances ( $σ_{1}^{2} = 1.5, σ_{2}^{2} = 1.6$ ) and are introduced only in specificity.

Finally, Figure 6 presents the results of the three models when there are no outlying and influential studies neither in Se nor in Sp, and the random-effects are generated from the bivariate normal distribution. Figure 6 depicts that the three models yielded similar estimates of sensitivity and specificity when there are no outlying and influential studies. This result further illustrates that our proposed BNL model is robust to outlying and influential studies. Similar results have been observed when the random-effects are generated from the bivariate Laplace distribution.

Figure 6. Point estimates of sensitivity (left) and specificity (right) of the three methods when data are generated from the bivariate normal distribution, $σ_{12} = - 0.2, n_{1} = 100, n_{2} = 200$ , outlying and influential studies are generated using large between-study variances ( $σ_{1}^{2} = 1.5, σ_{2}^{2} = 1.6$ ) and are introduced neither in sensitivity nor in specificity.

To understand the effect of smaller within-study sample sizes and varying the between-study covariance between test-sensitivity and test-specificity, we performed additional simulation study by considering “small” (n₁=30, n₂=50) and “medium” (n₁=50, n₂=100) within-study sample sizes, and zero (σ₁₂=0) and non-zero (σ₁₂=−0.5) between-study covariance in addition to the parameters considered in Section 4.1. We summarize our findings below and provide the detailed results with additional figures in Appendix E of the on-line supplementary material.

Figures 7 to 10 present our additional simulation study results when there is non-zero covariance between test-sensitivity and specificity, the within-study sample sizes are medium (n₁ = 50, n₂ = 100), outlying and influential studies are generated using large between-study variances ( $σ_{1}^{2}$ =1.5, $σ_{2}^{2}$ =1.6), and data are generated from the bivariate Laplace distribution. According to these figures, the absolute (relative) change in sensitivity between the BNL and the BNN model and between the BNL and the BBN model varies from 5% (9%) to 10% (16%) and from 1% (3%) to 8% (13%), respectively. In specificity, these differences correspond respectively to 7% (11%) to 11% (17%) and 5% (7%) to 8% (11%). When the random-effects are generated from the bivariate normal distribution, the absolute (relative) change in sensitivity between the proposed BNL and BNN model, and between the BNL and BBN model varies from 5% (10%) to 11% (27%), and from 1% (1%) to 6% (11%), respectively. On the other hand, the absolute (relative) change in specificity between the BNL and BNN model, and between the BNL and BBN model, respectively, ranges between 3% (5%) and 9% (14%), and between 2% (2%) and 6% (9%).

Figure 7.

Point estimates of sensitivity (left) and specificity (right) of the three methods when data are generated from the bivariate Laplace distribution, $σ_{12} = - 0.5, n_{1} = 50, n_{2} = 100$ , outlying and influential studies are generated using large between-study variances ( $σ_{1}^{2} = 1.5, σ_{2}^{2} = 1.6$ ) and are introduced both in sensitivity and specificity.

Figure 8. Point estimates of sensitivity (left) and specificity (right) of the three methods when data are generated from the bivariate Laplace distribution, $σ_{12} = - 0.5, n_{1} = 50, n_{2} = 100$ , outlying and influential studies are generated using large between-study variances ( $σ_{1}^{2} = 1.5, σ_{2}^{2} = 1.6$ ) and are introduced only in sensitivity.

Figure 9.

Figure 10.

In summary, the proposed BNL model produced less biased point estimates of the pooled sensitivity and specificity, precise corresponding confidence intervals, and is robust to model misspecification. Moreover, it was observed that the influence of outlying and influential studies on both the point and interval estimates of Se and Sp upsurges as both the within-study sample sizes and between-study covariance increase.

5 Illustrative examples

This section will discuss the results of the standard and proposed methods when applied to the two datasets introduced in section 2.

5.1 NEE data

Table 1 displays results of the three methods when fitted to the NEE dataset. The results indicate that the presence of three outlying and influential studies made a significant difference in the point and interval estimates of the pooled sensitivity and specificity among the models. A particularly significant difference is observed between the results of the standard BBN model and the proposed BNL model. Study 3, which was detected as influential in sensitivity, caused the BBN model to yield an overall Se which is larger than that of the proposed model by 4.4% (6.4%) in absolute (relative) value. On the other hand, the two influential studies in specificity resulted in a 5.6% and 6.6% drop, respectively, in the absolute and relative Sp of the BBN model compared to the BNL model. The impact of these influential studies in absolute (relative) value corresponded to a 1.6% (2.3%) increase in pooled Se and a 2.2% (2.6%) decrease in pooled Sp of the BNN model compared to the BNL model. Moreover, the proposed model produced precise estimates as it yielded narrower confidence intervals. Finally, both the AIC and BIC selected the proposed model as the best model to explain the NEE dataset.

Table 1.

Sensitivity (95% CI), specificity (95% CI) in percent and model comparison statistics for the NEE data.

Model	Se (95% CI)	Sp (95% CI)	AIC	BIC
BBN	73.2 (55.6, 85.6)	79.1 (39.4, 95.7)	158.2	162.6
BNN	70.4 (58.9, 79.7)	82.5 (60.4, 93.6)	68.6	73.1
BNL	68.8 (57.5, 78.3)	84.7 (80.6, 88.0)	68.2	72.6

BBN: bivariate binomial-normal; BNN: bivariate normal-normal; BNL: bivariate normal-Laplace.

5.2 MMSE data

We present the meta-analysis results of the methods when applied to the MMSE data in Table 2. As expected, the three models yielded similar point estimates in the overall sensitivity. However, the two influential studies in specificity caused both the BNN and BBN models to result in a respective 6.0% (7.4%) and 5.8% (7.3%) reduction in absolute (relative) Sp compared to the proposed BNL model. Additionally, our proposed model yielded precise estimates as observed from the CIs, and also chosen by the AIC and BIC to be the best model to describe the MMSE data.

Table 2.

Sensitivity (95% CI), specificity (95% CI) in percent and model comparison statistics for the MMSE data.

Model	Se (95% CI)	Sp (95% CI)	AIC	BIC
BBN	59.5 (46.9, 68.9)	74.3 (53.1, 88.1)	132.3	136.2
BNN	60.0 (46.7, 78.9)	74.2 (58.9, 85.2)	48.8	52.7
BNL	58.3 (57.5, 78.3)	80.1 (69.5, 87.7)	47.5	51.4

BBN: bivariate binomial-normal; BNN: bivariate normal-normal; BNL: bivariate normal-Laplace.

6 Discussion

This article aimed at addressing the need for a robust statistical method for meta-analysis of DTA data when there are outlying and/or influential studies. With the aim of making a robust inference, a new bivariate random-effects meta-analysis model that automatically downweights the impact of outlying and influential studies on the point and interval estimates of the overall sensitivity and specificity is proposed and compared against the well-established meta-analysis methods for DTA studies.

The proposed method is motivated by using two real-life meta-analyses and had also been demonstrated and validated using simulated data. The simulation study was designed to demonstrate the performance of our proposed BNL model when outlying and influential studies are present only in Se, only in Sp, both in Se and Sp, and neither in Se nor in Sp. Under each of the four settings, regardless of the other parameters varied, our proposed model produced robust point estimates in Se and Sp compared to the standard models. Of the standard models, the BBN model yielded relatively robust estimates compared to the BNN model.

In terms of confidence intervals for sensitivity and specificity, our proposed BNL model produced precise confidence intervals relative to the standard models regardless of the parameters varied in our simulation including the distributional assumption for the random-effects. Comparing the standard models, the BBN model yielded narrower confidence intervals than the BNN model. Results of the two published meta-analyses were also in line with our simulation study as the proposed BNL model produced robust and precise point and interval estimates of the overall sensitivity and specificity. Moreover, the proposed model was shown to have demonstrated a better fit according to the information criteria (AIC and BIC).

Although we have chosen to use parameter estimates of the BNN model as starting values to start the MCEM algorithm of our proposed BNL model, we observed that the algorithm converged to similar values when employing different starting values, showing that the proposed method is robust to the choice of starting values. Additionally, the proposed method converged for all the scenarios considered in our simulation study and real-data examples. Obviously, as for every EM algorithm, fitting our proposed BNL model takes longer time than the standard methods. However, due to its demonstrated superior performance to the standard models and the need for making an accurate inference to inform evidence-based decision-making, it is worth waiting the extra time to use the robust method we proposed here. In fact, the computational time can significantly be reduced by using parallel computing as we did for our study in this paper. We implemented our methodology using the R open-source statistical software, and the R code is available online as a supplementary material on the publisher’s website.

Our proposed BNL model has potential limitations. Although the EM algorithm is in general slower than other methods of finding MLE, the MCEM algorithm is even slower due to the required MC samples at each iteration of the algorithm. As a result, the BNL model takes longer computational time compared to the standard methods. Additionally, the BNL model does not discuss the impact of outlying and/or influential studies on the summary receiver operating characteristics (SROC) curve.

In summary, we proposed a flexible bivariate random-effects model for meta-analysis of DTA studies with the aim of accommodating outlying and influential studies, as a direct response to the study of Negeri and Beyene⁴ which called for such methods to be developed. Our methodology is based on the assumption that the within-and between-study variation in sensitivity and specificity can approximately be described by the bivariate normal and Laplace distribution, respectively. That led to a new bivariate random-effects model which has the same number of parameters as the standard models but which possesses flexibility (as it avoids removing potential outlying and influential studies) and robustness to the distributional assumption for the random-effects and presence of outlying and influential studies. Due to its flexibility and robustness properties, we believe that our proposed BNL method is a good addition to the meta-analytic literature and an alternative to the standard models for practitioners and researchers in general.

Supplemental Material

sj-zip-1-smm-10.1177_0962280220925840 - Supplemental material for Robust bivariate random-effects model for accommodating outlying and influential studies in meta-analysis of diagnostic test accuracy studies

Supplemental material, sj-zip-1-smm-10.1177_0962280220925840 for Robust bivariate random-effects model for accommodating outlying and influential studies in meta-analysis of diagnostic test accuracy studies by Zelalem F Negeri and Joseph Beyene in Statistical Methods in Medical Research

Supplemental Material

sj-pdf-2-smm-10.1177_0962280220925840 - Supplemental material for Robust bivariate random-effects model for accommodating outlying and influential studies in meta-analysis of diagnostic test accuracy studies

Supplemental material, sj-pdf-2-smm-10.1177_0962280220925840 for Robust bivariate random-effects model for accommodating outlying and influential studies in meta-analysis of diagnostic test accuracy studies by Zelalem F Negeri and Joseph Beyene in Statistical Methods in Medical Research

Footnotes

Acknowledgements

JB holds the John D. Cameron Endowed Chair in the Genetic Determinants of Chronic Diseases, McMaster University.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: JB would like to acknowledge funding from the Natural Sciences and Engineering Research Council of Canada (NSERC) and the Canadian Institutes of Health Research (CIHR).

ORCID iDs

Zelalem F Negeri

Joseph Beyene

Supplemental material

Supplementary material for this article is found on-line on the publisher’s website.

References

Rutter

Gatsonis

CA.

A hierarchical regression approach to meta-analysis of diagnostic test accuracy evaluations.

Stat Med 2001; 20: 2865–2884.

Reitsma

Glas

Rutjes

, et al. Bivariate analysis of sensitivity and specificity produces informative summary measures in diagnostic reviews. J Clin Epidemiol 2005; 58: 982–990.

Chu

Cole

SR.

Bivariate meta-analysis of sensitivity and specificity with sparse data: a generalized linear mixed model approach.

J Clin Epidemiol 2006; 59: 1331.

Negeri

Beyene

Statistical methods for detecting outlying and influential studies in meta-analysis of diagnostic test accuracy studies. Stat Meth Med Res 2020; 29: 1227--1242.

Baker

Jackson

A new approach to outliers in meta-analysis.

Health Care Manag Sci 2008; 11: 121–131.

Lee

Thompson

SG.

Flexible parametric models for random-effects distributions.

Stat Med 2008; 27: 418–434.

Viechtbauer

Cheung

MWL.

Outlier and influence diagnostics for meta-analysis.

Res Syn Meth 2010; 1: 112–125.

Gumedze

Jackson

A random effects variance shift model for detecting and accommodating outliers in meta-analysis.

BMC Med Res Methodol 2011; 11: 19.

Beath

KJ.

A finite mixture method for outlier detection and robustness in meta-analysis.

Res Syn Meth 2014; 5: 285–293.

10.

Baker

Jackson

New models for describing outliers in meta-analysis.

Res Syn Meth 2016; 7: 314–328.

11.

Jackson

White

IR.

When should meta-analysis avoid making hidden normality assumptions?

Biom J 2018; 60: 1040–1058.

12.

Negeri

Beyene

Skew-normal random-effects model for meta-analysis of diagnostic test accuracy (DTA) studies. Biom J 2020; DOI: 10.1002/bimj.201900184

13.

Kotz

Kozubowski

Podgorski

The Laplace distribution and generalizations: a revisit with applications to communications, economics, engineering, and finance. Boston, MA: Springer Science & Business Media, 2001.

14.

Vestergaard

Macaskill

Holt

, et al. Dermoscopy compared with naked eye examination for the diagnosis of primary melanoma: a meta-analysis of studies performed in a clinical setting. Br J Dermatol 2008; 159: 669–676.

15.

Arevalo-Rodriguez

Smailagic

Figuls

, et al. Mini-mental state examination (MMSE) for the detection of Alzheimer’s disease and other dementias in people with mild cognitive impairment (MCI). Cochrane Database Syst Rev 2015; DOI: 10.1002/14651858.CD010783.pub2

16.

Doebler

mada: meta-analysis of diagnostic accuracy. R package version 0.5. 7, 2015.

17.

Gasparrini

Armstrong

Kenward

MG.

Multivariate meta-analysis for non-linear and other multi-parameter associations.

Stat Med 2012; 31: 3821–3839.

18.

Chu

Guo

Zhou

Bivariate random effects meta-analysis of diagnostic studies using generalized linear mixed models.

Med Decis Making 2010; 30: 499–508.

19.

Huang

Metatron: meta-analysis for classification data and correction to imperfect reference. R package version 0.1-1. 2014, 2018.

20.

Paul

Riebler

Bachmann

, et al. Bayesian bivariate meta-analysis of diagnostic test studies using integrated nested Laplace approximations. Stat Med 2010; 29: 1325–1339.

21.

Eltoft

Kim

Lee

TW.

On the multivariate Laplace distribution. IEEE Signal Process Lett 2006; 13: 300–303.

22.

Wang

Zhang

Zhao

Multivariate Laplace filter: a heavy-tailed model for target tracking. In: 19th international conference on pattern recognition, Tampa, FL, USA, 8--11 December 2008, pp. 1–4. Piscataway, NJ: IEEE.

23.

Jose

Thomas

MM.

Multivariate normal-Laplace distribution and processes.

Statistica 2014; 74: 27–44.

24.

Dempster

Laird

Rubin

DB.

Maximum likelihood from incomplete data via the EM algorithm. J R Stat Soc Ser B 1977; 39: 1–38.

25.

Tanner

MA.

Tools for statistical inference: observed data and data augmentaion methods. Lecture Notes Stat 1993; 67: 1--117.

26.

McCulloch

CE.

Maximum likelihood algorithms for generalized linear mixed models. J Am Stat Assoc 1997; 92: 162–170.

27.

Booth

Hobert

JP.

Maximizing generalized linear mixed model likelihoods with an automated Monte Carlo EM algorithm. J R Stat Soc Ser B 1999; 61: 265–285.

28.

Levine

Casella

Implementations of the Monte Carlo EM algorithm. J Comput Graph Stat 2001; 10: 422–439.

29.

Levine

Fan

An automated (Markov chain) Monte Carlo EM algorithm. J Stat Comput Simul 2004; 74: 349–360.

30.

Neath

, et al. On convergence properties of the Monte Carlo EM algorithm. In: Jones G and Shen X (eds) Advances in modern statistical theory and applications: a festschrift in honor of Morris L. Eaton. Beachwood, OH: Institute of Mathematical Statistics, 2013. pp. 43–62.

31.

Robert

Casella

Introducing Monte Carlo methods with R. Vol. 18. New York, NY: Springer-Verlag, 2010.

32.

Oakes

Direct calculation of the information matrix via the EM. J R Stat Soc Ser B 1999; 61: 479–482.

33.

Robert

Casella

Monte Carlo statistical methods. New York, NY: Springer Science & Business Media, 2004.

34.

Chen

Zhang

Davidian

A Monte Carlo EM algorithm for generalized linear mixed models with flexible random effects distribution.

Biostatistics 2002; 3: 347–360.

35.

Hamza

Reitsma

Stijnen

Meta-analysis of diagnostic studies: a comparison of random intercept, normal-normal, and binomial-normal bivariate summary ROC approaches.

Med Decis Making 2008; 28: 639–649.

36.

Negeri

Shaikh

Beyene

Bivariate random-effects models for meta-analysis of diagnostic test accuracy studies using arcsine-based transformations. Biom J 2018; 60: 827–844.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.04 MB

0.00 MB

0.47 MB