In VNyaga/NMADAS: Arm-Based Network Meta-Analaysis of Diagnostic Accuracy Studies

Network meta-analysis (NMA) allow combining efficacy information from multiple comparisons from trials assessing different therapeutic interventions for a given disease and to estimate unobserved comparisons from a network of observed comparisons. Applying NMA on diagnostic accuracy studies is a statistical challenge given the inherent correlation of sensitivity and specificity.

In this tutorial we present two models for NMA of diagnostic accuracy studies. The first is a conceptually simple and novel hierarchical arm-based (AB) model which expresses the logit transformed sensitivity and specificity as sum of fixed effects for test, correlated study-effects and a random error associated with various tests evaluated in given study is proposed.

Despite the familiariry and popularity of the normal distribution, modeling the latent sensitivity and specificity using the normal distribution after transformation is neither natural nor computationally convenient.

Thus the second model uses a bivariate beta distribution to model sensitivity and specificity in their natural scale. Using the beta distribution in regression has the following advantages, that the probabilities are modeled in their proper scale rather than a monotonic transform of the probabilities. Secondly, the model is flexible as it allows for asymmetry often present in the distribution of bounded variables such as proportions which is the case with sparse data common in meta-analysis. Thirdly, the model provides parameters with direct meaningful interpretation.

Introduction

Network meta-analyses (NMA) have classically been used to extend conventional pairwise meta-analyses by combining and summarizing direct and indirect evidence on multiple `therapeutic' interventions for a given condition when the set of evaluated interventions/treatments differs among studies. By borrowing strength from the indirect evidence, there is a potential gain in precision of the estimates \cite{Higgins}. Furthermore, the estimates may be less biased and more robust. Such an approach uses the data efficiently and is line with the principle of intention-to-treat (ITT) \cite{Fisher} in randomized clinical trials which requires that all valid available data should be used even when a part of the data is missing.

In a diagnostic test accuracy study, an index test and possibly one or more comparator tests are administered to each tested subject. A standard or reference test or procedure is also applied to all the patients to classify them as having the target condition or not. The patients results are then categorized by the index and reference test as true positive, false positive, true negative and false negative. The diagnostic accuracy of the index test is represented as a bivariate outcome and is typically expressed as sensitivity and specificity at a defined test cutoff. Differences due to chance, design, conduct, patients/participants, interventions, tests and reference test imply there will be heterogeneity often in opposite direction for the two typical accuracy outcomes: sensitivity and specificity. While traditional meta-analyses allow for comparison between two tests, there are often multiple tests for the diagnosis of a particular disease outcome. To present the overall picture, inference about all the tests for the same condition and patient characteristics is therefore required. The simultaneous analysis of the variability in the accuracy of multiple tests within and between studies may be approached through a network meta-analysis.

In combining univariate summaries from studies where the set of tests differs among studies two types of linear mixed models have been proposed. The majority of network meta-analyses express treatment effects in each study as a contrast relative to a baseline treatment in the respective study \cite{Higgins, Lumley}. This is the so called contrast-based (CB) model. Inspired by the CB models developed for interventional studies, Menten and Lesaffre (2015) \cite{Menten} introduced a CB model for diagnostic test accuracy data to estimate the average log odds ratio for sensitivity and specificity of the index test relative to a baseline or comparator test.

The second type of models is the classical two-way ANOVA model with random effects for study and fixed effect for tests \cite{Senn, Whitehead, Piepho12}, the so called arm-based (AB) model. The AB model is based on the assumption that the missing arms or tests are missing at random. While the two types of models yield similar results for the contrasts with restricted maximum likelihood (REML) procedures, the CB model is generally not invariant to changes in the baseline test in a subset of studies and yields an odds ratio (OR) making it difficult to recover information on the absolute diagnostic accuracy (the marginal means), relative sensitivity or specificity of a test compared to another or differences in accuracy between tests, measures that are easily interpretable and often used in clinical epidemiology. It is common knowledge that the OR is only a good approximation of relative sensitivity/specificity when the outcome is rare but this is often not the case in diagnostic studies. Moreover, the AB model is simpler when the baseline/comparator treatment varies from one study to another or when the number of tests varies substantially among studies. By accommodating more complex variance-covariance structures AB models have been shown to be superior to CB models \cite{Zhang}.

We present two two-way ANOVA models in a diagnostic data setting. The first is a generalized linear mixed model is analogous to randomized trials with complete block designs or repeated measures in analysis of variance models where studies are equivalent to blocks. The model expresses the logit transformed sensitivity and specificity as sum of fixed effects for test, correlated study-effects and a random error associated with various tests evaluated in given study is proposed.

This generalized linear mixed model assumes that the transformed study-specific sensitivity and specificity (often logit transformed) have approximately a normal distribution with constant variance. For proportions such a model is structurally flawed because the variance depends on the underlying probability. The constant variance condition is rarely satisfied (indeed as the probability moves to 0 or 1, the variance moves towards zero and is highest when the probability is 0.5). Furthermore, the dependence between the mean and variance implies that the parameter space for the two parameters is constrained. This is in contrast with the functionally independent and unbounded variance parameters in the normal distribution. The assumption of normality for the random effects is not natural and difficult to defend empirically when there is a small number of studies. With a small number of studies, there could be problems estimating the between-study variability requiring informative priors \cite(Higgins et al. 2009) which are often difficult to specify. Alternatively, the often unrealistic assumption of homogeneity is applied \cite(Chung et al. 2013). Since the random-effects are latent, it is difficult to validate a model, yet the particular parametric model for the random-effects can possibly impact the conclusions of the meta-analysis \cite(Higgins et al. 2009) \cite(Aitkin 1999). Furthermore, the parameters have a study-specific interpretation and do not naturally describe the overall mean sensitivity and specificity and the heterogeneity thereof.

Since both sensitivity and specificity take values in the interval (0, 1), it is natural to model them using a bivariate beta without the need for a transformation. This allows for the asymmetry typically present in the distribution of variables with restricted range such as proportions and accommodates for overdispersion appropriately. The beta distribution is conjugate to the binomial distribution and therefore the random-effects are easily integrated out resulting in a beta-binomial marginal distribution. Such a modeling approach yields direct marginal effects, averaged over all the studies, as mostly of interest in meta-analyses.

The second model is built by marginal beta distributions for sensitivity and specificity being linked by a copula density. The resulting bivariate beta density describes the distribution of sensitivity and specificity jointly. The advantage of using the copula approach is that model estimation proceeds in two separate stages. First, the marginal distributions for sensitivity and specificity are estimated separately and then the dependence structure is estimated. Furthermore different copula densities can be used each resulting in a new different bivariate beta density.

This approach in the two models is efficient because the correlation structure allows the model to borrow information from the `imputed' missing data to obtain adjusted sensitivity and specificity estimates for all the tests. The two models mainly assume that results missing for some tests and studies are missing at random.

The NMADAS package

The NMADAS package is an extension of rstan [@rstan], the R interface to Stan [@stan] for network meta-analysis of diagnostic test accuracy data. Stan is a probabilistic programming language which has implemented Hamilton Monte Carlo(MHC) and uses No-U-Turn sampler (NUTS) [@Gelman]. The package facilitates easy application of complex models and their visualization within the Bayesian framework with much shorter run-times.

The NMADAS package is available via the Comprehensive R Archive Network (CRAN) at https://CRAN.R-project.org/package=NMADAS. With a working internet connection, the NMADAS package is installed and loaded in R with the following commands

#install.packages("NMADAS", dependencies = TRUE)
library(NMADAS) 

The NMADAS package provide functions to fit bivariate beta-binomial distributions constructed as a product of two beta marginal distributions and copula densities discussed earlier. The package also provides forest plots for a model with categorical covariates or with intercept only. Given the chosen copula function, a beta-binomial distribution is assembled up by the cdtamodel function which returns a cdtamodel object. The main function fit takes the nmadasmodel object and fits the model to the given dataset and returns a nmadasfit object for which print, summary and plot methods are provided for.

Demonstration dataset

To illustrate the use of the two models in network meta-analysis of diagnostic test accuracy data, we analyse data on a diversity of cytological or molecular tests to triage women with equivocal or mildly abnormal cervical cells \cite{Arbyn12, Arbyn13a, Arbyn13b, Roelens, Verdoodt}. A Pap smear is a screening test used to detect cervical precancer. When abnormalities in the Pap smear are not high grade, a triage test is needed to identify the women who need referral for further diagnostic work-up. There are several triage options, such as repetition of the Pap smear or HPV DNA or RNA assays. HPV is the virus causing cervical cancer \cite{Bosch}. Several other markers can be used for triage as well, such as p16 or the combinations of p16/Ki67 which are protein markers indicative for a transforming HPV infection~\cite{Roelens, Arbyn09} .

The data are derived from a comprehensive series of meta-analyses on the accuracy of triage with HPV assays, cervical cytology or molecular markers applied on cervical specimens in women with ASC-US (atypical squamous cells of unspecified significance) \cite{Arbyn12, Arbyn13a, Arbyn13b, Roelens, Verdoodt}. Studies were included in the analysis if they performed besides one or more triage test a verification with a reference standard based on colposcopy and biopsies.

In total, the accuracy of 10 tests in detecting intraepithelial neoplasia lesion of grade two or worse (CIN2+) were evaluated. These tests were: hrHPV DNA testing with HC2 (HC2), Conventional Cytology (CC), Liquid-Based Cytology (LBC), generic PCRs targeting hrHPV DNA (PCR) and commercially available PCR-based hrHPV DNA assays such as: Abbott RT PCR hrHPV, Linear Array, and Cobas-4800; assays detecting mRNA transcripts of five (HPV Proofer) or fourteen (APTIMA) HPV types HPV types; and protein marker (p16) identified by cytoimmunochestry p16 which is over-expressed as a consequence of HPV infection. 72 studies with at least one test and maximum of six tests were included allowing assessment of the accuracy of the ten triage tests.

This dataset is available within the package and the following commands

data(demodata)

loads the data into the R enviroment. The first six row of the dataset are

head(demodata)

study is the study identifier, TP is the number of true positives,FP is the number of false positives, TN is the number of true negatives, FN is the number of false negatives, and Test is the test identifier.

The network plot of this data is plotted as follows

networkplot.nmadas(data=demodata, S.ID="study", T.ID="Test")

The size of the nodes in figure~\ref{Fig:1} is proportional to the number of studies evaluating a test and thickness of the lines between the nodes is proportional to the number of direct comparisons between tests. The size of the node and the amount of information in a node consequently influence the standard errors of the marginal means and the relative measures. From the network plot, test 1 (HC2) and test 11 (APTIMA) were the most commonly assessed tests. The network in figure~\ref{Fig:1} is dis-connected because not all direct comparisons are present.

Statistical methods for meta-analysis

The hierarchical model

Suppose there are \textit{K} tests and \textit{I} studies. Studies assessing two tests ($ k = 2 $) are called two-arm' studies while those with $k > 2 $ aremulti-arm' studies. For a certain study \textit{i}, let $(Y_{i1k}, ~Y_{i2k})$ denote the true positives and true negatives, $(N_{i1k}, ~N_{i2k})$ the diseased and healthy individuals and $(\pi_{i1k}, ~\pi_{i2k})$ the `unobserved' sensitivity and specificity respectively with test \textit{k} in study \textit{i}. Given study-specific sensitivity and specificity, two independent binomial distributions describe the distribution of true positives and true negatives among the diseased and the healthy individuals as follows; $$ Y_{ijk} ~|~ \pi_{ijk}, ~x_i ~\sim~ bin(\pi_{ijk}, ~N_{ijk}), ~i ~=~ 1, ~\ldots I, ~j ~=~ 1, ~2, ~k~ = 1,~ \ldots ~K, $$
where $x_i$ generically denotes one or more covariates, possibly affecting $\pi_{ijk}$.

Consider a design where there is at least one test per study. The study serves as a block where all diagnostic accuracy tests are hypothetically evaluated of which some are missing. This modelling approach has potential gain in precision by borrowing strength from studies with single tests as well as multi-arm studies. The proposed single-factor design with repeated measures model is written as follows $$ logit(\pi_{ijk}) = \mu_{jk} + \eta_{ij} + \delta_{ijk} \nonumber\ \begin{pmatrix} \eta_{i1} \nonumber\ \eta_{i2} \end{pmatrix} \sim N \bigg (\begin{pmatrix} 0 \nonumber\ 0 \end{pmatrix}, \boldsymbol{\Sigma} \bigg ) \nonumber\ \boldsymbol{\Sigma} = \begin{bmatrix} \sigma^2_1 ~~~~ \rho\sigma_1\sigma_2 \nonumber\ \rho\sigma_1\sigma_2 ~~~~ \sigma^2_2 \end{bmatrix} \nonumber\ (\delta_{ij1}, \delta_{ij1}, \ldots \delta_{ijK}) \sim N(\textbf{0}, diag(\tau^2_j)) $$ where $\mu_{1k}$ and $\mu_{2k}$ are the mean sensitivity and specificity in a hypothetical study with random-effects equal to zero respectively. $\eta_{ij}$ is the study effect for healty individuals \textit{(j = 1)} or diseased individuals \textit{(j = 2)} and represents the deviation of a particular study \textit{i} from the mean sensitivity (j=1) or specificity (j=2), inducing between-study correlation. The study effects are assumed to be a random sample from a population of such effects. The between-study variability of sensitivity and specificity and the correlation thereof is captured by the parameters $\sigma_1^2$, $\sigma_2^2$, and $\rho$ respectively. $\delta_{ijk}$ is the error associated with the sensitivity (\textit{j=1}) or specificity (\textit{j=2}) of test \textit{k} in the $i^{th}$ study. Conditional on study \textit{i}, the repeated measurements are independent with variance constant across studies such that $\boldsymbol{\tau}j^2$ = $(\tau{j1}^2, \ldots, \tau_{jK}^2)$ is a \textit{K} dimensional vector of homogeneous variances.

This model is fitted by specifying marginals = "normal" in the nmamodel function as follows

modelf <-  nmamodel.nmadas(marginals = 'normal', fullcov = TRUE)

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

The model is then fitted to the data as follows

print(modelf)
str(modelf)

fitf <- fit.nmamodel(
  nma.model = modelf,
  data = demodata,
  S.ID = "study",
  T.ID = "Test")

In case $\tau_{jk}^2 = \tau_j^2$ (variances homogeneous across tests), the shared random element $\eta_{ij}$ within study \textit{i} induce a non-negative correlation between any two test results \textit{k} and $k\prime$ from healthy individuals \textit{(j = 1)} or from diseased individuals \textit{(j = 2)} equal to $\rho_j = \frac{\sigma_j^2}{\sigma_j^2 + \tau_j^2}$ (implying that a covariance matrix with compound symmetry). The compound symmetry variance-covariance is the defualt. The model with the reduced variance-covariance matrix is specified as follows

modelr <-  nmamodel(marginals = 'normal')

fitr <- fit.nmamodel(
  nma.model = modelr,
  data = demodata,
  S.ID = "study",
  T.ID = "Test")

While it might seem logical to expect and allow for similar correlation between any two sensitivities or specificities in a given study, the variances $\tau_{jk}^2$ of different sensitivities or specificities of the same study may be different. In such instances, the unstructured covariance matrix is more appropriate as it allows varying variances between the tests (in which case $\boldsymbol{\tau}j^2$ is a \textit{K} dimension vector of the unequal variances). The correlation between the $k^{th}$ and $k^{'th}$ test result is then equal to $\rho{jkk}' = \frac{\sigma_j^2}{\sqrt{\sigma_j^2 + \tau_{jk}^2~\times (\sigma_j^2 + \tau_{jk'}^2 )}}$. $\rho_j$ or $\rho_{jkk}'$ is called the intra-study correlation coefficient which also measures the proportion of the variability in $logit(\pi_{ijk})$ that is accounted for by the between study variability. It takes the value 0 when $\sigma_j^2 = 0$ (if study effects convey no information) and values close to 1 when $\sigma_j^2$ is large relative to $\tau_j^2$ and the studies are essentially all identical. When all components of $\boldsymbol{\tau}_j^2$ equal to zero, the model reduces to fitting separate bivariate random-effect meta-analysis (BRMA)~\cite{reitsma, chu} model for each test.

In essence, the model separates the variation in the studies into two components: the within-study variation $diag(\tau_j^2)$ referring to the variation in the repeated sampling of the study results if they were replicated, and the between-study variation $\boldsymbol{\Sigma}$ referring to variation in the studies true underlying effects.

Since $\pi_{ijk}$ lies within 0 and 1, a more natural, flexible and computationally more convenient distribution to model the latent variables is the beta distribution which is conjugate to the binomial distribution. The joint distribution of $(\pi_{i1k}, \pi_{i2k})$ can be easily constructed following the copula theory.

Definition of copula function

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

According to the theorem of [@Skylar], there exists for every bivariate (multivariate in extension) distribution a copula representation C which is unique for continuous random variables. If the joint cumulative distribution function and the two marginals are known, then the copula function can be written as $$C(u, ~v, ~\omega) = F(F_1^{-1} (u), ~F_2^{-1} (v), ~\omega),~ 0 \le~ u, ~v ~\le~ 1$$ A 2-dimensional copula is in fact simply a 2-dimensional cumulative function restricted to the unit square with standard uniform marginals. A comprehensive overview of copulas and their mathematical properties can be found in [@Nelsen]. To obtain the joint probability density, the joint cumulative distribution $F(x_1, x_2)$ should be differentiated to yield $$f(x_1, ~x_2) = f_1(x_1) ~f_2(x_2) ~c(F_1(x_1), ~F_2 (x_2), ~\omega)$$ where $f_1$ and $f_2$ denote the marginal density functions and c the copula density function corresponding to the copula cumulative distribution function C. Therefore from Equation~\ref{eq:3}, a bivariate probability density can be expressed using the marginal and the copula density, given that the copula function is absolutely continuous and twice differentiable.

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

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

Though the package has no capabilities to perform meta-regression it is worth mentioning an important implication when Sklar's theorem is extended to a meta-regression setting with covariates. According to [@Patton], it is important that the conditioning variable remains the same for both marginal distributions and the copula, as otherwise the joint distribution might not be properly defined. This implies that covariate information should be introduced in both the marginals and the association parameters of the model.

We now turn our attention to different bivariate beta distributions $f(\pi_{i1k},\pi_{i2k})$ constructed using different copula densities $c$. We discuss their implications and demonstrate their application in meta-analysis of diagnostic accuracy studies. The five copulas considered can model negative correlation and in some cases the positive correlation as well. The beta distributions $f(\pi_{ijk})$ are parameterized using $\mu_{jk}$, $\theta_j$ and $\delta_{jk}$ parameters as follows, $$ f_{j}(x_{ijk}) \sim beta(\mu_{jk}, \theta_j\delta_{jk}) $$ where $\mu_{jk}$ is the mean sensitivity $(j = 1)$ or specificity $(j = 2)$ of test $k$, $\theta_j$ captures the common overdispersion among the sensitivities $(j = 1)$ or specificities $(j = 2)$ in a given study as a result of repeated tests in a study, and $\delta_{jk}$ captures the test specific extra overdispersion in sensitivity $(j = 1)$ or specificity $(j = 2)$ of test $k$.

The natural $\alpha_{jk}$ and $beta_{jk}$ are recovered as follows $$ \alpha_{jk} = \mu_{jk}\times \bigg( \frac{1 - \theta_j\delta_{jk}}{\theta_j\delta_{jk}}\bigg) \ \beta_{jk} = (1 -\mu_{jk})\times \bigg( \frac{1 - \theta_j\delta_{jk}}{\theta_j\delta_{jk}}\bigg) $$

Frank copula

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

Since $\omega_k \in \mathbb{R}$, both positive and negative correlation can be modelled, making this one of the more comprehensive copulas. When $\omega_k$ is 0, sensitivity and specificity of test $k$ are independent . For $\omega_k > 0$, sensitivity and specificity exhibit positive quadrant dependence and negative quadrant dependence when $\omega_k< 0$.

The frank copula is the default and is specified and fitted as follows

frank <- nmamodel.nmadas()

fitfgm <- fit.nmamodel(
  nma.model = frank,
  data = demodata,
  S.ID = "study",
  T.ID = "Test")

Bivariate Gaussian copula

Given the density and the distribution function of the univariate and bivariate standard normal distribution with correlation parameter $\omega \in (-1, 1)$, the bivariate Gaussian copula function and density is expressed [@Meyer] as $$C(F(\pi_{i1k}), ~F(\pi_{i2k}), ~\omega_k) = \Phi_2(\Phi^{-1}(F(\pi_{i1k})),~ \Phi^{-1}(F(\pi_{i2k})),~ \omega_k), $$ $$c(F(\pi_{i1k}), ~F(\pi_{i2k}), ~\omega_k) =~ \frac{1}{\sqrt{1~-~\omega_k^2}} ~exp\bigg(\frac{2~\omega~\Phi^{-1}(F(\pi_{i1k})) ~\Phi^{-1}(F(\pi_{i2k})) - \omega_k^2~ (\Phi^{-1}(F(\pi_{i1k}))^2 + \Phi^{-1}(F(\pi_{i2k}))^2)}{2~(1 - \omega_k^2)}\bigg ) $$.

The copula is specified and fitted as follows

gauss <- nmamodel.nmadas(copula = "gaussian", marginals = "beta")

fitgauss <- fit.nmamodel(
  nma.model = gauss,
  data = demodata,
  S.ID = "study",
  T.ID = "Test")

Farlie-Gumbel-Morgenstern copula (FGM)

This popular copula by [@Farlie], [@Gumbel] and [@Morg] is defined as $$ C(F(\pi_{i1k} ), ~F(\pi_{i2k}),~\omega_k) ~= ~F(\pi_{i1k})~F(\pi_{i2k})[1 ~+~ \omega_k~(1 - F(\pi_{i1k}))~(1 ~-~ F(\pi_{i2k}))], \nonumber \
c(F(\pi_{i1k} ), ~F(\pi_{i2k}),~\omega_k) ~=~ [1 ~+~ \omega_k~(2~F(\pi_{i1k}) ~-~ 1)~(2~F(\pi_{i2k}) ~-~ 1)]. $$
for $\omega_k \in (-1, 1)$. This copula is appropriate for data with weak dependence since $|\rho_s| ~\leq~ 1/3$.

The copula is specified and fitted as follows

fgm <- nmamodel.nmadas(copula = "fgm", marginals = "beta")

fitfgm <- fit.nmamodel(
  nma.model = fgm,
  data = demodata,
  S.ID = "study",
  T.ID = "Test")

Clayton copula

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

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

The copulas are specified and fitted as follows

c90 <- nmamodel.nmadas(copula = "c90", marginals = "beta")

fitc90 <- fit.nmamodel(
  nma.model = c90,
  data = demodata,
  S.ID = "study",
  T.ID = "Test")

c270 <- nmamodel.nmadas(copula = "c270", marginals = "beta")

fitc270 <- fit.nmamodel(
  nma.model = c270,
  data = demodata,
  S.ID = "study",
  T.ID = "Test")

Of course other copula functions that allow for negative association can be chosen.

Ranking of the tests

While ranking of tests using rank probabilities and rankograms is an attractive feature of univariate NMA, it is still a challenge to rank competing diagnostic tests especially when a test does not outperform the others on both sensitivity and specificity. The superiority of a diagnostic test could be quantified using a superiority index introduced by Deutsch et al. \cite{Deutsch} expressed as $$ S_k = \frac{2a_k + c_k} {2b_k + c_k}, $$ where $a_k$ is the number of tests to which test \textit{k} is superior (higher sensitivity and specificity), $b_k$ is the number of tests to which test \textit{k} is inferior (lower sensitivity and specificity), and $c_k$ the number of tests with equal performance as test \textit{k} (equal sensitivity and specificity). \textit{S} ranges from 0 to $\infty$ with; $S$ tending to $\infty $ and $S $ tending to $0$ as the number of tests to which test \textit{k} is superior and inferior increases respectively, and $S$ tending to $1$ the more the tests are equal. Since the number of tests not comparable to test \textit{k} do not enter into the calculation of \textit{S} the index for different tests may be based on different sets of tests.

Model diagnostics

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

Model comparison and selection

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

library(httr)
mylink <- GET(url="https://www.dropbox.com/s/wm1wzehrglmqog9/NMADAS.RData?dl=1")
load(rawConnection(mylink$content))

Results

By default, chains = 3 and cores = 3 and need not be specified unless otherwise. From the code above, 28000 samples are drawn from each of the 3 chains, the first 1000 samples are discarded and thereafter every 30$^{th}$ draw kept such that each chain has 900 post-warm-up draws making a total of 2700 post-warm-up draws. The seed value, seed = 3, specifies a random number generator to allow reproducibility of the results and cores = 3 allows for parallel-processing of the chains by using 3 cores, one core for each chain. They were no initial values specified and in that case, the program randomly generates random values satisfying the parameter constraints.

References

VNyaga/NMADAS documentation built on May 6, 2019, 11:20 a.m.