Calculating confidence intervals for the Youden optimal value and for the predictor value that yields the optimal value


For binary response models, the Receiver Operating Characteristic (ROC) curve is an important tool for assessing the fit of the model. The ROC curve is a plot involving two statistics, sensitivity and specificity, which are associated with each distinct observed predicted event probability produced by the model. A very commonly used measure of fit associated with the ROC curve is the area under the ROC curve (AUC). A value of AUC near 1 indicates that the model provides good discrimination between events and nonevents (see SAS Note 22630). The ROC curve and the AUC are available in the LOGISTIC procedure.

It is often of interest to determine the value(s) of the predictor(s) in the model that can be used to decide whether a given observation should be predicted as an event or as a nonevent. Since the fitted model provides only a predicted event probability for each observation, this requires finding an optimal cutoff probability corresponding to a point on the ROC curve, which can be compared to the observation's predicted probability to allow its classification as an event or nonevent.

There are several criteria that can be used to select an optimal point on the ROC curve. Multiple criteria are available in SAS® Viya® with PROC LOGISTIC and in SAS® 9.4 with the ROCPLOT macro (SAS Note 25018). One of the most common is the Youden criterion. The Youden criterion is computed as sensitivity+specificity-1. It is the height of each point above the diagonal line representing an uninformative model in a plot of the ROC curve. The point that maximizes the height is the Youden optimal point, J. The Youden criterion can be computed for each distinct point on the ROC curve. The OUTROC= option in the MODEL statement of PROC LOGISTIC provides a data set of the distinct points on the ROC curve. Using the variables in OUTROC= data set, the Youden criterion is computed for each point as

Youden=_SENSIT_ - _1MSPEC_;

where _1MSPEC_ is equal to 1-specificity. The optimal Youden point is then the observation with the largest value of this variable, J = max(Youden). Note that when there are multiple predictors in the model, several observations with different combinations of predictor values can achieve the optimal value.

Since each observation in the OUTROC= data set is associated with a predicted probability (given in the _PROB_ variable), you can find optimal predictor value(s) by finding the observation(s) with the same predicted probability as is associated with J. This probability will be called the Youden cutpoint, pJ. The predicted probabilities for all of the input observations are available using the PREDICTED= (or P=) option in the OUTPUT statement in PROC LOGISTIC.

In addition to a point estimate of J, a confidence interval for it can be obtained. Two ways of doing this are shown. First, by fitting a saturated Poisson model to the cell counts from using the Youden cutpoint to classify the observations. This is done in the same way as shown in SAS Note 24170. Alternatively, a resampling approach can be taken by generating bootstrap samples and obtaining the ROC curve for each.

A confidence interval can also be obtained for the predictor value that yields J. This can be obtained by refitting the model in PROC PROBIT and using the INVERSECL option or again by using a resampling approach.

While the following shows how you can obtain confidence intervals related to the Youden criterion, the same approach can be applied when using other optimality criteria available from PROC LOGISTIC or the ROCPLOT macro.

Example

The following statements use the remission data in the example titled "Stepwise Logistic Regression and Predicted Values" in the PROC LOGISTIC documentation.

Plot the ROC Curve and Find the Youden Optimal Point, J

These statements fit a logistic model with LI as the predictor. The OUTROC= option displays the ROC curve and saves the data set of points on the ROC curve. The P= option in the OUTPUT statement saves a data set with the predicted probabilities for all of the input observations:

   proc logistic data=remission;
      model remiss(event='1') = li / outroc=ROCpts;
      output out=outp p=p;
      run;

The ROC plot from the fitted model is shown below:

The maximum Youden criterion value, J, and the value of LI that maximizes the Youden criterion, LIJ, can be displayed using the ROCPLOT macro as follows:

   %rocplot(inroc=ROCpts, inpred=outp, p=p,
            id=li _opty_, optcrit=youden)

The macro also plots the ROC curve and adds labels to some of the points, including a Y symbol indicating the Youden optimal point:

LIJ can also be found using the following statements. The OUTROC= data set is sorted in descending order of the Youden criterion so that the first observation is the optimal value, J. That observation can be used to find any matches in the OUTPUT OUT= data set that provide the associated predictor value, LIJ:

   data ROCpts; set ROCpts;
      Youden=_sensit_-_1mspec_;  
      run;
   proc sort data=ROCpts; 
      by descending Youden;
      run;
   data opt; 
      if _n_=1 then set ROCpts;
      set outp;
      if _prob_=p;
      run;
   proc print;
      id _prob_ Youden;
      var li _sensit_ _1mspec_;
      format _prob_ 12.10;
      title 'Optimal cutpoint and probability';
      run;

First, the optimal results (_PROB_ = pJ, Youden = J ) from the above statements are shown:

Next, the optimal results as displayed by the ROCPLOT macro:

Confidence Interval for J, the Optimal Youden Criterion

In the OUTROC= data set, the variables _POS_, _NEG_, _FALPOS_, and _FALNEG_ are the cell counts of the observed by predicted table that results from using the Youden cutpoint, pJ, to classify the observations as predicted events or nonevents. One way to obtain a confidence interval for J, is to arrange the cell counts as separate observations and then use PROC NLMIXED to fit a saturated Poisson model. Using the formulas for the sensitivity, specificity, and Youden criterion in ESTIMATE statements can then provide point and confidence interval estimates of those statistics. See SAS Note 24170 for more information.

   proc transpose data=or(obs=1) out=optcounts;
      var _pos_ _falpos_ _falneg_ _neg_;
      run;
   data optcounts; set optcounts;
      do test=1,0; do response=1,0;
        set optcounts; output;
      end; end;
      run;
   proc nlmixed data=optcounts df=1e8;
      mu=exp(t0r0*(test=0 and response=0) + t0r1*(test=0 and response=1) + 
             t1r0*(test=1 and response=0) + t1r1*(test=1 and response=1));
      model col1 ~ poisson(mu);
      /* Cell Counts and Margins */
         n11=exp(t1r1); n10=exp(t1r0); n01=exp(t0r1); n00=exp(t0r0);
      /* Test margins and Response margins */
         t1=n11+n10; t0=n01+n00;
         r1=n11+n01; r0=n10+n00;
         total=n11+n10+n01+n00;
      /* Compute sensitivity, specificity, and Youden criterion */
         sens=n11/r1; spec=n00/r0;
         estimate 'Sensitivity' sens;
         estimate 'Specificity' spec;
         estimate 'Youden' sens+spec-1;
      ods select additionalestimates;
      title 'Youden Optimality Criterion';
      run;

Below are the results from the ESTIMATE statements. Values of the sensitivity, specificity, and the Youden optimal value, J, agree with the results from the statements above. Additionally, a test and confidence interval for each are provided. The estimated J significantly exceeds zero (p<0.0001) and has confidence interval (0.5153, 0.9291):

Confidence Interval for LIJ, the Optimal Predictor Value

The Youden optimal point on the ROC curve, J, has an associated cutoff probability, pJ, that can be used to match the predicted probability of one or more observations in the input data set. As shown above, LIJ = 0.9. That is, an observation with LI=0.9 has predicted probability matching pJ = 0.23693, yielding the maximum Youden criterion, J = 0.7222. Given that value of J, you can obtain a confidence interval (referred to as a fiducial interval) on LIJ. This is done by the following statements that refit the logistic model in PROC PROBIT and include the INVERSECL(PROB=pJ) option:

   proc probit data=remission inversecl(prob=0.2369268106);
      model remiss(event='1') = li / d=logistic;
      run;

The INVERSECL option provides an estimate of LI that yields event probability 0.2369 matching the probability associated with the maximized Youden criterion (J = 0.7222). The estimate is LI = 0.9 with confidence interval (-0.1338, 1.2570):

Bootstrap Estimates and Confidence Limits for J and LIJ

An alternative approach uses bootstrap resampling to obtain distributions for J and LIJ. 95% confidence limits for each are then provided by the 2.5 and 97.5 quantiles of their respective distributions.

In the following statements, PROC SURVEYSELECT generates 1000 bootstrap samples, each with the same proportions of events and nonevents as in the original data. The SEED= option is used so that these results can be reproduced, but in general you would either omit this option or specify a unique value. The logistic model is then fit on each sample, saving the OUTROC= data set and the parameter estimates (OUTEST=). The Youden criterion is computed for each observation in the OUTROC= data set of each sample. PROC SUMMARY creates a data set containing J and pJ for each sample. The DATA OptCalc step uses the logistic model and parameter estimates in each sample to solve for LIJ using the concepts further discussed in SAS Note 69873. The logistic model is

logit(pJ) = log(pJ / (1-pJ)) = b0 + b1*LIJ .

Solving for LIJ,

LIJ = [log(pJ / (1-pJ)) - b0] / b1 .

With J and LIJ computed in all samples representing their distributions, PROC UNIVARIATE provides the mean, median, standard error, and the 2.5 and 97.5 quantiles for both statistics and these are displayed by the PROC PRINT steps:

   proc sort data=remission out=r; by remiss; run;
   proc surveyselect data=r out=bootsamps outhits method=urs seed=3943
                     samprate=1 reps=1000(repname=_sample_);
      strata remiss;
      run;
   proc sort data=bootsamps; by _sample_; run;
   ods exclude all; options nonotes;
   proc logistic data=bootsamps outest=logparms;
      by _sample_;
      model remiss(event='1') = li / outroc=ROCpts;
      run;
   options notes; ods select all;
   data ROCpts; set ROCpts;
      Youden=_sensit_-_1mspec_;  
      run;
   proc summary data=ROCpts nway;
      class _sample_; var Youden;
      output out=MaxY max=J maxid(youden(_prob_))=p_J; 
      run;
   data OptCalc; 
      merge MaxY logparms(keep=intercept li);
      LI_J = (log(p_J/(1-p_J))-intercept)/li;
      run;
   proc univariate data=OptCalc noprint;
      var J LI_J;
      output out=YoudenCI pctlpts=2.5 97.5 pctlpre=CL_J CL_LI_J 
             mean=MeanJ MeanLI_J median=MedianJ MedianLI_J std=SEJ SELI_J;
      run;
   proc print noobs label; 
      var MeanJ MedianJ SEJ CL_J:;
      label MeanJ="Mean J" MedianJ="Median J" SEJ="StdErr J"
            CL_J2_5="Lower" CL_J97_5="Upper";
      title "Max Youden and 95% Confidence Limits";
      run;
   proc print noobs label; 
      var MeanLI_J MedianLI_J SELI_J CL_LI_J:;
      label MeanLI_J="Mean LI_J" MedianLI_J="Median LI_J"
            SELI_J="StdErr LI_J" CL_LI_J2_5="Lower" CL_LI_J97_5="Upper";
      title "Optimal LI and 95% Confidence Limits";
      run;

Below are the results from bootstrapping. The point estimate (mean or median) and confidence interval for J are quite similar to the nonresampling method shown above. The point estimate and upper limit for LIJ are also similar to the above method, though the lower estimate is considerably larger suggesting that its distribution is heavily piled up at 0.9 with quantiles up to 50% at that value: