Providing software solutions since 1976

• Log in
• Create profile
• Welcome
  [Log out]
• Edit profile
• My SAS

• Print  |
• Feedback  |
• Share

  • E-mail
  

  E-mail This Page

  Your Name:*
  Your E-mail:*
  Recipient E-mail:*
  Send a copy to me
  Personal Message:
  
  

  * indicates a required field

PROC LOGISTIC automatically computes the proportional odds test when the response is ordinal and the default logit link is used. For such a response, several cumulative logits are simultaneously modeled while only a single logit is modeled for a binary response. Conceptually, you could fit a model that has a complete set of parameter estimates for each of the cumulative logits. The simpler model that PROC LOGISTIC fits constrains each predictor's parameter estimates to be the same across all of the logits. This means that the fitted surfaces for the logits are all parallel and they are only allowed to differ by a constant shift that necessitates the separate intercepts that you get when you fit an ordinal model. When the logit link is used, this parallelism assumption also implies that the effect of a given predictor is the same regardless of where you "cut" the ordinal response to dichotomize it. The proportional odds test in PROC LOGISTIC simply tests whether the parameters are the same across logits, simultaneously for all predictors. PROC GENMOD fits the same proportional odds model, but it does not provide a proportional odds test.

Peterson and Harrell (1990) note that this proportional odds/parallelism test is known to be liberal when the sample size is small, which means that the p -value for the test could be artificially too small in small samples, possibly leading to inappropriate rejection of the proportional odds assumption. Because of this, graphical assessment of the parallelism assumption is useful and is discussed and illustrated in this note. However, nonsignificant test results (large p -values) are not affected and still indicate insufficient evidence of nonparallelism. The book by Stokes, et. al. suggests doing crosstabulations of the response with each predictor involved in the model. If all cell counts are about five or larger, then the sample size should be adequate.

Alternative models for nonproportional odds

If the sample size is adequate and the test or the graphs convince you to reject the assumption of proportional odds, then there are two alternative models that are often used. The simplest to fit is the generalized logit model  which can be fit by simply specifying the LINK=GLOGIT option in the MODEL statement of PROC LOGISTIC. This model treats the response as nominal (unordered) rather than ordinal and has a full set of parameters for each generalized logit.

The second alternative is the partial proportional odds model . This model continues to treat the response as ordinal, but allows you to assume proportional odds for some predictors while not for others. To do this, the model contains separate parameters across the logits for model effects exhibiting a lack of proportionality. A single parameter is used for effects where proportionality holds. At the extreme, the fully nonproportional odds model can be fit in which separate parameters are used for every effect in the model. There are two ways to fit this model. The first uses the Generalized Estimating Equations (GEE) method via the REPEATED statement in PROC GENMOD. The second uses PROC NLMIXED. Since these two procedures use different model parameterizations and estimation methods, the results may not be identical.

The GENMOD approach requires reshaping your data to have multiple observations per original observation corresponding to the multiple cumulative logits. The resulting groups of multiple observations form the clusters of the GEE analysis. This approach is discussed and illustrated in detail in an example from Chapter 15 (beginning at the DATA DENT statement) of the text by Stokes et. al. and is only summarized here. In this example, proportional and partial proportional odds models are fit to data from a study of a new analgesic for dental pain relief. Increasing relief from pain was recorded on a 5-point scale (0-4). PROC LOGISTIC is used to fit the proportional odds model involving three categorical predictors: the research center (CENTER=1 or 2), the treatment dosage (TRT=ACL, TL, ACH, TH, or P), and baseline pain (BASELINE=0 or 1). The proportional odds test is significant and tables (produced with PROC FREQ) of each predictor against the response show generally adequate sample size. The DENT2 data set is created as described above and used in PROC GENMOD to fit a fully unconstrained (nonproportional odds) model. The interaction of LOGTYPE with a predictor enables the predictor's effect to vary across logits. The nonsignificant test for one of the interactions indicates that this predictor satisfies the proportional odds (equal parameters) assumption so that a final model constrains this predictor but leaves the others unconstrained.

Fitting the proportional odds model in PROC NLMIXED

The NLMIXED approach does not require data manipulation but the model must be written showing each parameter and statements must be included defining the log likelihood function. Since no CLASS statement is available, you must create dummy variables that will represent any categorical predictor.

The following statements illustrate fitting the proportional odds model to the dental data mentioned above using PROC LOGISTIC and PROC NLMIXED. In PROC LOGISTIC, the CLASS statement creates reference-coded dummy variables for each of the three categorical predictors. The ORDER=DATA option causes predictor levels to be ordered as they first appear in the data set. The DESCENDING response variable option allows you to model the probability of the higher response levels^NOTE. The OUTPUT statement saves the cumulative predicted probabilities and the individual predicted probabilities to data set LOG.

In PROC NLMIXED, the same reference-coded dummies (T1-T4, B, and C) are created . Four cumulative probabilities (CP1-CP4) are computed on the five levels of the response (RESP=0, 1, 2, 3, 4) . The model is contained in these statements. INT1-INT4 are the four intercepts of the proportional odds model. Note that the same name is used across the four cumulative probabilities for each of the other model parameters (BC1, BB0, BACL, BTL, BACH, and BTH). Initial values for the intercepts are provided in the previous PARMS statement. The increasing initial values help ensure properly increasing cumulative probabilities and avoid estimation problems. All other parameters are started at the default value of one. The individual probabilities are computed from the cumulative probabilities by taking differences. For instance, the probability of response 3 is obtained by subtracting the probability of response 4 from the probability of response 3 or 4. These statements result in modeling the probability of higher response levels as in the LOGISTIC model^NOTE. The next three statements ensure that the predicted probability values are valid and define the log likelihood to be maximized by the procedure. The ID and PREDICT statements save the cumulative probabilities and the predicted probability of the observed response to data set NLM.

      proc logistic data=dent;
         class center baseline trt / param=ref order=data;
         model resp(descending)=center baseline trt;
         output out=log predprobs=(c i);
         run;
         
      proc nlmixed data=dent;
         parms Int4 = -1, Int3 = 0, Int2 = 1, Int1 = 2;

      t1=(trt='ACL'); t2=(trt='TL'); t3=(trt='ACH'); t4=(trt='TH');
         b=(baseline=0); c=(center=1);
            
      cp4= 1/(1 + exp(-(Int4 + bc1*c + bb0*b + bACL*t1 + bTL*t2 + bACH*t3 + bTH*t4)));
         cp3= 1/(1 + exp(-(Int3 + bc1*c + bb0*b + bACL*t1 + bTL*t2 + bACH*t3 + bTH*t4)));
         cp2= 1/(1 + exp(-(Int2 + bc1*c + bb0*b + bACL*t1 + bTL*t2 + bACH*t3 + bTH*t4)));
         cp1= 1/(1 + exp(-(Int1 + bc1*c + bb0*b + bACL*t1 + bTL*t2 + bACH*t3 + bTH*t4)));
      
      if      resp=4 then ip = cp4;         /* CP4 is Pr(resp=4) */
         else if resp=3 then ip = cp3-cp4;     /* CP3 is Pr(resp=4 or 3) */
         else if resp=2 then ip = cp2-cp3;     /* CP2 is Pr(resp=4 or 3 or 2) */
         else if resp=1 then ip = cp1-cp2;     /* CP1 is Pr(resp=4 or 3 or 2 or 1) */
         else                ip = 1-cp1;       /* IP  is Pr(resp=observed level) */
         
      p = (ip>0 and ip<=1)*ip + (ip<=0)*1e-8 + (ip>1);
         loglik = log(p);
         model resp ~ general(loglik);
         
      id cp1-cp4;
         predict ip out=nlm;
         run;

The test of the proportional odds assumption in PROC LOGISTIC is significant (p =.0089) indicating that proportional odds does not hold and suggesting that separate parameters are needed across the logits for at least one predictor. A visual assessment of the assumption is provided by plotting the empirical logits.

The parameter estimates from PROC NLMIXED agree closely with those from PROC LOGISTIC. The standard errors differ slightly since LOGISTIC constructs Wald tests while NLMIXED uses t-tests based on the available degrees of freedom. The predicted probabilities produced by the two procedures (in data sets LOG and NLM) are also very similar.

Fitting the fully nonproportional odds model in PROC NLMIXED

These statements fit the fully nonproportional odds model which allows you to test for proportionality in each of the model effects and overall. Note that unique names are used across the cumulative probabilities for each model effect. For instance, instead of the common parameter BC1 for the CENTER effect in the proportional odds model above, there are now four parameters, BC11, BC12, BC13, and BC14. Similarly for the BASELINE and TRT effects. The four CONTRAST statements with "PO" in their labels provide proportional odds tests for each of the predictors as well as an overall test comparable to the test provided by PROC LOGISTIC. Preceding these, additional contrasts are done to test the effect of each predictor on logit 1 which contrasts response level 0 (no relief) with the higher response levels representing varying degrees of relief.

      proc nlmixed data=dent;
         parms Int4 = -1, Int3 = 0, Int2 = 1, Int1 = 2;
      
         t1=(trt='ACL'); t2=(trt='TL'); t3=(trt='ACH'); t4=(trt='TH');
         b=(baseline=0); c=(center=1);
      
      cp4= 1/(1 + exp(-(Int4 + bc14*c + bb04*b + bACL4*t1 + bTL4*t2 + bACH4*t3 + bTH4*t4)));
         cp3= 1/(1 + exp(-(Int3 + bc13*c + bb03*b + bACL3*t1 + bTL3*t2 + bACH3*t3 + bTH3*t4)));
         cp2= 1/(1 + exp(-(Int2 + bc12*c + bb02*b + bACL2*t1 + bTL2*t2 + bACH2*t3 + bTH2*t4)));
         cp1= 1/(1 + exp(-(Int1 + bc11*c + bb01*b + bACL1*t1 + bTL1*t2 + bACH1*t3 + bTH1*t4)));
      
         if      resp=4 then ip = cp4;
         else if resp=3 then ip = cp3-cp4;
         else if resp=2 then ip = cp2-cp3;
         else if resp=1 then ip = cp1-cp2;
         else                ip = 1-cp1;
         
         p = (ip>0 and ip<=1)*ip + (ip<=0)*1e-8 + (ip>1);
         loglik = log(p);
         model resp ~ general(loglik);
         
         id cp1-cp4;
         predict ip out=nlm;
         
         contrast 'center on logit1 (4321 vs 0)' bc11;
         contrast 'baseline on logit1 (4321 vs 0)' bb01;
         contrast 'trt on logit1 (4321 vs 0)' 
                       bACL1, bTL1, bACH1, bTH1;
      contrast 'center PO' 
                       bc11-bc14, bc12-bc14, bc13-bc14;
         contrast 'baseline PO' 
                       bb01-bb04, bb02-bb04, bb03-bb04;
         contrast 'trt PO' 
                       bACL1-bACL4, bACL2-bACL4, bACL3-bACL4,
                       bTL1-bTL4, bTL2-bTL4, bTL3-bTL4,
                       bACH1-bACH4, bACH2-bACH4, bACH3-bACH4,
                       bTH1-bTH4, bTH2-bTH4, bTH3-bTH4;
         contrast 'overall PO' 
                       bc11-bc14, bc12-bc14, bc13-bc14,
                       bb01-bb04, bb02-bb04, bb03-bb04,
                       bACL1-bACL4, bACL2-bACL4, bACL3-bACL4,
                       bTL1-bTL4, bTL2-bTL4, bTL3-bTL4,
                       bACH1-bACH4, bACH2-bACH4, bACH3-bACH4,
                       bTH1-bTH4, bTH2-bTH4, bTH3-bTH4;
         run;

Notice that the parameter estimates of the fully nonproportional odds model has 4 sets of parameter estimates compared to the single set in the proportional odds model. The logit to which the parameter estimate applies is given in the last character of its name. For instance, the estimates Int4, bc14, bb04, bACL4, bTL4, bACH4, and bTH4 are the estimates associated with the fourth logit which contrasts the highest level of relief against all other levels.

The overall test of proportional odds is similar to that from PROC LOGISTIC (p =.0059). Of the three predictors, only the test of the TRT effect suggests the possibility of nonproportionality (p =.0761). This seems consistent with the plots of the empirical logits. The smaller number of degrees of freedom available to the proportionality tests for the other two effects make it harder to detect nonproportionality. The large number of degrees of freedom for the overall test provides increased statistical power to find nonproportionality in the model overall.

Fitting a partial proportional odds model in PROC NLMIXED

If we are willing to assume that proportional odds holds for the BASELINE and CENTER predictors, then a partial proportional odds model can be fit which allows for nonproportionality only in the TRT predictor. A common parameter is used across the logits for BASELINE (BB0) and CENTER (BC1) while separate parameters are used for the TRT levels. The following statements fit this partial proportional odds model. The CONTRAST statements test the common effects of BASELINE and CENTER, the effect of TRT on no relief (logit 1), and the proportionality for TRT.

      proc nlmixed data=dent;
         parms Int4 = -1, Int3 = 0, Int2 = 1, Int1 = 2;
      
         t1=(trt='ACL'); t2=(trt='TL'); t3=(trt='ACH'); t4=(trt='TH');
         b=(baseline=0); c=(center=1);
      
         cp4= 1/(1 + exp(-(Int4 + bc1*c + bb0*b + bACL4*t1 + bTL4*t2 + bACH4*t3 + bTH4*t4)));
         cp3= 1/(1 + exp(-(Int3 + bc1*c + bb0*b + bACL3*t1 + bTL3*t2 + bACH3*t3 + bTH3*t4)));
         cp2= 1/(1 + exp(-(Int2 + bc1*c + bb0*b + bACL2*t1 + bTL2*t2 + bACH2*t3 + bTH2*t4)));
         cp1= 1/(1 + exp(-(Int1 + bc1*c + bb0*b + bACL1*t1 + bTL1*t2 + bACH1*t3 + bTH1*t4)));
      
         if      resp=4 then ip = cp4;
         else if resp=3 then ip = cp3-cp4;
         else if resp=2 then ip = cp2-cp3;
         else if resp=1 then ip = cp1-cp2;
         else                ip = 1-cp1;
         
         p = (ip>0 and ip<=1)*ip + (ip<=0)*1e-8 + (ip>1);
         loglik = log(p);
         model resp ~ general(loglik);
         
         id cp1-cp4;
         predict ip out=nlm;
         contrast 'center all logits' bc1;
         contrast 'baseline all logits' bb0;
         contrast 'trt on logit1 (4321 vs 0)' 
                       bACL1, bTL1, bACH1, bTH1;
         contrast 'trt PO' 
                       bACL1-bACL4, bACL2-bACL4, bACL3-bACL4,
                       bTL1-bTL4, bTL2-bTL4, bTL3-bTL4,
                       bACH1-bACH4, bACH2-bACH4, bACH3-bACH4,
                       bTH1-bTH4, bTH2-bTH4, bTH3-bTH4;
         run;

The table of parameter estimates shows the separate intercepts and separate TRT parameters for each logit, but a single parameter for each of BASELINE and CENTER.

The contrasts for BASELINE and CENTER essentially repeat the test of their single parameters above (though this time with F tests). The tests of TRT on logit 1 and for proportionality are largely unchanged from the fully nonproportional odds model.

__________

Peterson, B. and F. Harrell (1990), "Partial Proportional Odds Models for Ordinal Response Variables," Applied Statistics , 39:205-217.

Stokes, M.E., Davis, C.S., and Koch, G.G.

__________

NOTE: If you want to model the probability of lower, rather than higher, response levels omit the DESCENDING response variable option in the MODEL statement of PROC LOGISTIC. In PROC NLMIXED, you could simply use increasing values of RESP in the statements defining the individual probability. For instance:

         if      resp=0 then ip = cp4;
         else if resp=1 then ip = cp3-cp4;
         else if resp=2 then ip = cp2-cp3;
         else if resp=3 then ip = cp1-cp2;
         else                ip = 1-cp1;

But the results will be more meaningfully labeled if you also change the names of the parameters and cumulative probabilities as below for the proportional odds model:

     proc nlmixed data=dent;
         parms Int0 = -1, Int1 = 0, Int2 = 1, Int3 = 2;

         t1=(trt='ACL'); t2=(trt='TL'); t3=(trt='ACH'); t4=(trt='TH');
         b=(baseline=0); c=(center=1);

         cp0= 1/(1 + exp(-(Int0 + bc1*c + bb0*b + bACL*t1 + bTL*t2 + bACH*t3 + bTH*t4)));
         cp1= 1/(1 + exp(-(Int1 + bc1*c + bb0*b + bACL*t1 + bTL*t2 + bACH*t3 + bTH*t4)));
         cp2= 1/(1 + exp(-(Int2 + bc1*c + bb0*b + bACL*t1 + bTL*t2 + bACH*t3 + bTH*t4)));
         cp3= 1/(1 + exp(-(Int3 + bc1*c + bb0*b + bACL*t1 + bTL*t2 + bACH*t3 + bTH*t4)));

         if      resp=0 then ip = cp0;         /* CP0 is Pr(resp=0) */
         else if resp=1 then ip = cp1-cp0;     /* CP1 is Pr(resp=0 or 1) */
         else if resp=2 then ip = cp2-cp1;     /* CP2 is Pr(resp=0 or 1 or 2) */
         else if resp=3 then ip = cp3-cp2;     /* CP3 is Pr(resp=0 or 1 or 2 or 3) */
         else                ip = 1-cp3;       /* IP  is Pr(resp=observed level) */

         p = (ip>0 and ip<=1)*ip + (ip<=0)*1e-8 + (ip>1);
         loglik = log(p);
         model resp ~ general(loglik);

         id cp1-cp4;
         predict ip out=nlm;
         run;

Operating System and Release Information