meta4diag: Bayesian Bivariate Meta-analysis of Diagnostic Test Studies for Routine Practice

This paper introduces the \proglang{R} package \pkg{meta4diag} for implementing Bayesian bivariate meta-analyses of diagnostic test studies. Our package \pkg{meta4diag} is a purpose-built front end of the \proglang{R} package \pkg{INLA}. While \pkg{INLA} offers full Bayesian inference for the large set of latent Gaussian models using integrated nested Laplace approximations, \pkg{meta4diag} extracts the features needed for bivariate meta-analysis and presents them in an intuitive way. It allows the user a straightforward model-specification and offers user-specific prior distributions. Further, the newly proposed penalised complexity prior framework is supported, which builds on prior intuitions about the behaviours of the variance and correlation parameters. Accurate posterior marginal distributions for sensitivity and specificity as well as all hyperparameters, and covariates are directly obtained without Markov chain Monte Carlo sampling. Further, univariate estimates of interest, such as odds ratios, as well as the SROC curve and other common graphics are directly available for interpretation. An interactive graphical user interface provides the user with the full functionality of the package without requiring any \proglang{R} programming. The package is available through CRAN \url{https://cran.r-project.org/web/packages/meta4diag/} and its usage will be illustrated using three real data examples.


Introduction
A meta-analysis summarises the results from multiple studies with the purpose of finding a general trend across the studies. It plays a central role in several scientific areas, such as medicine, pharmacology, epidemiology, education, psychology, criminology and ecology (Borenstein et al., 2011). A bivariate meta-analysis of diagnostic test studies is a special type of meta-analysis that summarises the results from separately performed diagnostic test studies while keeping the two-dimensionality of the data (Van Houwelingen et al., 2002;Reitsma et al., 2005). Results of a diagnostic test study are commonly provided as a two-by-two table, which is a cross tabulation including four numbers: the number of patients tested positive that are indeed diseased (according to some gold standard), those tested positive that are not diseased, those tested negative that are however diseased and finally those tested negative that are indeed not diseased. Usually the table entries are referred to as true positives (TP), false positives (FP), false negatives (FN) and true negatives (TN), respectively. Those entries are used to compute pairs of sensitivity and specificity indicating the quality of the respective diagnostic test. The main goal of a bivariate meta-analysis is to derive summary estimates of sensitivity and specificity from several separately performed test studies. For this purpose pairs of sensitivity and specificity are jointly analysed and the inherent correlation between them is incorporated using a random effects approach (Reitsma et al., 2005;Chu and Cole, 2006). Related accuracy measures, such as likelihood ratios (LRs), which indicate the discriminatory performance of positive and negative tests, LR+ and LR− respectively, can be also derived. Further, frequently used estimates include diagnostics odds ratios (DORs) illustrating the effectiveness of the test or risk differences which are related to the Youden index (Altman, 1990;Youden, 1950). Reitsma et al. (2005) proposed to model logit sensitivity and logit specificity using a bivariate normal likelihood, whereby the mean vector itself is modelled using a bivariate normal distribution (normal-normal model). Our new package meta4diag follows the approach proposed by Chu and Cole (2006) and  using an exact binomial likelihood (binomial-normal model), which has been shown to outperform the approximate normal likelihood in terms of bias, mean-squared error (MSE) and coverage. Furthermore, it does not require a continuity correction for zero cells in the two-bytwo table (Harbord et al., 2007). Recently, Chen et al. (2011) and Kuss et al. (2014) proposed a third alternative, the beta-binomial model, where sensitivity and specificity are not modelled after the logit transformation but on the original scale using a beta distribution. The inherent correlation is then incorporated via copulas (Kuss et al., 2014). In the absence of covariates or in the case that all covariates affect both sensitivity and specificity (Harbord et al., 2007), the binomial-normal model can be reparameterised into the hierarchical summary receiver operating characteristic (HSROC) model (Rutter and Gatsonis, 2001;Harbord et al., 2007). In contrast to the binomial-normal model the HSROC model uses a scale parameter and an accuracy parameter, which are functions of sensitivity and specificity and defines an underlying hierarchical SROC (summary receiver operating characteristic) curve.
Different statistical software environments, such as SAS software (SAS Institute Inc., 2003), Stata (StataCorp, 2011) and R (R Core Team, 2015), have been used in the past ten years to conduct bivariate meta-analysis of diagnostic test studies. Within a frequentist setting the SAS PROC MIXED routine and PROC NLMIXED routine can be used to fit the normal-normal and binomial-normal model, see for example Van Houwelingen et al. (2002); Arends et al. (2008); Hamza et al. (2009). The SAS macro METADAS provides a user-friendly interface for the binomial-normal model and the HSROC model (Takwoingi and Deeks, 2008). Within Stata the module metandi fits the normal-normal model using an adaptive quadrature (Harbord and Whiting, 2009), while the module mvmeta performs maximum likelihood estimation of multivariate random-effects models using a Newton-Raphson procedure (White, 2009;Gasparrini et al., 2012). The R package mada (Doebler, 2015;Doebler and Holling, 2012), a specialised version of mvmeta, is specifically designed for the analysis of diagnostic accuracy. The package provides both univariate modelling of log odds ratios and bivariate binomial-normal modelling of sensitivity and specificity. A continuity correction is used for zero cells in the two-by-two tables.
Since the number of studies involved in a meta-analysis of diagnostic tests commonly is small, often less than 20 studies, and data within each two-by-two table can be sparse, the use of numerical algorithms for maximising the likelihood of the above complex bivariate model might be problematic and lead to non-convergence (Paul et al., 2010). Bayesian inference that introduces prior information for the variance and correlation parameters in the bivariate term is therefore attractive (Harbord, 2011). Markov chain Monte Carlo (MCMC) algorithms can be implemented through the generic frameworks WinBUGS (Lunn et al., 2000), OpenBUGS (Lunn et al., 2009) or JAGS (Plummer et al., 2003). There exist further specialised R-packages for analysing diagnostic test studies in Bayesian setting, such as bamdit or HSROC (Verde, 2015;Schiller and Dendukuri, 2015). Instead of modelling the link-transformed sensitivity and specificity directly, the package bamdit models the differences (D i ) and sums (S i ) of the link-transformed sensitivity and specificity jointly. The quantities D i and S i are roughly independent by using these linear transformations, so that Verde (2010) used a zero centered prior for the correlation of D i and S i to represent vague prior information. Consequently, JAGS is used for model estimation. In contrast, HSROC builds on the HSROC model to jointly analyse sensitivity and specificity with and without a gold standard reference test. Uniform priors on a restricted interval are thereby assumed for all the hyperparameters and model estimation is carried out using a Gibbs sampler (Chen and Peace, 2013, chapter 10). However, the use of Bayesian approaches is still limited in practice which might be partly caused by the fact that many applied scientists feel not comfortable with using MCMC sampling-based procedures (Harbord, 2011). Implementation needs to be performed carefully to ensure mixing and convergence. Furthermore, MCMC based methods are often time consuming, in particular, when interest lies in simulation studies which require several MCMC runs. Paul et al. (2010) proposed to perform full Bayesian inference using integrated nested Laplace approximations (INLA) which avoids MCMC entirely (Rue et al., 2009). The Rpackage INLA, see www.r-inla.org, implements Bayesian inference using INLA for the large set of latent Gaussian models. However, we understand that the range of options and the required knowledge of available features in INLA might be overwhelming for the applied user interested in only one specific model. Here, we present a new R package meta4diag which is a purpose-built package defined on top of INLA extracting only the features needed for bivariate meta regression. Our package meta4diag implements the binomial-normal model. Model definition is straightforward, and output statistics and graphics of interest are directly available. Therefore, users do not need to know the structure of the general INLA output object. Although its greatest strength, another criticism towards Bayesian inference is the choice of prior distributions. Our package meta4diag allows the user to specify prior distributions for the hyperparameters using intuitive statements based on the recently proposed framework of penalised complexity (PC) priors (Simpson et al., 2014). Alternatively, standard prior distributions or userspecific prior distributions can be used. Our package is appealing for routine use and applicable without any deep knowledge of the programming language R via the integrated Graphical User Interface (GUI) offering roll-down menus and dialog boxes implemented using the R package shiny (RStudio, Inc, 2013).
The rest of this paper is organised as follows. In Section 2 we introduce the binomialnormal model and discuss its estimation within a Bayesian inference setting. Here, specific emphasis is given on the definition of prior distributions. Section 3 illustrates the functionality of the package meta4diag. Model output and available graphics are described based on the previously analysed Telomerase (Glas et al., 2003), Scheidler (Scheidler et al., 1997) and Catheter (Chu et al., 2010) datasets. Further, the user-friendly GUI is presented. Finally, a conclusion is given in Section 4.
2 Introducing the statistical framework 2.1 Binomial-normal model for bivariate meta-analysis In a bivariate meta-analysis, each study presents the number of true positives (TP), false positives (FP), true negatives (TN), and false negatives (FN). Let Se = TP/(TP + FN) denote the true positive rate (TPR) which is known as sensitivity and Sp = TN/(TN+FP) the true negative rate (TNR) which is known as specificity. Chu and Cole (2006) proposed the following bivariate generalised linear mixed effects model to summarise the results of several diagnostic studies, i = 1, . . . , I, by modelling sensitivity and specificity jointly: Here, µ, ν denote the intercepts for logit(Se i ) and logit(Sp i ), respectively, and U i , V i study-level covariates vectors with corresponding coefficient parameters α and β. The covariance matrix of the random effects φ i and ψ i is parameterised using between-study variances σ 2 φ , σ 2 ψ and correlation ρ. The most-commonly-used logit link function can be replaced by other monotone link functions, such as the probit or the complementary log-log transformation. We assume that both sensitivity and specificity are modelled with the same link function. If desired, model (1) can easily be changed to model sensitivity and the false positive rate (1 − Sp), or the false negative rate (1−Se) and specificity, or 1−Se and 1−Sp, instead of sensitivity and specificity, causing the corresponding change in parameter estimates. Different model options are available through the argument model.type in the package meta4diag, see Section 3.3.

Specification of prior distributions
We specify prior distributions for all parameters, i.e., the three hyperparameters σ 2 φ , σ 2 ψ and ρ, as well as the fixed effects µ, ν, α and β. Per default a normal prior with zero mean and large variance is used for the fixed effects µ, ν, α and β. The user is free to specify any prior distribution for σ 2 φ , σ 2 ψ and ρ including the newly proposed penalised complexity (PC) priors, see Simpson et al. (2014) for details. One of the four principles underlying PC priors is Occam's razor. The idea is to see a certain model component as a flexible extension of a base model (commonly a simpler model) to which we would like to reduce if not otherwise indicated by the data. Thinking of a Gaussian random effect with mean zero and covariance matrix σ 2 I, the base model would be σ 2 = 0, i.e., the absence of the effect. A PC prior puts maximum density mass at the base model and decreasing mass with increasing distance away from the base model. The PC prior for the variance components σ 2 φ or σ 2 ψ is discussed in Simpson et al. (2014, Section 2.3) and corresponds to an exponential prior with parameter λ for the standard deviation σ φ or σ ψ , respectively. A simple choice to set λ is to provide (u, a) such that P (σ > u) = a leading to λ = − log(u)/a with u > 0 and 0 < a < 1. Figure 1 shows an example of the PC prior for the variance. In practice, the PC prior for the variance parameter in a diagnostic meta-analysis could be derived from the belief of the interval that sensitivities or specificities lie in. For example, choosing the contrast P(σ > 3) = 0.05 corresponds to believing that the sensitivities or specificities lie in the interval [0.5, 0.95] with probability 0.95 (Wakefield, 2007).  Figure 1: An example of the PC prior for the variance calibrated such that P (σ > 1) = 0.05. The black line is the prior density and the shaded area denotes the density weight a = 0.05 when the standard deviation is larger than u = 1.
For the correlation parameter ρ, Harbord (2011) proposed to use a stronger prior than the normal prior for the Fisher's z-transformed correlation, which was used in Paul et al. (2010). Motivated by the nature of diagnostic tests he proposed to use a prior which is not centred around zero but defined around some (negative) base value ρ 0 instead (Reitsma et al., 2005). Using the PC prior framework the above suggestions can be implemented directly. Simpson et al. (2014, Appendix A.3) derives the PC-prior for the correlation parameter in an autoregressive model of first order assuming the base model being defined at ρ 0 = 0 and identical statistical behaviour left and right of 0. Although slightly tedious, this derivation can be generalised to an arbitrary ρ 0 and asymmetrical behaviour to the left and right of ρ 0 (Guo et al., 2015). Within meta4diag() we offer three strategies to intuitively define a PC-prior for ρ given an arbitrary value of ρ 0 . Similar as for the variance, probability contrasts are used to define the prior intuitively.
Here, (ρ 0 , u 1 , a 1 , u 2 , a 2 ) are the hyperparameters needed to define the prior density. Figure 2 shows examples of the PC prior for the correlation using the three different strategies. The prior density used in Paul et al. (2010) is shown as the gray dashed lines for comparison. The parameters for the strategies are motivated based on the estimations results from Menke (2014), who analysed 50 independent bivariate meta-analyses which were selected randomly from the literature within a Bayesian setting, and Diaz (2015), who reported frequentist estimates based on a literature review of 61 bivariate metaanalyses of diagnostic accuracy published in 2010. According to these two publications, the distribution of the correlation seems asymmetric around zero. We find that around half of the correlation point estimates are negative, with a mode around −0.2. Only a small proportion are larger than 0.4 and values larger than 0.8 are rare. Based on these findings, we choose three differently behaved PC priors that are all defined around Defining the parameters of the prior distributions based on probability contrasts seems very intuitive. As illustrated it is straightforward to incorporate available prior knowledge into the prior distributions, while still have the option to define vague priors using less stringent probability contrasts. Although we recommend to specify priors for the variance and correlation components separately, our package also offers the option to use an inverse Wishart distribution as a prior for the entire covariance matrix.

Package overview
The meta4diag package provides functions for fitting bivariate meta-analyses within a full Bayesian setting as outlined in Section 2. The package is available via CRAN https: //cran.r-project.org/web/packages/meta4diag/ and can be directly installed in R by typing R> install.packages("meta4diag") (given a working internet connection and the appropriate access rights on the computer). Within this paper we use package version 2.0.5 and INLA version 0.0-1466099681. Of note, meta4diag requires INLA to be installed, which can be done using R> install.packages("INLA", repos = "http://www.math.ntnu.no/inla/R/testing").
Once the package and its dependencies are installed all analyses presented throughout this work are reproducible.
The meta4diag package consists of one major function called meta4diag(). This function estimates the Bayesian bivariate regression model for diagnostic test studies, assuming each study provides TP, FP, TN and FN. Several studies can be grouped according to a categorical variable. Posterior estimates for parameters of the bivariate model as well as common plots and summary statistics are directly available. Inference is thereby performed using INLA, which provides accurate deterministic approximations to all model parameters and linear summary estimates. Based on the output of meta4diag() different plots of interest can be generated and also non-linear summary estimates, for example the diagnostics odds ratio (DOR), are available based on Monte-Carlo estimation, whereby iid samples are generated from the approximated posterior distribution using a built-in function of INLA.
The package includes three data sets which will be used in the following subsections to illustrate the functionality of meta4diag. The data sets differ in their structure and the availability of covariates. The first data set, called Telomerase, was presented by Glas et al. (2003) and consists of 10 diagnostic test studies. There is no covariate information available. The low number of studies involved makes this data set challenging when using maximum likelihood procedures, see for example Riley et al. (2007); Paul et al. (2010). The second data set, called Scheidler, was presented in Scheidler et al. (1997) and combines three meta-analyses to compare the utility of three types of diagnostic imaging procedures to detect lymph node metastases in patients with cervical cancer. The third data set, called Catheter, consists of 33 studies from a diagnostic accuracy analysis presented by Chu et al. (2010) and provides disease prevalence as additional covariate.

General data structure required
The first argument data in the function meta4diag() is the data set. It should be given as a data frame with a minimum of 4 columns named TP, FP, TN and FN. If there is no column named studynames providing study names, the meta4diag() function will generate an additional column setting the study name indicators to study 1 , · · · , study n , where n is the number of studies in the meta-analysis. Further columns are considered to be covariates. The data set Telomerase can thus be defined using five columns, where the first column provides study name indicators and the remaining four provide values of TP, FP, TN and FN.

Analysing a standard meta-analysis without covariate information
Here, we show how to analyse the Telomerase data set which represents a meta-analysis of studies that use the telomerase maker for the analysis of bladder cancer. To analyse the dataset, we first load the INLA and the meta4diag package in R using R> library("INLA") R> library("meta4diag") We then call the function meta4diag() as follows: R> res = meta4diag(data = Telomerase, model.type = 1, + var.prior = "PC", var2.prior = "PC", cor.prior = "Normal", + var.par = c(3, 0.05), cor.par = c(0, 5), + link = "logit", nsample = 10000) The data set is transferred as the first argument followed by the argument model.type = 1, saying that we would like to model sensitivity and specificity jointly. Of note, the argument model.type can be any integer from 1 to 4 depending on which two accuracy measures are going to be modelled. When model.type = 1, sensitivity and specificity are modelled jointly. The sensitivity and (1-specificity), (1-sensitivity) and specificity and (1-sensitivity) and (1-specificity) will be jointly modelled when model.type = 2, model.type = 3 and model.type = 4, respectively. The argument var.prior is a character string to specify the prior distribution for the (transformed) variance component of the first accuracy measure, i.e., here the sensitivity. The options are "PC" for the PC prior, "Tnormal" for the truncated normal prior, "Hcauchy" for the half-Cauchy prior and "Unif" for the uniform prior, which are all defined on the standard deviation scale. Alternatively "Invgamma" for the inverse gamma prior or any user specified prior defined on the variance scale can be chosen. A user-specified prior for the variance is chosen by setting var.prior = " Table" and providing a 2-column data frame to var.par. The first column provides support points for the variance which should be in [0, ∞], and the second column provides the corresponding prior density at these points. Of note, the usage of " Table" prior in meta4diag is different from that in INLA. While INLA requires the user to define the " Table" prior on the internal parameterisation of the hyperparameter, the user of meta4diag can work on the original scale. The argument var2.prior is a character string to specify the prior distribution for the second variance component. The options are the same as for the argument var.prior.
The argument cor.prior is a character string defining the prior distribution for the (transformed) correlation parameter between the two accuracy measures. The options are "PC" for the PC prior defined on the correlation scale, "Normal" for the normal distribution defined on the Fisher's z-transformed correlation, "Beta" for the beta distribution defined on a suitable transformation, see documentation, and " Table" for an user specific prior defined on the correlation scale. The " Table" prior for the correlation should be provided as a 2-column data frame, where the first column provides suitable support points within [−1, 1], and the second column provides the corresponding density mass of those points. Alternatively, if at least one of three arguments var.prior, var2.prior and cor.prior is set to "Invwishart", an inverse Wishart distribution will be used for the covariance matrix ignoring any other prior definitions for the remaining arguments. The arguments var.par, var2.par, cor.par are numerical vectors specifying the hyperparameters for the priors for variance and correlation parameters. If the inverse Wishart prior is used the hyperparameters can be set in wishart.par. Prior definitions including parameterisations of the different options are given in the package documentation of meta4diag() or makePriors(). Of note, the arguments var.prior, var2.prior and cor.prior are not case sensitive, i.e. var.prior="pc" is valid if one uses it to indicate the PC prior for the first variance component.
Here, we use the logit link function by using link = "logit". Alternative options are "probit" for the probit link and "cloglog" for the complementary log-log transformation. The argument quantiles requires a numerical vector with values in [0, 1] defining which posterior quantiles should be returned. The default setting is c(0.025, 0.5, 0.975), and these three quantiles will always be returned. The argument nsample is an integer specifying the number of iid samples, generated from the approximated posterior distribution, which are used to compute any non-linear function of interest, such as DOR, LR+ or LR−.
To get summary information for all parameters of the model, we use the function summary() R> summary(res)  Here, also the time needed to fit the model as well as the estimated correlation between the two linear predictors, here µ and ν, is shown. This correlation is different from the hyperparameter correlation provided in cor, which corresponds toρ, i.e the posterior correlation between random effects.
To plot the posterior marginal distribution of σ 2 φ , say, we call the function plot() with argument var.type = "var1". When defining separate prior distributions for the variance and correlation parameters and setting overlay.prior = TRUE the prior distribution is shown in the same device. The posterior marginal distributions of σ 2 φ and σ 2 ψ together with their corresponding prior distribution are shown in Figure 3. Valid values of var.type are the names of the fixed effects (i.e., "mu" and "nu" for this dataset), "var1", "var2" or "rho". The argument save can be set to FALSE (default) to indicate that the resulting figures are not saved on the computer, or to a file name, (e.g., "posterior v1.pdf"), to indicate that the plot is saved as "./meta4diagPlot/posterior v1.pdf", where "./" denotes the current working directory and the directory meta4diagPlot is created automatically if it does not exist. Alternatively, the argument save can be set to TRUE to indicate that the plot is saved in the directory meta4diagPlot whereby the name is chosen according to var.type. Many standard R plotting arguments, such as xlab, ylab, xlim, ylim and col, can also be set in the plot() function.
R> fitted(res, accuracy.type = "TPR") A commonly used graphic to illustrate the results of a meta-analysis is the so-called forest plot (Lewis and Clarke, 2001). Figure 4 shows the forest plot including 95% credible intervals for the telomerase data set as obtained using the forest() function.
R> forest(res, accuracy.type = "sens", est.type = "mean", cut=c(0.4, 1), + nameShow = T, ciShow = T, dataShow = "center", + text.cex = 1.5, arrow.lwd = 1.5)  Forest plot for true positive rate (sensitivity) The arguments nameShow, dataShow, estShow require a logical value indicating whether the study names, the given observations (values of TP, FP, TN and FN) and values of credible intervals are displayed as texts in the forest plot, respectively. The corresponding texts are right aligned when the arguments are set to be TRUE. They could also be "left", "right" or "center" specifying the different alignments. The argument accuracy.type is defined as in the fitted() function. The argument est.type requires a character string specifying the summary estimate to be used. The options are "mean" (default) and "median". The arguments text.cex specifies the text size, while arrow.lwd specifies the line width for the credible lines.

Ito_1998
The two functions crosshair() and SROC() are available to study the result in the ROC space with sensitivity on the y-axis and 1-specificity on the x-axis. Figure 5 shows a crosshair plot displaying the individual studies in ROC space with paired confidence intervals representing sensitivity and specificity (Phillips et al., 2010). Figure 6 shows a summary receiver operating characteristic curve (SROC) which is only available when no separate covariates are included for the two model components, here sensitivity and specificity, as only then the bivariate meta-regression approach is equivalent to the HSROC approach (Rutter and Gatsonis, 2001). The corresponding commands are R> crosshair(res, est.type = "mean", col = c(1:10)) The argument dataShow specifies whether the original data are shown. The argument crShow and prShow are Boolean and indicate whether a credible region or prediction region, respectively, is shown. The argument sroc.type takes an integer value from 1 to 5. When sroc.type=1, the function used to define the SROC line corresponds to " The regression line 1" in Arends et al. (2008); Chappell et al. (2009). The values sroc.type=2, sroc.type=3, sroc.type=4 and sroc.type=5 correspond to "The major axis method", "The Moses and Littenberg's regression line", "The regression line 2" and "The Rutter and Gatsonis's SROC curve", respectively. Different choices may result in different SROC lines when the correlation for sensitivity and specificity is positive. We refer to Chappell et al. (2009) for more details and a comparison of the different formulations.

Incorporating additional sub-data stratification
The Scheidler data set contains the results of a meta-analysis conducted by Scheidler et al. (1997) to compare the utility of three types of diagnostic imaging, lymphangiography (LAG), computed tomography (CT) and magnetic resonance (MRI), to detect lymph node metastases in patients with cervical cancer. The dataset consists of a total of 44 studies: the first 17 for CT, the following 17 for LAG and the last 10 for MRI. The Scheidler data set is provided in the package as a data frame with length 44. It has a special column named "modality" that specifies to which imaging technology, namely CT, LAG or MRI, each study belongs to. The first five lines of the data set are given as R> data("Scheidler", package = "meta4diag") R> head (Scheidler) studynames modality TP FP FN TN 1 Grumbine_1981 CT 0 1 6 17 2 Walsh_1981 CT 12 3 3 7 3 Brenner_1982 CT 4 1 2 13 4 Villasanta_1983 CT 10 4 3 25 5 vanEngelshoven_1984 CT 3 1 4 12 6 Bandy_1985 CT 9 3 3 29 There are two obvious ways to analyse this data set. First, analyse the meta-analysis of each imaging technology separately, which gives each study its own estimates of the hyperparameters. Second, analyse all studies together and incorporate the stratification using a technology-specific intercept.
To analyse all subdata separately, we call the function meta4diag() three times assuming for each subset model (1)  Prior distributions as well as other model details, such as the link function, can be changed as described in Section 3.3.
To plot the results of all three analyses in one device, we can use the SROC() function with the argument add=TRUE, see Figure 7a.
From Figure 7a and Figure 7b, we can see that the estimated summary points are almost the same in both analyses. However, the credible regions change slightly using the different model formulations. More striking are the changes in the SROC curves, in particular for the LAG subset (blue). Looking at the data there is no obvious trend that sensitivity increases along with increasing 1−specificity. The estimated posterior correlationρ is 0.1809[−0.55, 0.79]. Chappell et al. (2009) stated that it is not appropriate to use SROC curves whenρ is close to zero or positive. Using separate analyses, we assume that each subdata has its own random effect properties. While using the full data set with a categorical covariate, we assume that all the subdata share the same covariance matrix. The choice of how to model the data is up to the user. However, when the argument covariates is used in the modelling, i.e., continuous covariates are included, the overall summary points, the confidence region and the prediction region are no longer available through the function SROC(), and only the study specified summary points can be obtained instead, see the example in Section 3.5.
The corresponding forest plot for this dataset is shown in Figure 8. The plot is automatically separated into three parts due to the column modality with three different levels.

Use of continuous covariate information
The Catheter Segment Culture data consists of 33 studies from a diagnostic accuracy analysis by Chu et al. (2010). The studies analysed semi-quantitative (19 studies) and quantitative (14 studies) catheter segment culture for the diagnosis of intravascular devicerelated blood stream infection. In the dataset a column with name type indicates whether a study is based on semi-quantitative or quantitative catheter segment culture. We consider type as a categorical covariate in the model so that it should be set to modality = "type". We choose this dataset as an example because it contains an additional column with name "prevalence" providing disease prevalence information which can be considered as a continuous covariate. To analyse the dataset, we first load the Catheter dataset R> data("Catheter", package = "meta4diag") R> head ( Forest plot for true positive rate (sensitivity) Figure 8: Forest plot for the Scheidler data. The plot is separated into three parts relating to the three sub-data sets.
Consider that we would like to use the model (3) That means we would like to model sensitivity and 1−specificity jointly as proposed by Chu et al. (2010) for this dataset. This can be done by setting model.type = 2. As the Catheter Segment Culture data contains one categorical covariate type and one continuous covariate prevalence, the argument modality is set to be "type" and argument covariates is set to be "prevalence". R> res = meta4diag(data = Catheter, model.type = 2, + var.prior = "PC", var2.prior = "PC", cor.prior = "PC", + var.par = c(3, 0.05), cor.par = c(1, -0.1, 0.5, -0.95, 0.05, NA, NA), + modality = "type", covariates = "prevalence", + quantiles = c(0.125, 0.875), nsample = 10000) Currently only one categorical covariate can be included in the model, whereas there is no limitation for the number of continuous covariates. In order to include more than one continuous covariate in the model, the user can provide a vector giving the names of all covariates to be included or the respective column numbers in the data frame.
Here, we choose a PC prior for all hyperparameters. The vector of parameters for the PC prior of the correlation parameter must always be of length 7 specifying strategy, ρ 0 , ω, u 1 , α 1 , u 2 , α 2 . However, u 2 and α 2 are not required when using strategy = 1, u 1 and α 1 are not required when strategy = 2 and there is no need to specify ω when strategy = 3, see Section 2.2. To obtain the 12.5% and 87.5% quantiles in addition to the default 2.5%, 50% and 97.5% quantiles we set quantiles = c(0.125, 0.875). Summary estimates are again obtained using the function summary() R> summary(res) A forest plot for the log diagnostic odds ratio is given in Figure 9. Here, 75% credible intervals are shown which is specified by setting the argument intervals = c(0.125, 0.875) within the function forest().

Estimates
Forest plot for log diagnostic odds ratio (ldor) Figure 9: Forest plot for the log diagnostic odds ratio (LDOR) of the Catheter data set. The study names, original dataset, estimated mean and corresponding 75% credible intervals are also shown.
Of note, when the argument covariates is available, the summary estimates cannot be returned through the function forest(). Similarly, the summary points, confidence region and prediction region in the SROC plot are not available. The SROC curve in contrast is still available. However, it does not depend on the choice of the argument sroc.type, but is computed according to Walter (2002) by fitting a regression equation respectively. After fitting the regression line, the equation of the SROC curve can be obtained as

Graphical user interface
To make Bayesian diagnostic meta-analysis easier to use for applied scientists, a crossplatform, interactive and user-friendly graphical user interface (GUI) has been implemented. The GUI can be used to load the data, set and graphically inspect the priors as the hyperparameters are manually changed by sliders (see Figure 12), and run the model. The results of the analysis are shown directly in the interface and can be saved for later use. The GUI only requires the basic knowledge of R required to start R, load the libraries and run the command that starts the GUI. Within the interface all options are visualised as buttons or drop-down menus, and help for each option is found as tooltips when the user moves the mouse over the option or the "Description area". The interface has been tested in the browsers "Internet Explorer", "Mozilla Firefox", "Google Chrome" and "Safari" on Linux, Mac and Windows 10 operating systems. The GUI is started by loading the libraries meta4diag and INLA and then calling the function meta4diagGUI() with R> library("meta4diag") R> library("INLA") R> meta4diagGUI() The start window of the GUI is shown in Figure 10 and is divided into three areas A, B and C. A contains the toolbar and has buttons for running INLA and writing the results to a text file, and buttons for starting the tutorial, saving the results to an R object for further study in R and for quitting the interface. B has 6 tabs that contain the various control panels, which are used to set up the analysis, such as the data control panel, the prior control panel, and the model control panel. The options within these three panels must be set before pressing the "RunINLA" button in A. The "Forest" control panel and "SROC" control panel in B are used for choosing plotting settings, and can be used both before and after running the model, and the "Fitted" panel allows the user to inspect the estimates for different choices of accuracy types and can also be set after running the model. Lastly, C has 10 tabs where the first is a welcome page and the rest are used to view the data and the results. Figure 11, which contains a screenshot of the "Prior Control Panel" on the left-hand side and the "Model Control Panel" on the right-hand side, gives an example of how the user can set the model and the prior. The description of the options in each panel is integrated in the GUI through tooltips, but can also be found in the package documentation (see meta4diag() for details). The left-hand side screenshot only shows how the user can showing different pages (welcome message, data set, summary results, graphics). In particular, the "Data Control Panel" is shown in the "tool control panel" area. Users could upload their own datasets for analysing or choose a example dataset for understanding the package. The "Welcome" page is shown in the "view area". The basic information for modelling and the description of bivariate meta-analysis of diagnostic test studies are shown in this page. set the prior distribution for the first variance component, but the panel also contains options for setting the priors on the second variance component and the correlation in the bivariate model. Figure 12 illustrates how the user can explore different settings of the hyperparameters interactively by sliding the sliders corresponding to each parameter. When the PC prior is selected for the correlation parameter, the user may use either of the specification strategies described in Section 2.2 to set the hyperparameters. The "Model Control Panel" shown on the right-hand side of Figure 11 is used to specify the model type, link function, quantiles of interest, and more.
After setting the options in the first three control panels and clicking the "RunINLA" button, the chosen dataset will be loaded and analysed. The results of the analysis will be shown in the view area (C) and, for example, the SROC plot can be viewed in the SROC tab of the view area (C) as shown in Figure 13. The other tabs can be used to view summary estimates, study-specific accuracy estimates, posterior marginal plots and forest plot, and in each case the R code that generated the figure or text is also shown. If the data, the model or prior settings are changed, the user must push the "RunINLA" button again to update the results. Figure 11: Details for two tool panels. Left: "Prior Control Panel". In this panel, users can set the prior distributions for the first variance component, second variance component and the correlation. In particular, the specification of the PC-prior for the first variance component is shown. A "Description area" is shown to explain what the prior is. The red "Invalid!" indicates that the given value for the hyper parameter α is not valid. The interval of the valid values can be seen from the tooltips of the indicator "Invalid". Right: "Model Control Panel".

Conclusion
The present paper demonstrates the usage of the novel R-package meta4diag for analysing bivariate meta-analyses of diagnostic test studies with R, and illustrates its usage using three examples from the literature. The package is built on top of the R-package INLA and thus provides full Bayesian inference without the need for Markov Chain Monte Carlo techniques. This is especially important when several or complex meta-analyses are studied, or simulation studies shall be performed, as then the time speed-up becomes obvious. The model can be easily specified, whereby the user does not need to know any INLA-specific details. Quantities relevant in the field of diagnostic meta-regression are internally computed and returned directly without requiring the user to work with the general and complex INLA output.
One of the biggest advantages, besides of being fast, compared to other software packages for Bayesian inference is the flexible and at the same time intuitive prior specification framework. In particular the newly proposed PC priors (Simpson et al., 2014) are supported. Here, the user has the possibility to incorporate expert knowledge in the form of probability contrasts. Guo et al. (2015) compared the performance of different PC priors with previous proposed priors in the bivariate model through an intensive simulation study and a real data set. Both informative and less informative PC priors were studied, and results indicated that the PC priors perform at least as good as previously used priors.
A graphical user interface makes the package also attractive for users who prefer to work with interactive windows offering selection menus. The GUI provides the full functionality of the package. In addition the user can inspect the priors directly and change them interactively. By offering fast inference within a Bayesian framework, intuitive choice of prior distributions and the GUI we feel that this package has great potential for routine practice. As a future research direction, we would like to expand the functionality of this package to a three-variate model analysing sensitivity, specificity and disease prevalence jointly. Further, we would like to investigate how to extend meta4diag when the assumption of a perfect reference test is not given.

Acknowledgments
We are grateful for helpful discussions with Håvard Rue and Geir-Arne Fuglstad.