A Bayesian semiparametric approach to correlated ROC surfaces with stochastic order constraints

Diagnostic accuracy in the context of an ordinal true disease status has been studied by Nakas and Yiannoutsos (2004) , Xiong et al. (2006) , Chi and Zhou (2008) , Inacio et al. (2011) , and Kang and Tian (2013) , among others. In the example of 3-class disease status, let Y i = Y i 1 , … , Y i n i be measurements obtained from the i th population, where i = 1 , 2 , 3 and n i is the number of participants in the i th population. In the motivating example, these three populations are no disease (ND), mild disease (MD), and severe disease (SD). Assume further that participants in population 3 tend to have higher measurements than participants from population 2, who tend to have higher measurements than participants in population 1. Each participant is classified into one of the three populations based on the following decision rule ( Nakas and Yiannoutsos, 2004 ): for two ordered thresholds c 1 < c 2 , if Y ≤ c 1 then decision is “population 1”, if c 1 < Y ≤ c 2 then decision is “population 2”, otherwise decision is “population 3”. For a pair of thresholds c 1 , c 2 , the probability p i of correct classification into the i th population can be computed as follows: p 1 = P r Y 1 j ≤ c 1 = F 1 c 1 , j = 1 , … , n 1 , p 2 = P r c 1 ≤ Y 2 j ≤ c 2 = F 2 c 2 - F 2 c 1 , j = 1 , … , n 2 , p 3 = P r Y 3 j ≥ c 2 = 1 - F 3 c 2 , j = 1 , … , n 3 , where F i is the distribution of measurements for participants from population i , i = 1 , 2 , 3 . The ROC surface is defined as a function of p 1 and p 3 for all cutoff points ( c 1 , c 2 ) , with c 1 < c 2 , ROCS p 1 , p 3 = F 2 F 3 - 1 1 - p 3 - F 2 F 1 - 1 p 1 , if F 1 - 1 p 1 ≤ F 3 - 1 1 - p 3 , 0 , otherwise. To measure the overall diagnostic accuracy of the test, the volume under the ROC surface (VUS) is defined as V U S = ∫ 0 1 ∫ 0 1 R O C S p 1 , p 3 d p 1 d p 3 . Mossman (1999) has shown that V U S = P r Y 1 < Y 2 < Y 3 , meaning that the V U S is equal to the probability that three randomly selected measurements, one from each population, is classified correctly. The range of VUS is 1/6 to 1, with the two extremes attainable by an uninformative and a perfect test, respectively. In this paper, we take a Bayesian nonparametric approach to estimating F i , i = 1 , 2 , 3 . In particular, we use the Dirichlet process mixtures (DPM) ( Antoniak, 1974 ): (1) F = F y ; G = ∫ K y ; θ d G θ , G ~ D P α G 0 , where K ( y ; θ ) is a kernel with parameter θ , α and G 0 are precision and base parameters in a Dirichlet process (DP). Ferguson (1973) established the Dirichlet process and its existence, and also showed that draws from a DP are discrete with a probability of one. Sethuraman (1994) provided an alternative representation which makes the discreteness property of a DP explicit via a stick-breaking construction: if a random probability measure G is distributed according to a DP with α and base measure G 0 , then (2) G = ∑ k = 1 ∞ ϕ k δ θ k where the elements θ 1 , θ 2 , … are i i d realizations from G 0 , δ θ k is a point mass at θ k , and the ϕ k are weights constructed through the “stick-breaking” process ϕ k = ν k ∏ s < k 1 - ν s , with v k ~ Beta ( 1 , α ) . The use of DPM in equation (2) avoids the discreteness of DP by adding an additional convolution with a continuous kernel K ( y ; θ ) . This nonparametric approach allows flexibility in the distributions of the measurements that may exhibit skewness and multi-modality, and is motivated by the PRS data; moreover, as explained in the next subsection, DPM provides a convenient framework to incorporate the a priori constraints of interest, both in the center and in the variability of the distributions. In passing we note that DPM has been successfully applied to many statistical problems, including the estimation of ROC curves; see, for example, Erkanli et al. (2006) . Motivated by the PRS, we considered ordered multiple populations and multiple settings. Let Y i j = Y i j 1 , … , Y i j n i j T be a ( n i j × 1 ) response vector measured on i th population at j th setting, i = 1 , … , I and j = 1 , … , J . In the PRS, I = 3 for the populations of ND, MD, and SD, and J = 2 for settings 1 and 2, respectively. We developed scale and location mixture models for the distribution functions of each population, F i j ( y ) ≡ ∫ ∫ K ( y ; θ , σ ) d H i ( θ ) d G j ( σ ) , where { K ( y ; θ , σ ) : θ ∈ ℛ , σ ∈ ℛ + is a parametric family of distributions with a location parameter θ and a scale parameter σ . We assume that a gold standard exists for population classification in this section and will relax this assumption in Section 5 . We consider two order constraints simultaneously: stochastic order and variability order. First, for the j th setting, j = 1 , … , J , the response vector is stochastically ordered, i.e., Y 1 j ≤ s t Y 2 j ≤ s t … ≤ s t Y I j . That is, participants in the ND population is more likely to have lower values of Y than participants in the MD population, who tend to have lower measurements than participants in the SD population. Second, for the i th population, i = 1 , … , I , the response vector has a variability order between settings, i.e., Y i 1 ≥ v a r Y i 2 ≥ v a r … ≥ var Y i J , in the sense that participants are more likely to have smaller variability in measurement in setting J than J - 1 , in setting J - 1 than J - 2 and so on. For detailed definitions of stochastic and variability orders, see Gelfand and Kottas (2001) and Kottas and Gelfand (2002) . The DPM offers a natural framework for incorporating both the stochastic and variability order constraints when considering ROC curves ( Gelfand and Kottas, 2001 ; Hwang and Chen, 2015 ; Kottas and Gelfand, 2002 ). We now extend the result to the case of ROC surfaces. Theorem 1. Let S - ( h ) = s u p S - h u 1 , … , h u m denote the number of sign changes of a function h in U ∈ R where the supremum is taken with respect to all sets u 1 , … , u m , with u i ∈ U , m arbitrary but finite, and S - t 1 , … , t m is the number of sign changes of the sequence t 1 , … , t m exclusive of zero terms. Consider a parametric family of distributions K ( y ; θ , σ ) : θ ∈ ℛ , σ ∈ ℛ + , where θ is a location parameter, σ is a scale parameter, and ℛ + is the positive real line. Suppose kernel K ( y ; θ , σ ) is differentiable and strictly decreasing in θ and also differentiable in σ , such that for any σ 1 < σ 2 , S - K y ; θ , σ 2 - K y ; θ , σ 1 = 1 , the sign sequence being +,−, where the crossing point y 0 of K y ; θ , σ 2 and K y ; θ , σ 1 is the same for all σ 1 and σ 2 . For distribution H i on ℛ , and distribution G j on ℛ + , define the mixtures F i j y ; H i , G j ≡ ∫ ℛ + ∫ ℛ K ( y ; θ , σ ) d H i ( θ ) d G j ( σ ) , i = 1 , … , I , j = 1 , … , J . Then if H i ≤ s t H i + 1 , we have F i j y ; H i , G j ≤ s t F i + 1 , j y ; H i + 1 , G j for i = 1 , … , ( I - 1 ) , j = 1 , … , J , and if G j ≤ s t G j + 1 , we have F i j y ; H i , G j ≤ v a r F i , j + 1 y ; H i , G j + 1 for i = 1 , … , I , j = 1 , … , ( J - 1 ) . A proof of this theorem is provided in Supplementary Materials . Theorem 1 illustrates a way to construct multiple distributions satisfying both stochastic and variability orders through mixing. In particular, imposing stochastic orders with respect to the location and scale parameters of the kernel induces the stochastic and variability ordered mixture distributions. In practice, the Gaussian kernel K ( y ; θ , σ ) = Φ ( y ; θ , σ ) satisfies the conditions in Theorem 1 and hence can be used. Following Gelfand and Kottas (2001) , we choose product measures for H and G to construct mixing distributions that are stochastically ordered. In particular, we let H 1 = P 1 , H 2 = P 1 P 2 , H 3 = P 1 P 2 P 3 for the location parameter and G 1 = Q 1 Q 2 , G 2 = Q 1 for the scale parameter. This product measure construction enforces a sufficient but not necessary condition for stochastic ordering among multiple distributions and induces the space 𝒫 = H 1 , H 2 , H 3 , G 1 , G 2 : H 1 ≤ s t H 2 ≤ s t H 3 and G 2 ≤ s t G 1 where the desired stochastic and variability orders are achieved. Although the multiple populations at each setting can be assumed independent, the multiple settings within a given population are possibly correlated. This can happen if repeated measurements are taken from the same participants at different settings, as in the PRS data. In an early analysis of the PRS data, this correlation was not modeled explicitly ( Hwang and Chen, 2015 ). Many different approaches are available in the literature to induce correlations in DP, including ANOVA dependent DP (DDP) ( De Iorio et al., 2004 ), hierarchical DP ( Teh et al., 2006 ), and nested DP ( Rodriguez et al., 2008 ), among others. In this paper, we use the “common-weight” DDP introduced by MacEachern (1999) that, under the stick-breaking representation of DP in (2), allows locations of the point masses θ to be independent realizations from a stochastic process on covariates, and hence induces correlations. More specifically for the PRS data, we model the correlation by letting P j ~ D P α j P j 0 , for the j th population, j = 1 , … , 3 , where P j 0 = d B V N μ j , Σ j , a bivariate normal distribution with mean μ j and an arbitrary variance-covariance matrix Σ j . Non-zero estimates of the off-diagonal elements of Σ j allows correlated measurements in population j between the two settings in consideration.

In this section, we illustrate our approach through an application to data from the Physician Reliability Study (PRS) ( Schliep et al., 2012 ) that investigated measurement agreement among physicians in OB/GYN in diagnosing endometriosis. Although the PRS invited 12 physicians to conduct the diagnosis in 4 settings, we focus on the 4 regional experts (REs) in the first two settings. Moreover, we classify participants into 3 disease populations, defined as ND, MD, and SD; see Section 1 . This use of three instead of five ordinal stages eases computational burden and facilitates graphical presentation via 3-D surfaces. Substantively, grouping endometriosis stages I and II and separately III and IV is a common practice in clinical setting, as the latter two stages are usually indistinguishable under MRI ( Schliep et al., 2017 ). We account for two order constraints in the PRS. First, participants in the ND population are more likely to have lower values of rASRM score than participants from the MD population, who tend to have lower values than participants in the SD population. To account for this a priori information, we incorporate a stochastic order constraint for the rASRM score across the three populations at each of the two settings. Second, the rASRM score at setting 2 is more likely to have smaller variability than setting 1, as setting 2 provides more clinical information than setting 1. In PRS data, the stochastic order constraint is satisfied, and the variability order constraint is partially satisfied. In particular, both disease populations have smaller variability at setting 2 that 1, while the ND population has higher variability in setting 2. In the diagnosis of endometriosis, there does not exist a gold standard in practice that is universally accepted ( Hwang and Chen, 2015 ). To circumvent this situation, we follow Hwang and Chen to utilize the diagnoses of the four international experts (IEs) to construct an imperfect reference standard in the analysis. Details are provided in Supplementary Materials . We apply the four models specified in Section 5.2 of Supplementary Materials to the PRS data. Figure 1 shows the estimated predictive densities of logarithm of the rASRM scores from these models (95% credible bands are provided Figures S5 – S8 in Supplementary Materials ) This figure indicates strong evidence suggesting non-standard features in the distributions of the rASRM score after transformation, particularly the skewness and multimodality. Although it might be customary to model the logarithm transformed test scores as normally distributed, this figure suggests that such a practice might not be adequate here. For example, the rASRM scores after log-transformation are generally skewed to the right in the ND and MD populations and to the left in the SD population. Moreover, there might be evidence in supporting bimodality in these transformed scores, especially in the ND population at setting 1 and the MD population at setting 2. There is no substantial difference in the fitted densities among the four models. This is likely a consequence of the PRS data being mainly in concordance with the a priori constraints. The ROC surfaces are produced in Figures 2 and 3 for the NO and SO models and VO and JO models, respectively. The surface evaluates the diagnostic accuracy of the test for the three disease populations at all thresholds. The three coordinates correspond to the probabilities of correct classification into the three disease populations (see Section 2.1 ). Table 3 presents the estimated volume under the surface (VUS) measures associated with the estimated ROC surfaces for each of the four models at each of the two settings. As VUS ranges from 1/6 to 1, the estimates of all four models in Table 3 are generally large, indicating reasonably good performance of the REs in diagnosing endometriosis in both settings. Moreover, examining surgeons’ notes in addition to the intra-uteral images (Setting 2) can improve the diagnostic power considerably upon using images only (Setting 1). We note that the VUS estimates are not substantially different among the four models within each setting. This is not surprising since the PRS data are mainly in concordance with the order constraints. Stochastic order constraint does not have an impact on the estimates because the direction of the actual data agree perfectly with the stochastic order. The results are also not affected by variability order constraint because two of the disease populations in PRS (MD and SD) are in concordance with the variability order. While the rARSM scores from the ND population do indicate an opposite direction in variability between the two settings compared to the a priori variability constraint, its effect might be dominated by that from the other two diseased populations. The net effect of the two forces results in VUSs that are always larger in setting 2 than setting 1. Finally, we observe that the JO model has narrower credible intervals compared to the other three models. This might suggest that the interplay between SO and VO can produce large efficiency gains. We note that the VUS estimates are quite similar (see Table S1 in Supplementary Materials ) when we vary values of a l and b l , l = 1 , 2 , of the inverse gamma distribution (notations defined in Section S3 in Supplementary Materials ), suggesting robustness of the estimates in Table 3 to hyperparameters. In sensitivity analysis for no gold standard, we use various γ 0 , γ 1 , λ 0 , and λ 1 (notations defined in Section S8 in Supplementary Materials ) over a wide range of reasonable values: γ 0 = λ 0 = - 2 , - 3 , - 4.6 , - 6 , - 7 , and γ 1 = λ 1 = 6 , 8 , 9 , 10 , 12 . The estimates of VUS measures are found to be close to those from our primary analysis.

Compared to existing literature, our paper advances the development of order constrained analysis from ROC curves to surfaces. This is important as in many diagnostic accuracy studies, the true disease status takes an ordinal form, and the practice of simply combining groups to construct a dichotomous standard might result in a loss of information and even distort the underlying distributions. We also adequately model the correlation between test scores when the tests are conducted on the same study participants repeated over time, as is the case in the PRS. In addition to the two types of constraints considered in this paper, there might be other constraints that suit a particular application in hand. For example, one can consider a weaker constraint in the form of stochastic precedence (see, e.g., Chen and Dunson (2004) and Kottas (2011) ). Moreover, the constraints considered so far are in terms of the distributions of the test scores. It is possible that some a priori constraints are more appropriate in terms of the ROC surfaces or the associated VUS. How to account for these constraints is a current research topic of the authors. Model selection is an integrated and important component of statistical data analysis. Yet care has to be taken in the context of modeling a priori constraints. In this paper, we found that the diagnostic accuracy measures were all estimated similarly for the four models considered. Questions can be asked on which of the four models we should select for this particular dataset. To answer this, we need to pause to think about the original reasons for conducting order constrained statistical analysis. The fact that a constraint is believed to be true a priori indicates that the belief originates subjectively, independent of the empirical data. As such, it would be an inconsistent practice if we subject the a priori belief to the empirical evidence and use a single realized dataset from the underlying stochastic process to test its validity. In general methodological development involving a priori constraints, the use of data to guide model selection in relation to these constraints should not be advised. Philosophical considerations aside, even when data agree with the constraints so that all four models return similar point estimates of model parameters, the model with the believed a priori constraints will generally produce narrower credible intervals of the estimates, a well-known feature of the constrained statistical analysis. When one firmly believes in the order constraints, the appropriate approach is to fit the data with a model with the constraints, regardless whether the data agree with the constraints. Methods have been proposed when the constraints are thought to be “soft” in the sense that the prior probability that the constraint is true is less than one; see Danaher et al. (2012) . In general it is of interest to compare the proposed approach with some standard ones, such as those implemented in R package DiagTest3Grp ( Xiong et al., 2006 ; Luo and Xiong, 2012 ). However, the inability of the standard approaches to accommodate stochastic order constraints and the lack of a gold standard in the Physician Reliability Study data make this comparison difficult. Instead, we answer two related questions to provide some insights. One, what is the gain in using the proposed approach compared to one that does not accommodate the constraints? And two, what is the gain in using DPM compared to a simpler parametric approach? In both the simulation studies and real data analysis, we saw that incorporation of constraints can improve statistical efficiency of the VUS estimates when constraints agree with data. In cases where constraints do not agree with data, the use of them can change the point estimates dramatically, as our simulation studies demonstrated. To answer the second question, we applied a Tri-Normal model using the DiagTest3Grp package to simulation data that are generated according to scenario IX in Table S2 in the Supplementary Materials and compared the results to those from the proposed approach under the NO (no-constraint) model. Results in Table S3 in Supplementary Materials suggest that the posterior estimates from these two approaches are generally close to each other and to the true values. The take-home messages from answers to these two questions are as follows. First, if researchers do not possess any a priori order constraints, both the proposed and standard approaches will produce similar results in VUS estimates. However, in situations where a priori order constraints are available, standard approaches will produce undesired results, either in point estimates that do not reflect the extra prior information or in loss in statistical efficiency. Although it was illustrated with the endometriosis example, the proposed approach can be readily applied to other medical research situations. For example, in obstetrics, it is common to monitor fetal development by collecting repeated anthropometric ultrasound measurements of fetus at different weeks of gestation. The predictive power of these ultrasound measurements in relation to some birth outcomes, e.g., birth size categorized by small for gestational age (SGA, birth weight less than 10th percentile), normal for gestational age (NGA, birth weight in between 10th and 90th percentiles), and large for gestational age (LGA, birth weight greater than 90th percentile), are of considerable interest ( Liu and Albert, 2014 ). In this context, it is reasonable to assume a priori that the anthropometric ultrasound measures in the earlier weeks are less accurate, therefore more variable, than those in later weeks, given that fetuses at early gestation ages are small and hence difficult to take accurate measures. See Foster et al. (2017) for a visualization of this example. By formularizing this assumption into a stochastic variability constraint and incorporating it into the estimation of the predictive powers of ultrasound at different weeks of gestation, we force the estimated variability to conform with the a priori restriction when data do not exhibit such agreement or gain statistical efficiency when they do.

To evaluate the performance of our proposed approach, we simulated data that mimic the PRS with a 3-class true disease status and 2 settings. We considered four different variability order directions between the two settings, varying from totally in concordance with the a priori variability constraint to totally in discordance with it. In all these cases, the data were generated to always agree with the a priori stochastic order constraint. For each order direction, we considered two different distributions for the test scores (normal and a mixture of normals) and two sample sizes (120 and 250). The resulting eight scenarios are presented in Table 1 . For every scenario, 100 datasets were generated, with each fit by one of the four models: (i) No order constraint (NO), (ii) Stochastic order only (SO), (iii) Variability order only (VO), and (iv) Joint order (JO) where both stochastic and variability orders are considered. The results are reported in Figures S1 – S4 in Supplementary Materials for predictive density estimates and in Table 2 for VUS estimates. Supplementary Materials provide more details of the simulation studies, including data generation, estimations, and results interpretations. Given that there is no substantial difference in estimates in Table 2 between the two distributional assumptions, we concentrate on the results based on the mixture of normal cases (scenarios V to VIII). Overall, we see similar patterns of results whether the sample size is 120 or 250, although a slight efficiency gain can be noticed in the larger sample size cases. VUS estimates from each of the four models depend on the degree of concordance between the simulated data and the order constraints. As the NO model does not pose any constraint, and the SO model is always satisfied in the generated data, the resulting estimates are close to each other and are close to the true VUSs. The VO and JO models also produce VUS estimates that are close to the truth when there is no or low discordance between data and the variability constraint (scenarios V and VIII, respectively). However, when there is reasonable discordance (scenarios VII and VI), the VUS estimates deviate significantly from the truth. This is expected and is a consequence of the spreading-shrinking process in the predictive density estimations in response to the variability order constraint. To see an example, consider scenario VI with N = 120 where all three populations have larger variability in setting 2 than 1, with the implied VUSs under the data simulation specifications at 0.65 and 0.55 for settings 1 and 2, respectively. When the VO and JO models are applied to the data generated in this scenario, the variability order constraint forces the distributions in all three populations to satisfy the constraint. As a net result of this, the VUS estimates under the VO model now become 0.54 and 0.62 for settings 1 and 2, respectively, and 0.49 and 0.58 under the JO model. We also note that compared to those from NO and SO models, the VUS estimates from VO and JO tend to have slightly narrower 95% credible intervals in scenarios I and V, suggesting efficiency gains when data agree with constraints. More pronounced efficiency gains can be obtained when we consider a scenario similar to V but with smaller true VUS’s. The data generating mechanism and corresponding estimation results are presented in Tables S2 (scenario X) and S4 in the Supplementary Materials , respectively. The numbers in Table S4 demonstrate a clear improvement in estimation efficiency as more constraints are considered. For example, for setting 1 under sample size of 120, the VUS estimate from the JO model has a 95% credible interval width of 0.17, compared to 0.23 from the NO model. This corresponds to a 25% efficiency gain in VUS estimate when both stochastic and variability order constraints are considered compared to when none is.

The corresponding likelihood function is given in Supplementary Materials where we also provide details on how an equivalent hierarchical model specification can be constructed to alleviate the computational challenge in evaluating the multiple integrals. The proposed model can be fitted using the MCMC algorithm with Gibbs sampling ( Casella and George, 1992 ) and Metropolis-Hastings steps ( Chib and Greenberg, 1995 ). We use a truncation approximation for the DP priors ( Ishwaran and Zarepour, 2000 ) and let P i = ∑ h i = 1 M p h i δ θ h i * , Q j = ∑ l j = 1 L q l j δ σ l j * , i = 1 , 2 , 3 , j = 1 , 2 , where truncation levels M and L can be determined a priori at a moderately large number (e.g., 30). Although M and L can be specified alternatively, e.g., by examining the moments of the random weights, pre-specifying them at moderately large values provides a simple and satisfactory approximation. More details regarding our sampling algorithm including full conditional posterior distributions can be found in Supplementary Materials . At each MCMC iteration, we obtain draws of model unknowns p , q , θ , σ , μ , Σ and compute measures P i s and Q j s. These quantities are then used to obtain the distribution functions F i j ( y ) , i = 1 , 2 , 3 , j = 1 , 2 . ROC surfaces can be fitted by R O C S j = F 1 j c 1 , F 2 j c 2 - F 2 j c 1 , 1 - F 3 j c 2 , - ∞ < c 1 < c 2 < ∞ , where ( c 1 , c 2 ) are pre-specified values, usually over a grid spanning the range of Y (e.g., 0 to 5 by 0.05 for the example in Section 5 where logarithm transformation is taken on Y ). At each iteration, we also obtain the volume under the surface (VUS) measure using V U S = P r Y 1 < Y 2 < Y 3 where Y 1 , Y 2 , Y 3 are sampled from predictive densities of the corresponding populations at each iteration ( Mossman, 1999 ). With T iterations, we obtain a sample of size T from the posterior distribution which we can then use to make inference using, e.g., posterior means and credible intervals.

In many epidemiological and biomedical studies, there is interest in examining whether a specific test can diagnose a disease status better than competing tests. The receiver’s operating characteristic (ROC) technique provides a very useful way to visualize and summarize data evidence when the test scores are continuous. In these applications, it is possible that some a priori information exists regarding the test score distributions, either between different disease populations for a single test or between multiple correlated tests. There is little work in the literature that considers constrained diagnostic accuracy analysis when the true disease status is binary, with the exception of Gelfand and Kottas (2001) , Kottas and Gelfand (2002) , and Hanson et al. (2008) . The lack of existing work is more pronounced when the disease status is ordinal, as no such prior work can be found in the literature to the best knowledge of the authors. We will present an interesting example in this section to motivate this research, and will discuss other applications in Section 6 . Endometriosis is a disease that affects women in which tissue that normally lines the inside uterus grows outside into other parts of the body. Symptoms of endometriosis include pain, sometimes severe, excessive bleeding, and infertility. Large variation exists in endometriosis diagnosis, which in turn leads to under-diagnosis and under-treatment of this disease. To study the accuracy of endometriosis diagnosis, the Physician Reliability Study (PRS) invited physicians in OB/GYN to review clinical information of about 150 women ( Schliep et al., 2012 ). These physicians were presented with different combinations of clinical information and were asked to provide the revised American Society of Reproductive Medicine (rASRM) score, a continuous quantity that is algorithmically computed based on multiple observations of women’s uteruses, including depth of endometriosis in peritoneum, degree of obliteration of posterior cul de sac, and appearance of adhesion, among others. A higher value of rASRM is indicative of higher likelihood of endometriosis or more advanced stages. The physicians conducted the diagnoses in four successive settings: setting 1 had participants’ intra-uteral digital images as the sole clinical information, setting 2 had both the participants’ intra-uteral digital images and surgeon’s notes taken in an early laparoscopy session, setting 3 had the clinical information in setting 2 plus MRI results of the participants, and setting 4 had the clinical information in setting 3 plus the participants’ histopathology reports. With the aim of estimating diagnostic accuracy of the physicians under different settings, the PRS design presents a natural order of sequentially increased amounts of clinical information available to the physicians for diagnoses, and is believed to pose an a priori constraint where the resulting rASRM score should have smaller variability in a later setting. Such information, if incorporated adequately, can potentially enhance statistical efficiency when estimating the diagnostic accuracy measures. An early work analyzed the PRS data assuming a dichotomous true disease status (No Disease/Disease), focusing on the data from four regional experts in the first two settings ( Hwang and Chen, 2015 ). While of interest itself, consideration of the dichotomous true disease status, which can be regarded as a collapsed version of the ordinal endometriosis disease staging (no disease, stages I through IV), may not describe the full picture of the underlying diagnostic performance and can be suboptimal due to a loss of information. Furthermore, when considering a priori constraints, as we propose to do in this paper, this data collapsing practice can lead to potentially distorted representation of the rASRM distributions. For example, when the dichotomous true disease status is considered, the PRS data present a larger variability in setting 2 than setting 1 in both the No Disease and Disease populations, contrary to the a priori variability constraint that, as more clinical information becomes available, the rASRM scores should have smaller variability. When a 3-class ordinal true disease status—No Disease (ND), Mild Disease (MD, stages I and II) and Severe Disease (SD, stages III and IV)—is instead considered, the data actually agree with the a priori variability constraint in two of the three populations. This suggests that we may obtain a more realistic picture of the diagnostic performance by using the ordinal rather than the dichotomous true disease status. For these reasons, it is desirable to approach the diagnostic accuracy question in the context of an ordinal true disease status. The overarching goal of this paper is to develop a framework for estimating diagnostic accuracy measures of the OB/GYNs for diagnosing endometriosis under several particular features of the PRS data: (i) the true disease status takes the form of a 3-class ordinal variable, (ii) there is a stochastic ordering relationship in the rASRM scores between different disease populations within a setting, in the sense that the rASRM scores tend to be larger in the MD population compared to the ND population and larger in the SD population than in the MD population, (iii) there is also a variability order between different settings within a disease population, in the sense that the rASRM scores tend to have smaller variability in setting 2 compared to setting 1, and (iv) the rASRM scores between different settings are potentially correlated as they are from the same participant. In the PRS data, the overall sample correlation in rASRM scores between settings 1 and 2 is 0.81. To accommodate these features, we propose to estimate the receiver operating characteristic (ROC) surfaces ( Mossman, 1999 ; Nakas and Yiannoutsos, 2004 ; Xiong et al., 2006 ) and the associated volume under the surface (VUS), a summary measure similar to the area under the curve measure for ROC curves. This will be achieved by using Dirichlet process mixtures to model the rASRM score distributions flexibly, including their skewness and multi-modality characteristics. To account for the between-setting correlations in rASRM scores, we model the pair of distributions (one pair for each disease population) by considering dependent Dirichlet process mixtures that have bivariate distributions with non-zero correlations as their bases. Finally we consider the nice property of stochastic ordering in the context of mixture distributions ( Gelfand and Kottas, 2001 ; Kottas and Gelfand, 2002 ) to accommodate ordering constraints. Compared to existing literature, this paper makes several contributions. First, we extend the methodological framework of Hwang and Chen (2015) to the case of an ordinal true disease status. In particular, we prove that the ordering constraints can be achieved by product measures in the framework of Dirichlet process mixtures when ROC surfaces are of interest. We are not aware of any other work in the literature that have addressed order constrained analysis in estimating ROC surfaces or VUSs. Second, as the rASRM scores between settings are from the same physicians for the same study participants, and hence are likely correlated, we directly and explicitly model this correlation in our proposed framework. Substantively, we show that the ordinal true disease status is a more reasonable standard to consider in PRS data than its dichotomous counterpart. This paper is organized as follows. Section 2 will outline the methodological framework, including the ROC surfaces and the incorporation of order constraints through Bayesian Dirichlet process mixtures. Correlations between multiple ROC surfaces are also considered in this section. Section 3 specifies the computational procedures and delineates the detailed steps in the Markov chain Monte Carlo algorithm. In Section 4 , we conduct a comprehensive simulation analysis to evaluate the performance of our proposed approach. Section 5 analyzes the PRS data in detail and interprets the estimated results. Finally, Section 6 concludes with summaries and future research directions.

Additional supporting information may be found online in the Supporting Information section at the end of the article. Web Appendices referenced in Sections 2.2 , 3.1 , 4 , 5.2 , 5.3 , and 6 and R code to implement the method are available with this article at the Biometrics website on Wiley Online Library.