How the multiple intercepts in an ordinal model are used in computing predicted probabilities in PROC PROBIT


Before illustrating how the intercepts of an ordinal model are involved in the computation of predicted values, first note that you do not need to compute these yourself either by hand or by writing a program. The P= option in the OUTPUT statement does this for you. You can compute predicted values for new observations (also called scoring) as well as for the original observations that are used to fit the model. However, if you need to compute predicted values without the aid of SAS, the following information shows you how the PROBIT procedure computes predicted values for ordinal models.

Using the cheese example in the PROC LOGISTIC chapter of the SAS/STAT User's Guide, the following statements fit an ordinal probit model with the single predictor ADDQUANT. The ODS OUTPUT statement creates a data set that contains the parameter estimates from the fitted model. PROC SORT and the ORDER=DATA option are used in order to model the higher levels (most liked) of the taste-rating response variable, Y. The P= option in the OUTPUT statement requests an output data set that contains the predicted cumulative probabilities, Pr(Y ≤ y( i ) ), where y( i ) represents the i th ordered value of the response.

data Cheese;
   do addquant = 1 to 4;
      do y = 1 to 9;
         input freq @@;
         output;
      end;
   end;
   label y='Taste Rating';
   datalines;
0  0  1  7  8  8 19  8  1
6  9 12 11  7  6  1  0  0
1  1  6  8 23  7  5  1  0
0  0  0  1  3  7 14 16 11
;
ods output parameterestimates=pe;
proc sort data=Cheese;
   by descending y;
   run;
proc probit data=Cheese order=data;
   weight freq;
   class y;
   model y = addquant;
   output out=out p=cp;
   run;

The ordinal model that PROC PROBIT always fits is as follows:

Pr(Y ≤ y( i ) ) = F(α + αi + β•addquant) , i = 2, 3, ... , k-1

where F is usually the normal or logistic cumulative distribution function (for probit or logistic models, respectively), k is the number of response levels, and i refers to the ordered values in the Response Profile table. If you use PROC SORT and the ORDER=DATA option to put the response values in descending order, as in this example, the lower ordered values are associated with the higher response values so that y( 1 ) = 9, y( 2 ) = 8, etc.:

PROC PROBIT always models the probabilities of the lower ordered values, which means that you are modeling probabilities of the higher response levels when you use PROC SORT and the ORDER=DATA option to put the response values in descending order. Because the values of the response are the integers 1, 2, ... , 9 in this example, the logistic and probit models are, respectively, as follows:

Pr(Y ≥ i ) = logistic(α + αi + β•addquant)
Pr(Y ≥ i ) = probnorm(α + αi + β•addquant)

The PROBNORM function evaluates the normal cumulative distribution function, often denoted as Φ. The LOGISTIC function evaluates the logistic distribution function. Equivalently, use the CDF function with 'NORMAL' or 'LOGISTIC' as the first argument, designating the distribution. See SAS Language Reference: Dictionary for details.

The estimates of the probit model's parameters α, αi, and β are displayed to 10 decimal places by the following statements:

proc print data=pe;
   var parameter estimate;
   format estimate 15.10;
   run;

Using the parameters of the fitted probit model, the estimated cumulative probabilities are computed as follows. It is important to use many significant digits of the parameter estimates to avoid rounding errors.

Pr(Y ≥ 9 ) = probnorm(-2.3624248815 + 0.2757618404*addquant)
Pr(Y ≥ 8 ) = probnorm(-2.3624248815 + 0.7150632656 + 0.2757618404*addquant)
Pr(Y ≥ 7 ) = probnorm(-2.3624248815 + 1.3283837985 + 0.2757618404*addquant)
Pr(Y ≥ 6 ) = probnorm(-2.3624248815 + 1.6826398116 + 0.2757618404*addquant)
Pr(Y ≥ 5 ) = probnorm(-2.3624248815 + 2.2169522954 + 0.2757618404*addquant)
Pr(Y ≥ 4 ) = probnorm(-2.3624248815 + 2.6570246437 + 0.2757618404*addquant)
Pr(Y ≥ 3 ) = probnorm(-2.3624248815 + 3.1187431200 + 0.2757618404*addquant)
Pr(Y ≥ 2 ) = probnorm(-2.3624248815 + 3.5635390005 + 0.2757618404*addquant)
Pr(Y ≥ 1 ) = 1

These are the values that are produced by the P= option in the OUTPUT statement. Note that the predicted probability of an individual response value is that value's cumulative probability minus the cumulative probability of the next-lower ordered value. When PROC SORT and the ORDER=DATA option are used to put the response values in descending order, the next-lower ordered value might be a logically higher value. These individual values are not computed by PROC PROBIT, but they can be produced by the PREDPROBS=I option in the OUTPUT statement of PROC LOGISTIC, or they can be computed as shown below.

The following statements compute both the cumulative and individual predicted values for each level of the predictor variable.

data recomp;
   do addquant=1 to 4;
   /* cumulative probabilities */
      cp9=probnorm(-2.3624248815 +                0.2757618404*addquant);
      cp8=probnorm(-2.3624248815 + 0.7150632656 + 0.2757618404*addquant);
      cp7=probnorm(-2.3624248815 + 1.3283837985 + 0.2757618404*addquant);
      cp6=probnorm(-2.3624248815 + 1.6826398116 + 0.2757618404*addquant);
      cp5=probnorm(-2.3624248815 + 2.2169522954 + 0.2757618404*addquant);
      cp4=probnorm(-2.3624248815 + 2.6570246437 + 0.2757618404*addquant);
      cp3=probnorm(-2.3624248815 + 3.1187431200 + 0.2757618404*addquant);
      cp2=probnorm(-2.3624248815 + 3.5635390005 + 0.2757618404*addquant);
      cp1=1;
   /* individual probabilities */
      ip9=cp9-0;
      ip8=cp8-cp9;
      ip7=cp7-cp8;
      ip6=cp6-cp7;
      ip5=cp5-cp6;
      ip4=cp4-cp5;
      ip3=cp3-cp4;
      ip2=cp2-cp3;
      ip1=cp1-cp2;
      output;
   end;
   run;

The following tables show that the previous program reproduces the values from the P= option in PROC PROBIT. First, the cumulative values that are produced by the P= option are displayed, followed by the cumulative values that are computed and stored in the data set RECOMP:

proc print data=out noobs;
   where y = 1;
   var addquant _level_ cp;
   run;

  proc print data=recomp noobs;
     var addquant cp:;
     run;

Next, the individual values that are computed and stored in the data set RECOMP are displayed:

    proc print data=recomp noobs;
       var addquant ip:;
       run;