§ INTRODUCTION

Before a diagnostic test is approved for being routinely used in practice, its ability to distinguish say, diseased from nondiseased individuals, must be narrowly evaluated. Throughout we assume that the true disease status of the individuals is known and the task is, compared to the truth, to quantify how accurate the test being investigated is. Before proceeding, it is worth noting that although our focus is on medical diagnosis, the problem of binary classification is a much wider one, finding applications in fields as diverse as biology and finance, to name only two.

The receiver operating characteristic (ROC) curve [Metz, 1978] is, unarguably, the most popular tool used for evaluating the discriminatory ability of continuous-outcome diagnostic tests. The ROC curve displays the false positive fraction (FPF) against the true positive fraction (TPF) for all possible threshold values used to dichotomise the test result. The ROC curve thus provides a global description of the trade-off between the FPF and the TPF of the test as the threshold changes. Plenty of parametric and semi/nonparametric methods are available for estimating ROC curves, either from frequentist or Bayesian viewpoints and we refer the interested reader to <cit.>, <cit.>, and [Gonçalves et al., 2014], and references therein.

However, it is known that in many situations, a test's outcome and, by extension, its discriminatory capacity, can be affected by additional information (covariates); <cit.> provides several examples of covariates that can affect the result of a diagnostic test. For instance, patient characteristics, such as age and gender, are important covariates to be considered as diagnostic accuracy is likely to vary according to them. In these cases pooling the test outcomes regardless of their covariate values may lead to erroneous or, at least, oversimplified conclusions and decisions. Interest should therefore be focused on assessing the accuracy of the test, but taking into account covariate information. Two different ROC-based measures that incorporate covariate information have been proposed: the covariate-specific or conditional ROC curve <cit.> and the covariate-adjusted ROC curve [Janes and Pepe, 2009]. The formal definition of both measures is given in Section <ref>. In brief, a covariate-specific ROC curve is an ROC curve that conditions on a specific covariate value, thus describing the accuracy of the test in the `subpopulation' defined by that covariate value. On the other hand, the covariate-adjusted ROC curve is a weighted average of covariate-specific ROC curves. Regarding estimation, since the seminal paper of [Pepe, 1998], a plethora of methods have been proposed in the literature for the estimation of the covariate-specific ROC curve and associated summary measures. Without being exhaustive, we mention the work of [Faraggi, 2003], [Rodríguez-Álvarez et al., 2011, Rodríguez-Álvarez et al., 2011], [Inácio de Carvalho et al., 2013], and [Inácio de Carvalho et al., 2017]. A detailed review can be found in [Rodríguez-Álvarez et al., 2011] and [Pardo-Fernández et al., 2014]. With respect to the covariate-adjusted ROC curve, estimation has been discussed in [Janes and Pepe, 2009], [Rodríguez-Álvarez et al., 2011], [Guan et al., 2012], and [Inácio de Carvalho and Rodríguez-Álvarez, 2018].

In a slightly different context, in the machine learning community, the topic of covariate-dependent classification has received much attention recently, being related to the concept of fairness <cit.>. As an example, [Buolamwini, 2017] reports that the evaluation of four gender classifiers revealed that a significant gap exists when comparing gender classification accuracies of females versus males.

There are a few R [R Core Team, 2020] packages for ROC curve analysis available on the Comprehensive R Archive Network (CRAN) and, as far as we are aware, all of them implementing frequentist approaches. The package sROC [Wang, 2012] contains functions to perform nonparametric, kernel-based, estimation of ROC curves, while pROC [Robin et al., 2011] offers a set of tools to visualise, smooth, and compare ROC curves, but covariate information cannot be explicitly taken into account in any of these packages. Packages ROCRegression (available at <https://bitbucket.org/mxrodriguez/rocregression>) and npROCRegression [Rodriguez-Alvarez and Roca-Pardinas, 2017] provide routines to estimate semiparametrically and nonparametrically, under a frequentist framework, the covariate-specific ROC curve. We also mention  OptimalCutpoints [López-Ratón et al., 2014] and  ThresholdROC [Perez Jaume et al., 2017] that provide a collection of functions for point and interval estimation of optimal thresholds for continuous diagnostic tests. To the best of our knowledge, there is no statistical software package implementing Bayesian inference for ROC curves and associated summary indices and optimal thresholds.

To close this gap, in this paper we introduce the ROCnReg package that allows conducting Bayesian inference for the (pooled or marginal) ROC curve, the covariate-specific ROC curve, and the covariate-adjusted ROC curve. For the sake of generality, frequentist approaches are also implemented. Specifically, in what concerns estimation of the pooled ROC curve, ROCnReg implements the frequentist empirical estimator described in [Hsieh and Turnbull, 1996], the kernel-based approach proposed of [Zou et al., 1997], the Bayesian Bootstrap method of [Gu et al., 2008], and the Bayesian nonparametric method based on a Dirichlet process mixture of normal distributions model proposed by [Erkanli et al., 2006]. Regarding the covariate-specific ROC curve, ROCnReg implements the frequentist normal method of [Faraggi, 2003] and its semiparametric counterpart as described in [Pepe, 1998], the kernel-based approach of [Rodríguez-Álvarez et al., 2011], and the Bayesian nonparametric model, based on a single-weights dependent Dirichlet process mixture of normal distributions, proposed by [Inácio de Carvalho et al., 2013]. As for the covariate-adjusted ROC curve, the ROCnReg package allows estimation using the frequentist semiparametric approach of [Janes and Pepe, 2009], the frequentist nonparametric method discussed in [Rodríguez-Álvarez et al., 2011], and the recently proposed Bayesian nonparametric estimator of [Inácio de Carvalho and Rodríguez-Álvarez, 2018]. Table <ref> shows a summary of all methods implemented in the package. In addition, ROCnReg also provides functions to obtain ROC-based optimal thresholds to perform the classification/diagnosis using two different criteria, namely, the Youden index and the criterion that sets a target value for the false positive fraction. These are implemented for both the ROC curve, the covariate-specific and the covariate-adjusted ROC curve. A detailed description of the methods is presented in Section <ref>.

Method     Description
2lPooled ROC curve
emp     (Frequentist) empirical estimator [Hsieh and Turnbull, 1996].

kernel     (Frequentist) kernel-based approach [Zou et al., 1997].

BB      Bayesian bootstrap method [Gu et al., 2008].

dpm      Nonparametric Bayesian approach based on Dirichlet process mixtures of normal distributions [Erkanli et al., 2006].
2lCovariate-specific ROC curve
sp     (Frequentist) parametric and semiparametric induced ROC regression approach [Pepe, 1998, Faraggi, 2003]

kernel     Nonparametric (kernel-based) induced ROC regression approach [Rodríguez-Álvarez et al., 2011].

bnp     Nonparametric Bayesian model based on a single-weights dependent Dirichlet process mixture of normal distributions [Inácio de Carvalho et al., 2013].
2lCovariate-adjusted ROC curve
sp     (Frequentist) semiparametric method [Janes and Pepe, 2009].

kernel      Nonparametric (kernel-based) induced ROC regression approach [Rodríguez-Álvarez et al., 2011].

bnp     Nonparametric Bayesian model based on a single-weights dependent Dirichlet process mixture of normal distributions and the Bayesian bootstrap [Inácio de Carvalho and Rodríguez-Álvarez, 2018].

Overview of ROC estimation methods included in the ROCnReg package.

The remainder of the paper is organised as follows. In Section <ref> we formally introduce the (pooled or marginal) ROC curve, the covariate-specific ROC curve, and the covariate-adjusted ROC curve. The description of the estimation methods implemented in the ROCnReg package is given in Section <ref>. In Section <ref> the usage of the main functions and capabilities of ROCnReg are described and illustrated using a real example. The paper concludes with a discussion in Section <ref>.

§ NOTATION AND DEFINITIONS

This section sets out the formal definition of the pooled or marginal ROC curve, the covariate-specific ROC curve, and the covariate-adjusted ROC curve. Also, it describes the most commonly used summary measures of accuracy, namely, the area under the ROC curve, the partial area under the ROC curve, and the Youden Index. For conciseness, we intentionally avoid giving too many details and refer the interested reader to [Pepe, 2003] (and references therein) for an extensive account of many aspects of ROC curves with and without covariates.

In what follows, we denote as $Y$ the outcome of the diagnostic test and as $D$ the binary variable indicating the presence ($D = 1$) or absence ($D = 0$) of disease. We also assume that along with $Y$ and the true disease status $D$, a covariate vector $\mathbf{X}$ is also available, and that it may encompass both continuous and categorical covariates. For ease of notation, the covariate vector $\mathbf{X}$ is assumed to be the same in both the diseased ($D = 1$) and nondiseased ($D = 0$) populations, although this is not necessarily the case in practice (e.g., disease stage is, obviously, a disease-specific covariate). By a slight abuse of notation, we use the subscripts $D$ and $\bar{D}$ to denote (random) quantities conditional on, respectively, $D = 1$ and $D = 0$. For example, $Y_{D}$ and $Y_{\bar{D}}$ denote the test outcomes in the diseased and nondiseased populations.

§.§ Pooled ROC curve

In the case of a continuous-outcome diagnostic test, the classification is usually made by comparing the test result $Y$ against a threshold $c$. If the outcome is equal or above the threshold, $Y \geq c$, the subject will be considered as diseased. On the other hand, if the test result is below the threshold, $Y < c$, he/she will be classified as nondiseased. The ROC curve is then defined as the set of all possible false positive fractions, $\text{FPF}\left(c\right) = \Prob(Y \geq c \mid D = 0) = \Prob(Y_{\bar{D}} \geq c)$, and true positive fractions, $\text{TPF}\left(c\right) = \Prob(Y \geq c \mid D = 1) = \Prob(Y_D \geq c)$, which can be obtained by varying the threshold value $c$, i.e.,

\[
\left\{\left(\text{FPF}\left(c\right), \text{TPF}\left(c\right)\right): c \in \mathbb{R} \right\}.
\]

It is common to represent the ROC curve as $\left\{\left(p, \text{ROC}(p)\right): p \in [0,1] \right\}$, where

\begin{equation}
\text{ROC}(p) = 1 - F_D\left\{F_{\bar{D}}^{-1}(1-p)\right\},
\label{ROC}
\end{equation}

with $F_{D}\left(y\right) = \Prob(Y_{D} \leq y)$ and $F_{\bar{D}}\left(y\right) = \Prob(Y_{\bar{D}} \leq y)$ denoting the cumulative distribution function (CDF) of $Y$ in the nondiseased and diseased groups, respectively. Several indices can be used as global summary measures of the accuracy of a test. The most widely used is the area under the ROC curve (AUC), defined as

\begin{equation}
\text{AUC} = \int_{0}^{1}\text{ROC}\left(p\right)\text{d}p.
\label{AUC1}
\end{equation}

In addition to its geometric definition, the AUC has also a probabilistic interpretation <cit.>

\begin{equation}
\text{AUC} = \Prob\left(Y_{D} \geq Y_{\bar{D}}\right),
\label{AUC2}
\end{equation}

that is, the AUC is the probability that a randomly selected diseased subject has a higher test outcome than that of a randomly selected nondiseased subject. The AUC takes values between $0.5$, in the case of an uninformative test that classifies individuals no better than chance, and $1.0$ for a perfect test. We note that an AUC below $0.5$ simply means that the classification rule should be reversed. As it is clear from its definition, the AUC integrates the ROC curve over the whole range of FFPs. Depending on the circumstances, however, interest might lie only on a relevant interval of FPFs or TPFs, which leads to the notion of partial area under the ROC curve (pAUC). The pAUC over a range of FPFs $\left(0, u_1\right)$, where $u_1$ is typically low and represents the largest acceptable FPF, is defined as

\begin{equation}
\text{pAUC}\left(u_1\right) = \int_{0}^{u_1}\text{ROC}\left(p\right)\text{d}p.
\label{pAUC}
\end{equation}

On the other hand, the pAUC over a range of TPFs $(v_1,1)$, where $v_1$ is typically large and represents the lowest acceptable TPF, is defined as

\begin{equation}
\text{pAUC}_{\text{TPF}}\left(v_1\right) = \int_{v_1}^{1}\text{ROC}_{\text{TNF}}\left(p\right)\text{d}p,
\label{pAUC_TPF}
\end{equation}

where $\text{ROC}_{\text{TNF}}$ is a $270^\circ$ rotation of the ROC curve, which can be expressed as

\begin{equation}
\text{ROC}_{\text{TNF}}(p) = F_{\bar{D}}\{F_{D}^{-1}(1-p)\}.
\label{ROC_TNF}
\end{equation}

The curve (<ref>) is referred to as the true negative fraction (TNF) ROC curve, since TNF ( $= 1 - \text{FPF}$) is plotted on the $y$-axis.
We shall highlight that the argument $p$ in the ROC curve stands for a false positive fraction, whereas in the $\text{ROC}_{\text{TNF}}$ curve it stands for a true positive fraction.

Another summary index of diagnostic accuracy is the Youden Index [Shapiro, 1999, Youden, 1950]

\begin{eqnarray}
\text{YI} & = & \max_{c}\left\{\text{TPF}(c)-\text{FPF}(c)\right\} \label{YI1}\\
& = & \max_{c}\left\{F_{\bar{D}}\left(c\right) - F_{D}\left(c\right)\right\} \label{YI2}\\
& = & \max_{p}\left\{\text{ROC}(p)-p\right\} \label{YI3}.
\end{eqnarray}

The YI ranges from $0$ to $1$, taking the value of $0$ in the case of an uninformative test and $1$ for a perfect test. As for the AUC, a YI below $0$ means that the classification rule should be reversed. The value $c^{*}$ which maximises Equation (<ref>) (or, equivalently, Equation (<ref>)) is frequently used in practice to classify subjects as diseased or nondiseased. It should be noted that this index is equivalent to the Kolmogorov–Smirnov measure of distance between the distributions of $Y_{D}$ and $Y_{\bar{D}}$ <cit.>.

§.§ Covariate-specific ROC curve

The conditional or covariate-specific ROC curve, given a covariate value $\mathbf{x}$, is defined as

\begin{eqnarray}
\text{ROC}(p\mid\mathbf{x}) = 1-F_{D}\{F_{\bar{D}}^{-1}(1-p\mid\mathbf{x})\mid\mathbf{x}\},
\label{ROCConditional}
\end{eqnarray}

where $F_{D}(y\mid\mathbf{x}) = \Prob(Y_{D}\leq y\mid \mathbf{X}_{D}= \mathbf{x})$ and $F_{\bar{D}}(y\mid\mathbf{x}) = \Prob(Y_{\bar{D}}\leq y\mid\mathbf{X}_{\bar{D}}=\mathbf{x})$ are the conditional CDFs of the test in the diseased and nondiseased groups, respectively. In this case, a number of possibly different ROC curves (and therefore accuracies) can be obtained for different values of $\mathbf{x}$. Thus, the covariate-specific ROC curve is an important tool that helps to understand and determine the optimal and suboptimal populations where to apply the tests on. Similarly to the unconditional case, the covariate-specific TNF-ROC curve is given by

\begin{eqnarray}
\text{ROC}_{\text{TNF}}(p\mid\mathbf{x}) = F_{\bar{D}}\{F_{D}^{-1}(1-p\mid\mathbf{x})\mid\mathbf{x}\},
\label{ROCConditional_TNF}
\end{eqnarray}

and the covariate-specific AUC, pAUC, and Youden index are

\begin{eqnarray}
\text{AUC}(\mathbf{x}) & = & \int_{0}^{1}\text{ROC}(p\mid\mathbf{x})\text{d}p, \label{AUCConditional}\\
\text{pAUC}(u_1 \mid \mathbf{x}) & = & \int_{0}^{u_1}\text{ROC}(p\mid\mathbf{x})\text{d}p, \label{pAUCConditional}\\
\text{pAUC}_{\text{TPF}}(v_1 \mid \mathbf{x}) & = & \int_{v_1}^{1}\text{ROC}_{\text{TNF}}(p\mid\mathbf{x})\text{d}p, \label{pAUCConditional2}\\
\text{YI}(\mathbf{x}) & = & \max_{c} \lvert \text{TPF}(c\mid\mathbf{x}) - \text{FPF}(c \mid \mathbf{x}) \rvert \label{YIConditional1}\\
& = &  \max_{c}\lvert F_{\bar{D}}(c\mid\mathbf{x}) - F_{D}(c\mid\mathbf{x}) \rvert \label{YIConditional2}\\
& = &  \max_{p}\lvert \text{ROC}(p\mid\mathbf{x})- p \rvert \label{YIConditional3}.
\end{eqnarray}

The value $c^{*}_{\mathbf{x}}$ that achieves the maximum in (<ref>) (or (<ref>)) is called the optimal covariate-specific YI threshold and can be used to classify a subject, with covariate value $\mathbf{x}$, as diseased or nondiseased.

§.§ Covariate-adjusted ROC curve

The covariate-specific ROC curve and associated AUC, pAUCs, and YI described in Section <ref> depict the accuracy of the test for specific covariate values. However, it would be undoubtedly useful to have a global summary measure that also takes covariate information into account. Such summary measure was developed by [Janes and Pepe, 2009], who proposed the covariate-adjusted ROC (AROC) curve, defined as

\begin{equation}
\text{AROC}(p)=\int \text{ROC}(p\mid\mathbf{x})\text{d}H_{D}(\mathbf{x}),
\label{aroc1}
\end{equation}

where $H_{D}(\mathbf{x}) = \Prob(\mathbf{X}_{D}\leq \mathbf{x})$ is the CDF of $\mathbf{X}_{D}$. That is, the AROC curve is a weighted average of covariate-specific ROC curves, weighted according to the distribution of the covariates in the diseased group. Equivalently, as shown by [Janes and Pepe, 2009], the AROC curve can also be expressed as

\begin{align}\label{aroc2}
\text{AROC}(p) & = \Prob\{Y_{D}>F_{\bar{D}}^{-1}(1-p\mid \mathbf{X}_{D})\} \nonumber \\
& = \Prob\{1-F_{\bar{D}}(Y_D\mid\mathbf{X}_{D})\leq p\}.
\end{align}

As will be seen in Section <ref>, Expression (<ref>) is very convenient when it comes to estimating the AROC curve. Also, it emphasises that the AROC curve at a FPF of $p$ is the overall TPF when the thresholds used for defining a positive test result are covariate-specific and chosen to ensure that the FPF is $p$ in each subpopulation defined by the covariate values.

In contrast to the pooled ROC curve (see Expression (<ref>)) and the covariate-specific ROC curve (see Expression (<ref>)), the AROC curve cannot be expressed in terms of the (conditional) CDFs of the test in each group. This does not, however, preclude the possibility of defining AROC-based summary accuracy measures, yet more care is needed. Thus, for the AROC curve, the area under the AROC, as well as the partial areas and YI are expressed as follows

\begin{eqnarray}
\text{AAUC} & = & \int_{0}^{1}\text{AROC}(p)\text{d}p, \label{AAUCConditional}\\
\text{pAAUC}(u_1) & = & \int_{0}^{u_1}\text{AROC}(p)\text{d}p, \label{pAAUCConditional}\\
\text{pAUC}_{\text{TPF}}(v_1) & = & \int_{\text{AROC}^{-1}(v_1)}^{1}\text{AROC}(p)\text{d}p - \{1-\text{AROC}^{-1}(v_1)\}v_1, \label{pAUCConditional3}\\
\text{YI}_{\text{AROC}} & = &  \max_{p}\left\{\text{AROC}(p)-p\right\} \label{YIAROC}.
\end{eqnarray}

Note, in particular, that the expressions for both the partial area under the AROC curve over a range of TPFs and for the YI are defined in terms of the AROC curve. For the YI, once the value that achieves the maximum in (<ref>) is calculated, say $p^{*}$, covariate-specific threshold values can be obtained as follows

\[
c^{*}_{\mathbf{x}} = F_{\bar{D}}^{-1}(1-p^{*}\mid \mathbf{X}_{D} = \mathbf{x}).
\]

Note that, by construction, these threshold values will ensure that the FPF is $p^{*}$ in each subpopulation defined by the covariate values; however, the TPF may vary with the covariate values, i.e.,

\[
\text{TPF}\left(c^{*}_{\mathbf{x}}\right) = 1 - F_{D}\left(c^{*}_{\mathbf{x}} \mid \mathbf{X}_{D} = \mathbf{x} \right).
\]

To finish this part, we mention that when the accuracy of a test is not affected by covariates, this does not necessarily means that the covariate-specific ROC curve (which in this case is the same for all covariate values) coincides with the pooled ROC curve. It does coincide, however, with the AROC curve <cit.>. As such, in all cases where covariates affect the results of the test, even though they might not affect its discriminatory capacity, inferences based on the pooled ROC curve might be misleading. In such cases, the AROC curve should be used instead. This also applies to the selection of (optimal) threshold values, which might be covariate-specific (i.e., possibly different for different covariate values).

§ METHODS

In this section we describe the different methods for ROC inference (with and without covariate information) implemented in the ROCnReg package.

§.§ Pooled ROC curve

Let $\{y_{\bar{D}i}\}_{i=1}^{n_{\bar{D}}}$ and $\{y_{Dj}\}_{j=1}^{n_D}$ be two independent random samples of test outcomes from the nondiseased and diseased groups of size $n_{\bar{D}}$ and $n_D$, respectively.

§.§.§ Empirical estimator

The function pooledROC.emp estimates the pooled ROC curve using the empirical estimator proposed by [Hsieh and Turnbull, 1996], which consists in estimating the CDFs of the test in each group by its empirical counterpart, that is,

\begin{equation*}
\widehat{F}_{\bar{D}}(y) = \frac{1}{n_{\bar{D}}}\sum_{i=1}^{n_{\bar{D}}}I(y_{\bar{D}i}\leq y),\qquad \widehat{F}_{D}(y) = \frac{1}{n_D}\sum_{j=1}^{n_D}I(y_{Dj}\leq y).
\end{equation*}

These empirical estimates are then plugged into Equations (<ref>) and (<ref>) to obtain, respectively, an estimate of the ROC and $\text{ROC}_{\text{TNF}}$ curves.

In what concerns estimation of the AUC (Expression (<ref>)), pAUC, and  $\text{pAUC}_{\text{TNF}}$ (Expressions (<ref>) and (<ref>), respectively) these are computed empirically by means of the Mann–Whitney U statistic. With respect to the Youden Index (and associated threshold value), it is obtained by maximising, over a grid of possible threshold values, the expression in (<ref>), with $F_D$ and $F_{\bar{D}}$ being replaced by their empirical estimators.

§.§.§ Kernel estimator

The function pooledROC.kernel estimates the pooled ROC curve using the kernel-based estimator proposed by [Zou et al., 1997] and [Zou et al., 1998], which is based on estimating the CDFs of the test as follows

\begin{equation*}
\widehat{F}_{\bar{D}}(y)=\frac{1}{n_{\bar{D}}}\sum_{i=1}^{n_{\bar{D}}}\Phi\left(\frac{y-y_{\bar{D}i}}{h_{\bar{D}}}\right),\qquad \widehat{F}_{D}(y)=\frac{1}{n_D}\sum_{j=1}^{n_D}\Phi\left(\frac{y-y_{Dj}}{h_D}\right),
\end{equation*}

where $\Phi(y)$ stands for the standard normal distribution evaluated at $y$. For the bandwidths, $h_{\bar{D}}$ and $h_D$, which control the amount of smoothing, two options are popular. Silverman's rule of thumb <cit.>, which sets the bandwidth as

\begin{equation*}
h_d=0.9\min\{\text{SD}(\textbf{y}_d), \text{IQR}(\textbf{y}_d)/1.34\}n_d^{-0.2},\quad d\in\{\bar{D},D\},
\end{equation*}

where $\text{SD}(\textbf{y}_d)$ and $\text{IQR}(\textbf{y}_d)$ are the standard deviation and interquantile range, respectively, of $\textbf{y}_d=(y_{d1},\ldots,y_{dn_{d}})$. Another alternative criterion is to select the bandwidth by using least squares cross-validation <cit.>.

Here, both the AUC, pAUC, and  $\text{pAUC}_{\text{TNF}}$ (Expressions (<ref>), (<ref>), and (<ref>)) are computed numerically using Simpson's rule. Regarding the Youden Index (and associated threshold value), it is obtained by maximising, over a grid of possible threshold values, Expression (<ref>), with $F_D$ and $F_{\bar{D}}$ being replaced by their kernel estimators.

Uncertainty estimation for both the empirical and kernel estimators is conducted through bootstrap resampling.

§.§.§ Bayesian bootstrap estimator

The function pooledROC.bb implements the Bayesian bootstrap (BB) approach proposed by [Gu et al., 2008]. Their estimator relies on the notion of placement value <cit.>, which is simply a standardisation of the test outcomes with respect to a reference group. Specifically, $U_D= 1-F_{\bar{D}}(Y_D)$ is to be interpreted as a standardisation of a diseased test outcome with respect to the distribution of test results in the nondiseased population. The ROC curve can be regarded as the CDF of $U_D$

\begin{equation}\label{rocbb}
\Prob(U_D\leq p) = \Prob\{1-F_{\bar{D}}(Y_D)\leq p\}=1-F_{D}\{F_{\bar{D}}^{-1}(1-p)\}=\text{ROC}(p),\quad 0\leq p \leq 1.
\end{equation}

This representation in (<ref>) of the ROC curve provided the rationale for the two-step algorithm of [Gu et al., 2008], which can be described as follows. Let $S$ be the number of iterations.

Step 1: Computation of the placement value  based on the BB.

For $s=1,\ldots,S$, let

\begin{equation*}
U_{Dj}^{(s)}=\sum_{i=1}^{n_{\bar{D}}}q_{1i}^{(s)}I\left(y_{\bar{D}i}\geq y_{Dj}\right), \quad j=1,\ldots,n_{D},
\end{equation*}

where $\left(q_{11}^{(s)},\ldots,q_{1n_{\bar{D}}}^{(s)}\right)\sim\text{Dirichlet}(n_{\bar{D}};1,\ldots,1)$.

Step 2: Generate a realisation of the ROC curve.
Based on (<ref>), generate a realisation of $\text{ROC}^{(s)}(p)$, the cumulative distribution function of $(U_{D1}^{(s)},\ldots,U_{Dn_D}^{(s)})$, where

\begin{equation*}
\text{ROC}^{(s)}(p)=\sum_{j=1}^{n_D}q_{2j}^{(s)}I\left(U_{Dj}^{(s)}\leq p\right),\quad \left(q_{21}^{(s)},\ldots,q_{2n_{D}}^{(s)}\right)\sim\text{Dirichlet}(n_D;1,\ldots,1).
\end{equation*}

The $S$ posterior samples give rise to an ensemble of ROC curves $\{\text{ROC}^{(1)}(p),\ldots,\text{ROC}^{(S)}(p)\}$ from which the posterior mean (or median) can be computed, e.g.,

\begin{equation*}
\widehat{\text{ROC}}^{\text{BB}}(p)=\frac{1}{S}\sum_{s=1}^{S}\text{ROC}^{(s)}(p),
\end{equation*}

and a $95\%$ pointwise credible band can be obtained from the $2.5\%$ and $97.5\%$ percentiles of the same ensemble.

The Bayesian bootstrap estimator leads to closed-form expressions for the AUC and pAUC, which are, respectively, given by

\begin{align*}
\text{AUC}^{(s)} &= \int_0^1\text{ROC}^{(s)}(p)\text{d}p = 1-\sum_{j=1}^{n_D}q_{2j}^{(s)}U_{Dj}^{(s)},\\
\text{pAUC}^{(s)}(u_1)&=\int_{0}^{u_1}\text{ROC}^{(s)}(p)\text{d}p = u_1-\sum_{j=1}^{n_D}q_{2j}^{(s)}\min\left\{u_1,U_{Dj}^{(s)}\right\}.
\end{align*}

It is easy to show that

\[
\text{pAUC}^{(s)}_{\text{TPF}}(v_1) = \int_{v_1}^{1}\text{ROC}^{(s)}_{\text{TNF}}\left(p\right)\text{d}p = \sum_{i=1}^{n_{\bar{D}}}q_{1i}^{(s)}\max\left\{v_1,U_{\bar{D}i}^{(s)}\right\} - v_1,
\]

\begin{equation*}
U_{\bar{D}i}^{(s)}=\sum_{j=1}^{n_D}q_{2j}^{(s)}I\left(y_{Dj}\geq y_{\bar{D}i}\right), \quad i=1,\ldots,n_{\bar{D}},
\end{equation*}

and it is also easy to demonstrate that the $\text{ROC}_{\text{TNF}}$ curve is the survival function of the placement value $U_{\bar{D}}=1-F_{D}(Y_{\bar{D}})$. With respect to the Youden Index, it is obtained by maximising, over a grid of possible threshold values, the following expression

\[
\text{YI}^{(s)} = \max_{c}\left\{F^{(s)}_{\bar{D}}\left(c\right) - F^{(s)}_{D}\left(c\right)\right\},
\]

\[
F^{(s)}_{\bar{D}}\left(c\right) = \sum_{i=1}^{n_{\bar{D}}}q_{1i}^{(s)}I\left(y_{\bar{D}i} \leq c\right) \;\;\; \mbox{and} \;\;\; F^{(s)}_{D}\left(c\right) = \sum_{j=1}^{n_{D}}q_{2j}^{(s)}I\left(y_{Dj} \leq c\right).
\]

As for the ROC curve, point estimates for the AUC, pAUC, $\text{pAUC}_{\text{TPF}}$, YI, and $c^{*}$ can be obtained by averaging over the respective ensembles of $S$ realisations, with credible bands derived from the percentiles of the same ensembles.

§.§.§ Dirichlet process mixture of normal distributions estimator

The Bayesian nonparametric approach, based on a Dirichlet process mixture (DPM) of normal distributions, for estimating the pooled ROC curve [Erkanli et al., 2006] is implemented in the pooledROC.dpm function. In this case, as implicit by the name, the CDFs of the test in each group are estimated via  a Dirichlet process mixture of normal distributions, that is, it is assumed that the CDF, say in the diseased group (the one in the nondiseased group, $\bar{D}$, follows analogously), is of the form

\begin{equation}\label{cdfdpm1}
F_{D}(y)=\int \Phi(y\mid \mu,\sigma^2)\text{d}G_{D}(\mu,\sigma^2), \qquad G_D\sim\text{DP}(\alpha_D,G_D^{*}(\mu,\sigma^2)),
\end{equation}

where $\Phi(y\mid \mu,\sigma^2)$ denotes the CDF of the normal distribution with mean $\mu$ and variance $\sigma^2$ evaluated at $y$. Here $G_D\sim\text{DP}(\alpha_D,G_D^{*})$ is used to denote that the mixing distribution $G_D$ follows a Dirichlet process (DP) [Ferguson, 1973] with centring distribution $G_D^{*}$, for which $E(G_D)=G_D^{*}$, and precision parameter $\alpha_D$. The centring distribution $G_D^{*}$ encapsulates any prior knowledge that might be known about $G_D$. Larger values of $\alpha_D$ lead to realisations of $G_D$ closer to $G_D^{*}$, while smaller values lead to realisations of $G_D$ with substantial variation around $G_D^{*}$. Usually, due to conjugacy reasons, $G_D^{*}(\mu,\sigma^2)\equiv \text{N}(\mu\mid m_{0D},S_{0D})\Gamma(\sigma^{-2}\mid a_D,b_D)$.

For ease of posterior simulation and because it provides a highly accurate approximation, we make use of the truncated stick-breaking representation of the DP [Ishwaran and James, 2001], according to which $G_D$ can be written as

\begin{equation*}
\end{equation*}

where $(\mu_{Dl},\sigma^2_{Dl})\overset{\text{iid}}\sim G_D^{*}(\mu,\sigma^2)$, for $l=1,\ldots,L_D$, and the weights follow the so-called (truncated) stick-breaking construction: $\omega_{D1}=v_{D1}$, $\omega_{Dl}=v_{Dl}\prod_{r<l}(1-v_{Dr})$, $l=2,\ldots,L_D$, and $v_{D1},\ldots,v_{D,L_{D}-1}\overset{\text{iid}}\sim\text{Beta}(1,\alpha_D)$. Further, one must set $v_{DL_{D}}=1$ in order to ensure that the weights add up to one. Further, one must set $v_{DL_{D}}=1$ in order to ensure that the weights add up to one. With regard to the parameter $\alpha_D$, in ROCnReg a prior distribution is placed on it. In particular, and due to conjugacy reasons, a gamma distribution is considered, i.e., $\alpha_D\sim\Gamma(a_{\alpha_{D}},b_{\alpha_{D}})$. The CDF in (<ref>) can therefore be written as

\begin{equation*}
\end{equation*}

where we shall note that $L_D$ is not the exact number of components expected to be observed, but rather an upper bound on it. Some comments are in order regarding how to specify $\alpha_D$ (or $a_{\alpha_{D}}$ and $b_{\alpha_{D}}$) and how to set $L_D$. Firstly, note that in ROCnReg a prior distribution is placed on  $\alpha_D$. In particular, and due to conjugacy reasons, a gamma distribution is considered, i.e., $\alpha_D\sim\Gamma(a_{\alpha_{D}},b_{\alpha_{D}})$.  This parameter is intrinsically related to the number of occupied mixture components, say $L_D^{*}$. Using the results shown by [Liu, 1996] we have that for moderate to large sample sizes, the conditional prior mean and variance of the number of occupied components given a fixed $\alpha_D$ and sample size $n_D$ are, respectively,

\begin{equation*}
E(L_{D}^*\mid\alpha_D)=\alpha_D\log\left(\frac{\alpha_D+n_D}{\alpha_D}\right),\quad \text{var}(L_{D}^*\mid\alpha_D)=\alpha_D\left\{\log\left(\frac{\alpha_D+n_D}{\alpha_D}\right)-1\right\}.
\end{equation*}

These expressions can be averaged over the $\Gamma(a_{\alpha_{D}},b_{\alpha_{D}})$ prior for $\alpha_D$ to obtain $E(L_{D}^*)$ and $\text{var}(L_{D}^*)$, thus selecting $a_{\alpha_{D}}$ and $b_{\alpha_{D}}$ to agree with a prior guess at the expected and variance number of distinct mixture components.  In practice, and in the absence of knowledge about the number of occupied components, it is common to use hyperparameters' values that encourage, a priori, a small number of occupied components (e.g., $a_{\alpha_{D}}=1$ and $b_{\alpha_{D}}=1$ or $a_{\alpha_{D}}=2$ and $b_{\alpha_{D}}=2$). The results derived by [Liu, 1996] can also be used to guide the selection of $L_D$. It might be reasonable to set $L_D>E(L_{D}^*\mid\alpha_D)+2\sqrt{\text{var}(L_{D}^*\mid\alpha_D)}$.

Because the full conditional distributions for all model parameters are available in closed-form, posterior simulation can be easily conducted through Gibbs sampler (see the details, for instance, in Ishwaran and James, 2002). At iteration $s$ of the Gibbs sampler procedure, the ROC curve is computed as

\begin{equation*}
\text{ROC}^{(s)}(p)=1-F_D^{(s)}\left\{F_{\bar{D}}^{-1(s)}(1-p)\right\},  \quad s=1,\ldots, S,
\end{equation*}

\begin{equation} \label{cdfdpm_2}
\end{equation}

and where the inversion is performed numerically. There is a closed-form expression for the AUC [Erkanli et al., 2006] given by

\begin{equation*}
\text{AUC}^{(s)}=\sum_{k=1}^{L_{\bar{D}}}\sum_{l=1}^{L_D}\omega_{\bar{D}k}^{(s)}\omega_{Dl}^{(s)}\Phi\left(\frac{b_{kl}^{(s)}}{\sqrt{1+a_{kl}^{2(s)}}}\right),\quad b_{kl}^{(s)}=\frac{\mu_{Dl}^{(s)}-\mu_{\bar{D}k}^{(s)}}{\sigma_{Dl}^{(s)}},\quad a_{kl}^{(s)}=\frac{\sigma_{\bar{D}k}^{(s)}}{\sigma_{Dl}^{(s)}}.
\end{equation*}

Also, when $L_D = L_{\bar{D}}$ = 1, there are closed-form expressions for the pAUC and $\text{pAUC}_{\text{TPF}}$ which are used in the package <cit.>. For the pAUC/$\text{pAUC}_{\text{TPF}}$, when $L_D > 1$ or $L_{\bar{D}} > 1$, the integrals are approximated numerically using Simpson's rule. The Youden index/optimal threshold is computed as in the Bayesian bootstrap method, with the obvious difference that here the CDFs are expressed as in (<ref>). At the end of the sampling procedure, we have an ensemble of $S$ ROC curves and AUCs/pAUCs/$\text{pAUC}_{\text{TPF}}\text{s}$/YIs/optimal thresholds which, as before, allows obtaining point and interval estimates.

§.§ Covariate-specific ROC curve

We now let $\{(\mathbf{x}_{\bar{D}i},y_{\bar{D}i})\}_{i=1}^{n_{\bar{D}}}$ and $\{(\mathbf{x}_{Dj},y_{Dj})\}_{j=1}^{n_D}$ be two independent random samples of test outcomes and covariates from the nondiseased and diseased groups of size $n_{\bar{D}}$ and $n_D$, respectively. Further, for all $i = 1,\ldots,n_{\bar{D}}$ and $j = 1,\ldots,n_D$, let $\mathbf{x}_{\bar{D}i}=(x_{\bar{D}i,1},\ldots, x_{\bar{D}i,q})^{\top}$ and $\mathbf{x}_{Dj}=(x_{Dj,1},\ldots, x_{Dj,q})^{\top}$ be $q$-dimensional vectors of covariates, which can be either continuous or categorical.

§.§.§ Induced semiparametric linear model

The function cROC.sp implements the induced ROC approaches proposed by [Faraggi, 2003] and [Pepe, 1998]. Both authors assume a location-scale regression model of the following form for the test outcomes in each group

\begin{equation}
Y_{\bar{D}}=\tilde{\mathbf{X}}_{\bar{D}}^{\top}\boldsymbol{\beta}_{\bar{D}}+\sigma_{\bar{D}}\varepsilon_{\bar{D}},\qquad Y_D=\tilde{\mathbf{X}}_D^{\top}\boldsymbol{\beta}_D+\sigma_D\varepsilon_D,
\label{loc-scale-sp}
\end{equation}

where $\tilde{\mathbf{X}}_{\bar{D}}^{\top}=(1,\mathbf{X}_{\bar{D}}^{\top})$ and $\boldsymbol{\beta}_{\bar{D}}=(\beta_{\bar{D}0},\ldots,\boldsymbol{\beta}_{\bar{D}q})^{\top}$ is a $(q+1)$-dimensional vector of (unknown) regression coefficients; $\tilde{\mathbf{X}}_{D}$ and $\boldsymbol{\beta}_{D}$ are analogously defined. The error terms $\varepsilon_{\bar{D}}$ and $\varepsilon_{D}$ have mean zero, variance one, are independent of each other and of the covariate, and have distribution functions given by $F_{\varepsilon_{\bar{D}}}$ and $F_{\varepsilon_{D}}$, respectively. Under these assumptions, we have

\begin{equation}
F_{\bar{D}}(y\mid\mathbf{x}) = F_{\varepsilon_{\bar{D}}}\left(\frac{y - \tilde{\mathbf{x}}^{\top}\boldsymbol{\beta}_{\bar{D}}}{\sigma_{\bar{D}}}\right)\;\;\;\mbox{and}\;\;\;F_{D}(y\mid\mathbf{x}) = F_{\varepsilon_{D}}\left(\frac{y - \tilde{\mathbf{x}}^{\top}\boldsymbol{\beta}_{D}}{\sigma_D}\right),
\label{cROC_YI_sp}
\end{equation}

with $\tilde{\mathbf{x}}^{\top}=(1,\mathbf{x}^{\top})$.

The approaches of [Faraggi, 2003] and [Pepe, 1998] differ in the assumptions made about the error terms. More concretely, [Faraggi, 2003]'s method assumes that the error term in both groups follows a standard normal distribution, i.e., $F_{\varepsilon_{\bar{D}}}(y)=F_{\varepsilon_{D}}(y)=\Phi(y)$, and can be summarised by the following three steps:

* Estimate the regression coefficients $\boldsymbol{\beta}_{\bar{D}}$ and $\boldsymbol{\beta}_D$ by ordinary least squares, on the basis of the samples $\{(\mathbf{x}_{\bar{D}i},y_{\bar{D}i})\}_{i=1}^{n_{\bar{D}}}$ and $\{(\mathbf{x}_{Dj},y_{Dj})\}_{j=1}^{n_D}$, respectively.

* Estimate $\widehat{\sigma}_{D}^{2}$ as

\begin{equation*}
\widehat{\sigma}_{D}^{2}=\frac{\sum_{j=1}^{n_D}\left(y_{Dj}-\tilde{\mathbf{x}}_{Dj}^{\top}\widehat{\boldsymbol{\beta}}_D\right)^{2}}{n_D-q-1},
\end{equation*}

with $\widehat{\sigma}_{\bar{D}}^{2}$ similarly estimated.

* For a given covariate vector $\mathbf{x}$, compute the covariate-specific ROC curve as follows

\begin{equation}
\widehat{\text{ROC}}(p\mid\mathbf{x})=1-\Phi\left\{a(\mathbf{x})+b\Phi^{-1}(1-p)\right\},
\label{cROC-a-b-sp}
\end{equation}

\begin{equation}
a(\mathbf{x}) =\tilde{\mathbf{x}}^{\top}\frac{(\widehat{\boldsymbol{\beta}}_{\bar{D}}-\widehat{\boldsymbol{\beta}}_D)}{\widehat{\sigma}_D},\quad \text{and}\quad b =\frac{\widehat{\sigma}_{\bar{D}}}{\widehat{\sigma}_D}.
\label{a-b-sp}
\end{equation}

Regarding the covariate-specific AUC, pAUC, and $\text{pAUC}_{\text{TPF}}$ (Expressions (<ref>), (<ref>), and (<ref>)) they admit closed-form expressions <cit.>.

As an alternative, [Pepe, 1998] suggests to estimate the CDF of the errors in each group by the corresponding empirical CDF of the estimated standardised residuals. Therefore, the first two steps of the estimation procedure remain the same, but now we have the following extra step

\begin{equation*}
\widehat{F}_{\varepsilon_{D}}(y)=\frac{1}{n_D}\sum_{j=1}^{n_D}I(\widehat{\varepsilon}_{Dj}\leq y), \qquad \widehat{\varepsilon}_{Dj}=\frac{y_{Dj}-\tilde{\mathbf{x}}_{Dj}^{\top}\widehat{\boldsymbol{\beta}}_D}{\widehat{\sigma}_D}.
\end{equation*}

The empirical CDF of the standardised residuals in the nondiseased group is estimated in a similar fashion. The covariate-specific ROC curve is finally computed in an analogous way as for the method of [Faraggi, 2003] as

\begin{equation*}
\widehat{\text{ROC}}(p\mid\mathbf{x})=1-\widehat{F}_{\varepsilon_{D}}\left\{a(\mathbf{x})+b\widehat{F}_{\varepsilon_{\bar{D}}}^{-1}(1-p)\right\}.
\end{equation*}

Here, the covariate-specific AUC and pAUC (Expressions (<ref>) and (<ref>)) also admit closed forms. However, especially for large datasets, their calculation can be very time-consuming. As a consequence, in ROCnReg they are computed numerically using Simpson's rule; in our experience Simpson's rule provides almost identical results to the ones obtained using the closed-form expressions. In what concerns the covariate-specific $\text{pAUC}_{\text{TNF}}$, it is interesting to note that

\[
\widehat{\text{ROC}}_{\text{TNF}}(p\mid\mathbf{x}) = 1-\widehat{F}_{\varepsilon_{D}}\left\{a^{*}(\mathbf{x})+b^{*}\widehat{F}_{\varepsilon_{\bar{D}}}^{-1}(1-p)\right\},
\]

\[
a^{*}(\mathbf{x}) =\tilde{\mathbf{x}}^{\top}\frac{(\widehat{\boldsymbol{\beta}}_{D}-\widehat{\boldsymbol{\beta}}_{\bar{D}})}{\widehat{\sigma}_{\bar{D}}}\;\;\;\mbox{and}\;\;\; b^{*} =\frac{\widehat{\sigma}_D}{\widehat{\sigma}_{\bar{D}}}.
\]

The covariate-specific $\text{pAUC}_{\text{TNF}}$ (Expression (<ref>)) is then computed numerically using Simpson's rule based on the previous expressions.

Finally, in pretty much the same way as for the pooled ROC curve, the covariate-specific Youden Index (and associated threshold value) is obtained by maximising, over a grid of possible threshold values, the expression in (<ref>), making use of result (<ref>).

§.§.§ Induced kernel-based approach

The kernel-based approach of [González-Manteiga et al., 2011] and [Rodríguez-Álvarez et al., 2011] is available in the cROC.kernel. Differently to all the other estimating approaches for the covariate-specific ROC curve presented in this section, it can only deal with one continuous covariate. Similarly to the approaches of [Pepe, 1998] and [Faraggi, 2003], it also assumes a location-scale regression model for the test outcomes in each group, but the effect of the covariate is not assumed to be linear and the variance is allowed to depend on the covariate. Specifically, the models postulated in each group are as follows

\begin{equation}\label{locscaleker}
Y_{\bar{D}}=\mu_{\bar{D}}(X_{\bar{D}})+\sigma_{\bar{D}}(X_{\bar{D}})\varepsilon_{\bar{D}},\qquad Y_D=\mu_{D}(X_{D})+\sigma_D(X_{D})\varepsilon_D,
\end{equation}

where $\mu_{D}(x)=\E(Y_D\mid X_D=x)$ and $\sigma_{D}^2(x) = \VAR(Y_D\mid X_D=x)$ are the regression and variance functions, respectively, with $\mu_{\bar{D}}(x)$ and $\sigma_{\bar{D}}^2(x)$ being analogously defined. The error terms $\varepsilon_{\bar{D}}$ and $\varepsilon_{D}$ have mean zero, variance one, are independent of each other and of the covariate, and have distribution functions given by $F_{\varepsilon_{\bar{D}}}$ and $F_{\varepsilon_{D}}$, respectively.

Both the regression and variance functions are estimated using local polynomial kernel smoothers [Fan and Gijbels, 1996]. In particular, local constant (Nadaraya–Watson) or local linear estimators are employed for the regression function, whereas for the variance function only local constant estimators are used. Estimation in ROCnReg makes use of the R package np by [Hayfield and Racine, 2008]. We note that estimation proceeds in a sequential manner: 1) the regression function, say in the diseased group and denoted by $\widehat{\mu}_{D}$, is estimated first on the basis of $\{(x_{Dj},y_{Dj})\}_{j=1}^{n_{D}}$, and 2) the variance function is estimated next on the basis of the sample $\{(x_{Dj},[y_{Dj}-\widehat{\mu}_{D}(x_{Dj})]^{2})\}$. Both steps involve the selection of a bandwidth parameter which is done via cross-validation. As in the model of [Pepe, 1998], the CDFs $F_{\varepsilon_{D}}$ and $F_{\varepsilon_{\bar{D}}}$ are estimated via the empirical CDF of the standardised residuals, that is,

\begin{equation*}
\widehat{F}_{\varepsilon_{D}}(y)=\frac{1}{n_D}\sum_{j=1}^{n_D}I(\widehat{\varepsilon}_{Dj}\leq y), \qquad \widehat{\varepsilon}_{Dj}=\frac{y_{Dj}-\widehat{\mu}_D(x_{Dj})}{\widehat{\sigma}_D(x_{Dj})}.
\end{equation*}

with the empirical CDF of the standardised residuals in the nondiseased group estimated analogously. Finally, the covariate-specific ROC curve is computed in an analogous way as before as

\begin{equation*}
\widehat{\text{ROC}}(p\mid x)=1-\widehat{F}_{\varepsilon_{D}}\left\{a(x)+b(x)\widehat{F}_{\varepsilon_{\bar{D}}}^{-1}(1-p)\right\},
\end{equation*}

\[
a(x) = \frac{\widehat{\mu}_{\bar{D}}(x) - \widehat{\mu}_{D}(x)}{\widehat{\sigma}_{D}(x)}\;\;\;\mbox{and}\;\;\; b(x) = \frac{\widehat{\sigma}_{\bar{D}}(x)}{\widehat{\sigma}_{D}(x)}.
\]

Estimation of the covariate-specific AUC, pAAUC, $\text{pAUC}_{\text{TNF}}$, and YI follows a similar reasoning as the one described previously for the induced semiparametric linear model when no assumptions are made regarding the distribution of the error terms.

For both the induced semiparametric linear model and the induced kernel approach, uncertainty quantification is done through a bootstrap of the residuals. For further details see, for instance, [Rodríguez-Álvarez et al., 2011].

§.§.§ Bayesian nonparametric approach based on a dependent Dirichlet process mixture of normal distributions

The Bayesian nonparametric approach for conducting inference about the covariate-specific ROC curve of [Inácio de Carvalho et al., 2013], which is based on a single-weights dependent Dirichlet process mixture of normal distributions, is implemented in the function cROC.bnp. By opposition to the previously described approaches to ROC regression, this method rests on directly modelling the CDF of test outcomes separately in the diseased and nondiseased groups. In a single-weights dependent Dirichlet process mixture of normals model [De Iorio et al., 2009], the conditional CDF in the diseased group is modelled as follows

\begin{equation*}
F_{D}(y_{Dj}\mid\mathbf{X}_D=\mathbf{x}_{Dj})=\int \Phi(y_{Dj}\mid \mu_{D}(\mathbf{x}_{Dj},\boldsymbol{\beta}),\sigma^2)
\text{d}G_D(\boldsymbol{\beta},\sigma^2),\qquad G_D\sim\text{DP}(\alpha_D,G_{D}^{*}(\boldsymbol{\beta},\sigma^2)),
\end{equation*}

with the conditional CDF in the nondiseased group, $\bar{D}$, following in an analogous manner. As in the no-covariate case, by making use of Sethuraman's truncated representation of the DP, we can write the conditional CDF as

\begin{align*}
&F_{D}(y_{Dj}\mid\mathbf{x}_{Dj}) = \sum_{l=1}^{L_D}\omega_{Dl}\Phi(y_{Dj}\mid \mu_{Dl}(\mathbf{x}_{Dj}),\sigma_{Dl}^2),\\
&\omega_{D1}=v_{D1},\quad \omega_{Dl}=v_{Dl}\prod_{r<l}(1-v_{Dr}),\quad l=2,\ldots, L_{D},\\
&v_{Dl}\overset{\text{iid}}\sim\text{Beta}(1,\alpha_D),\quad l = 1,\ldots,L_{D}-1,\quad v_{DL_{D}}=1.
\end{align*}

It is worth mentioning that although the variance of each component does not depend on covariates, the overall variance of the mixture does depend on covariates  as it can be written as

\begin{equation*}
\text{var}(y_{Dj}\mid\mathbf{x}_{Dj})=\sum_{l=1}^{L_D}\omega_{Dl}\sigma_{Dl}^{2}+\sum_{l=1}^{L_D}\omega_{Dl}\left\{\mu_{Dl}^{2}(\mathbf{x}_{Dj})-\left(\sum_{l=1}^{L_D}\omega_{Dl}\mu_{Dl}(\mathbf{x}_{Dj})\right)^2\right\}.
\end{equation*}

Note that by assuming that the weights, $w_{Dl}$, do not vary with covariates, the model might has limited flexibility in practice [MacEachern, 2000]. This issue can, however, be largely mitigated by using a flexible formulation for $\mu_{Dl}(\mathbf{x}_{Dj})$, which is needed not only for the model to be able to recover nonlinear trends, but also to recover flexible shapes that might arise due to a dependence of the weights on the covariates. As such, the ROCnReg package allows to model the mean function of each component using an additive smooth structure

\begin{equation}\label{additive}
\mu_{Dl}(\mathbf{x}_{Dj})=\beta_{Dl0}+f_{l1}(x_{Dj,1})+\ldots+f_{lq}(x_{Dj,q}),\qquad l=1,\ldots,L_D,
\end{equation}

where the smooth functions, $f_{lm}$ $(m = 1, \ldots, q)$, are approximated using a linear combination of B-splines basis functions. To avoid notational burden we have assumed that all $q$ covariates are continuous and modelled in a flexible way. However, the function cROC.bnp cal also deal with categorical covariates, linear effects of continuous covariates, as well as, interactions. For the reasons mentioned before, we recommend that all continuous covariates are modelled as in (<ref>). Nonetheless, posterior predictive checks, as illustrated in Section sec:illustration, can also be used to informally validate the fitted model. We write

\begin{equation}\label{mu_ddp}
\mu_{Dl}(\mathbf{x}_{Dj})=\mu_{D}(\mathbf{x}_{Dj},\boldsymbol{\beta}_{Dl})=\mathbf{z}_{Dj}^{\top}\boldsymbol{\beta}_{Dl}, \quad l=1,\ldots,L_D,\quad j =1,\ldots,n_D,
\end{equation}

where $\mathbf{z}_{Dj}$ is the $j$th column of the design matrix that contains the intercept, the continuous covariates that are modelled in a linear way (if any), the cubic B-splines basis representation for those modelled in a flexible way, the categorical covariates (if any), and their interaction(s) (if believed to exist). Also, $\boldsymbol{\beta}_{Dl}$ collects, for the $l$th component, the regression coefficients associated with the aforementioned covariates. For the covariate effects modelled using B-splines, an important issue is the selection of the number and location of the knots at which to anchor the basis functions, as this has the potential to impact inferences, more so for the former rather than the latter. The selection of the number of knots can be assisted by a model selection criterion, for example, (the adaptation to the case of mixture models of) the deviance information criterion (DIC) [Celeux et al., 2006], the log pseudo marginal likelihood (LPML) [Geisser and Eddy, 1979], and the widely applicable information criterion (WAIC) [Gelman et al., 2014]. In turn, for the location of the interior knots themselves we follow [Rosenberg, 1995] and use the quantiles of the covariate values.

The regression coefficients and variances associated with each of the $L_D$ components are sampled from the conjugate centring distribution $(\boldsymbol{\beta}_{Dl},\sigma_{Dl}^{-2})\overset{\text{iid}}\sim\text{N}(\mathbf{m}_D,\mathbf{S}_D)\Gamma(a_D,b_D)$, with conjugate hyperpriors $\mathbf{m}_D\sim\text{N}(\mathbf{m}_{D0},\mathbf{S}_{D0})$ and $\mathbf{S}_D^{-1}\sim\text{Wishart}(\nu_D,(\nu_D\Psi_D)^{-1})$ (a Wishart distribution with degrees of freedom $\nu_D$ and expectation $\Psi_D^{-1}$). Hyperparameters $\mathbf{m}_{D0}$ and $\Psi_D$ must be chosen to represent the prior belief in the regression coefficients and in their covariance matrix, whereas $\mathbf{S}_{D0}$ and $\nu_D$ are chosen to represent the confidence in the prior belief of $\mathbf{m}_{D0}$ and $\Psi_D$, respectively. The values of $a_D$ and $b_D$ on the prior for the components' variances can be chosen to represent belief in the variance of the regression model. With regard to the specification of $\alpha_D$ and $L_D$, an analogous reason to the DPM model (no-covariate case) can be followed and, as before, we consider $\alpha_D\sim\Gamma(a_{\alpha_{D}},b_{\alpha_{D}})$. The blocked Gibbs sampler is used to simulate draws from the posterior distribution and details about it can be found, for instance, in the Supplementary Materials of [Inácio de Carvalho et al., 2017].

Similarly to the analogous model for the no covariate case, at iteration $s$ of the Gibbs sampler procedure, the covariate-specific ROC curve is computed as

\begin{equation*}
\text{ROC}^{(s)}(p\mid\mathbf{x})=1-F_D^{(s)}\left\{F_{\bar{D}}^{-1(s)}(1-p\mid\mathbf{x})\mid\mathbf{x}\right\},  \quad s=1,\ldots, S,
\end{equation*}

\begin{equation}\label{cdfdpm}
F_D^{(s)}(y\mid\mathbf{x})=\sum_{l=1}^{L_D}\omega_{Dl}^{(s)}\Phi\left(y\mid \mathbf{z}^{\top}\boldsymbol{\beta}_{Dl}^{(s)},\sigma_{Dl}^{2(s)}\right),\quad
\end{equation}

and where the inversion is performed numerically. A point estimate for $\text{ROC}(p\mid\mathbf{x})$ can be obtained by computing the mean (or the median) of the ensemble $\{\text{ROC}^{(1)}(p\mid\mathbf{x}),\ldots,\text{ROC}^{(S)}(p\mid\mathbf{x})\}$, with pointwise credible bands derived from the percentiles of the ensemble.

Although the results presented in [Erkanli et al., 2006] can be extended to derive a closed-form expression for the covariate-specific AUC, for computational reasons, in ROCnReg the integral in (<ref>) is approximated using Simpson's rule, and the same applies for the partial areas.
Conditionally on a specific covariate value, the computation of the Youden index and of the optimal threshold proceeds in a similar way as in the DPM model <cit.>. As for the covariate-specific ROC curve, point and interval estimates can be obtained from the corresponding covariate-specific ensemble of each summary measure.

§.§ Covariate-adjusted ROC curve

All estimators for the covariate-adjusted ROC curve make use of Equation (<ref>) and rely on the following three steps

* Estimation of the conditional distribution of test outcomes in the nondiseased group, $F_{\bar{D}}(y_{\bar{D}i}\mid\mathbf{x}_{\bar{D}i})$.

* Computation of the placement value $U_D=1-F_{\bar{D}}(Y_D\mid\mathbf{X}_D)$ where, by a slight abuse of notation, we are designating it by the same letter used for the unconditional case.

* Estimation of the cumulative distribution function of $U_D$.

§.§.§ Frequentist approaches

The approaches used for estimation of the AROC curve proposed by [Janes and Pepe, 2009] and [Rodríguez-Álvarez et al., 2011] only differ in Step 1. Specifically, once one has an estimate of the conditional CDF in the nondiseased group, say $\widehat{F}_{\bar{D}}(\cdot\mid\mathbf{x})$, Step 2 in the two approaches consists of trivially computing the diseased placement values as

\begin{equation*}
\widehat{U}_{Dj}=1-\widehat{F}_{\bar{D}}(y_{Dj}\mid\mathbf{x}_{Dj}),\quad j = 1,\ldots, n_D.
\end{equation*}

Next, in Step 3, the AROC curve at a false positive fraction of $p$ is estimated via the empirical distribution
function of the placement values calculated in the previous step, $\{\widehat{U}_{Dj}\}_{j=1}^{n_D}$, that is,

\begin{equation*}
\widehat{\text{AROC}}(p)=\frac{1}{n_D}\sum_{j=1}^{n_D}I(\widehat{U}_{Dj}\leq p).
\end{equation*}

With regard to Step 1, both authors assume a location-scale regression model for the test outcomes in the nondiseased group and, as such and as explained in the previous section, the conditional CDF of the test results can be written using the CDF of the regression errors, i.e.,

\begin{equation*}
\end{equation*}

While [Janes and Pepe, 2009] assume a location-scale model of the form of (<ref>), [Rodríguez-Álvarez et al., 2011] rely on specification (<ref>). The estimation of the mean and variance functions follow exactly the same procedures as those described in the induced semiparametric linear model (for Janes and Pepe, 2009) and induced kernel-based approach (for Rodríguez-Álvarez et al., 2011) for the covariate specific ROC curve (the only difference being that here we only need to perform the estimation for the nondiseased group). At last, and also as in the estimators for the covariate-specific ROC curve, $F_{\varepsilon_{\bar{D}}}$ can be either assumed to be the standard normal distribution or left unspecified and estimated empirically on the basis of the standardised residuals. In both cases, the AAUC and pAAUC can be computed as follows

\begin{align*}
\widehat{\text{AAUC}} &= \int_{0}^{1}\text{AROC}(p)\text{d}p=1-\frac{1}{n_D}\sum_{j=1}^{n_D}\widehat{U}_{Dj},\\
\widehat{\text{pAAUC}(u_1)} &= \int_{0}^{u_1}\text{AROC}(p)\text{d}p=u_1-\frac{1}{n_D}\sum_{j=1}^{n_D}\min\left\{u_1,\widehat{U}_{Dj}\right\},
\end{align*}

whereas the $\text{pAAUC}_{\text{TNF}}$ is computed as in Equation (<ref>) using numerical integration methods (function integrate in R package stats).

§.§.§ Bayesian nonparametric approach

Recently, [Inácio de Carvalho and Rodríguez-Álvarez, 2018] proposed a Bayesian nonparametric estimator for the AROC curve that combines a dependent Dirichlet process mixture model and the Bayesian bootstrap. As in the Bayesian nonparametric approach for estimating the covariate-specific ROC curve, in Step 1, rather than assuming a location-scale regression model for the test outcomes in the nondiseased group, the entire conditional distribution is modelled using a single-weights dependent Dirichlet process mixture of normal distributions. Again, here, we also recommend to use cubic B-splines transformations of all continuous covariates. Using the same notation as before, we model the conditional density as

\begin{equation*}
\end{equation*}

Once Step 1 has been completed and given a posterior sample from the parameters of interest, the
corresponding realisation of the placement value of a diseased subject in the nondiseased population is easily computed as

\begin{equation*}
U_{Dj}^{(s)}=1-F_{\bar{D}}^{(s)}(y_{Dj}\mid\mathbf{x}_{Dj})=\sum_{l=1}^{L_{\bar{D}}}\omega_{\bar{D}l}^{(s)}\Phi\left(y_{Dj}\mid\mathbf{z}_{Dj}^{\top}\boldsymbol{\beta}_{\bar{D}l}^{(s)},\sigma_{\bar{D}l}^{2(s)}\right), \quad j =1,\ldots,n_D,\quad s=1,\ldots, S.
\end{equation*}

Finally, in Step 3, the cumulative distribution function of $\{U_{Dj}^{(s)}\}_{j=1}^{n_D}$ is estimated through the Bayesian bootstrap

\begin{equation*}
\text{AROC}^{(s)}(p)=\sum_{j=1}^{n_D}q_{j}^{(s)}I\left(U_{Dj}^{(s)}\leq p\right),\quad (q_{1}^{(s)},\ldots,q_{n_{D}}^{(s)})\sim\text{Dirichlet}(n_D; 1,\ldots,1).
\end{equation*}

As before, closed-form expressions do exist for the AAUC and pAAUC

\begin{align*}
\text{AAUC}^{(s)}&=\int_{0}^{1}\text{AROC}^{(s)}(p)\text{d}p=1-\sum_{j=1}^{n_D}q_j^{(s)}U_{Dj}^{(s)},\\
\text{pAAUC}^{(s)}(u_1)&=\int_{0}^{u_1}\text{AROC}^{(s)}(p)\text{d}p=u_1-\sum_{j=1}^{n_D}q_j^{(s)}\min\left\{u_1,U_{Dj}^{(s)}\right\},
\end{align*}

and the $\text{pAAUC}_{\text{TNF}}$ (Equation (<ref>)) is computed using numerical integration methods. With regard to the YI, it is obtained by directly plugging in $\text{AROC}^{(s)}(p)$ in Expression (<ref>).

A point estimate for $\text{AROC}(p)$ can be obtained by computing the mean (or the median) of the ensemble $\{\text{AROC}^{(1)}(p),\ldots,\text{AROC}^{(S)}(p)\}$, that is,

\begin{equation*}
\widehat{\text{AROC}}(p)=\frac{1}{S}\sum_{s=1}^{S}\text{AROC}^{(s)}(p), %\quad 0\leq p \leq 1,
\end{equation*}

and the percentiles of the ensemble can be used to provide pointwise credible bands/credible intervals. The same applies for the AAUC and pAAUC.

§ PACKAGE PRESENTATION AND ILLUSTRATION

This section describes the main functions in the ROCnReg package and illustrates their usage using endocrine data from a cross-sectional study carried out by the Galician Endocrinology and Nutrition Foundation (FENGA). A detailed description of this dataset can be found in [Tomé Martínez de Rituerto et al., 2009]. It has also been previously analysed in [Rodríguez-Álvarez et al., 2011, Rodríguez-Álvarez et al., 2011] and [Inácio de Carvalho and Rodríguez-Álvarez, 2018]. For confidentiality reasons, the data used in this paper correspond to a synthetic dataset that was obtained by mimicking the original one and can be found in the ROCnReg package under the name endosyn. A summary of the data follows.

R> library("ROCnReg")
R> data("endosyn")
R> summary(endosyn)

cvd_idf            age          gender          bmi
Min.   :0.0000   Min.   :18.25   Men  :1317   Min.   :12.60
1st Qu.:0.0000   1st Qu.:29.57   Women:1523   1st Qu.:23.19
Median :0.0000   Median :39.28                Median :26.24
Mean   :0.2433   Mean   :41.43                Mean   :26.69
3rd Qu.:0.0000   3rd Qu.:50.84                3rd Qu.:29.74
Max.   :1.0000   Max.   :84.66                Max.   :46.20

The dataset is comprised of $2840$ individuals ($1317$ men and $1523$ women, variable gender), with an age range between $18$ and $85$ years old. Variable bmi contains the body mass index (BMI) values, and  cvd_idf is the variable that indicates the presence (1) or absence (0) of two or more cardiovascular disease (CVD) risk factors. Following previous studies, the CVD risk factors considered include raised triglycerides, reduced HDL-cholesterol, raised blood pressure, and raised fasting plasma glucose. Note that from the $2840$ individuals, about $24\%$ present two or more CVD risk factors.

Using the ROCnReg package, in the subsequent sections we will illustrate how to ascertain, through the pooled ROC curve, the discriminatory capacity of the BMI (our diagnostic test) in differentiating individuals with two or more CVD risk factors (those belonging to the diseased class $D$) from those having none or just one CVD risk factor (and that therefore belong to the nondiseased group $\bar{D}$). Also, in Section <ref> we will show how to evaluate, through the covariate-specific ROC curve, the possible modifying effect of age and gender on the discriminatory capacity of the BMI. Finally, Section <ref> focuses on the covariate-adjusted ROC curve as a global summary measure of the BMI discriminatory ability, when taking the age and gender effects into account. We mention that in the Appendices we show the usage of the package for those methods not described in the main text.

§.§ Pooled ROC curve

The ROCnReg package allows estimating the pooled ROC curve by means of the four methods described in Section <ref>. Recall that function pooledROC.emp implements the empirical estimator, pooledROC.kernel the kernel-based approach, and pooledROC.BB and pooledROC.dpm correspond, respectively, to the Bayesian bootstrap estimator and the approach based on Dirichlet process mixtures (of normal distributions). The input arguments in the functions are method-specific (details can be found in the manual accompanying the package), but in all cases numerical and graphical summaries can be obtained by calling the functions print.pooledROC, summary.pooledROC, and plot.pooledROC, which can be abbreviated by print, summary, and plot.

By way of example, we present here the syntax using the pooledROC.dpm function. Recall that our aim is to ascertain, using the endosyn dataset, the discriminatory capacity of the BMI in differentiating individuals with two or more CVD risk factors from those having none or just one CVD risk factor.

R> set.seed(123, "L'Ecuyer-CMRG") # for reproducibility
R> pROC_dpm <- pooledROC.dpm(marker = "bmi", group = "cvd_idf",
+    tag.h = 0, data = endosyn, standardise = TRUE, p = seq(0, 1, l = 101),
+    compute.lpml = TRUE, compute.WAIC = TRUE, compute.DIC = TRUE,
+    pauc = pauccontrol(compute = TRUE, focus = "FPF", value = 0.1),
+    density = densitycontrol(compute = TRUE), prior.h = priorcontrol.dpm(),
+    prior.d = priorcontrol.dpm(),
+    mcmc = mcmccontrol(nsave = 8000, nburn = 2000, nskip = 1),
+    parallel = "snow", ncpus = 2, cl = NULL)

Before describing in detail the previous call, we first present the control functions that are used. In particular

pauccontrol(compute = FALSE, focus = c("FPF", "TPF"), value = 1)

can be used to indicate whether the pAUC should be computed (by default it is not computed), and in case it is computed (i.e., compute = TRUE ), whether the focus should be placed on restricted FPFs (pAUC, see (<ref>)) or on restricted TPFs ($\text{pAUC}_{\text{TPF}}$, see (<ref>)). In both cases, the upper bound $u_1$ (if focus is the FPF) or the lower bound $v_1$ (if focus is the TPF) should be indicated in value. In addition to the pooled ROC curve, AUC, and pAUC (if required), the function pooledROC.dpm also allows computing the probability density function (PDF) of the test outcomes in both the diseased and nondiseased groups. In order to do so, we use

densitycontrol(compute = FALSE, grid.h = NA, grid.d = NA)

By default, PDFs are not returned by the function pooledROC.dpm, but this can be changed by setting compute = TRUE, and through grid.h and grid.d the user can specify a grid of test results where the PDFs are to be evaluated in, respectively, the nondiseased and diseased groups. Value NA signals auto initialisation, with default a vector of length $200$ in the range of the test results. Regarding the hyperparameters for the Dirichlet process mixture of normals model (used for the estimation of the PDFs/CDFs of the test outcomes in each group), they can be controlled using

priorcontrol.dpm(m0 = NA, S0 = NA, a = 2, b = NA, aalpha = 2, balpha = 2,
L = 10)

A detailed description of these hyperparameters is found in Section <ref>. Finally, to set the various parameters controlling the Markov chain Monte Carlo (MCMC) procedure (which in our case is simply a Gibbs sampler) we use

mcmccontrol(nsave = 8000, nburn = 2000, nskip = 1)

Here nsave is an integer value with the total number of scans to be saved, nburn is the number of burn-in scans, and nthin is the thinning interval. Unless due to memory usage reasons, we recommend not thinning and instead monitoring the effective sample size of the MCMC chain.

Coming back to the pooledROC.dpm function, through marker the user specifies the name of the variable containing the test results; in our case, these are the values of the BMI. The name of the variable that distinguishes diseased (two or more CVD risk factors, $D$) from nondiseased individuals (none or one CVD risk factor, $\bar{D}$) is represented by the argument group, and the value codifying nondiseased individuals in group is specified by tag.h. The data argument is a data frame containing the data and all needed variables. Setting standardise = TRUE (the default) will standardise (i.e., subtract the mean and divide by the standard deviation) the test outcomes, which may help improving the MCMC mixing. The set of FPFs at which to estimate the pooled ROC curve is specified in the argument p. The LPML, WAIC and the DIC are computed by setting, respectively, the arguments compute.lpml, compute.WAIC, and compute.DIC to TRUE. Argument pauc is an (optional) list of values to replace the default values returned by the function pauccontrol. Here, we ask for the pAUC to be computed, with focus on restricted FPFs and upper bound $u_1 = 0.1$. Similarly, the argument density is an (optional) list of values to replace the default values returned by the function densitycontrol, as it is the argument mcmc. Through prior.h and prior.d arguments we specify the hyperparameters in the nondiseased and diseased classes, respectively. Again, both arguments are (optional) lists of values to replace the default values returned by the function priorcontrol.dpm. Different hyperparameters' default values are setted depending on whether test outcomes are standardised or not. Finally, arguments parallel, ncpus and cl allow to perform parallel computations (based on the R-package parallel). In particular, through parallel the user specifies the type of parallel operation: either "no" (default), "multicore" (not available on Windows Operating Systems) or "snow". Argument ncpus is used to indicate the number of processes to be used in parallel operation (when parallel = "multicore" or parallel = "snow"), and cl is an optional parallel or snow cluster to be used when parallel = "snow". If cl is not supplied (as in our example), a cluster on the local machine is created for the duration of the call.

A numerical summary of the fitted model can be obtained by calling the function summary, that provides, among other information, the estimated AUC (posterior mean) and $95\%$ credible interval and, if required, the LPML, WAIC, and DIC, separately, in the nondiseased (denoted here as Group H) and diseased (Group D) classes.

R> summary(pROC_dpm)

pooledROC.dpm(marker = "bmi", group = "cvd_idf", tag.h = 0, data = endosyn,
standardise = TRUE, p = seq(0, 1, l = 101), compute.lpml = TRUE,
compute.WAIC = TRUE, compute.DIC = TRUE, pauc = pauccontrol(compute = TRUE,
focus = "FPF", value = 0.1), density = densitycontrol(compute = TRUE),
prior.h = priorcontrol.dpm(), prior.d = priorcontrol.dpm(),
mcmc = mcmccontrol(nsave = 8000, nburn = 2000, nskip = 1),
parallel = "snow", ncpus = 2)

Approach: Pooled ROC curve - Bayesian DPM
Area under the pooled ROC curve: 0.758 (0.739, 0.777)
Partial area under the pooled ROC curve (FPF = 0.1): 0.168 (0.139, 0.199)

Model selection criteria:
Group H       Group D
WAIC               12491.338      4016.333
WAIC (Penalty)         6.351         4.087
LPML               -6245.669     -2008.167
DIC                12491.208      4016.253
DIC (Penalty)          6.286         4.047

Sample sizes:
Group H     Group D
Number of observations        2149         691
Number of missing data           0           0

To complement these numerical results, the ROCnReg package furnishes graphical results that can be used to further explore the fitted model. Specifically, the function plot depicts the estimated pooled ROC curve and AUC (posterior means), jointly with $95\%$ credible intervals

R> plot(pROC_dpm, cex.main = 1.5, cex.lab = 1.5, cex.axis = 1.5, cex = 1.5)

The result of the above code is shown in Figure <ref>.

Graphical results as provided by the plot.pooledROC function for an object of class pooledROC.dpm. Posterior mean and 95% pointwise credible band for the pooled ROC curve and corresponding posterior mean and 95% credible interval for the AUC.

By means of density = densitycontrol(compute = TRUE) in the call to the function, the estimates of the PDFs of the BMI in both classes are to be returned. This information can be accessed through component dens in object pROC_dpm (i.e., pROC_dpm$dens), which is a list with elements h and d associated with the nodiseased and diseased groups, respectively. Each of the two elements is itself a list of two components: (1) grid, a vector that contains the grid of test results at which the PDFs have been evaluated (estimated); and (2) dens, a matrix with the PDFs at each iteration of the MCMC procedure. We can use these results to plot, e.g, the posterior mean (and 95% pointwise credible bands) of the PDF of the BMI in the healthy and diseased populations (see Figure <ref>, obtained using the R package ggplot2 by Wickham, 2016). As can be observed, the estimated densities obtained under the DPM method follow very closely the histograms of the data. Further, the estimated densities available in dens can be used, as advised by [Gelman et al., 2013] (p. 553), to monitor convergence of the MCMC chains. The well-known label switching problem often leads to poor mixing of the chains of the component-specific parameters, but this may not impact convergence and mixing of the induced density/distribution of interest. For instance, Figure <ref> shows trace plots of the MCMC iterations (after burn-in) of the PDFs of the BMI in the two groups, for different (and randomly selected) values of the BMI, and Figure <ref> depicts the corresponding effective sample sizes (obtained using the R package coda by Plummer et al., 2006). Note that both plots give evidence of a good mixing and do not suggest lack of convergence. For conciseness, the R-code for producing Figures <ref>, <ref> and <ref> is not provided here, but in the replication code that accompanies this paper (<https://bitbucket.org/mxrodriguez/rocnreg>).

[DPM model with 10 mixture components in each group]

[Normal model in each group]

Histogram of the (observed) BMI and posterior mean jointly along with 95% pointwise credible bands (red lines) of the PDF of the BMI obtained using (a) a Dirichlet process mixture of normals model (object pROC_dpm); and (b) a normal model (object pROC_normal). Left: Healthy individuals (none or one CVD risk factor). Right: Diseased individuals (two or more CVD risk factors).

Trace plots of the MCMC draws (after burn-in) of the PDFs of the BMI based on model pROC_dpm. Results are shown separately for the healthy and diseased populations and for different values of the BMI.

Effective sample size of the MCMC chains (after burn-in) of the PDFs of the BMI based on model pROC_dpm in the healthy and diseased population. In both cases, results are shown along BMI values.

It is worth noting that the function pooledROC.dpm also allows fitting a normal distribution in each group; this is just a particular case (for which $L_D=L_{\bar{D}}=1$) of the more general DPM model. In order to fit such model, one simply needs to set $L=1$ in the prior.d and prior.h arguments. The code follows.

R> pROC_normal <- pooledROC.dpm(marker = "bmi", group = "cvd_idf",
+    tag.h = 0, data = endosyn,
+    standardise = TRUE, p = seq(0, 1, l = 101),
+    compute.lpml = TRUE, compute.WAIC = TRUE, compute.DIC = TRUE,
+    pauc = pauccontrol(compute = TRUE, focus = "FPF", value = 0.1),
+    density = densitycontrol(compute = TRUE),
+    prior.h = priorcontrol.dpm(L = 1), prior.d = priorcontrol.dpm(L = 1),
+    mcmc = mcmccontrol(nsave = 8000, nburn = 2000, nskip = 1),
+    parallel = "snow", ncpus = 2, cl = NULL)

For the sake of space we omit from the summary the call to the function

R> summary(pROC_normal)

Call: [...]

Approach: Pooled ROC curve - Bayesian DPM
Area under the pooled ROC curve: 0.748 (0.727, 0.768)
Partial area under the pooled ROC curve (FPF = 0.1): 0.224 (0.194, 0.253)

Model selection criteria:
Group H       Group D
WAIC               12639.946      4048.987
WAIC (Penalty)         2.426         2.252
LPML               -6319.973     -2024.493
DIC                12639.502      4048.703
DIC (Penalty)          1.983         1.979

Sample sizes:
Group H     Group D
Number of observations        2149         691
Number of missing data           0           0

The fit of the DPM and normal models, in each group, can be compared on the basis of the WAIC, DIC, and/or the LPML. Remember that for the LPML, the higher its value, the better the model fit, while for the WAIC and DIC it is the other way around. By comparing these values, provided in the summary of each fitted model, we can conclude that the three criteria favour, in both the diseased and (especially in the) nondiseased groups, the more general DPM model. This is also corroborated by the plot of the fitted densities in each group shown in Figure <ref>.

We now estimate the pooled ROC curve using the empirical estimator (function pooledROC.emp), and comparisons with the results obtained using the DPM approach are provided.

R> pROC_emp <- pooledROC.emp(marker = "bmi", group = "cvd_idf",
+    tag.h = 0, data = endosyn, p = seq(0, 1, l = 101),
+    pauc = pauccontrol(compute = TRUE, focus = "FPF", value = 0.1), B = 500,
+    parallel = "snow", ncpus = 2)

R> summary(pROC_emp)

Call: [...]

Approach: Pooled ROC curve - Empirical
Area under the pooled ROC curve: 0.76 (0.743, 0.778)
Partial area under the pooled ROC curve (FPF = 0.1): 0.169 (0.14, 0.201)

Sample sizes:
Group H     Group D
Number of observations        2149         691
Number of missing data           0           0

Note that the posterior means for the AUC and pAUC obtained using the DPM method ($0.758$ and $0.168$, respectively) are very similar to the point estimates using the empirical approach ($0.760$ and $0.169$). This can also be observed when plotting the estimated ROC curves under the two methods (Figure <ref>).

Estimated ROC curve using the empirical approach (black line) and using the DPM method (red line–posterior mean).

We finish this section by showing how to use ROCnReg to obtain an (optimal) threshold value which could be further used to `diagnose' an individual as diseased (two or more CVD risk factors) or healthy/nondiseased (none or only one CVD risk factor). To that aim, and for pooledROC objects (i.e., those obtained using functions pooledROC.dpm, pooledROC.BB, pooledROC.emp, and pooledROC.kernel), we use the function compute.threshold.pooledROC, which allows obtaining (optimal) threshold values using two criteria: the YI and the one that sets a target value for the FPF. For illustration, we show here the results using the YI criterion.

R> th_pROC_dmp <- compute.threshold.pooledROC(pROC_dpm, criterion = "YI",
+    parallel = "snow", ncpus = 2)
R> th_pROC_dmp

compute.threshold.pooledROC(object = pROC_dpm, criterion = "YI",
parallel = "snow", ncpus = 2)

est       ql       qh
26.57403 26.25361 26.91658

est        ql        qh
0.4043355 0.3716829 0.4359299

est        ql        qh
0.3710136 0.3406197 0.4021607

est        ql        qh
0.7753490 0.7450765 0.8041786
\end{CodeOutput}
\end{CodeChunk}

As can be observed, the function returns the posterior mean (\code{est}) and $95%$ credible interval (lower bound: \code{ql}, upper bound: \code{qh}) for the YI and associated threshold value, as well as for the FPF and TPF linked to this cutoff value. For our example, the (posterior mean of the) YI is $0.40$ and the YI-based threshold value is a BMI value of $26.6$, which falls in the nutritional status defined as pre-obesity by the World Health Organization. By using this YI-based threshold value, we would have a FPF of $0.37$ and a TPF of $0.78$.
\subsection{Covariate-specific ROC curve}\label{sec:croc_ilu}
We now turn our attention to the inclusion of covariates into the ROC analysis. As shown in Table \ref{tab:overview} and Section \ref{sec:methods}, with \pkg{ROCnReg} the user can estimate the covariate-specific ROC curve by means of three approaches: function \code{cROC.sp} implements the frequentist parametric and semiparametric induced ROC regression estimator, \code{cROC.kernel} corresponds to the nonparametric, kernel-based, counterpart of \code{cROC.sp}, and \code{cROC.bnp} stands for the Bayesian approach based on a single-weights dependent Dirichlet process mixture of normal distributions. As for the functions in \pkg{ROCnReg} for estimating the pooled ROC curve, the input arguments are method-specific, and we refer the reader to the manual for details. For all methods, numerical and graphical summaries are obtained using functions \code{print.cROC}, \code{summary.cROC} and \code{plot.cROC}. Also, for objects of class \code{cROC.bnp}, \pkg{ROCnReg} provides the function \code{predictive.ckecks}, which implements tools for assessing model fit via posterior predictive checks.

Recall that, when including covariate information into the ROC analysis, interest resides in evaluating if and how the discriminatory capacity of the test varies with such covariates. In particular, in our endocrine study we aim at evaluating the possible effect of both \code{age} and \code{gender} in the discriminatory capacity of the BMI. In what follows, we do that using the \code{cROC.bnp} function, and two different models are fitted. One which considers a normal distribution in each group and that incorporates the \code{age} effect in a linear way, and a second one which caps the maximum number of mixture components in each group at 10 (i.e., $L_D=L_D̅=10$) and that models the \code{age} effect using cubic B-splines (and thus allows for a nonlinear effect of age). Following [Rodr{\'i}guez-{\'A}lvarez \emph{et~al.}, 2011, Rodr{\'i}guez-{\'A}lvarez \emph{et~al.}, 2011], both models consider the interaction between \code{age} and \code{gender}. For clarity, we first focus on the code when modelling \code{age} effect in a linear way, and use it to describe in detail the different arguments of the \code{cROC.bnp} function.

\begin{Schunk}
\begin{Sinput}
R> # Dataframe for predictions
R> agep <- seq(22, 80, l = 30)
R> endopred <- data.frame(age = rep(agep,2),
+    gender = factor(rep(c("Women", "Men"), each = length(agep))))

R> set.seed(123, "L'Ecuyer-CMRG") # for reproducibility
R> cROC_bp <- cROC.bnp(formula.h = bmi ~ gender*age,
+    formula.d = bmi ~ gender*age, group = "cvd_idf", tag.h = 0,
+    data = endosyn, newdata = endopred,
+    standardise = TRUE, p = seq(0, 1, l = 101),
+    compute.lpml = TRUE, compute.WAIC = TRUE, compute.DIC = TRUE,
+    pauc = pauccontrol(compute = FALSE),
+    prior.h = priorcontrol.bnp(L = 1), prior.d = priorcontrol.bnp(L = 1),
+    density = densitycontrol(compute = TRUE),
+    mcmc = mcmccontrol(nsave = 8000, nburn = 2000, nskip = 1),
+    parallel = "snow", ncpus = 2)
\end{Sinput}
\end{Schunk}

As can be seen, many arguments coincide with those of the function \code{pooledROC.dpm} (described in Section \ref{sec:pooledroc_ilu}). We thus focus here on those that are specific to \code{cROC.bnp}. The arguments \code{formula.h} and \code{formula.d} are \code{formula} objects specifying the model for the regression functions (see Equation \eqref{mu_ddp}) in, respectively, the nondiseased and diseased classes. They are similar to the formula used with the \code{glm} function, except that nonlinear functions (modelled by means of cubic B-splines) can be added using function \code{f} (an example will follow). Note that in both cases, the left-hand side of the formulas should include the name of the test/marker (in our case \code{bmi}). In our application, and for both groups, the model for the component's means includes, in addition to the linear effect of \code{age} and \code{gender}, the (linear) interaction between the two (i.e., \code{gender*age} $≡$ \code{gender + age + gender:age}). Through \code{newdata} the user can specify a new data frame containing the values of the covariates at which the covariate-specific ROC curve and AUC (and also pAUC and PDFs, if required) are to be computed. Finally, \code{prior.h} (the same holds for \code{prior.d}) is an (optional) list of values to replace the defaults returned by \code{priorcontrol.bnp}
\begin{Schunk}
\begin{Sinput}
priorcontrol.bnp(m0 = NA, S0 = NA, nu = NA, Psi = NA, a = 2, b = NA,
aalpha = 2, balpha = 2, L = 10)
\end{Sinput}
\end{Schunk}
which allows setting the hyperparameters for the dependent Dirichlet process mixture of normals model (see Section \ref{sec:methods}). In our example, we only modified the upper bound for the number of components in the mixture model (by default $L = 10$) and set it to $1$. With this configuration, the model for the covariate-specific ROC curve can be regarded as a Bayesian counterpart of the induced ROC approach proposed by [Faraggi, 2003] and we denote it as the Bayesian normal linear model (for the test outcomes).

In this case, the \code{summary} of the fitted model provides the following information. %For the sake of page space, in what follows we omit from the summary the call to the function.
\begin{Schunk}
\begin{Sinput}
R> summary(cROC_bp)
\end{Sinput}
\begin{Soutput}
Call: [...]

Approach: Conditional ROC curve - Bayesian nonparametric

Parametric coefficients
Group H:
Post. mean    Post. quantile 2.5%    Post. quantile 97.5%
(Intercept)           26.1454                25.8769                 26.4154
genderWomen           -0.9163                -1.2721                 -0.5643
age                    1.1953                 0.9192                  1.4655
genderWomen:age        1.1940                 0.8420                  1.5413

Group D:
Post. mean    Post. quantile 2.5%    Post. quantile 97.5%
(Intercept)           29.1876                28.7677                 29.6117
genderWomen            2.0807                 1.3839                  2.7680
age                    0.6592                 0.2217                  1.1075
genderWomen:age       -0.7741                -1.4798                 -0.0868

ROC curve:
Post. mean    Post. quantile 2.5%    Post. quantile 97.5%
(Intercept)           -0.6946                -0.8187                 -0.5767
genderWomen           -0.6843                -0.8665                 -0.5013
age                    0.1224                 0.0017                  0.2396
genderWomen:age        0.4494                 0.2737                  0.6313
b                      0.9377                 0.8810                  0.9957

Model selection criteria:
Group H       Group D
WAIC               12174.999      4008.025
WAIC (Penalty)         6.280         5.643
LPML               -6087.499     -2004.013
DIC                12173.684      4007.397
DIC (Penalty)          4.999         5.076

Sample sizes:
Group H     Group D
Number of observations        2149         691
Number of missing data           0           0
\end{Soutput}
\end{Schunk}
The first thing to note is that, in this case, the \code{summary} function does not provide the estimated AUC as there is one (possibly different) AUC for each combination of covariate values. Also, given that: (1) only one component has been considered for modelling the CDFs of test results in the diseased and nodiseased groups, and (2) covariate effects have been modelled in a linear way, the \code{summary} function provides the posterior mean (and quantiles) of the (parametric) coefficients associated with the regression functions (see Equation \eqref{loc-scale-sp}) and with the covariate-specific ROC curve (see Equation \eqref{a-b-sp}). We note that since in the call to the function we have specified \code{standardise = TRUE} (and consequently both the test outcomes and covariates are standardised), the regression coefficients are on the scale of the standardised covariates. If we focus on the coefficients for the covariate-specific ROC curve, it seems that the discriminatory capacity of the BMI decreases with age, with the decrease being more pronounced in women (note that the expression of the covariate-specific ROC curve in \eqref{cROC-a-b-sp} implies that positive coefficients correspond with a decrease in discriminatory capacity). These results are possibly better judged by plotting the estimated covariate-specific ROC curves and associated AUCs. This can be done using the \code{plot} function. For the covariate-specific ROC curve, the depicted graphics will depend on the number and nature of the covariates included in the analyses. In particular, for our application, we obtain, separately for men and women, the covariate-specific ROC curves (and AUCs) along age. These are shown in Figure \ref{cROC_bp_plot}, obtained using the code
\begin{Schunk}
\begin{Sinput}
R> op <- par(mfrow = c(2,2))
R> plot(cROC_sp, ask = FALSE)
R> par(op)
\end{Sinput}
\end{Schunk}

\begin{figure}[h!]
\begin{center}
\includegraphics[width=14cm]{cROC_bp.pdf}
\end{center}
\caption{Graphical results as provided by the \code{plot.cROC} function for an object of class \code{cROC.bnp}. Results for the model that includes the linear interaction between \code{age} and \code{gender} and one mixture component. Top row: Posterior mean of the covariate-specific ROC curve along age, separately for men and women. Bottom row: Posterior mean and 95\% pointwise credible band for the covariate-specific AUC along age, separately for men and women.}
\label{cROC_bp_plot}
\end{figure}

In this example we have modelled the \code{age} effect linearly and only one mixture component was considered. However, \pkg{ROCnReg} also allows for modelling the effect of continuous covariates in a nonlinear way, either using cubic B-spline basis expansions (through the function \code{cROC.bnp}) or  kernel-based smoothers (via the function \code{cROC.kernel}). Also, as noted before, using only one mixture component for the dependent Dirichlet process mixture of normals model (function \code{cROC.bnp}) is equivalent to consider a (Bayesian) normal model, which might be too restrictive for most data applications. In what follows, we provide more flexibility to the model for the covariate-specific ROC curve by means of (a) increasing the number of mixture components, and (b) modelling the \code{age} effect in a nonlinear way (recall our considerations in Section 4 about the lack of flexibility of the single-weights dependent Dirichlet process mixture of normals model when covariates effects on the components' means are modelled linearly). The former is done by modifying the value of \code{L} in the arguments \code{prior.h} and \code{prior.d}, with $10$ being the default value. Regarding the latter, this is done by making use of the function \code{f} when specifying the component's mean functions through \code{formula.h} and \code{formula.d}. In particular, in our application we are interested in modelling the factor-by-curve interaction between \code{age} and \code{gender} (i.e., we model \code{age} effect ``separately'' for men and women). This is done using, e.g., \code{bmi ~ gender + f(age, by = gender, K = c(3,5))}. Through the argument \code{K} we indicate the number of internal knots for constructing the cubic B-spline basis that is used to approximate the nonlinear effect of \code{age} (with the quantiles of \code{age} used to anchor the knots). Note that we can specify a different number of internal knots for men and women (\code{K = c(3,5)}), where the order of vector \code{K} should match the ordering of the levels of the factor \code{gender}. When using a cubic B-spline basis, one must choose the number of interior knots (in \pkg{ROCnReg} the location is always based on the quantiles of the corresponding covariates). Here, the task of selecting the number of interior knots is assisted by the WAIC, DIC, and/or LPML, i.e., we fit the model for different number of internal knots and consider the model that provided the lowest WAIC/DIC or the highest LPML (this is done in both the healthy and diseased populations and we remark that the number of knots does not need to be the same in the two). The final model is shown below.

\begin{Schunk}
\begin{Sinput}
R> # Levels of gender, and its ordering.
R> # Needed if we want to specify different
R> # number of knots for men and women
R> levels(endosyn$gender)
\end{Sinput}
\begin{Soutput}
[1] "Men"   "Women"
\end{Soutput}
\begin{Sinput}
R> set.seed(123, "L'Ecuyer-CMRG") # for reproducibility
R> cROC_bnp <- cROC.bnp(
+    formula.h = bmi ~ gender + f(age, by = gender, K = c(0,0))
+    formula.d = bmi ~ gender + f(age, by = gender, K = c(3,3)),
+    group = "cvd_idf", tag.h = 0, data = endosyn, newdata = endopred,
+    standardise = TRUE, p = seq(0, 1, l = 101),
+    compute.lpml = TRUE, compute.WAIC = TRUE, compute.DIC = TRUE,
+    pauc = pauccontrol(compute = FALSE),
+    prior.h = priorcontrol.bnp(L = 10), prior.d = priorcontrol.bnp(L = 10),
+    density = densitycontrol(compute = TRUE),
+    mcmc = mcmccontrol(nsave = 20000, nburn = 8000, nskip = 1),
+    parallel = "snow", ncpus = 2)

R> summary(cROC_bnp)
\end{Sinput}
\begin{Soutput}
Call: [...]

Approach: Conditional ROC curve - Bayesian nonparametric

Model selection criteria:
Group H       Group D
WAIC               11833.926      3915.618
WAIC (Penalty)        28.215        28.501
LPML               -5917.101     -1958.074
DIC                11831.484      3913.118
DIC (Penalty)         26.994        27.251

Sample sizes:
Group H     Group D
Number of observations        2149         691
Number of missing data           0           0
\end{Soutput}
\end{Schunk}
\begin{Schunk}
\begin{Sinput}
R> op <- par(mfrow = c(2,2))
R> plot(cROC_sp, ask = FALSE)
R> par(op)
\end{Sinput}
\end{Schunk}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=14cm]{cROC_bnp.pdf}
\end{center}
\caption{Graphical results as provided by the \code{plot.cROC} function for an object of class \code{cROC.bnp}. Results for the model that includes the factor-by-curve interaction between \code{age} and \code{gender} and 10 mixture components. Top row: Posterior mean of the covariate-specific ROC curve along age, separately for men and women. Bottom row: Posterior mean and 95\% pointwise credible band for the covariate-specific AUC along age, separately for men and women.}
\label{cROC_bnp_plot}
\end{figure}

The graphical results are shown in Figure \ref{cROC_bnp_plot}. Note that, especially for women, age displays a marked nonlinear effect. Recall that for objects of class \code{cROC.bnp}, and if required in the call to the function, the \code{summary} function provides, separately for the diseased and nondiseased/healthy groups, the WAIC, LPML, and DIC. Note that, in both cases, the three criteria support the use of the more flexible model that uses B-splines and 10 mixture components for modelling the distribution of the BMI (model \code{cROC_bnp}) over the more restrictive Bayesian normal linear model (model \code{cROC_bp}). Since WAIC, LPML, and DIC are relative criteria, posterior predictive checks are also available in \pkg{ROCnReg} through the function \code{predictive.checks}. Specifically, the function generates replicated datasets from the posterior predictive distribution in both class $D$ and $\bar{D}$ and compares them to the test values (BMI values in our application) using specific statistics. For the choice of such statistics we follow [Gabry \emph{et~al.}, 2019], who suggest choosing statistics that are `orthogonal' to the model parameters. Since we are using a location-scale mixture of normals model for the test outcomes, we use here the skewness and kurtosis and check how well the posterior predictive distribution captures these two quantities.
\begin{Schunk}
\begin{Sinput}
R> op <- par(mfrow = c(2,3))
R> pc_cROC_bp <- predictive.checks(cROC_bp,
+    statistics = c("kurtosis", "skewness"), devnew = FALSE)
R> par(op)

R> op <- par(mfrow = c(2,3))
R> pc_cROC_bnp <- predictive.checks(cROC_bnp,
+    statistics = c("kurtosis", "skewness"), devnew = FALSE)
R> par(op)
\end{Sinput}
\end{Schunk}
Results are shown in Figure \ref{cROC_bnp_pred_checks}. As can be seen, the model that includes the factor-by-curve interaction between \code{age} and \code{gender} and 10 mixture components performs quite well in capturing both quantities, while the Bayesian normal linear model fails to do so. Also shown in Figure \ref{cROC_bnp_pred_checks} (and provided by function \code{predictive.checks}) are the kernel density estimates of $500$ randomly selected datasets drawn from the posterior predictive distribution, in each group, compared to the kernel density estimate of the observed BMI (in each group). Again, the more flexible model, as opposed to the Bayesian normal linear model, is able to simulate data that are very much similar to the observed BMI values.
\begin{figure}[h!]
\begin{center}
\subfigure[Model including the linear interaction between \code{age} and \code{gender} and one component]{\includegraphics[width=15cm]{cROC_bp_predchecks.pdf}}
\subfigure[Model including the factor-by-curve interaction between \code{age} and \code{gender} and 10 mixture components]{\includegraphics[width=15cm]{cROC_bnp_predchecks.pdf}}
\end{center}
\caption{Graphical results as provided by the \code{predictive.checks} function for an object of class \code{cROC.bnp}. Histograms of the statistics \textit{skewness} and \textit{kurtosis} computed from $8000$ draws from the posterior predictive distribution in the diseased and nondiseased group. The red line is the estimated statistic from the observed BMI values. The right-hand side plots show the kernel density estimate of the observed BMI (solid black line), jointly with the kernel density estimates for $500$ simulated datasets drawn from the posterior predictive distributions.}
\label{cROC_bnp_pred_checks}
\end{figure}

As for the pooled ROC curve (Section \ref{sec:pooledroc_ilu}), \pkg{ROCnReg} also provides a function that allows obtaining (optimal) threshold values for the covariate-specific ROC curve. For illustration, instead of the threshold values based on the Youden index, we now use the criterion that sets a target value for the FPF. The code for model \code{cROC_bnp}, when setting the $\text{FPF} = 0.3$, is as follows.
\begin{Schunk}
\begin{Sinput}
R> th_fpf_cROC_bnp <- compute.threshold.cROC(cROC_bnp,
+    criterion = "FPF", FPF = 0.3, newdata = endopred,
+    parallel = "snow", ncpus = 2)
R> names(th_fpf_cROC_bnp)
\end{Sinput}
\begin{Soutput}
[1] "newdata"    "thresholds" "TPF"    "FPF"    "call"
\end{Soutput}
\end{Schunk}
In addition to the data frame \code{newdata} containing the covariate values at which the thresholds are computed, function \code{compute.threshold.cROC} returns the covariate-specific \code{thresholds} corresponding to the specified FPF (posterior mean and 95\% pointwise credible intervals) as well as the covariate-specific \code{TPF} (posterior mean and 95\% pointwise credible intervals) attached to these thresholds. Although \pkg{ROCnReg} does not provide a function for plotting the results obtained using \code{compute.threshold.cROC}, graphical results can be easily obtained. For simplicity, we only show here the code for the covariate-specific threshold values (\code{thresholds}), but a similar code can be used to plot the covariate-specific TPFs (\code{TPF}). Both plots are shown in Figure \ref{cROC_bnp_th}. As can be observed, for a FPF of $0.3$, the BMI age-specific thresholds tend to increase with age both for men and women, although for the latter there is a slight decrease after an age of about 70 years old. The age-specific TPFs corresponding to the thresholds for which the FPF is $0.3$ show a nonlinear behaviour and these are in general higher for women than for men (of the same age).
\begin{Schunk}
\begin{Sinput}
df <- data.frame(age = th_fpf_cROC_bnp$newdata$age,
+    gender = th_fpf_cROC_bnp$newdata$gender,
+    y = th_fpf_cROC_bnp$thresholds[[1]][,"est"],
+    ql = th_fpf_cROC_bnp$thresholds[[1]][,"ql"]),
+    qh = th_fpf_cROC_bnp$thresholds[[1]][,"qh"]))

R> g0 <- ggplot(df, aes(x = age, y = y, ymin = ql, ymax = qh)) +
+    geom_line() +
+    geom_ribbon(alpha = 0.2) +
+    labs(title = "Covariate-specific thresholds for a FPF = 0.3",
+    x = "Age (years)", y = "BMI") +
+    theme(strip.text.x = element_text(size = 20),
+    plot.title = element_text(hjust = 0.5, size = 20),
+    axis.text = element_text(size = 20),
+    axis.title = element_text(size = 20)) +
+    facet_wrap(~gender)
R> print(g0)
\end{Sinput}
\end{Schunk}
\begin{figure}[h!]
\begin{center}
\subfigure[BMI threshold values]{\includegraphics[width=12cm]{cROC_bnp_fpf_thresholds.pdf}}
\subfigure[True positive fractions]{\includegraphics[width=12cm]{cROC_bnp_fpf_TPF.pdf}}
\end{center}
\caption{Top row: Posterior mean (solid black line) and 95\% pointwise credible band for the BMI threshold values, along age, corresponding to a FPF of 0.3. Bottom row: Posterior mean (solid black line) and 95\% pointwise credible band of the TPFs, along age, corresponding to the BMI threshold values for which $\text{FPF}=0.3$.}
\label{cROC_bnp_th}
\end{figure}

Finally, we mention that for conciseness we have not shown here how to perform convergence diagnostics of the MCMC chains for models fitted using \code{cROC.bnp}. In very much the same way as shown in Section \ref{sec:pooledroc_ilu} for the object \code{pROC_dpm}, using the information contained in component \code{dens} in the list of returned values (if required), one can produce trace plots of the conditional densities at some sampled values, as well as, obtain the corresponding effective sample sizes. Some results are provided in Appendix \ref{AppB}, and the associated code can be found in the \proglang{R} replication code that accompanies this paper.
\subsection{Covariate-adjusted ROC curve}\label{sec:aroc_ilu}
In this last section we illustrate how to conduct inference about the covariate-adjusted ROC curve using  \pkg{ROCnReg}. Similarly to the covariate-specific ROC curve, three approaches are available, namely, function \code{AROC.sp} implements the frequentist approaches that postulate that test outcomes in the nondiseased group follow a linear model with the CDF of the error term being either a standard normal distribution or estimated via the empirical CDF of the standardised residuals, \code{AROC.kernel} corresponds to the kernel-based counterpart, and \code{AROC.bnp} implements the Bayesian (nonparametric) approach based on a single-weights dependent Dirichlet process mixture of normal distributions and the Bayesian bootstrap.

Recall that the AROC curve is a global summary measure of diagnostic accuracy that takes covariate information into account. In the context of our endocrine application we seek to study the overall discriminatory capacity of the BMI for detecting the presence of CVD risk factors when adjusting for age and gender. Here we focus on how to estimate the AROC curve using the \code{AROC.bnp} function. The function syntax is exactly similar to the one of \code{cROC.bnp}, with the only difference being that we only need to specify the formula for the components' means in the nondiseased population. The code and respective summary follow.
\begin{Schunk}
\begin{Sinput}
R> AROC_bnp <- AROC.bnp(
+    formula.h = bmi ~ gender + f(age, by = gender, K = c(0,0)),
+    group = "cvd_idf", tag.h = 0, data = endosyn, standardise = TRUE,
+    p = seq(0, 1, l = 101),
+    compute.lpml = TRUE, compute.WAIC = TRUE, compute.DIC = TRUE,
+    pauc = pauccontrol(compute = FALSE), prior = priorcontrol.bnp(L = 10),
+    density = densitycontrol.aroc(compute = FALSE),
+    mcmc = mcmccontrol(nsave = 8000, nburn = 2000, nskip = 1),
+    parallel = "snow", ncpus = 2)

R> summary(AROC_bnp)
\end{Sinput}
\begin{Soutput}
Call: [...]

Approach: AROC Bayesian nonparametric
Area under the covariate-adjusted ROC curve: 0.653 (0.625, 0.682)

Model selection criteria:
Group H
WAIC               11833.926
WAIC (Penalty)        28.215
LPML               -5917.101
DIC                11831.484
DIC (Penalty)         26.994

Sample sizes:
Group H     Group D
Number of observations        2149         691
Number of missing data           0           0
\end{Soutput}
\end{Schunk}
The area under the AROC curve (posterior mean) is $0.653$ (95\% credible interval: $(0.625, 0.682)$) thus revealing a reasonable good ability of the BMI to detect the presence of CVD risk factors when teasing out the age and gender effects. As for the pooled ROC curve and the covariate-specific ROC curve, a \code{plot} function is also available (result in Figure \ref{AROC_bnp}).
\begin{Schunk}
\begin{Sinput}
R> plot(AROC_bnp, cex.main = 1.5, cex.lab = 1.5, cex.axis = 1.5, cex = 1.3)
\end{Sinput}
\end{Schunk}
We finish with a comparison of the AROC curve with the pooled ROC curve that was obtained earlier by using a DPM model with 10 components in each group. In Figure \ref{AROC_bnp_pooled} we show the plots of the two curves and, as can be noticed, the pooled ROC curve lies well above the AROC curve, thus evidencing the need for incorporating covariate information into the analysis.
\begin{Schunk}
\begin{Sinput}
R> plot(AROC_bnp$p, AROC_bnp$ROC[,1],
+    type = "l", xlim = c(0,1), ylim = c(0,1),
+    xlab = "FPF", ylab = "TPF",
+    main = "Pooled ROC curve vs AROC curve",
+    cex.main = 1.5, cex.lab = 1.5, cex.axis = 1.5, cex = 1.5)
R> lines(AROC_bnp$p, AROC_bnp$ROC[,2], col = 1, lty = 2)
R> lines(AROC_bnp$p, AROC_bnp$ROC[,3], col = 1, lty = 2)
R> lines(pROC_dpm$p, pROC_dpm$ROC[,1], col = 2)
R> lines(pROC_dpm$p, pROC_dpm$ROC[,2], col = 2, lty = 2)
R> lines(pROC_dpm$p, pROC_dpm$ROC[,3], col = 2, lty = 2)
R> abline(0, 1, col = "grey", lty = 2)
\end{Sinput}
\end{Schunk}

\begin{figure}[h!]
\begin{center}
\subfigure[AROC curve]{\includegraphics[width=7.5cm]{AROC_bnp.pdf}\label{AROC_bnp}}
\subfigure[AROC vs Poooled ROC curve]{\includegraphics[width=7.5cm]{pooled_ROC_AROC.pdf}\label{AROC_bnp_pooled}}
\end{center}
\caption{(a) Age/gender-adjusted ROC curve: posterior mean and 95\% pointwise credible band. (b) Age/gender-adjusted ROC curve (in black) and pooled ROC curve (estimated using a DPM of normals model) (in red). Solid lines represent the posterior means and dashed lines the  95\% pointwise credible bands.}
\label{AROC_bnp_complete}
\end{figure}
%We end remarking that in the Supplementary Materials accompanying this paper we show the usage of the package for those methods not described in the main manuscript.
\section{Summary and future plans} \label{sec:summary}
In this paper we have introduced the capabilities of the \proglang{R} package \pkg{ROCnReg} for conducting inference about the pooled ROC curve, the covariate-specific ROC curve, and the covariate-adjusted ROC curve and their associated summary indices. As we have illustrated, the current version of the package provides several options for estimating ROC curves, both under frequentist and Bayesian paradigms, either parametrically, semiparametrically, or nonparametrically. To the best of our knowledge, this is the first software package implementing Bayesian inference for ROC curves. Several additions/extensions are planned in the future and these, among others, include:
\begin{itemize}
\item Make the MCMC algorithm faster by implementing its time-consuming parts in \proglang{C++}.
\item Incorporate methods for non-binary disease status (e.g., no disease, mild disease, severe disease). That is, implement ROC surface models.
\item Implement new (optimal) threshold criteria (e.g., YI including costs).
\end{itemize}
\section*{Computational details}
The results in this paper were obtained using \proglang{R}~4.0.2 with the \pkg{ROCnReg}~1.0-4 package. The \pkg{ROCnReg} package has multiple dependencies: \pkg{graphics}, \pkg{grDevices}, \pkg{parallel}, \pkg{splines}, \pkg{stats}, \pkg{moments} [Komsta and Novomestky, 2015], \pkg{nor1mix} [Maechler, 2019], \pkg{Matrix} [Bates and Maechler, 2019], \pkg{spatstat} [Baddeley and Turner, 2005], \pkg{np} [Hayfield and Racine, 2008], \pkg{lattice} [Sarkar, 2008], \pkg{MASS} [Venables and Ripley, 2002] and \pkg{pbivnorm} [Genz and Kenkel, 2015]. \proglang{R} itself and all packages used are available from the Comprehensive \proglang{R} Archive Network (CRAN) at \url{https://CRAN.R-project.org/}.
\section*{Acknowledgments}
MX Rodr\'iguez-\'Alvarez was funded by project MTM2017-82379-R (AEI/FEDER, UE), by the Basque Government through the BERC 2018-2021 program and Elkartek project 3KIA (KK-2020/00049), and by the Spanish Ministry of Science, Innovation, and Universities (BCAM Severo Ochoa accreditation SEV-2017-0718). The work of V In\'acio was partially supported by FCT (Funda\c c\~ao para a Ci\^encia e a Tecnologia, Portugal), through the projects PTDC/MAT-STA/28649/2017 and \linebreak UID/MAT/00006/2020.

\begin{thebibliography}{59}
\newcommand{\enquote}[1]{``#1''}
\providecommand{\natexlab}[1]{#1}
\providecommand{\url}[1]{\texttt{#1}}
\providecommand{\urlprefix}{URL }
\expandafter\ifx\csname urlstyle\endcsname\relax
\providecommand{\doi}[1]{doi:\discretionary{}{}{}#1}\else
\providecommand{\doi}{doi:\discretionary{}{}{}\begingroup
\urlstyle{rm}\Url}\fi
\providecommand{\eprint}[2][]{\url{#2}}

[Baddeley and Turner, 2005]
Baddeley A, Turner R (2005).
\newblock \enquote{\pkg{spatstat}: An \proglang{R} Package for Analyzing
Spatial Point Patterns.}
\newblock \emph{Journal of Statistical Software}, \textbf{12}(6), 1--42.
\newblock \urlprefix\url{http://www.jstatsoft.org/v12/i06/}.

[Bates and Maechler, 2019]
Bates D, Maechler M (2019).
\newblock \emph{\pkg{Matrix}: Sparse and Dense Matrix Classes and Methods}.
\newblock \proglang{R}~package version 1.2-17,
\urlprefix\url{https://CRAN.R-project.org/package=Matrix}.

[Buolamwini, 2017]
Buolamwini JA (2017).
\newblock \emph{Gender {S}hades: {I}ntersectional {P}henotypic and
{D}emographic {E}valuation of {F}ace {D}atasets and {G}ender {C}lassifiers}.
\newblock Ph.D. thesis, Massachusetts Institute of Technology, School of
Architecture and Planning, Program in Media Arts and Sciences.

[Celeux \emph{et~al.}, 2006]
Celeux G, Forbes F, Robert CP, Titterington DM (2006).
\newblock \enquote{Deviance {I}nformation {C}riteria for {M}issing {D}ata
\newblock \emph{Bayesian Analysis}, \textbf{1}(4), 651--673.

[De~Iorio \emph{et~al.}, 2009]
De~Iorio M, Johnson WO, M{\"u}ller P, Rosner GL (2009).
\newblock \enquote{Bayesian {N}onparametric {N}onproportional {H}azards
{S}urvival {M}odeling.}
\newblock \emph{Biometrics}, \textbf{65}(3), 762--771.

[Dunn and Smyth, 1996]
Dunn PK, Smyth GK (1996).
\newblock \enquote{{R}andomized {Q}uantile {R}esiduals.}
\newblock \emph{Journal of Computational and Graphical Statistics},
\textbf{5}(3), 236--244.

[Erkanli \emph{et~al.}, 2006]
Erkanli A, Sung M, Jane~Costello E, Angold A (2006).
\newblock \enquote{Bayesian {S}emi-{P}arametric {ROC} {A}nalysis.}
\newblock \emph{Statistics in Medicine}, \textbf{25}(22), 3905--3928.

[Fan and Gijbels, 1996]
Fan J, Gijbels I (1996).
\newblock \emph{Local Polynomial Modelling and Its Applications}.
\newblock Monographs on {S}tatistics and {A}pplied {P}robability. Chapman \&

[Faraggi, 2003]
Faraggi D (2003).
\newblock \enquote{Adjusting {R}eceiver {O}perating {C}haracteristic {C}urves
and {R}elated {I}ndices for {C}ovariates.}
\newblock \emph{Journal of the Royal Statistical Society D}, \textbf{52}(2),

[Ferguson, 1973]
Ferguson TS (1973).
\newblock \enquote{A {B}ayesian {A}nalysis of {S}ome {N}onparametric
\newblock \emph{The Annals of Statistics}, \textbf{1}(2), 209--230.

[Gabry \emph{et~al.}, 2019]
Gabry J, Simpson D, Vehtari A, Betancourt M, Gelman A (2019).
\newblock \enquote{Visualization in {B}ayesian {W}orkflow.}
\newblock \emph{Journal of the Royal Statistical Society A}, \textbf{182}(2),

[Geisser and Eddy, 1979]
Geisser S, Eddy WF (1979).
\newblock \enquote{A {P}redictive {A}pproach to {M}odel {S}election.}
\newblock \emph{Journal of the American Statistical Association},
\textbf{74}(365), 153--160.

[Gelman \emph{et~al.}, 2013]
Gelman A, Carlin JB, Stern HS, Dunson DB, Vehtari A, Rubin DB (2013).
\newblock \emph{Bayesian Data Analysis}.
\newblock Chapman and Hall/CRC.

[Gelman \emph{et~al.}, 2014]
Gelman A, Hwang J, Vehtari A (2014).
\newblock \enquote{Understanding {P}redictive {I}nformation {C}riteria for
{B}ayesian {M}odels.}
\newblock \emph{Statistics and Computing}, \textbf{24}(6), 997--1016.

[Genz and Kenkel, 2015]
Genz A, Kenkel B (2015).
\newblock \emph{\pkg{pbivnorm}: Vectorized Bivariate Normal CDF}.
\newblock \proglang{R} package version 0.6.0,
\urlprefix\url{https://CRAN.R-project.org/package=pbivnorm}.

[Gon{\c{c}}alves \emph{et~al.}, 2014]
Gon{\c{c}}alves L, Subtil A, Oliveira MR, Bermudez P (2014).
\newblock \enquote{{ROC} {C}urve {E}stimation: {A}n {O}verview.}
\newblock \emph{REVSTAT--Statistical Journal}, \textbf{12}(1), 1--20.

[Gonz{\'a}lez-Manteiga \emph{et~al.}, 2011]
Gonz{\'a}lez-Manteiga W, Pardo-Fern{\'a}ndez JC, van Keilegom I (2011).
\newblock \enquote{ROC {C}urves in {N}on-{P}arametric {L}ocation-{S}cale
{R}egression {M}odels.}
\newblock \emph{Scandinavian Journal of Statistics}, \textbf{38}(1), 169--184.

[Gu \emph{et~al.}, 2008]
Gu J, Ghosal S, Roy A (2008).
\newblock \enquote{Bayesian {B}ootstrap {E}stimation of {ROC} {C}urve.}
\newblock \emph{Statistics in Medicine}, \textbf{27}, 5407--5420.

[Guan \emph{et~al.}, 2012]
Guan Z, Qin J, Zhang B (2012).
\newblock \enquote{Information {B}orrowing {M}ethods for {C}ovariate-{A}djusted
{ROC} {C}urve.}
\newblock \emph{Canadian Journal of Statistics}, \textbf{40}(3), 569--587.

[Hayfield and Racine, 2008]
Hayfield T, Racine JS (2008).
\newblock \enquote{Nonparametric {E}conometrics: {T}he \pkg{np} {P}ackage.}
\newblock \emph{Journal of Statistical Software}, \textbf{27}(5).
\newblock \urlprefix\url{http://www.jstatsoft.org/v27/i05/}.

[Hillis and Metz, 2012]
Hillis SL, Metz CE (2012).
\newblock \enquote{An {A}nalytic {E}xpression for the {B}inormal {P}artial
{A}rea under the {ROC} {C}urve.}
\newblock \emph{Academic Radiology}, \textbf{19}(12), 1491 -- 1498.

[Hsieh and Turnbull, 1996]
Hsieh F, Turnbull BW (1996).
\newblock \enquote{Nonparametric and {S}emiparametric {E}stimation of the
{R}eceiver {O}perating {C}haracteristic {C}urve.}
\newblock \emph{The Annals of Statistics}, \textbf{24}(1), 25--40.

[Hutchinson and Mitchell, 2019]
Hutchinson B, Mitchell M (2019).
\newblock \enquote{50 {Y}ears of {T}est ({U}n){F}airness: {L}essons for
{M}achine {L}earning.}
\newblock In \emph{Proceedings of the Conference on Fairness, Accountability,
and Transparency}, FAT* '19, pp. 49--58. ACM, New York, NY, USA.
\newblock ISBN 978-1-4503-6125-5.

[In{\'a}cio~de Carvalho \emph{et~al.}, 2017]
In{\'a}cio~de Carvalho V, de~Carvalho M, Branscum AJ (2017).
\newblock \enquote{Nonparametric {B}ayesian {C}ovariate-{A}djusted {E}stimation
of the {Y}ouden {I}ndex.}
\newblock \emph{Biometrics}, \textbf{73}(4), 1279--1288.

[In{\'a}cio~de Carvalho \emph{et~al.}, 2013]
In{\'a}cio~de Carvalho V, Jara A, Hanson TE, de~Carvalho M (2013).
\newblock \enquote{{B}ayesian {N}onparametric {ROC} {R}egression {M}odeling.}
\newblock \emph{Bayesian Analysis}, \textbf{8}(3), 623--646.

[In{\'a}cio~de Carvalho and Rodr{\'i}guez-{\'A}lvarez, 2018]
In{\'a}cio~de Carvalho V, Rodr{\'i}guez-{\'A}lvarez MX (2018).
\newblock \enquote{Bayesian {N}onparametric {I}nference for the
{C}ovariate-{A}djusted {ROC} {C}urve.}
\newblock \emph{arXiv:1806.00473 [stat.ME]}.

[Ishwaran and James, 2001]
Ishwaran H, James LF (2001).
\newblock \enquote{Gibbs {S}ampling {M}ethods for {S}tick-{B}reaking {P}riors.}
\newblock \emph{Journal of the American Statistical Association},
\textbf{96}(453), 161--173.

[Ishwaran and James, 2002]
Ishwaran H, James LF (2002).
\newblock \enquote{Approximate {D}irichlet {P}rocess {C}omputing in {F}inite
{N}ormal {M}ixtures: {S}moothing and {P}rior {I}nformation.}
\newblock \emph{Journal of Computational and Graphical Statistics},
\textbf{11}(3), 508--532.

[Janes and Pepe, 2009]
Janes H, Pepe MS (2009).
\newblock \enquote{{A}djusting for {C}ovariate {E}ffects on {C}lassification
{A}ccuracy using the {C}ovariate-{A}djusted {R}eceiver {O}perating
{C}haracteristic {C}urve.}
\newblock \emph{Biometrika}, \textbf{96}(2), 371--382.

[Komsta and Novomestky, 2015]
Komsta L, Novomestky F (2015).
\newblock \emph{\pkg{moments}: Moments, cumulants, skewness, kurtosis and
related tests}.
\newblock \proglang{R}~package version 0.14,
\urlprefix\url{https://CRAN.R-project.org/package=moments}.

[Liu, 1996]
Liu JS (1996).
\newblock \enquote{Nonparametric {H}ierarchical {B}ayes via {S}equential
\newblock \emph{The Annals of Statistics}, \textbf{24}(3), 911--930.

[L{\'o}pez-Rat{\'o}n \emph{et~al.}, 2014]
L{\'o}pez-Rat{\'o}n M, Rodr{\'\i}guez-{\'A}lvarez MX, Cadarso-Su{\'a}rez C,
Gude-Sampedro F (2014).
\newblock \enquote{\pkg{OptimalCutpoints}: an \proglang{R} {P}ackage for
{S}electing {O}ptimal {C}utpoints in {D}iagnostic {T}ests.}
\newblock \emph{Journal of Statistical Software}, \textbf{61}(8), 1--36.

[MacEachern, 2000]
MacEachern SN (2000).
\newblock \enquote{Dependent {D}irichlet {P}rocesses.}
\newblock \emph{Unpublished manuscript, Department of Statistics, The Ohio
State University}, pp. 1--40.

[Maechler, 2019]
Maechler M (2019).
\newblock \emph{\pkg{nor1mix}: Normal aka Gaussian (1-d) Mixture Models (S3
Classes and Methods)}.
\newblock \proglang{R}~package version 1.3-0,
\urlprefix\url{https://CRAN.R-project.org/package=nor1mix}.

[Metz, 1978]
Metz CE (1978).
\newblock \enquote{{B}asic {P}rinciples of {ROC} {A}nalysis.}
\newblock \emph{Seminars in Nuclear Medicine}, \textbf{8}(4), 283 -- 298.

[Pardo-Fern{\'a}ndez \emph{et~al.}, 2014]
Pardo-Fern{\'a}ndez JC, Rodr{\'i}guez-{\'A}lvarez MX, van Keilegom I (2014).
\newblock \enquote{A {R}eview on {ROC} {C}urves in the {P}resence of
\newblock \emph{{REVSTAT}--Statistical Journal}, \textbf{12}(1), 21--41.

[Pepe, 1998]
Pepe MS (1998).
\newblock \enquote{Three {A}pproaches to {R}egression {A}nalysis of {R}eceiver
{O}perating {C}haracteristic {C}urves for {C}ontinuous {T}est {R}esults.}
\newblock \emph{Biometrics}, \textbf{54}(1), 124--135.

[Pepe, 2003]
Pepe MS (2003).
\newblock \emph{The {S}tatistical {E}valuation of {M}edical {T}ests for
{C}lassification and {P}rediction}.
\newblock Oxford Statistical Sciences Series. Oxford University Press, New

[Perez~Jaume \emph{et~al.}, 2017]
Perez~Jaume S, Skaltsa K, Pallar{\`e}s N, Carrasco~Jordan JL (2017).
\newblock \enquote{\pkg{ThresholdROC}: {O}ptimum {T}hreshold {E}stimation
{T}ools for {C}ontinuous {D}iagnostic {T}ests in \proglang{R}.}
\newblock \emph{Journal of Statistical Software}, \textbf{82}(4), 1--21.

[Plummer \emph{et~al.}, 2006]
Plummer M, Best N, Cowles K, Vines K (2006).
\newblock \enquote{CODA: Convergence Diagnosis and Output Analysis for MCMC.}
\newblock \emph{R News}, \textbf{6}(1), 7--11.
\newblock \urlprefix\url{https://journal.r-project.org/archive/}.

[\proglang{R} Core Team}, 2020]
{\proglang{R} Core Team} (2020).
\newblock \emph{\proglang{R}: {A} Language and Environment for Statistical
\newblock \proglang{R} Foundation for Statistical Computing, Vienna, Austria.
\newblock \urlprefix\url{https://www.R-project.org/}.

[Robin \emph{et~al.}, 2011]
Robin X, Turck N, Hainard A, Tiberti N, Lisacek F, Sanchez JC, M{\"u}ller M
\newblock \enquote{\pkg{pROC}: {A}n {O}pen-{S}ource {P}ackage for \proglang{R}
and \proglang{S+} to {A}nalyze and {C}ompare {ROC} {C}urves.}
\newblock \emph{BMC Bioinformatics}, \textbf{12}, 77.

[Rodriguez-Alvarez and Roca-Pardinas, 2017]
Rodriguez-Alvarez MX, Roca-Pardinas J (2017).
\newblock \emph{\pkg{npROCRegression}: Kernel-Based Nonparametric ROC
Regression Modelling}.
\newblock \proglang{R} package version 1.0-5,
\urlprefix\url{https://CRAN.R-project.org/package=npROCRegression}.

[Rodr{\'i}guez-{\'A}lvarez \emph{et~al.}, 2011]
Rodr{\'i}guez-{\'A}lvarez MX, Roca-Pardi{\~{n}}as J, Cadarso-Su{\'a}rez C
\newblock \enquote{A {N}ew {F}lexible {D}irect {ROC} {R}egression {M}odel:
{A}pplication to the {D}etection of {C}ardiovascular {R}isk {F}actors by
{A}nthropometric {M}easures.}
\newblock \emph{Computational Statistics \& Data Analysis}, \textbf{55}(12),

[Rodr{\'i}guez-{\'A}lvarez \emph{et~al.}, 2011]
Rodr{\'i}guez-{\'A}lvarez MX, Roca-Pardi{\~{n}}as J, Cadarso-Su{\'a}rez C
\newblock \enquote{{ROC} curve and {C}ovariates: {E}xtending {I}nduced
{M}ethodology to the {N}on-{P}arametric {F}ramework.}
\newblock \emph{Statistics and Computing}, \textbf{21}(4), 483--499.

[Rodr{\'i}guez-{\'A}lvarez \emph{et~al.}, 2011]
Rodr{\'i}guez-{\'A}lvarez MX, Tahoces PG, Cadarso-Su{\'a}rez C, Lado MJ
\newblock \enquote{Comparative {S}tudy of {ROC} {R}egression {T}echniques --
{A}pplications for the {C}omputer-{A}ided {D}iagnostic {S}ystem in {B}reast
{C}ancer {D}etection.}
\newblock \emph{Computational Statistics \& Data Analysis}, \textbf{55}(1),

[Rosenberg, 1995]
Rosenberg PS (1995).
\newblock \enquote{Hazard {F}unction {E}stimation {U}sing B-{S}plines.}
\newblock \emph{Biometrics}, \textbf{51}(3), 874--887.

[Sarkar, 2008]
Sarkar D (2008).
\newblock \emph{\pkg{Lattice}: Multivariate Data Visualization with
\proglang{R}}.
\newblock Springer-Verlag, New York.
\newblock ISBN 978-0-387-75968-5,
\urlprefix\url{http://lmdvr.r-forge.r-project.org}.

[Shapiro, 1999]
Shapiro DE (1999).
\newblock \enquote{The {I}nterpretation of {D}iagnostic {T}ests.}
\newblock \emph{Statistical Methods in Medical Research}, \textbf{8}(2),

[Silverman, 1986]
Silverman BW (1986).
\newblock \emph{Density Estimation for Statistics and Data Analysis}.
\newblock Chapman \& Hall/CRC.

[Tom{\'e} Mart{\'i}nez~de Rituerto \emph{et~al.}, 2009]
Tom{\'e} Mart{\'i}nez~de Rituerto MA, Botana MA, Cadarso-Su{\'a}rez C,
Rego-Iraeta A, Fern{\'a}ndez-Mari{\~{n}}o A, Mato JA, Solache I,
Perez-Fernandez R (2009).
\newblock \enquote{{P}revalence of {M}etabolic {S}yndrome in {G}alicia (NW
Spain) on {F}our {A}lternative {D}efinitions and {A}ssociation with {I}nsulin
\newblock \emph{Journal of Endocrinological Investigation}, \textbf{32}(6),

[Venables and Ripley, 2002]
Venables WN, Ripley BD (2002).
\newblock \emph{Modern Applied Statistics with S}.
\newblock Fourth edition. Springer, New York.
\newblock ISBN 0-387-95457-0,
\urlprefix\url{http://www.stats.ox.ac.uk/pub/MASS4}.

[Wand and Jones, 1994]
Wand MP, Jones MC (1994).
\newblock \emph{Kernel Smoothing}.
\newblock Chapman \& Hall/CRC.

[Wang, 2012]
Wang XF (2012).
\newblock \emph{\pkg{sROC}: Nonparametric Smooth ROC Curves for Continuous
\newblock \proglang{R} package version 0.1-2,
\urlprefix\url{https://CRAN.R-project.org/package=sROC}.

[Wickham, 2016]
Wickham H (2016).
\newblock \emph{\pkg{ggplot2}: Elegant Graphics for Data Analysis}.
\newblock Springer-Verlag New York.
\newblock ISBN 978-3-319-24277-4.
\newblock \urlprefix\url{https://ggplot2.tidyverse.org}.

[Youden, 1950]
Youden WJ (1950).
\newblock \enquote{Index for {R}ating {D}iagnostic {T}ests.}
\newblock \emph{Cancer}, \textbf{3}(1), 32--35.

[Zhou \emph{et~al.}, 2011]
Zhou XH, McClish DK, Obuchowski NA (2011).
\newblock \emph{Statistical Methods in Diagnostic Medicine}.
\newblock John Wiley \& Sons, New York.

[Zou \emph{et~al.}, 1997]
Zou KH, Hall WJ, Shapiro DE (1997).
\newblock \enquote{Smooth {N}on-{P}arametric {R}eceiver {O}perating
{C}haracteristic (ROC) {C}urves for {C}ontinuous {D}iagnostic {T}ests.}
\newblock \emph{Statistics in Medicine}, \textbf{16}(19), 2143--2156.

[Zou \emph{et~al.}, 1998]
Zou KH, Tempany CM, Fielding JR, Silverman SG (1998).
\newblock \enquote{Original {S}mooth {R}eceiver {O}perating {C}haracteristic
{C}urve {E}stimation from {C}ontinuous {D}ata: {S}tatistical {M}ethods for
{A}nalyzing the {P}redictive {V}alue of {S}piral {CT} of {U}reteral
\newblock \emph{Academic Radiology}, \textbf{5}(10), 680--687.

\end{thebibliography}

\newpage
\begin{appendix}
\section{Further computational tools for pooled ROC curve}\label{AppA}
\subsection{Quantile residuals}
We start by illustrating a further model check that can be helpful for models fitted using the function \code{pooledROC.dpm}, for instance, in deciding if $L_{D}=1$ or $L_{D}>1$ (the same obviously applies to the nondiseased group). It is well-known that for a continuous random variable, say $Y_D$, with CDF given by $F_D$, that $F_D(Y_D)\sim\text{U}(0,1)$. As a consequence, quantile residuals defined by $\widehat{r}_{Dj}=\Phi^{-1}\{\widehat{F}_{D}(y_{Dj})\}$, for $j=1,\ldots,n_D$, should follow, approximately, a standard normal distribution if a correct model has been specified [Dunn and Smyth, 1996]. A quantile-quantile (QQ) plot can then be used to determine deviations of the quantile residuals from the standard normal distribution. Below we provide the code to construct, for the diseased population, such QQ plot using output from the object \code{pROC\_dpm} (that assumed $L_{D}=L_{\bar{D}}=10$) obtained using the function \code{pooledROC.dpm}. The code for the nondiseased population follows in a similar manner, and is provided in the \proglang{R} replication code that accompanies this paper.
\begin{Schunk}
\begin{Sinput}
R> library("nor1mix")
R> traj <- matrix(0, nrow = pROC_dpm$mcmc$nsave,
+    ncol = length(pROC_dpm$marker$d))
R> lgrid <- length(pROC_dpm$marker$d)
R> grid <- qnorm(ppoints(lgrid))
R> for (l in 1:pROC_dpm$mcmc$nsave) {
+    aux <- norMix(mu = pROC_dpm$fit$d$Mu[l,],
+    sigma = sqrt(pROC_dpm$fit$d$Sigma2[l,]),
+    w = pROC_dpm$fit$d$P[l,])
+    traj[l, ] <- quantile(qnorm(pnorMix(pROC_dpm$marker$d, aux)),
+    ppoints(lgrid), type = 2)
+ }
R> l.band <- apply(traj, 2, quantile, prob = 0.025)
R> trajhat <- apply(traj, 2, mean)
R> u.band <- apply(traj, 2, quantile, prob = 0.975)

R> op <- par(pty = "s")
R> plot(grid, trajhat, xlab = "Theoretical Quantiles",
+    ylab = "Sample Quantiles", main = "Healthy population",
+    cex.main = 2, cex.lab = 1.5, cex.axis = 1.5)
R> lines(grid, l.band, lty = 2, lwd = 2)
R> lines(grid, u.band, lty = 2, lwd = 2)
R> abline(a = 0, b = 1, col ="red", lwd = 2)
R> par(op)
\end{Sinput}
\end{Schunk}
\end{appendix}
The resulting QQ plots are shown in Figure \ref{qqplots_dpm} and show virtually no deviations from the standard normal distribution quantiles, thus revealing a good fit of the DPM model that assumes 10 mixture components in both the diseased and nondiseased groups. In contrast, those QQ plots obtained when fitting a normal model in each group (i.e., $L_{D}=L_{\bar{D}} = 1$; model \code{pROC\_normal} in the main manuscript), clearly show some deviations from the assumed normal distribution quantiles (see Figure \ref{qqplots_dpm_normal}). The code used follows.
\begin{Schunk}
\begin{Sinput}
R> traj_normal <- matrix(0, nrow = pROC_normal$mcmc$nsave,
+    ncol = length(pROC_normal$marker$d))
R> lgrid_normal <- length(pROC_normal$marker$d)
R> grid_normal <- qnorm(ppoints(lgrid_normal))
R> for (l in 1:pROC_normal$mcmc$nsave) {
+    traj_normal[l, ] <- quantile(qnorm(pnorm(pROC_normal$marker$d,
+    pROC_normal$fit$d$Mu[l], sqrt(pROC_normal$fit$d$Sigma2[l]))),
+    ppoints(lgrid_normal), type = 2)
+ }
R> l.band_normal <- apply(traj_normal, 2, quantile, prob = 0.025)
R> trajhat_normal <- apply(traj_normal, 2, mean)
R> u.band_normal <- apply(traj_normal, 2, quantile, prob = 0.975)
\end{Sinput}
\end{Schunk}

\begin{figure}[h!]
\begin{center}
\subfigure[DPM model with 10 mixture components in each group]{\includegraphics[width=7cm]{rqrd.pdf}
\includegraphics[width=7cm]{rqrh.pdf}\label{qqplots_dpm}}
\subfigure[Normal model in each group]{\includegraphics[width=7cm]{rqr_normal_d.pdf}
\includegraphics[width=7cm]{rqr_normal_h.pdf}\label{qqplots_dpm_normal}}
\caption{Quantile residuals of the BMI data. Posterior mean quantile residuals versus the theoretical quantiles of the standard normal distribution. The circles represent the posterior mean quantiles over all posterior samples, while the dashed lines represent the corresponding 95\% credible bands. Top row: DPM model with 10 components in each group (model \code{pROC\_dpm} in the main manuscript). Bottom row: normal model in each group (model \code{pROC\_normal} in the main manuscript).}
\label{qqplots_dpm_global}
\end{center}
\end{figure}

\subsection{Bayesian bootstrap and kernel estimators of the pooled ROC curve}
The following code is used to fit the Bayesian bootstrap ROC curve estimator for the pooled ROC curve (function \code{pooledROC.BB}). The number of iterations considered is \code{B = 5000} and, for the sake of illustration, we also compute the partial area under the curve corresponding to true positive fractions (TPF) or sensitivities between $0.8$ and $1$.
\begin{Schunk}
\begin{Sinput}
R> set.seed(123, "L'Ecuyer-CMRG") # for reproducibility
R> pROC_BB <- pooledROC.BB(marker = "bmi", group = "cvd_idf",
+    tag.h = 0, data = endosyn, p = seq(0, 1, l = 101),
+    pauc = pauccontrol(compute = TRUE, focus = "TPF", value = 0.8),
+    B = 5000, parallel = "snow", ncpus = 2)

R> summary(pROC_BB)
\end{Sinput}
\begin{Soutput}
Call: [...]

Approach: Pooled ROC curve - Bayesian bootstrap
Area under the pooled ROC curve: 0.76 (0.74, 0.779)
Partial area under the specificity pooled ROC curve (Se = 0.8): 0.423
(0.384, 0.462)

Sample sizes:
Group H     Group D
Number of observations        2149         691
Number of missing data           0           0
\end{Soutput}
\end{Schunk}
Note that all partial areas' values returned are normalised and, as such, what is being reported, in this case, is
\[
\text{pAUC}_{\text{TPF}}(0.8)/(1-0.8).
\]
The estimated pooled ROC curve and AUC (posterior means), jointly with 95\% credible intervals, are obtained as follows
\begin{Schunk}
\begin{Sinput}
R> plot(pROC_BB, cex.main = 1.5, cex.lab = 1.5, cex.axis = 1.5, cex = 1.5)
\end{Sinput}
\end{Schunk}
The result is shown in Figure \ref{pROC_BB_plot}.

\begin{figure}[h!]
\begin{center}
\includegraphics[width=8cm]{pooled_ROC_BB.pdf}
\end{center}
\caption{Graphical results as provided by the \code{plot.pooledROC} function for an object of class \code{pooledROC.BB}. Posterior mean and 95\% pointwise credible band of the pooled ROC curve and corresponding AUC (posterior mean and 95\% credible interval).}
\label{pROC_BB_plot}
\end{figure}

\noindent We shall present now the syntax associated to the kernel estimator of the pooled ROC curve (function \code{pooledROC.kernel}). In terms of arguments, \code{bw} specifies how the bandwidth should be computed, with \code{SRT} standing for Silverman's rule of thumb and \code{UCV} for least squares cross-validation. Additionally, here \code{B} stands for the number of bootstrap replications used to compute the confidence intervals/bands. The syntax follows.
\begin{Schunk}
\begin{Sinput}
R> set.seed(123, "L'Ecuyer-CMRG") # for reproducibility
R> pROC_kernel <- pooledROC.kernel(marker = "bmi", group = "cvd_idf",
+    tag.h = 0, data = endosyn, p = seq(0, 1, l = 101), bw = "SRT",
+    B = 500, method = "coutcome",
+    pauc = pauccontrol(compute = TRUE, focus = "TPF", value = 0.8),
+    parallel = "snow", ncpus = 2)

R> summary(pROC_kernel)
\end{Sinput}
\begin{Soutput}
Call: [...]

Approach: Pooled ROC curve - Kernel-based
Area under the pooled ROC curve: 0.755 (0.737, 0.774)
Partial area under the specificity pooled ROC curve (Se = 0.8): 0.408
(0.373, 0.446)

Group H     Group D
Bandwidths:       0.867       1.019

Bandwidth Selection Method: Silverman's rule-of-thumb

Sample sizes:
Group H     Group D
Number of observations        2149         691
Number of missing data           0           0
\end{Soutput}
\end{Schunk}
Graphical results are obtained with the following code, and present in Figure \ref{pROC_kernel_plot}
\begin{Schunk}
\begin{Sinput}
R> plot(pROC_kernel, cex.main = 1.5, cex.lab = 1.5, cex.axis = 1.5,
+    cex = 1.5)
\end{Sinput}
\end{Schunk}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=8cm]{pooled_ROC_kernel.pdf}
\end{center}
\caption{Graphical results as provided by the \code{plot.pooledROC} function for an object of class \code{pooledROC.kernel}. Estimate and 95\% pointwise bootstrap confidence interval of the pooled ROC curve and corresponding AUC.}
\label{pROC_kernel_plot}
\end{figure}
We finish this section with a comparison of the estimated pooled ROC curves obtained using all methods incorporated in \pkg{ROCnReg} (Figure \ref{pROC_kernel_comp}).
\begin{Schunk}
\begin{Sinput}
R> plot(pROC_emp$p, pROC_emp$ROC[,1], type = "s", xlim = c(0,1),
+    ylim = c(0,1), xlab = "FPF", ylab = "TPF",
+    main = "Pooled ROC curve - Different approaches",
+    cex.main = 1.5, cex.lab = 1.5, cex.axis = 1.5)
+    lines(pROC_dpm$p, pROC_dpm$ROC[,1], col = 2)
+    lines(pROC_BB$p, pROC_BB$ROC[,1], col = 3)
+    lines(pROC_kernel$p, pROC_kernel$ROC[,1], col = 4)
+    abline(0, 1, col = "grey", lty = 2)
+    legend("topleft", legend = c("Empirical", "DPM", "BB", "Kernel"),
+    lty = 1, col = 1:4, bty = "n", lwd = 2, cex = 1.5)
\end{Sinput}
\end{Schunk}

\begin{figure}[h!]
\begin{center}