An Integrated Bayesian Nonparametric Approach for Stochastic and Variability Orders in ROC Curve Estimation: An Application to Endometriosis Diagnosis

The receiver operating characteristic (ROC) curve is a graphical tool to assess the performance of continuous diagnostic tests in distinguishing between diseased and healthy populations. The theoretical ROC curve is a plot of all possible pairs of true positive rate (sensitivity) versus false positive rate (1-specificity) across all threshold values. A useful index of diagnostic accuracy is the area under the ROC curve (AUC), which can be interpreted as the probability that a randomly selected diseased subject has a higher value than a randomly selected healthy subject. Models for ROC curves are abundant, and range from parametric ( Hanley, 1996 ; Pepe, 1998 ; Faraggi and Reiser, 2002 ; Sorribas, March and Trujillano, 2002 ; Pepe, 2003 ) to nonparametric ( Hsieh and Turnbull, 1996 ; Lloyd, 1998 ; Pepe, 2000 ; Pepe, 2003 ; Peng and Zhou, 2004 ; Cai and Moskowitz, 2004 ), to Bayesian ( Peng and Hall, 1996 ; Gu, Ghosal and Roy, 2008 ; Gu and Ghosal, 2009 ), and cover both the case with and without a gold standard for the disease ( Jafarzadeh et al., 2010 ; Tang et al., 2013 ). More recently, nonparametric Bayesian approaches have also been developed ( Branscum, Gardner and Johnson, 2005 ; Erkanli et al., 2006 ; Branscum et al., 2008 ; Hanson, Branscum, and Gardner, 2008 ; Cheng et al., 2012 ; De Carvalho et al., 2013 ; Rodriguez and Martinez, 2014 ). In the application of ROC curves, it’s often the case that some a priori constraints may exist. It has been well established in the literature that when such constraints are available, incorporating them in the estimation can improve the statistical efficiency of the estimates and even result in different estimates that conform with biological theory or expert consensus ( Robertson, Wright and Dykstra, 1988 ). One commonly encountered constraint is that the diseased population tends to have larger test scores on average than the healthy population. Hanson, Kottas and Branscum (2008) modeled this constraint by a stochastic order in the healthy and diseased distributions. A less studied constraint may occur in multiple tests in that one test score tends to have larger variability than another. While these two orders have been studied separately in ROC curves estimation ( Gelfand and Kottas, 2001 ; Kottas and Gelfand, 2002 ; Dunson and Peddada, 2008 ; Hanson, Kottas and Branscum, 2008 ), joint modeling of them does not exist to our best knowledge. This paper proposes an integrated approach to jointly model stochastic and variability orders. This research was motivated by the Physician Reliability Study (PRS) ( Schliep et al., 2012 ) that examined the diagnostic accuracy of diagnosing endometriosis, a women’s health disorder that occurs when cells from the lining of the uterus grow in other parts of the body. With 12 physicians in obstetrics and gynecology (OB/GYN) invited to review clinical information of 156 participating women, the PRS provides the revised American Society for Reproductive Medicine (rASRM) score that ranges from 1 to 159. Higher rASRM score is indicative of higher likelihood of an endometriosis disease. The 12 physicians consist of three groups of obstetricians and gynecologists: the 4 international experts (IEs) are physicians with international reputation in the field and have considerably extensive experience in laparoscopy and endometriosis diagnosis, the 4 regional experts (REs) are physicians affiliated with Utah Medical Center with extensive experience in the field, and the 4 residents (RDs) are new graduates from medical schools conducting guided practice at Utah Medical Center. The reviewing and diagnosis were conducted at four different settings, with setting 1 having clinical information consisting of only the digital intrauterus image taken during a laparoscopy procedure. At setting 2, the clinical information was augmented to include surgeon’s notes. MRI and histopathology reports were further successively added to the clinical information pool at settings 3 and 4, respectively. In ordering the four settings this way, the PRS attempts to follow the usual practice in the field: a laparoscopy is usually done and the surgeon provides a diagnosis of endometriosis upon seeing the inside of the uterus (settings 1 and 2). When ambiguity still exists, other tests such as MRI (setting 3) and histopathology (setting 4) are then conducted to aid the diagnosis. One research question of interest is to estimate the different diagnostic accuracy of the group of 4 regional experts at settings 1 and 2. These REs are of interest because they represent everyday physicians in practice; we focus on settings 1 and 2 as they are the most commonly used steps in diagnosing endometriosis and as only a small number of women have MRI and histopathology results. In trying to estimate the two ROC curves of the REs at settings 1 and 2, two constraints were brought up by our expert collaborators. First, within each setting, the rASRM scores should be larger in the diseased population than in the healthy population. To reflect this feature, a stochastic order constraint is assumed between the two populations. Secondly, with the study design that more clinical information is available at setting 2 than 1, the variability in the rASRM scores within each population should be smaller at setting 2 than 1. According to our collaborators, this variability order constraint is sensible because the physicians are able to make fewer guesses in checking the various components of the rASRM scoring sheet that produced the rASRM score when these physicians have more clinical information. In this paper, we use the variability order notion proposed by Kottas and Gelfand (2002) due to its nice connection to the Bayesian nonparametric framework. Our interest is to account for these two a priori constraints simultaneously when we estimate the ROC curves for the group of REs. There are other features of the PRS data that make the analysis challenging. First, the rASRM scores can be non-normal, even after taking the logarithm transformation (see Figure 1 ). This warrants some flexible models for the scores that can accommodate the skewness and possible bimodality. To handle this feature, we propose to model these scores using a Bayesian nonparametric approach. More specifically, we will utilize the Dirichlet process mixture (DPM) for the rASRM score on the logarithm scale so that it affords the needed flexibility in the distributions. The DPM is also desirable as it enables us to extend the work of Gelfand and Kottas (2001) and Kottas and Gelfand (2002) to jointly model the stochastic and variability orders. Another challenge in the analysis is that there does not exist a gold standard for en-dometriosis in the PRS data. To address this issue, we exploit the fact that four of the 12 physicians are international experts (IEs) who tend to have greater experience in diagnosing endometriosis. The diagnostic results of these 4 IEs can be exploited to create an imperfect reference standard. We will propose a sensitivity analysis scheme and study the robustness of this mechanism by varying the sensitivity parameters in a scientifically reasonable way. This is an extension of the sensitivity approach introduced by Zhang, Chen and Albert (2012) . The rest of the paper is organized as follows. In Section 2, we take an initial look at the PRS data. In Section 3, we introduce the ROC curve, order constraints and Dirichlet process mixture models, and propose our integrated model of stochastic and variability orders using nonparametric Bayesian methods. Section 4 specifies priors and presents posterior computations of models developed in Section 3 using Markov chain Monte Carlo (MCMC) methodology. In Section 5, we apply our model to the PRS data and discuss the substantive results. In this section, we also consider the sensitivity analysis approach that exploits IEs’ imperfect reference standard to address the issue of not having a gold standard. In Section 6, we conduct simulation studies to illustrate the proposed model and to evaluate their performance. We conclude with a discussion of our proposed approach and future work in Section 7.

In the PRS, the 12 invited physicians consist of 4 international experts (IEs), 4 regional experts (REs) and 4 residents. In this paper, we focus on the diagnostic accuracy of the REs in diagnosing endometriosis at settings 1 and 2, separately, using their average rASRM scores. Averaging REs’ scores is reasonable here because the 4 REs are thought to be exchangeable as they all have similar years of experience and have been affiliated with the same institute. Out of the 156 women participants in the PRS, 113 at setting 1 and 126 at setting 2 have non-missing data in the average REs’ rASRM scores and form our working datasets. Table 1 , Figures 1 and 2 provide descriptive statistics of the rASRM scores based on these data. Several observations can be made: The rASRM scores within each setting are generally skewed and possibly have multimodality. The logarithm transformation does not help in achieving symmetry or uni-modality ( Figure 1 ); These scores are still non-symmetric and multi-modal within each setting by disease status ( Figure 2 ), where disease status is estimated by posterior modes from models specified in Section 5; Healthy and diseased populations combined, the scores are lower and have smaller variance in setting 1 than 2 ( Table 1 ); Stratified by disease status, the scores are higher in their means in the diseased than the healthy population in both settings; within each population, the variances are higher in setting 2 than 1 ( Table 1 ). Although there is no gold standard in the PRS for diagnosing endometriosis, the 4 IEs’ diagnosis can be exploited to provide an imperfect reference standard. In Section 5, we propose a sensitivity analysis approach to make use of the IEs’ diagnostic results to aid the estimation of the REs’ ROC curves. The 4 IEs made 1/0 (presence/absence) endometriosis diagnosis at each setting for each woman and we use their proportion of positive diagnosis in the sensitivity analysis (see Section 5 for details). The relationship between REs’ rASRM scores and the IEs’ proportion of positive diagnosis is presented in Figure 3 for settings 1 and 2, separately and combined. We observed that in general, higher REs’ scores are associated with higher proportion of IEs’ positive diagnosis, as expected. As an initial look, we fit both a binormal model ( Pepe, 2003 ) and a bi-DPM model where scores from the healthy and diseased populations are each modeled using a Dirichlet process mixture ( Erkanli et al., 2006 ). The “gold standard” for this initial analysis is constructed by IEs’ proportion of positive diagnosis in settings 1 and 2 combined as follows: a woman is assigned a presence of endometriosis if the proportion is greater than 0.5, and absence if less than 0.5. In the case of a proportion of 0.5, a one-time random draw from Bernoulli (0.5) is made and considered as her disease status. The results are demonstrated in Figure 4 , where both ROC curves and their corresponding AUCs are presented. From the initial analysis, we saw that setting 1 has a higher AUC than setting 2 in both the binormal and bi-DPM models. AUCs in the bi-DPM models tend to be smaller than their counterparts in the binormal models, but this is likely due to the different distributional assumptions used to model the rASRM scores. This initial analysis does not incorporate the constraints and assumes a known gold standard. The following two sections propose methodologies to address these two issues.

Let F H and F D denote the cumulative distribution functions (CDFs) of healthy and diseased populations, respectively. The ROC curve is defined by ROC ( u ) = 1 - F D { F H - 1 ( 1 - u ) } , for 0 ≤ u ≤ 1. A useful index of diagnostic accuracy is the area under the ROC curve (AUC), which has the definition of AUC = ∫ 0 1 ROC ( u ) d u . Bamber (1975) has shown that AUC = Pr ( Y D > Y H ), where Y D and Y H denote the measurements for diseased and healthy populations, respectively. A binormal ROC curve is a widely used parametric method, and assumes that test results follow normal distributions for the healthy and diseased populations after an unspecified transformation S: S ( Y H ) ~ N (0, 1) and S ( Y D ) ~ N ( μ , σ 2 ). Then, the functional form of the binormal ROC curve is ROC ( u ) = Φ( a + b Φ −1 ( u )), where a = μ / σ , b = 1/ σ and Φ denotes the standard normal cumulative distribution function. Also, the AUC for the binormal ROC curve is AUC = Φ ( a / 1 + b 2 ) ( Pepe, 2003 ). However, this binormal model can be overly restrictive in practice. If the data are not normally distributed, it is desirable to use more flexible class of distributions. The Dirichlet process (DP) is a commonly used prior in a nonparametric Bayesian framework when the distribution exhibits more flexibility in location and shape, which may not be well captured by parametric models. Suppose a real-valued continuous random variable Y i , i = 1, …, n , is sampled from an unknown distribution F . We could place a DP prior on the distribution F ( Ferguson, 1973 ). The DP with precision α and base parametric distribution F 0 , denoted by F ~ DP ( αF 0 ), can be described by the stick-breaking representation ( Sethuraman, 1994 ): F = ∑ h = 1 ∞ p h δ θ h ∗ , p h = V h ∏ l = 1 h - 1 ( 1 - V l ) , V h ~ Beta ( 1 , α ) , θ h ∗ ~ F 0 , where δ θ * is a point mass at θ * and all V h ’s and θ h ∗ ’s are independent. The θ h ∗ ’s are random atoms which represent all possible values of Y i . The corrersponding probability weights ( p 1 , p 2 , …) are generated from a stick-breaking procedure using independent Beta distributed V h ( h = 1, 2, …). The base measure F 0 corresponds to our best guess for unknown distribution F . The precision parameter α > 0 controls the variability in this guess for base measure, with F approaching F 0 as α → ∞. The unknown distribution F can be interpreted as a discrete distribution comprised of a countably infinite number of distinct values according to random weights. One simple extension to remove the discreteness restriction is a Dirichlet process mixture (DPM) model ( Antoniak, 1974 ). The DP can be employed as a prior for the mixture distribution, e.g., F ( y; F 1 ) = ∫ K ( y; θ ) dF 1 ( θ ), F 1 ~ DP ( αF 0 ). Then, the discretness constraint of DP can be solved by choosing a continuous distribution for the kernel, K ( y; θ ). In the application of ROC curves, a bi-DPM model can be constructed if we assume the healthy and diseased populations each follow a separate DPM ( Erkanli et al., 2006 ). In comparing two populations, let F 1 and F 2 be two population distributions on θ ∈ . The stochastic order of two populations is defined as follows: F 2 is stochastically larger than F 1 if F 1 ( c ) ≥ F 2 ( c ) for all c , denoted by F 1 ≤ st F 2 . Gelfand and Kottas (2001) showed how stochastic order can be retained through mixture distributions. Consider a parametric family of distributions { K ( y; θ ) : θ ∈ }, with K ( y; θ ) differentiable and strictly decreasing in θ , and mixture distribution F ( y; F i ) = ∫ K ( y; θ ) dF i ( θ ), i = 1, 2. Then, if F 1 ≤ st F 2 , F ( y; F 1 ) ≤ st F ( y; F 2 ). A convenient choice of F 1 and F 2 satisfying F 1 ≤ st F 2 is F 1 = G 1 and F 2 = G 1 G 2 , where G 1 and G 2 are distribution functions on ( Gelfand and Kottas, 2001 ; Kottas and Gelfand, 2002 ; Hanson, Kottas and Branscum, 2008 ). Then, the space = {( F 1 , F 2 ) : F 1 ≤ st F 2 } can be induced from the subspace = {( F 1 , F 2 ) : F 1 = G 1 , F 2 = G 1 G 2 }. Here, F 1 can be considered as the distribution of θ and F 2 as the distribution of max( θ , δ ), where θ ~ G 1 and independently δ ~ G 2 . The stochastically ordered mixture distribution is then given as F ( y ; F 1 ) = ∫ K ( y ; θ ) d G 1 ( θ ) , F ( y ; F 2 ) = ∫ ∫ K ( y ; max ( θ , δ ) ) d G 1 ( θ ) d G 2 ( δ ) , with F ( y; F 1 ) ≤ st F ( y; F 2 ). Independent DP priors can be used for G 1 and G 2 : G 1 ~ DP ( α 1 G 10 ), G 2 ~ DP ( α 2 G 20 ). Thus the DPM model turned out to be a useful framework to account for stochastic order restriction ( Gelfand and Kottas, 2001 ; Hanson, Kottas and Branscum, 2008 ). A similar approach can be taken to account for ordering in variability ( Kottas and Gelfand, 2002 ). Let X and Y be two random variables with distribution functions F 1 and F 2 . Following the notion of sign changes of a function, the number of sign changes of a function h is given by S − ( h ) = sup S − ( h ( u 1 ), … , h ( u m )) for u 1 ≤ ··· ≤ u m , u i ∈ U ⊆ , where m is arbitrary, but finite and S − (·) is the number of sign changes of the indicated sequence. Then, Y is said to be larger than X in variability order, denoted by F 1 ≤ var F 2 if F 1 and F 2 have the same location and S − F 2 − F 1 ) =1, the sign sequence being +, −. To see how variability order can be preserved through the mixture distribution, consider a parametric family of distributions { K ( y; θ ) : θ ∈ }, with K ( y; θ ) differentiable in θ , such that for any θ 1 < θ 2 , S − ( K ( y; θ 2 ) − K ( y; θ 1 )) = 1, the sign sequence being +, −. Then, if F 1 ≤ st F 2 , F ( y; F 1 ) ≤ var F ( y; F 2 ), where F ( y; F i ) = ∫ K ( y; θ ) dF i ( θ ), i = 1, 2. As in stochastic order, we can let F 1 = G 1 and F 2 = G 1 G 2 , where G 1 ~ DP ( β 1 G 10 ), G 2 ~ DP ( β 2 G 20 ), and choose a symmetric scale family kernel with a location parameter μ for the mixtures. It follows F ( y ; μ , F 1 ) = ∫ K ( y ; μ , θ ) d G 1 ( θ ) , F ( y ; μ , F 2 ) = ∫ ∫ K ( y ; μ , max ( θ , δ ) ) d G 1 ( θ ) d G 2 ( δ ) , with F ( y; μ , F 1 ) ≤ var F ( y; μ , F 2 ). Let Y ij = ( Y ij 1 , …, Y ijn ij ) T be a ( n ij × 1) response vector measured on i th population, i = 1 (healthy), 2 (diseased) at j th setting, j = 1, 2. We develop scale and location mixture models for the distribution functions of each population, F ij ( y ) ≡ ∫∫ K ( y; θ , σ ) dH i ( θ ) dG j ( σ ), i = 1, 2 and j = 1, 2. For simplicity in illustration, we assume for now that there exists a gold standard for disease and will relax this assumption in Section 5. It is possible that Y i 1 and Y i 2 , i = 1, 2 are correlated; direct modeling of this correlation is not straightforward in the Bayesian nonparametric framework. For this reason, we chose not to explicitly build a correlation structure, although the constraints across settings will induce some dependence. We account for two order constraints simultaneously: stochastic order and variability order. First, for the j th setting, j = 1, 2, the response vector ( Y 1 j , Y 2 j ) has a stochastic order between the healthy and diseased populations. That is, the diseased population is more likely to have higher values of Y than the healthy population. Second, for the i th population, i = 1, 2, the response vector ( Y i 1 , Y i 2 ) has a variability order between two settings. It implies that setting 2 is more likely to have smaller variability of Y than setting 1. The following result illustrates how stochastic and variability orders can be preserved through scale and location mixture distributions, which extends the work of Gelfand and Kottas (2001) and Kottas and Gelfand (2002) that have separately considered these two constraints. The proof of this result can be found in the supplementary materials . Consider a parametric family of distributions { K ( y; θ , σ ) : θ ∈ , σ ∈ }, where θ is a location parameter, σ is a scale parameter and is the positive real line. Suppose 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 distributions H 1 and H 2 on , and distributions G 1 and G 2 on , define the mixtures F ij ( y; H i , G j ) ≡ K ( y; θ , σ ) dHi ( θ ) dG j ( σ ), i , j = 1, 2. Then if H 1 ≤ st H 2 , F 1 j ( y; H 1 , G j ) ≤ st F 2 j ( y; H 2 , G j ) for j = 1; 2, and if G 1 ≤ st G 2 , F i 1 ( y; H i , G 1 ) ≤ var F i 2 ( y; H i , G 2 ) for i = 1, 2. In practice, we use K ( y; θ , σ ) = Φ( y; θ , σ ) as it satisfies the conditions in Result 1. We let H 1 = P 1 , H 2 = P 1 P 2 for the location parameter and G 1 = Q 1 Q 2 , G 2 = Q 1 for the scale parameter, and introduce latent variables, θ ~ P 1 , δ ~ P 2 and σ ~ Q 1 , ψ ~ Q 2 . The likelihood of the proposed model under jointly stochastic and variability order constraints is given by (1) L = ∏ i = 1 2 ∏ j = 1 2 ∏ k = 1 n i j F i j ( y ijk ; H i , G j ) , where the corresponding mixture distributions are as follows: F 11 ( y ; H 1 , G 1 ) = ∫ ∫ ∫ Φ ( y ; θ , max ( σ , ψ ) ) d P 1 ( θ ) d Q 1 ( σ ) d Q 2 ( ψ ) ; F 21 ( y ; H 2 , G 1 ) = ∫ ∫ ∫ ∫ Φ ( y ; max ( θ , δ ) , max ( σ , ψ ) ) d P 1 ( θ ) d P 2 ( δ ) d Q 1 ( σ ) d Q 2 ( ψ ) ; F 12 ( y ; H 1 , G 2 ) = ∫ ∫ Φ ( y ; θ , σ ) d P 1 ( θ ) d Q 1 ( σ ) ; F 22 ( y ; H 2 , G 2 ) = ∫ ∫ ∫ Φ ( y ; max ( θ , δ ) , σ ) d P 1 ( θ ) d P 2 ( δ ) d Q 1 ( σ ) .

We take independent DP priors with normal base distributions for P and inverse-gamma base distributions for Q to obtain the following hierarchical model: (2) Y 11 k 1 ∣ · ∼ ind N ( θ k 1 , max ( σ n 12 + n 22 + k 1 , ψ k 1 ) ) , k 1 = 1 , … , n 11 , Y 21 k 2 ∣ · ∼ ind N ( max ( θ n 11 + k 2 , δ k 2 ) , max ( σ n 11 + n 12 + n 22 + k 2 , ψ n 11 + k 2 ) ) , k 2 = 1 , … , n 21 , Y 12 k 3 ∣ · ∼ ind N ( θ n 11 + n 21 + k 3 , σ k 3 ) , k 3 = 1 , … , n 12 , Y 22 k 4 ∣ · ∼ ind N ( max ( θ n 11 + n 21 + n 12 + k 4 , δ n 21 + k 4 ) , σ n 12 + k 4 ) , k 4 = 1 , … , n 22 , where The hyperparameters μ i , τ i 2 , a i , b i , i = 1, 2 are all assumed fixed such that μ i is roughly located in the middle of the data and the rest exhibit large variance. The precision parameters ( α 1 , α 2 , β 1 , β 2 ) of the DP play an important role in deciding the number of distinct values in our model. For instance, for relatively small α 1 , each random weight is likely to be near one and as a result, many θ k 1 ’s will share the same distinct value. Ideally, we would like to assign a hyperprior to α and β . However, for simpler implementation, we conduct a sensitivity analysis for fixed α and β . If the results differ substantially with α and β , choice of α and β could be determined by comparing predictive and empirical distributions. As a result of the hierarchical presentation in ( 2 ), the likelihood function ( 1 ) can be rewritten as: (3) L = ∏ k 1 = 1 n 11 N ( y 11 k 1 ; θ k 1 , max ( σ n 12 + n 22 + k 1 , ψ k 1 ) ) × ∏ k 2 = 1 n 21 N ( y 21 k 2 ; max ( θ n 11 + k 2 , δ k 2 ) , max ( σ n 11 + n 12 + n 22 + k 2 , ψ n 11 + k 2 ) ) × ∏ k 3 = 1 n 12 N ( y 12 k 3 ; θ n 11 + n 21 + k 3 , σ k 3 ) × ∏ k 4 = 1 n 22 N ( y 22 k 4 ; max ( θ n 11 + n 21 + n 12 + k 4 , δ n 21 + k 4 ) , σ n 12 + k 4 ) Posterior computation proceeds with a hybrid MCMC algorithm consisting of Gibbs sampling ( Casella and George, 1992 ) and Metropolis-Hastings steps ( Chib and Greenberg, 1995 ). We employ a truncation approximation for the stick-breaking representation of the DP priors for P 1 , P 2 , Q 1 and Q 2 , proposed by Ishwaran and Zarepour (2000) . This method samples from the posteriors using MCMC algorithms, resulting in rapid mixing of the Markov chain and is easy to apply to normal mean mixture models. We let Appropriate truncation levels M and L can be determined by examining the moments of the random weights. As a simpler approximation, we recommend using a moderately large numbers (e.g., M =30) and sensitivity to the approximation could be then assessed by using a larger M and comparing results. The predictive density of Y ijk , i , j = 1, 2 can be obtained at each iteration of the MCMC: p ( Y 11 k 1 ∣ · ) ∝ ∑ h 1 = 1 M ∑ l 1 = 1 L ∑ l 2 = 1 L p h 1 q l 1 q l 2 N ( Y 11 k 1 ; θ h 1 ∗ , max ( σ l 1 ∗ , ψ l 2 ∗ ) ) , p ( Y 21 k 2 ∣ · ) ∝ ∑ h 1 = 1 M ∑ h 2 = 1 M ∑ l 1 = 1 L ∑ l 2 = 1 L p h 1 p h 2 q l 1 q l 2 N ( Y 21 k 2 ; max ( θ h 1 ∗ , δ h 2 ∗ ) , max ( σ l 1 ∗ , ψ l 2 ∗ ) ) , p ( Y 12 k 3 ∣ · ) ∝ ∑ h 1 = 1 M ∑ l 1 = 1 L p h 1 q l 1 N ( Y 12 k 3 ; θ h 1 ∗ , σ l 1 ∗ ) , p ( Y 22 k 4 ∣ · ) ∝ ∑ h 1 = 1 M ∑ h 2 = 1 M ∑ l 1 = 1 L p h 1 p h 2 q l 1 N ( Y 22 k 4 ; max ( θ h 1 ∗ , δ h 2 ∗ ) , σ l 1 ∗ ) . More details regarding our sampling algorithm including full conditional posterior distributions can be found in the supplementary materials . The estimation of ROC curves is obtained using predictive distribution function of Y ij , i , j = 1, 2. That is, RÔC ( u ) = 1− F̂ D { F̂ H −1 (1− u )}, for 0 ≤ u ≤ 1, where F̂ H (·) and F̂ D (·) are posterior distributions of Y in healthy and diseased populations, respectively, and s is a threshold value. To assess diagnostic accuracy, AUC is estimated by the posterior mean of the AUCs obtained using AUC = Pr ( Y D > Y H ) at each iteration of the MCMC, where Y H and Y D are sampled from predictive densities of the healthy and diseased populations, respectively at each iteration ( Bamber, 1975 ).

In this section, we build on the initial analysis in Section 2 and re-analyze the PRS data with two objectives in mind. First, we incorporate constraints in the estimation and second, we address the issue of no gold standard. As described earlier, the design of PRS involves successively augmented clinical information across settings for the physicians to use to diagnose endometriosis. In particular, the physicians were only provided with the women’s intrauterus digital images at setting 1. At setting 2, surgeons’ notes were added to the clinical file that was then made available to the physicians. As a result, the physicians were seeing a larger collection of clinical information for diagnosing endometriosis at setting 2 than at setting 1. In our collaborations with OB/GYN experts two order constraints were suggested: C1 : For a given physician at a given setting, the rASRM score should have a stochastic order between the diseased and healthy populations. It means that the diseased population is more likely to have higher rASRM score than the healthy population. C2 : For a given woman and a given physician, the rASRM score in setting 2 is more likely to have smaller variability than setting 1. While the stochastic order is consistent with what was observed in the data, the actual rASRM score has larger variance in setting 2 than 1 ( Table 1 ). It is of interest to see how our method models the data under the two order constraints, one of which has the opposite direction to the real data. In the initial analysis and model development sections, we have assumed a known gold standard for the disease. To handle the situation that no gold standard exists in the PRS, we construct an imperfect reference standard based on the diagnoses of the four IEs. Let D k be the true disease status of the k th subject, k = 1, …, N . We denote T jk as the sum of the IEs’ positive diagnoses and J jk as the number of non-missing IEs’ diagnoses on k th subject at j th setting, k = 1, …, N , j = 1, 2. We model D k as a binary latent variable by using the following logit model, P ( D k = 1 ∣ T k , J k ) = exp ( γ 0 + γ 1 [ 1 - T k J k ] ) 1 + exp ( γ 0 + γ 1 [ 1 - T k J k ] ) , where T k = T 1 k + T 2 k and J k = J 1 k + J 2 k . Ideally, we would like to take hyperpriors for γ 0 and γ 1 in the model, but that will result in weak identifiability. We instead conducted a sensitivity analysis in which we vary γ 0 and γ 1 in a scientifically reasonable way and assess robustness of the results. In particular, we chose γ 0 = 4.5 and γ 1 =−9 so that P ( D k = 1| T k / J k = 1) is close to 1 and P ( D k = 1| T k / J k = 0) is close to 0. We examine model robustness using various γ 0 and γ 1 in the PRS data. With the inclusion of the sensitivity analysis component, the likelihood function ( 3 ) now has the following form: L = ∏ k = 1 N [ N ( y 1 k ; θ k , max ( σ N + k , ψ k ) ) N ( y 2 k ; θ N + k , σ k ) ] 1 - D k × [ N ( y 1 k ; max ( θ k , δ k ) , max ( σ N + k , ψ k ) ) N ( y 2 k ; max ( θ N + k , δ N + k ) , σ k ) ] D k × [ 1 - exp ( γ 0 + γ 1 [ 1 - T k J k ] ) 1 + exp ( γ 0 + γ 1 [ 1 - T k J k ] ) ] 1 - D k [ exp ( γ 0 + γ 1 [ 1 - T k J k ] ) 1 + exp ( γ 0 + γ 1 [ 1 - T k J k ] ) ] D k . In the re-analysis, we fit a sequence of models, from no constraint to all constraints, using DPM for each population to allow distributional flexibility and the sensitivity approach in Section 5.2 to address no gold standard problem. These models are: No-constraint (NO) model : No constraint is imposed on the 4 distributions of the response variable ( Y ), except each follows a different DPM. This amounts to separate ROC estimation between the 2 settings; Stochastic order constraint only (SO) model : On top of (i), we impose the stochastic order constraint (C1) on the distributions within each setting, while leaving the variability unconstrained across settings; Variability order constraint only (VO) model : On top of (i), we impose the variability order constraint (C2) on the distributions within each population, while letting the distributional center unconstrained across populations; Joint (JO) model : On top of (i), we impose the stochastic order (C1) within each setting and variability order (C2) within each population. More details regarding the four model specifications can be found in the supplementary materials . To make the comparison between the four models fair, the DPM in models (i) to (iii) is implemented with respect to both the location and scale parameters, same as in model (iv). We chose normal mixture kernels K ( y; θ , σ ) = Φ( y; θ , σ ) for all models, and took the DP priors with normal base distribution N (0, 10 2 ) for the mean and inverse-Gamma base distribution IG (2, 1) for the variance. The precision parameters of DP priors were fixed at 1, i.e., α 1 = α 2 = β 1 = β 2 = 1, since we noticed the parameters were not very sensitive after comparing predicted and empirical distributions. For the truncation approximation of DP, we used M = L = 30 as a truncation level, which was large enough based on our empirical examination. We obtained results by running the MCMC algorithm for 50,000 iterations with 25,000 iteration burn-ins. Particularly, the multi-chain MCMC algorithm with various starting points was implemented in the joint model to produce satisfactory convergence results while saving on running times. Figure 5 contains the estimated predictive densities (and histograms) of log rASRM scores from models (i) – (iv). This figure suggests that the nonparametric feature of our models identifies very well the non-standard shape of the distributions, such as skewness in the healthy populations and bimodality in the diseased populations. This may not be possible if parametric models were taken to model these scores. The main results of applying the proposed methods to the PRS data are provided in Figure 6 and Table 2 , where the estimated ROC curves and corresponding posterior estimates of AUCs and their 95% credible intervals are presented for each of the four models at each of the two settings. While we observe that the REs are generally good at diagnosing endometriosis at both setting 1 (AUC ranges from 0.80 to 0.84) and setting 2 (AUC ranges from 0.77 to 0.85), we also notice the impact of the order constraints. The AUC estimates hardly changed from the no-constraint (NO) model to the SO model; this is expected as the PRS data already agree with the stochastic order constraint. Going from the NO to the VO model reveals some subtle changes. First, the direction of the estimated AUCs changed. While the AUC at setting 1 is larger than setting 2 in the NO model, it is smaller once the variability order is considered. This illustrates that although the data may suggest one directional order of the variability, it is the order constraint that dominates the inference when it is in disagreement with the data. We also notice the shrinkage of the 95% CI width in the VO model compared to the NO model, an indication of efficiency gain, although marginally. These two observations are more pronounced once the stochastic and variability orders are jointly considered in the model. As in the VO model, the AUC in setting 2 is larger than setting 1 in the JO model; however, the efficiency gain over the NO model is more substantial, with the width of 95% CIs at 0.11 for the JO model versus 0.23 for the NO model in setting 1 and 0.11 versus 0.19 in setting 2, respectively. This large efficiency gain in AUC estimates is likely a result of the interplay of the stochastic and variability orders. In conducting the sensitivity analysis, we varied γ 0 and γ 1 over a wide range of reasonable values: γ 0 = 2, 3, 4.5, 6, 7 and γ 1 =−15, −12, −9, −8, −7. When we used different values of γ 0 ’s with γ 1 =−9 and different values of γ 1 ’s with γ 0 = 4.5, respectively, we found the estimates of ROC curves and AUCs in all scenarios to be very close to those from our primary analysis.

We conduct simulation studies under four scenarios to assess the performance of our proposed approach. These scenarios consider different specifications in distributions, variability order directions, and sample sizes. We fit the same four models as in Section 5.3 to the simulated datasets and evaluate their performances in terms of bias and efficiency. The results of the simulation studies corroborate most of what was observed in the real data analyses: Our approach can adequately identify the distributions of the scores; There are no discernable differences in biases between the two sample sizes and the two distribution assumptions; There is efficiency gain when constraints are considered; Having the variability order in the opposite direction of the data results in estimated AUCs with reversed direction. Detailed simulation set up and rersults are reported in Section 5 of the supplementary materials .

In this paper, we have proposed an integrated framework for jointly incorporating order constraints, one on the center and the other on the variability of the distributions, and applied it to estimating multiple ROCs. The work advances early literature which addressed these constraints separately. This integrated framework was made possible by utilizing a Bayesian nonparametric technique, i.e., the Dirichlet process mixture, which also provided the needed distributional flexibility in the test scores. Substantively, we found that a regional OB/GYN physician can better diagnose endometriosis using both the intrauterus image and surgeon’s notes (setting 2) than using the image only (setting 1) after both order constraints are incorporated. This represents a reversal of the direction of AUCs between settings 1 and 2 when the order constraints were not considered. A closer look reveals that the change of direction comes as a result of the variability order, not of the stochastic order. That the stochastic order is not related to the AUC reversal is not surprising, as the PRS data already agree with this order constraint. In modeling the REs’ ROCs, we took the arithmatic average of the 4 individual rASRM scores in the RE group. This step is reasonable as we assumed that the 4 REs are exchangeable and as the research goal is on the group-specific, not physician-specific diagnostic accuracy. While alternative ways exist (e.g., obtain physician-specific ROCs and then average them), our approach represents our first step of studying REs’ diagnostic accuracy. We feel that if the accuracy obtained from the average score is not great, then that from individual scores will be worse. Using the average scores also enabled us to retain those women with some scores missing, which in turn allowed us to have reasonable sample sizes. An alternative to our imperfect reference standard approach in the PRS data is to model the diagnostic data with neither a gold standard nor an imperfect reference standard. Such an approach essentially amounts to a deconvolution of a two-component mixture; see Branscum et al. (2008) . However, as others and Branscum et al. have pointed out, this deconvolution approach inherently lacks identifiability unless strong prior information is used. When aided with no useful prior information, such an approach can result in slow mixing and poor convergence in the Markov chain Monte Carlo iterations. Indeed, when we ran a simulation study using no imperfect reference standard, we observed that under weak prior, the MCMC chain has poor mixing and converges to an AUC that is far away from the truth (see supplementary materials ). While it is common to estimate ROCs and the associated AUCs through modeling the distributions of the observed scores, alternatives have been suggested in the literature that directly models the ROCs. Of particular interest is the idea of using placement values (see e.g., Pepe (2000) and Alonzo and Pepe (2002) ). In this framework, a priori constraints on the ROCs, whether in the entire or partial range of specificity, can then be incorporated. This is a current research interest of the authors. In our proposed models, we have specified independent mixture models for the location ( θ ) and scale ( σ ) of the DPM kernel and used the idea of product measures to induce order constraints, i.e., H 1 = P 1 , H 2 = P 1 P 2 for stochastic order and G 1 = Q 1 Q 2 , G 2 = Q 1 for variability order (see the last paragraph of Section 3.5). However, the independent mixture can potentially lead to unsubstatiated clusters. Alternatively, one can specify a dependent mixture: for i = 1, 2, j = 1, 2 and k = 1, …, n ij , Y ijk ~ N ( θ i k , σ j k ) , ( θ i k , σ j k ) ~ H i j , where H ij ~ DP ( αH ij 0 ), say, with H ij 0 a bivariate distribution of θ and σ . Under this dependent mixture specification, it is not clear how the product measure idea can be used to induce the desired constraints. One possibility is to work with the multivariate distribution of the vector ξ k = ({ θ ik }, { σ jk }) and follow Hoff (2003) to use partially ordered latent observations to model stochastically ordered marginals, e.g., F ( θ 1 k ) and F ( θ 2 k ) or F ( σ 1 k ) and F ( σ 2 k ) (see also Dunson and Peddada (2008 )). This is an active research project of the authors. To gauge the impact of the two mixture specifications, we implemented the dependent mixture idea in the no-constraint model, and found quite similar results (see supplementary materials ). While it is convenient to use the product measures, i.e., P 1 , P 2 and Q 1 , Q 2 in Section 3, to induce the ordering, such a construction also results in a reduced space for the desired order constraints. It is therefore of interest to examine the impact of such a limitation and, if needed, to explore more general ways of incorporating these constraints in a Bayesian nonparametric framework. It is also of interest to consider other order constraints, such as the stochastic precedence on the distributional centers ( Chen and Dunson, 2004 ). In the PRS, the REs also gave ordinal diagnosis on the participants. Methodologically, it is of interest to explore possible extension of the current framework to the case with ordinal data.