Parametric modeling of receiver operating characteristics curves

1. Introduction

Receiver operating characteristics (ROC) analysis is extensively used in machine vision and machine learning on issues of detecting the presence or absence of a target in the field of view of a sensor (Bradley, 1997; Chang, 2010; Provost & Fawcett, 2001). Its wide use in clinical diagnostics is also well established (Obuchowski & Bullen, 2018; Swets, 1988; Vickers, 2008). Regardless of whether the interest is in machine vision or clinical diagnostics, availability of appropriate statistical fits for the ROC curves are essential. These fits will allow the user to identify the optimal operating points and draw appropriate conclusions regarding the sensor (machine vision) or a clinical test (medical diagnostics). They will also provide a means to test signal processing algorithms (combining data from multiple sensors in machine vision or combining data from multiple clinical observations) for performance enhancements before actual implementation (Brown & Davis, 2006; Hanley, 1996; Metz, 1989; Walsh, 1999). The goodness-of-fit (GOF) studies of the ROC curves undertaken in terms of binormal and bigamma fits have been around for a while, and interest has been directed at exploring other fits (Dorfman et al., 1997; Kester & Buntix, 2000; Hoyer and Kuss, 2020; Pundir & Amala, 2014; Rota & Antolini, 2014). Recent efforts based on fully parametrized bi-distributional models for ROC have examined several statistical models using a meta-analysis with pairs of sensitivity and specificity (Hoyer & Kuss, 2020; Kester & Buntix, 2000; Pundir & Amala, 2014). None of these reported models and accompanying analyses examined the impact of such models on the critical threshold values (Neyman-Person, Youden’s index and Bayes’ risk) for the interpretation of ROC curves (Chang, 1989; Lipovetsky & Conklin, 2020; Provost & Fawcett, 2001; Zweig & Campbell, 1993), and estimation of sensitivity, specificity, positive predictive values, etc. In other words, the efficacy of these statistical models to produce results that match those provided by the data (empirical results) has not been studied.

In this work, a methodology is proposed and tested to establish a statistical fit for the data under the hypothesis H 0 (disease absent in a clinical scenario or target absent in the case of machine vision) and another statistical fit for the data under the hypothesis H 1 (disease present in a clinical scenario or target present in the case of machine vision). These two fits are used to obtain a continuous fit (goodness-of-fit) for the ROC curve. Multiple hypothesis testing relying on chi square testing is used for both hypotheses individually to ascertain the best statistical fit in each case (Barton, 1956; Cochran, 1954; Watson, 1958). The densities studied include gamma, Nakagami, Rician, Rayleigh, lognormal, Weibull, inverse Gaussian, generalized gamma (Stacy’s distribution), Burr distribution, and generalized Pareto distribution. These densities were chosen considering various statistical models used in modeling sensor (radio frequency, wireless, sonar) and medical imaging (X-ray and ultrasound) data (Destrempes & Cloutier, 2010; Simon & Alouini, 2001). The Burr and Pareto distribution find extensive application in econometrics (Hogg & Klugman, 1983). Use of different densities for the data sets from the two hypotheses makes the effort in this research different from the existing body of research where same density is used for modeling data from both hypotheses (Hoyer & Kuss, 2020; Pundir & Amala, 2014).

A two-step methodology is followed here to carry out the analysis and interpretation of the ROC curves. In the first step, the area under the ROC curve (AUC) is examined in terms of its mean and standard deviation (Siva et al., 2023). In other words, AUC obtained from the data (empirical) is compared to the AUC calculated using the continuous fit for the ROC curve. In this regard, three different fits are chosen, a bigamma fit used by many researchers, a bi-generalized gamma fit and a biBest fit using the two fits obtained from hypothesis testing (Demidenko, 2012; Dorfman et al., 1997; Hanley, 1988; Rota & Antolini, 2014). A bi-generalized gamma fit is a general form of a bigamma fit because a three parameter Stacy’s distribution or the generalized gamma density becomes a gamma density when one of the parameters of the Stacy’s distribution is equal to 1. AUC comparison relies on both parametric and non-parametric bootstrapping (Davison et al., 2003). Because the sensor data values are always positive, binormal fits were not examined. In the second step, the ROC curves (empirical and continuous fits) are studied in terms of three thresholds used by researchers (Schisterman et al., 2005); distance to the top left corner of the ROC curve (Neyman-Pearson), the difference between probability of detection or sensitivity and probability of false alarm (Youden’s index J) and the accuracy (Bayes’ risk). A parametric hypothesis test (two-sided t -test) is conducted with the null hypothesis supporting the notion that in terms of the three thresholds, there is no statistically significant difference between empirical ROC and each of these different fits along the ROC curve (Rohatgi & Saleh, 2000).

The manuscript starts with the attributes of the data sets used in this study. The section on statistical methods provides details of the densities used. The section on methodology offers insights into the approaches for the analysis of the ROC fits using hypothesis (chi square) tests, bootstrapping, and t -tests. This is followed by results and discussion.

2. Data sets

The author has been active in research on statistical modeling of wireless channels and backscattered ultrasonic signals (medical ultrasound). During the past few years, the author incorporated data analytics in the undergraduate probability course for engineers where every student was provided with a data set containing samples of sensor outputs (target absent and target present). These randomly generated data sets were rigorously tested to ensure they had attributes of data that might fit both sensor outputs in machine vision and statistics of medical imaging data. Three data sets (#1 – #3) used in this study are displayed in Table 1. Each set has 70 values representing the Target Absent (hypothesis H₀) cohort and 60 values representing the Target Present cohort. The Target Present (hypothesis H₁) cohort has a mean value that is higher than the mean value of the Target Absent cohort.

Each data set implies existence of 131 samples of threshold values (130 samples of data along with the sample value of 0 allowing for the lowest possible measured value) for analysis that take the ROC curve from the bottom left corner [0,0] to the top right corner [1,1] of the ROC plot. Note that the ROC curve is drawn with the probability of false alarm (1-specificity) along the X-axis and the sensitivity or probability of detection along the Y-axis. For each threshold, it is possible to obtain a confusion matrix in terms of false and missed counts or a transition matrix in terms of the probabilities of false alarm and miss (Ltifi et al., 2020; Provost & Fawcett, 2001). If the probability density functions associated with the two hypotheses are represented as f X ( x | H k ) , k = 0,1, the probability of false alarm ( P F ) and the probability of detection ( P D ) become

$P_F={\displaystyle \underset{T}{\overset{\mathrm{\infty }}{\int }}}f(x|H_0)dx=\text{Prob}(x|H_0>T)=S_X(T|H_0),$(1)

$P_D={\displaystyle \underset{T}{\overset{\mathrm{\infty }}{\int }}}f(x|H_1)dx=\text{Prob}(x|H_1>T)=S_X(T|H_1).$(2)

In Eqs (1) and (2), 𝑻 is the threshold and S (.) is the survival function (Rohatgi & Saleh, 2000). The probability of miss ( P M ) is

P M = 1 - P D = ∫ 0 T f ( x | H 1 ) d x = Prob ( x | H 1 ⩽ T ) = F X ( T | H 1 ) .(3)

It should be noted that

F X ( T | H 1 ) = 1 - S X ( T | H 1 ) .(4)

In Eq. (3), F (.) is the cumulative distribution function (CDF). Notice the existence of the equality sign within the parenthesis of in Eq. (3) reflecting the definition of the CDF of a random variable (Rohatgi & Saleh, 2000). The ROC curve is the plot of P F versus P D as the threshold varies from 0 to the largest value of the pooled ( H 0 and H 1 ) data.

If N 0 and N 1 are the sample counts in the respective cohorts, a priori probabilities of the two hypotheses become,

P ⁢ ( H 0 ) = N 0 N 0 + N 1 , P ⁢ ( H 1 ) = N 1 N 0 + N 1 .(5)

If N F represents the false counts (number of samples from the target absent cohort that exceed the threshold 𝑻 ) and N C represents the correct counts (number of samples from the target present cohort that exceed the threshold 𝑻 ), the probably of false alarm P F and probability of detection P D become,

P F = N F N 0 , P D = N C N 1 .(6)

For every threshold N F and N C are available and the confusion matrix ( C X ) becomes (Ltifi et al., 2020; Provost & Fawcett, 2001),

(7)

While the elements of rows sum up to the target absent (first row) and target present (second row) counts, the elements of the columns sum up to target not detected (first column) and target detected (second column) respectively. The transition matrix ( T X ) becomes (Ltifi et al., 2020; Provost & Fawcett, 2001),

(8)

The transition matrix has the property that elements of each column add to 1. Furthermore, the probabilities of detecting the target, P 𝐷𝑒𝑡 and not detecting the target, P 𝑛𝑜𝑡𝐷𝑒𝑡 can be expressed as

[ P 𝑛𝑜𝑡𝐷𝑒𝑡 P 𝐷𝑒𝑡 ] = T X ⁢ [ P ⁢ ( H 0 ) P ⁢ ( H 1 ) ] .(9)

The transition or confusion matrix completely defines the properties of the sensor or a decision maker that uses the data. These matrices rely on the threshold, demonstrating the importance of the choice of an appropriate threshold.

Table 1b and Table 1c provide two additional data sets (#2 and #3) that were used in this study. They have attributes like those of set #1.

While bigamma fits (Dorfman et al., 1997) have been used to model the ROC curves specifically when the data contains only positive values, in this work other models, one for the target absent cohort and another one for the target present cohort, are explored. A specific set of models where both hypotheses (target absent and target present) modeled in terms of identical densities was reported earlier (Pundir & Amala, 2014). The probability densities considered in this work are listed below (in all case, x ⩾ 0 and H k , k = 0,1 correspond to the two hypotheses):

Gamma

$f_X(x|H_k)={\displaystyle \frac{x^{a_k-1}}{b_k^{a_k}\mathrm{\Gamma }\left(a_k\right)}}\exp (-{\displaystyle \frac{x}{b_k}}),a_k,b_k>0;$(10)

Generalized gamma (Stacy’s distribution expressed in a slightly different format)

$f_X(x|H_k)={\displaystyle \frac{c_k{x}^{c_k{a}_k-1}}{b_k^{a_k}\mathrm{\Gamma }\left(a_k\right)}}\exp (-{\displaystyle \frac{x^{c_k}}{b_k}}),a_k,b_k,c_k>0;$(11)

Weibull

$f_X(x|H_k)=\left({\displaystyle \frac{b_k}{a_k}}\right){\left({\displaystyle \frac{x}{a_k}}\right)}^{b_k-1}\exp (-{\left[{\displaystyle \frac{x}{a_k}}\right]}^{b_k}),a_k,b_k>0;$(12)

Nakagami

$f_X(x|H_k)=2\left({\displaystyle \frac{m_k}{{\mathrm{\Omega }}_k}}\right){\displaystyle \frac{1}{\mathrm{\Gamma }\left(m_k\right)}}{x}^{2m_k-1}\exp (-{\displaystyle \frac{m_k}{{\mathrm{\Omega }}_k}}{x}^2),{\mathrm{\Omega }}_k>0,m_k⩾{\displaystyle \frac{1}{2}};$(13)

Inverse Gaussian

$f_X(x|H_k)=\sqrt{{\displaystyle \frac{{\lambda }_k}{2\pi x^3}}}\exp (-{\displaystyle \frac{{\lambda }_k}{2{\mu }_k^{2}x}}{(x-{\mu }_k)}^2),{\lambda }_k,{\mu }_k>0;$(14)

Lognormal

$f_X(x|H_k)={\displaystyle \frac{1}{\sqrt{2\pi {\sigma }_k^2{x}^2}}}\exp (-{\displaystyle \frac{1}{2{\sigma }_k^2}}{[{\log }_{e}x-{\mu }_k]}^2),{\sigma }_k>0;$(15)

Rayleigh

$f_X(x|H_k)={\displaystyle \frac{x}{b_k^2}}\exp (-{\displaystyle \frac{x^2}{2b_k^2}}),b_k>0;$(16)

Rician

f X ( x | H k ) = x σ k 2 exp ( - x 2 + s k 2 2 ⁢ σ k 2 ) I 0 ( 𝑥𝑠 k σ k 2 ) , σ k , s k ⩾ 0 ;(17)

Burr

$f_X(x|H_k)={[1+{\left({\displaystyle \frac{x}{\alpha }}\right)}^c]}^{-\left(k+1\right)}{\displaystyle \frac{\text{𝑘𝑐}}{\alpha }}{\left({\displaystyle \frac{x}{\alpha }}\right)}^{c-1},\alpha >0,c>0,k>0;$(18)

Generalized Pareto

$f_X(x|H_k)={\displaystyle \frac{1}{\sigma }}{[1+k{\displaystyle \frac{\left(x-\theta \right)}{\sigma }}]}^{-\left(1+\frac{1}{k}\right)},x>\theta ,k>0,\sigma >0.$(19)

In Eq. (19), k and θ are real.

In the expressions for densities, Γ (.) is the gamma function and I 0 (.) is the 0^th order modified Bessel function of the first kind (Simon & Alouini, 2001). The generalized gamma density becomes the gamma density when the parameter c k = 1. The Rician density is included in the models so that existence of line-of-sight (LOS) component in wireless or radar communications can be considered. Rician density becomes Rayleigh when the LOS component, s k , becomes zero. Similarly, Nakagami density becomes Rayleigh when the Nakagami parameter m k becomes 1 (Simon & Alouini, 2001). Burr and generalized Pareto distributions are used in econometrics and under special conditions are related to each other, and some of the other densities (Hogg & Klugman, 1983). Normal densities (and hence binormal fits) are not included in this study considering the positive values of the data collected from the sensors.

The first step in analyzing the ROC curve through a continuous density model is to establish a statistical fit for the two cohorts of the data or the two hypotheses, H 0 and H 1 . A multi-hypothesis goodness-of-fit test is implemented using the chi square test. For each cohort, the best fit (the criterion yet to be defined) will be based on the results of chi square tests (Barton, 1956; Rohatgi & Saleh, 2000).

Chi square test is undertaken to estimate the test statistic and the p -value. The null hypothesis is not rejected if the p -value exceeds 0.05 (level of significance of 5%). This means that for each density listed in Section 3, a chi square test is conducted. The process is carried out for each cohort separately. The density with the highest value of p (exceeding 0.05) is taken as the best fit. It is understood that any one of the densities with p -value exceeding 0.05 might be acceptable, keeping in mind that p -value itself is random (Rohatgi & Saleh, 2000). This composite hypothesis testing results in one density as the best fit for the hypothesis H 0 and another density for the hypothesis H 1 , creating the pair identified as the bi-Best. It is possible that the two best fits may turn out to be the bigamma (gamma densities for both hypotheses). The pair may also be the bi generalized gamma (generalized gamma fits for both hypotheses). With the availability of the best fits, a fit for ROC curve is obtained, and a thorough analysis of the ROC curves can be carried out through bootstrapping, namely parametric bootstrapping to establish the appropriateness of the fits and non-parametric or standard bootstrapping for further comparison fits to the empirical set (Davison et al., 2003). Furthermore, t -tests are conducted to check whether the thresholds obtained from the data and those obtained from the fits match statistically (Rohatgi & Saleh, 2000). For this purpose, three different criteria are used; the Neyman Pearson, Youden’s index, and the Bayes’ (Lipovetsky & Conklin, 2020; Provost & Fawcett, 2001; Schisterman et al., 2005; Tilbury et al. 2000; Xu, 2012). The Neyman-Pearson criterion relies on the distance to the top left corner of the ROC curve ( P F = 0, P D = 1; the ideal operating point) and the optimum threshold is chosen based on the shortest distance. The distance to the top left corner of the ROC curve, 𝑫 , is,

D = P F 2 + ( 1 - P D ) 2 = P F 2 + P M 2 .(20)

The Youden’s index criterion is based maximizing the difference between P D and P F with the index 𝒀 defined as (Schisterman et al., 2005; Zweig & Campbell, 1993)

Y = P D - P F .(21)

Bayes’ criterion relies on cost factors involved with decisions. In the absence of cost considerations, Bayes’ criterion leads to a threshold that maximizes the accuracy of the system 𝑩 defined as

B = ( 1 - P F ) ⁢ P ⁢ ( H 0 ) + P D ⁢ P ⁢ ( H 1 ) .(22)

The continuous fits can be compared to the empirical ROC by conducting a paired t -test (2-sided) for all three thresholds. As stated earlier, with 130 samples of data, there will be 131 threshold values including the lowest possible value of the data, namely, a zero. Thus, t -tests use 131 pairs of values.

Matlab (www.mathworks.com) was used for the computational analysis relying primarily on the Statistics and Machine Learning Toolbox. All the densities other than generalized gamma are available in Matlab (version R2019a). Matlab also provides maximum likelihood estimates (MLE) of the parameters (of the densities), random numbers, probability density and cumulative distribution plots. The command chi2gof(.) performs the chi square tests providing various outputs (test statistic, degrees of freedom and p -value of the test).

Most of the densities are 2-parameter ones except for the Rayleigh, Burr, generalized gamma, and generalized Pareto. Rayleigh is a single parameter density while the generalized gamma, Burr, and generalized Pareto are 3-parameter densities. Of these, the generalized gamma is not available in Matlab. The generalized gamma variable (GG) is obtained from a gamma variable ( G ) as

𝐺𝐺 = G 1 c k .(23)

While fitdist(.) command in Matlab provides the MLE estimates of the parameters of the densities except the generalized gamma density. Matlab however allows MLE of parameters of any density using a ‘custom’ provision, making it possible to estimate the parameters of the generalized gamma density. While random numbers can be generated using the command, random(.), densities are obtained using the pdf(.) command and cumulative distribution functions obtained using the cdf(.) command. Continuous plots of the density are obtained using the command ksdensity(.), and histogram(.) provides the histogram with provisions for various normalized formats. Paired t -testing is accomplished through the command ttest(.). Equation (23) can be used to obtain the generalized gamma and its parameters, pdf and cdf. Traditional bootstrapping is carried out through bootstrp(.) command. Parametric bootstrapping is carried out through random number generation.

Data set #1 was analyzed first, followed by the other two sets. For all the densities, chi square tests were carried out in Matlab. For the generalized gamma, custom MLE routine in Matlab was used to obtain the parameters. Tables 2, 3, and 4 display the results of chi square tests, bootstrapping, and t -tests respectively for all three data sets. For each cohort (hypothesis), the densities are sorted in descending order of the p -values (Table 2).

OPEN IN VIEWER

Table 2a summarizes the results of chi square tests for set #1. The generalized Pareto and generalized gamma are taken as the best fits for the target absent and target present cohorts respectively. Several other densities also result in p -values exceeding 0.05. Figure 1 [left] displays the histogram of the data cohorts and the corresponding best fits and Fig. 1 [right] displays the empirical ROC curve (EMP) along with the bigamma fit (biG), bi-generalized gamma fit (bigG) and the best fit (biBest). The fits were obtained using the respective cumulative distribution functions as the threshold is varied as explained in Section 2. For set #1, all three ROC fits are sufficiently close to the empirical ROC curve. The similarities among these curves (three fits and the empirical) were examined through bootstrapping and paired hypothesis tests based on student’s t -tests (Rohatgi & Saleh, 2000).

The comparison of the ROC curves relied on both non-parametric and parametric bootstrapping (Davison et al. 2003). Traditional (non-parametric) bootstrapping (500 times) was used with the data set obtaining the 95% confidence interval, mean ( μ ) and standard deviation ( σ ) of AUC. Parametric bootstrapping was used to obtain the statistics of AUC obtained from the theoretical fits starting with the bigamma fits. This required generation of 70 random samples belonging to gamma density with parameters estimated from the data (H₀) and 60 samples belonging to gamma density with parameters estimated from the data (H₁). The ROC area with this set was calculated and the process repeated 500 times. This procedure leads to parametric bootstrapping for the bigamma case. This procedure was repeated with bi-generalized gamma fit and followed by the biBest fit. The results are displayed in Table 3a showing remarkable agreement with regards to the 95% confidence interval (CI) of the area (95% CI), its mean ( μ ) and standard deviation ( σ ). The parametric bootstrapping (500 boot samples) results are shown as EMP (empirical). The AUC values obtained directly from the ROC plots seen in Fig. 1 are close to the corresponding mean values obtained from bootstrapping.

While bootstrapping provides comparison of AUC values, there is also a need to examine the fits in terms of applicability to the estimation of the probabilities of false alarm and miss, and the positive predictive value through the determination of an optimum threshold (Liu, 2012; Schisterman et al., 2005). This step was accomplished through paired 2-sided t -tests. With the availability of 131 thresholds (130 value of the data and the lowest possible value of the data, namely 0), the empirical ROC can be paired with the theoretical ROC from the bigamma, bi-generalized gamma and the biBest fits based on three criteria, distance to the top left corner of the ROC plot (Neyman-Pearson) in Eq. (20), the Youden’s index in Eq. (21) and the Bayes’ criterion in Eq. (22). The results are tabulated in Table 4a. With the high p -values displayed, the 2-sided t -tests suggest that there are no statistically significant differences in the paired approaches, further demonstrating the strength of the theoretical fits. In other words, any interpretation of the ROC curve requiring a specific threshold (based on any one of the three criteria) is unaffected by a continuous fit.

OPEN IN VIEWER

Analysis was now carried out with set #2 (Table 1b). The summary of the chi square results is displayed in Table 2b. They show that the best fits for both data cohorts are Weibull densities. Several other densities also are acceptable fits based on the criterion of p -value exceeding 0.05. Figure 2 [left] displays the histogram and the best fits for each cohort. Figure 2 [right] shows the ROC curves (empirical and the three fits) along with the corresponding values of the area under the ROC curves (AUC).

Bootstrapping results are provided in Table 3b. While the best fit ROC results in slightly lower value of AUC, all the AUC values and the standard deviation are reasonably close. The mean values also are close to the AUC values displayed in Fig. 2b. Table 4b displays the results of the paired t -tests to examine whether there exist any statistically significant differences among the ROC fits (and the empirical curve) based on the three criteria. As before, results show that no statistically significant differences exist. The only exception is the threshold behavior associated with the bigamma fit for the case of Neyman Pearson criterion. The Burr and generalized Pareto distributions were not pursued because of the difficulties in estimation of their parameters.

OPEN IN VIEWER

Analysis was carried out with set #3 (Table 1c). The summary of the chi square results is displayed in Table 2c. They show that the best fits for the cohorts of data are gamma and Nakagami densities respectively. It is also seen that all the densities studied are acceptable fits for the target present cohort. Figure 3 [left] displays the histogram and the best fits for each cohort. Figure 3 [right] shows the ROC curves (empirical and the three fits). Table 3c summarizes the results of bootstrapping. All the AUC values and the standard deviation are reasonably close. Table 4c displays the results of the paired t -tests to examine whether there exist any statistically significant differences among the ROC fits (and the empirical curve) based on the three criteria. Except with the case of Neyman-Pearson criterion with regards to the bigamma fit, other fits show no statistically significant differences with the empirical ROC.

ROC analysis is an essential step in applications involving machine vision and medical diagnostics, and data analytics in general. The analysis of the data and subsequent decision-making process based on the interpretation of data are simplified with the availability of statistical fit for the ROC curves. This assumes that the ROC fits have been rigorously tested using established statistical tools. With fits assured for each cohort of the data sets (target absent and target present), the models and results presented here offer a very general approach to theoretical modeling of ROC curves (Kochański, 2022; Rota & Antolini, 2014).

All elements of a typical ROC curve such as the area under the ROC curve (AUC) and the three important thresholds (Neyman-Pearson, Youden’s J and Bayes’) are examined to ensure that the modeling is proper, and the model is appropriate. Combining hypothesis testing relying on chi square tests followed by bootstrapping and t -tests, it is possible to demonstrate the strength and validity of the statistical model. While the research reported here was conducted using simulated data, the approaches demonstrated here relying on established statistical methodologies, can still be applied to all data sets collected in experiments. Additional data sets and information on Matlab scripts will be available from the author upon request.