Skip to content

Latest commit

 

History

History
333 lines (223 loc) · 55.4 KB

22954 - The PROC LOGISTIC proportional odds test and fitting a partial proportional odds model.md

File metadata and controls

333 lines (223 loc) · 55.4 KB

本文由 简悦 SimpRead 转码, 原文地址 support.sas.com

PROC LOGISTIC automatically computes a test of the proportional odds assumption when the response is ......

PROC LOGISTIC automatically computes a test of the proportional odds assumption 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 reliably indicate adequacy of the proportional odds/parallel assumption. Stokes, et. al. (2012) suggests doing cross-tabulations 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 first 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. Beginning in SAS® 9.3 TS1M2, you can fit this model using the UNEQUALSLOPES option in the MODEL statement of PROC LOGISTIC. In earlier releases, the model can also be fit using PROC NLMIXED as shown in the addendum at the end of this note.

In the following sections, proportional, fully nonproportional, 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). Initially, PROC LOGISTIC is used to fit a 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 (not shown) of each predictor against the response show generally adequate sample size. The alternative models are then fit.

Derr (2013) also discusses the testing and assessment of the proportional odds assumption. When the proportional odds assumption is rejected, he illustrates fitting the partial proportional odds model.

Fitting the proportional odds model

The following statements fit the proportional odds model to the dental data using PROC LOGISTIC. 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 levelsNOTE. The OUTPUT statement saves the cumulative predicted probabilities and the individual predicted probabilities to data set LOG.

      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;


The test of the proportional odds assumption in PROC LOGISTIC is significant (p=0.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.

Score Test for the Proportional
Odds Assumption
Chi-SquareDFPr > ChiSq
35.2185180.0089
 
Analysis of Maximum Likelihood Estimates
Parameter DFEstimateStandard
Error
Wald
Chi-Square
Pr > ChiSq
Intercept41-2.80170.481433.8722<.0001
Intercept31-1.07320.44725.76000.0164
Intercept21-0.04340.44360.00960.9221
Intercept110.93300.44764.34400.0371
CENTER11-1.93880.255557.5778<.0001
BASELINE01-0.61370.33613.33420.0679
TRTACL11.53090.397814.81290.0001
TRTTL11.50450.399314.19400.0002
TRTACH11.82040.400220.6867<.0001
TRTTH12.02040.401025.3863<.0001

Fitting the fully nonproportional odds model

The addition of the UNEQUALSLOPES option causes PROC LOGISTIC to fit the fully nonproportional odds model which includes separate slope parameters for each effect on each of the logits. If you also include the EQUALSLOPES option (available beginning in SAS 9.4 TS1M2), the model includes both a common, reference slope (slope on last logit) for each effect as well as parameters representing the differences between the reference slope and the separate slopes. The results include joint tests of these sets of difference parameters in the "Type 3 Analysis of Effects" table. Each set is labeled U__effect_ which is a test of the proportional odds assumption for effect. These statements fit the fully nonproportional odds model, parametrized to test the proportional odds assumption in each model effect.

      proc logistic data=dent;
         class center baseline trt / param=ref order=data;
         model resp(descending) = center baseline trt / unequalslopes equalslopes;
         run;


The nonsignificant tests for the CENTER and BASELINE effects (p=0.1782 and p=0.1727, respectively) suggest that the proportional odds assumption is reasonable for those effects. The test of proportional odds for the TRT effect (p=0.0688) is marginally significant. This seems consistent with the plots of the empirical logits.

Type 3 Analysis of Effects
EffectDFWald
Chi-Square
Pr > ChiSq
CENTER142.7602<.0001
U_CENTER34.91380.1782
BASELINE12.67080.1022
U_BASELINE34.98710.1727
TRT432.1236<.0001
U_TRT1219.91220.0688

In the Parameter Estimates table, notice that the parameter estimates for each effect include a common, reference slope parameter (slope on RESP=1 logit) followed by three difference parameters for a total of four degrees of freedom. The three difference parameters for CENTER and BASELINE are all nonsignificant, as expected, but the difference parameters for TRT="ACH" show some significance which explains the marginally significant proportional odds test for the TRT effect.

Analysis of Maximum Likelihood Estimates
Parameter RESPDFEstimateStandard
Error
Wald
Chi-Square
Pr > ChiSq
Intercept 41-1.91430.79845.74870.0165
Intercept 31-0.87060.57982.25490.1332
Intercept 21-0.66240.53211.54950.2132
Intercept 111.19110.57244.33070.0374
CENTER1 1-2.23910.342442.7602<.0001
U_CENTER1410.75020.60001.56320.2112
U_CENTER1310.15330.39060.15420.6946
U_CENTER1210.46730.28992.59840.1070
U_CENTER1100...
BASELINE0 1-0.77440.47392.67080.1022
U_BASELINE041-0.35490.70360.25440.6140
U_BASELINE0310.25380.51070.24710.6191
U_BASELINE0210.68920.42482.63130.1048
U_BASELINE0100...
TRTACL 11.41900.46399.35570.0022
TRTTL 11.29230.47237.48480.0062
TRTACH 12.69770.533425.5815<.0001
TRTTH 12.28220.510220.0106<.0001
U_TRTACL41-0.20710.90750.05210.8195
U_TRTACL310.06440.56470.01300.9093
U_TRTACL210.16930.41800.16400.6855
U_TRTACL100...
U_TRTTL41-0.65491.02320.40970.5221
U_TRTTL310.31760.56390.31710.5733
U_TRTTL210.55260.39191.98790.1586
U_TRTTL100...
U_TRTACH41-1.73391.03342.81490.0934
U_TRTACH31-1.44940.65504.89590.0269
U_TRTACH21-1.01500.53553.59310.0580
U_TRTACH100...
U_TRTTH41-0.82381.00520.67170.4125
U_TRTTH31-0.69840.61561.28720.2566
U_TRTTH210.01290.46550.00080.9779
U_TRTTH100...

Note that you can also test for proportionality in each of the model effects and overall by omitting the EQUALSLOPES option and including appropriate TEST statements. When the UNEQUALSLOPES option appears without the EQUALSLOPES option, the parameters for each effect are the separate slopes on the logits rather than differences from the reference slope. The following four labeled TEST statements following the MODEL statement provide tests of proportionality. The first three TEST statements test for proportional odds in each of the three predictors. In the CENTER proportionality test (labeled CENTER_PO), the equality of the CENTER parameters across the four logits is tested. Similarly for the BASELINE proportionality test. For the TRT test, equality across the logits is tested simultaneously within each level of TRT. Parameter names used in the TEST statement are specified as described in "Parameter Names in the OUTEST= Data Set" in the Details section of the LOGISTIC documentation. The final TEST statement provides an overall test of proportional odds, similar to the test provided when the UNEQUALSLOPES option is not used. It combines the previous three tests into a single, joint test.

      proc logistic data=dent;
         class center baseline trt / param=ref order=data;
         model resp(descending) = center baseline trt / unequalslopes;
         CENTER_PO:test center1_4=center1_3=center1_2=center1_1;
         BASELINE_PO:test baseline0_4=baseline0_3=baseline0_2=baseline0_1;
         TRT_PO:test trtacl_4=trtacl_3=trtacl_2=trtacl_1,
                     trttl_4=trttl_3=trttl_2=trttl_1,
                     trtach_4=trtach_3=trtach_2=trtach_1,
                     trtth_4=trtth_3=trtth_2=trtth_1;
         OVERALL_PO:test center1_4=center1_3=center1_2=center1_1,
                         baseline0_4=baseline0_3=baseline0_2=baseline0_1,
                         trtacl_4=trtacl_3=trtacl_2=trtacl_1,
                         trttl_4=trttl_3=trttl_2=trttl_1,
                         trtach_4=trtach_3=trtach_2=trtach_1,
                         trtth_4=trtth_3=trtth_2=trtth_1;
         run;


In this parameterization of the fully nonproportional odds model, the separate slopes of each effect on each logit are given. The logit to which each parameter estimate applies is given in the RESP column. The value in the RESP column is the level after which the ordered response levels are cut to dichotomize the response. For instance, the value 4 indicates the logit comparing level 4 with levels 3, 2, and 1. The value 3 indicates the logit comparing levels 4 and 3 with levels 2 and 1, and so on.

Analysis of Maximum Likelihood Estimates
Parameter RESPDFEstimateStandard
Error
Wald
Chi-Square
Pr > ChiSq
Intercept 41-1.91430.79845.74870.0165
Intercept 31-0.87060.57982.25490.1332
Intercept 21-0.66240.53211.54950.2132
Intercept 111.19110.57244.33070.0374
CENTER141-1.48900.54667.42050.0064
CENTER131-2.08580.337238.2664<.0001
CENTER121-1.77190.290337.2503<.0001
CENTER111-2.23910.342442.7602<.0001
BASELINE041-1.12930.60613.47230.0624
BASELINE031-0.52060.42111.52800.2164
BASELINE021-0.08530.41040.04320.8354
BASELINE011-0.77440.47392.67080.1022
TRTACL411.21190.87381.92330.1655
TRTACL311.48330.55927.03590.0080
TRTACL211.58830.487410.62010.0011
TRTACL111.41900.46399.35570.0022
TRTTL410.63740.98530.41840.5177
TRTTL311.60980.55648.37160.0038
TRTTL211.84480.500813.56780.0002
TRTTL111.29230.47237.48480.0062
TRTACH410.96380.95551.01750.3131
TRTACH311.24830.56104.95120.0261
TRTACH211.68270.490011.79020.0006
TRTACH112.69770.533425.5815<.0001
TRTTH411.45830.94012.40620.1209
TRTTH311.58370.54758.36830.0038
TRTTH212.29510.496921.3320<.0001
TRTTH112.28220.510220.0106<.0001

The overall test of proportional odds from the fourth TEST statement (OVERALL_PO) is a Wald test and gives results (p=0.0038) similar to the score test provided with the proportional odds model above. As before, only the test of the TRT effect (TRT_PO) suggests the possibility of nonproportionality.

Linear Hypotheses Testing Results
LabelWald
Chi-Square
DFPr > ChiSq
CENTER_PO4.913830.1782
BASELINE_PO4.987130.1727
TRT_PO19.9122120.0688
OVERALL_PO38.0570180.0038

Fitting a partial proportional odds model

If you 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 by specifying the UNEQUALSLOPES=TRT option. A common parameter is used across the logits for BASELINE and CENTER while separate parameters are used for the TRT levels. The following statements fit this partial proportional odds model.

Additionally, suppose you want to test whether the TRT levels differ. Normally, this can be done using the DIFF option in an LSMEANS statement. However, this statement is not available when the UNEQUALSLOPES or EQUALSLOPES option is specified. Instead, TEST statements are used. The first four labeled TEST statements test for treatment differences in each of the four logits. The final TEST statement (TRT_PO) tests for proportional odds for TRT in this simplified model.

      proc logistic data=dent;
         class center baseline trt / param=ref order=data;
         model resp(descending)=center baseline trt / unequalslopes=trt;
         TRT_RESP4:test trtacl_4,trttl_4,trtach_4,trtth_4;
         TRT_RESP3:test trtacl_3,trttl_3,trtach_3,trtth_3;
         TRT_RESP2:test trtacl_2,trttl_2,trtach_2,trtth_2;
         TRT_RESP1:test trtacl_1,trttl_1,trtach_1,trtth_1;
         TRT_PO:test trtacl_4=trtacl_3=trtacl_2=trtacl_1,
                     trttl_4=trttl_3=trttl_2=trttl_1,
                     trtach_4=trtach_3=trtach_2=trtach_1,
                     trtth_4=trtth_3=trtth_2=trtth_1;
         run;


The table of parameter estimates shows the separate intercepts and separate TRT parameters for each logit, and the single parameter for each of BASELINE and CENTER in this partial proportional odds model.

Analysis of Maximum Likelihood Estimates
Parameter RESPDFEstimateStandard
Error
Wald
Chi-Square
Pr > ChiSq
Intercept 41-2.17280.78597.64410.0057
Intercept 31-0.89790.54692.69590.1006
Intercept 21-0.20450.49680.16940.6807
Intercept 110.87960.47233.46870.0625
CENTER1 1-1.97880.259558.1542<.0001
BASELINE0 1-0.56060.35372.51200.1130
TRTACL411.27630.85692.21850.1364
TRTACL311.46140.55656.89630.0086
TRTACL211.66510.483911.83820.0006
TRTACL111.33530.45018.80010.0030
TRTTL410.49600.94630.27470.6002
TRTTL311.66430.55249.07630.0026
TRTTL211.87280.488714.68870.0001
TRTTL111.17110.45196.71440.0096
TRTACH410.83200.90260.84970.3566
TRTACH311.31120.55945.49430.0191
TRTACH211.71050.486812.34900.0004
TRTACH112.57270.519524.5251<.0001
TRTTH411.28440.85492.25710.1330
TRTTH311.60810.54538.69720.0032
TRTTH212.33830.491822.6063<.0001
TRTTH112.21750.494220.1352<.0001

The tests for treatment differences suggest that the treatments differ on logits 1, 2, and 3, but not on logit 4 (p=0.4925). The test for proportional odds in the treatment effect (TEST_PO) continues to indicate possible lack of proportionality (p=0.0729).

Linear Hypotheses Testing Results
LabelWald
Chi-Square
DFPr > ChiSq
TRT_RESP43.404640.4925
TRT_RESP310.783440.0291
TRT_RESP224.38584<.0001
TRT_RESP131.94924<.0001
TRT_PO19.7037120.0729

Using SELECTION= to choose a partial proportional odds model

You can use the SELECTION= option to have PROC LOGISTIC determine which effects exhibit nonproportional odds and choose a final model. The following statements use stepwise effect selection to select a final model from a set of candidate effects that includes all equal and unequal slope parameters. Beginning in SAS 9.4 TS1M2, the EQUALSLOPES option makes the equal slope parameters available. Used with the UNEQUALSLOPES option, all equal and unequal slope parameters are available for effect selection. The INCLUDE=EQUALSLOPES option restricts the selection by requiring every model in the selection process to include the equal slope parameters. The selection process tests the unequal slope parameters for an effect and includes and retains them if significant at the 0.1 level. The DETAILS option shows the tests of the candidate effects at each step.

      proc logistic data=dent;
         class center baseline trt / param=ref order=data;
         model resp(descending)=center baseline trt / selection=stepwise slentry=0.1 slstay=0.1
            equalslopes unequalslopes include=equalslopes details;
         run;


The model in the first step includes the intercepts and equal slope parameters. This is the same proportional odds model as shown above. The analysis at this step includes tests of the remaining candidate effects, which are the unequal slope parameters for each of the three predictors. The U_ prefix indicates that these are the unequal slope parameters for each predictor. The score tests give similar results to the Wald tests for proportional odds shown above.

Analysis of Effects Eligible for Entry
EffectDFScore
Chi-Square
Pr > ChiSq
U_CENTER36.10380.1067
U_BASELINE35.61950.1317
U_TRT1221.86320.0391

Since the unequal slopes effect for TRT is significant at the 0.1 level, it is included in the next step of the model selection process. This results in the same partial proportional odds model shown above but is parameterized slightly differently with a set of equal slope parameters followed by sets of unequal slope parameters. For each TRT level, the equal slope parameter is used as the slope for the last logit (RESP=1), and the unequal slope parameters (U_TRT) are deviations from the equal slope parameter.

Analysis of Maximum Likelihood Estimates
Parameter RESPDFEstimateStandard
Error
Wald
Chi-Square
Pr > ChiSq
Intercept 41-2.17280.78597.64410.0057
Intercept 31-0.89790.54692.69590.1006
Intercept 21-0.20450.49680.16940.6807
Intercept 110.87960.47233.46870.0625
CENTER1 1-1.97880.259558.1542<.0001
BASELINE0 1-0.56060.35372.51200.1130
TRTACL 11.33530.45018.80010.0030
TRTTL 11.17110.45196.71440.0096
TRTACH 12.57270.519524.5251<.0001
TRTTH 12.21750.494220.1352<.0001
U_TRTACL41-0.05900.88640.00440.9469
U_TRTACL310.12610.55530.05160.8204
U_TRTACL210.32980.40850.65190.4194
U_TRTACL100...
U_TRTTL41-0.67510.97570.47870.4890
U_TRTTL310.49320.53800.84030.3593
U_TRTTL210.70180.37443.51360.0609
U_TRTTL100...
U_TRTACH41-1.74070.97613.18010.0745
U_TRTACH31-1.26150.64113.87250.0491
U_TRTACH21-0.86220.52482.69860.1004
U_TRTACH100...
U_TRTTH41-0.93310.91711.03520.3089
U_TRTTH31-0.60940.60311.02100.3123
U_TRTTH210.12080.45090.07170.7888
U_TRTTH100...

Since tests of the remaining nonproportional odds effects are not significant, the selection process ends.

Analysis of Effects Eligible for Entry
EffectDFScore
Chi-Square
Pr > ChiSq
U_CENTER35.64870.1300
U_BASELINE35.92370.1154

An additional example of fitting nonproportional odds models and using model selection can be found in the example titled "Partial Proportional Odds Model" in the LOGISTIC documentation.

References

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.

Addendum: Using NLMIXED to fit ordinal logistic models

The following section fits the same models as above using PROC NLMIXED.

Proportional odds model

In PROC NLMIXED, reference-coded dummies (T1-T4, B, and C) are created in the same way as done by the CLASS statement in PROC LOGISTIC. 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 modelNOTE. 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 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= logistic(Int4 + bc1*c + bb0*b + bACL*t1 + bTL*t2 + bACH*t3 + bTH*t4);
         cp3= logistic(Int3 + bc1*c + bb0*b + bACL*t1 + bTL*t2 + bACH*t3 + bTH*t4);
         cp2= logistic(Int2 + bc1*c + bb0*b + bACL*t1 + bTL*t2 + bACH*t3 + bTH*t4);
         cp1= logistic(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 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.

Fully nonproportional odds model

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.

The parameter estimates of the fully nonproportional odds model has four complete 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 parameter estimates associated with the fourth logit which contrasts the highest level of relief against all other levels.

      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= logistic(Int4 + bc14*c + bb04*b + bACL4*t1 + bTL4*t2 + bACH4*t3 + bTH4*t4);
         cp3= logistic(Int3 + bc13*c + bb03*b + bACL3*t1 + bTL3*t2 + bACH3*t3 + bTH3*t4);
         cp2= logistic(Int2 + bc12*c + bb02*b + bACL2*t1 + bTL2*t2 + bACH2*t3 + bTH2*t4);
         cp1= logistic(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;


Partial proportional odds model

A partial proportional odds model is fit with a common parameter across the logits for BASELINE (BB0) and CENTER (BC1) and separate parameters for the TRT levels. The following statements fit this partial proportional odds model in PROC NLMIXED.

      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= logistic(Int4 + bc1*c + bb0*b + bACL4*t1 + bTL4*t2 + bACH4*t3 + bTH4*t4);
         cp3= logistic(Int3 + bc1*c + bb0*b + bACL3*t1 + bTL3*t2 + bACH3*t3 + bTH3*t4);
         cp2= logistic(Int2 + bc1*c + bb0*b + bACL2*t1 + bTL2*t2 + bACH2*t3 + bTH2*t4);
         cp1= logistic(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;
         run;



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= logistic(Int0 + bc1*c + bb0*b + bACL*t1 + bTL*t2 + bACH*t3 + bTH*t4);
         cp1= logistic(Int1 + bc1*c + bb0*b + bACL*t1 + bTL*t2 + bACH*t3 + bTH*t4);
         cp2= logistic(Int2 + bc1*c + bb0*b + bACL*t1 + bTL*t2 + bACH*t3 + bTH*t4);
         cp3= logistic(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;


data dent; input patient center trt $ baseline ldose resp @@; datalines; 2 1 ACL 0 5.29832 0 131 1 TL 0 3.91202 2 1 1 ACH 1 5.99146 1 3 1 TH 0 4.60517 1 130 1 P 0 0.00000 0 132 1 P 0 0.00000 0 4 1 P 0 0.00000 0 133 1 P 0 0.00000 0 5 1 P 0 0.00000 0 134 2 ACH 0 5.99146 4 6 1 TL 1 3.91202 2 135 2 ACL 0 5.29832 4 7 1 ACH 0 5.99146 1 136 2 TH 0 4.60517 3 8 1 ACL 0 5.29832 0 137 2 ACL 0 5.29832 4 9 1 TL 1 3.91202 0 138 2 TL 0 3.91202 3 10 1 TL 1 3.91202 4 139 2 P 0 0.00000 4 11 1 ACL 0 5.29832 2 140 2 TL 0 3.91202 3 12 1 ACH 0 5.99146 0 141 2 TL 0 3.91202 3 13 1 P 0 0.00000 0 142 2 ACL 1 5.29832 3 14 1 TL 0 3.91202 0 143 2 ACH 0 5.99146 1 15 1 P 0 0.00000 0 144 2 ACH 0 5.99146 3 16 1 TH 1 4.60517 4 145 2 P 0 0.00000 1 17 1 TH 0 4.60517 2 146 2 P 0 0.00000 0 18 1 ACH 0 5.99146 1 147 2 ACH 0 5.99146 4 19 1 ACL 0 5.29832 0 148 2 TL 0 3.91202 2 20 1 TH 0 4.60517 3 149 2 TH 0 4.60517 3 21 1 P 1 0.00000 1 150 2 P 0 0.00000 0 22 1 TH 1 4.60517 0 151 2 TH 0 4.60517 2 23 1 TL 0 3.91202 2 152 2 ACL 1 5.29832 3 24 1 ACL 1 5.29832 0 153 2 TH 0 4.60517 2 25 1 P 0 0.00000 0 154 2 ACH 0 5.99146 3 26 1 ACH 0 5.99146 0 155 2 ACL 0 5.29832 1 27 1 ACL 0 5.29832 0 156 2 ACL 0 5.29832 0 28 1 P 0 0.00000 0 157 2 ACL 0 5.29832 3 29 1 ACH 0 5.99146 2 158 2 TH 0 4.60517 3 30 1 TL 0 3.91202 0 159 2 ACL 0 5.29832 1 31 1 P 0 0.00000 0 160 2 TL 0 3.91202 3 32 1 ACH 0 5.99146 1 161 2 P 0 0.00000 2 33 1 TL 0 3.91202 0 162 2 TH 0 4.60517 3 34 1 TH 1 4.60517 4 163 2 TH 0 4.60517 4 35 1 TL 0 3.91202 2 164 2 ACH 0 5.99146 3 36 1 ACH 0 5.99146 0 165 2 TH 0 4.60517 2 37 1 ACL 1 5.29832 1 166 2 P 0 0.00000 3 38 1 ACL 0 5.29832 3 167 2 ACH 0 5.99146 1 39 1 TH 0 4.60517 2 168 2 P 0 0.00000 2 40 1 TH 0 4.60517 0 169 2 TL 1 3.91202 2 41 1 ACL 0 5.29832 2 170 2 P 0 0.00000 0 42 1 TL 0 3.91202 3 171 2 TL 0 3.91202 0 43 1 ACL 0 5.29832 2 172 2 TL 0 3.91202 3 44 1 ACL 0 5.29832 0 173 2 ACH 1 5.99146 2 45 1 TH 0 4.60517 2 174 2 ACL 0 5.29832 3 46 1 ACH 0 5.99146 0 175 2 P 0 0.00000 4 47 1 P 0 0.00000 0 176 2 TL 1 3.91202 0 48 1 ACL 0 5.29832 0 177 2 ACH 0 5.99146 1 49 1 TL 0 3.91202 0 178 2 ACH 1 5.99146 2 50 1 TL 0 3.91202 2 179 2 ACL 0 5.29832 2 51 1 TL 0 3.91202 0 180 2 ACL 0 5.29832 2 52 1 ACH 0 5.99146 4 181 2 TH 0 4.60517 2 53 1 TH 0 4.60517 4 182 2 TL 0 3.91202 3 54 1 TH 0 4.60517 0 183 2 TH 0 4.60517 3 55 1 P 0 0.00000 0 184 2 ACH 0 5.99146 1 56 1 ACH 0 5.99146 3 185 2 P 0 0.00000 1 57 1 ACH 0 5.99146 2 186 2 ACL 0 5.29832 1 58 1 P 0 0.00000 0 187 2 TH 0 4.60517 2 59 1 TH 0 4.60517 0 188 2 ACH 1 5.99146 3 60 1 P 0 0.00000 0 189 2 TH 0 4.60517 2 61 1 TL 0 3.91202 1 190 2 TL 0 3.91202 3 62 1 P 0 0.00000 0 191 2 P 0 0.00000 1 63 1 TH 1 4.60517 1 192 2 TL 0 3.91202 3 64 1 TL 0 3.91202 0 193 2 P 0 0.00000 0 65 1 ACH 0 5.99146 2 194 2 TH 0 4.60517 2 66 1 ACL 0 5.29832 2 195 2 ACH 0 5.99146 4 67 1 P 0 0.00000 2 196 2 ACH 0 5.99146 2 68 1 TH 0 4.60517 1 197 2 ACL 0 5.29832 3 69 1 ACH 0 5.99146 0 198 2 P 0 0.00000 0 70 1 P 0 0.00000 0 199 2 P 0 0.00000 3 71 1 TL 0 3.91202 0 200 2 ACL 0 5.29832 0 72 1 ACH 0 5.99146 2 201 2 ACL 1 5.29832 4 73 1 P 0 0.00000 0 202 2 TH 0 4.60517 3 74 1 TL 0 3.91202 2 203 2 P 0 0.00000 1 75 1 TH 0 4.60517 2 204 2 TH 0 4.60517 1 76 1 ACL 1 5.29832 0 205 2 TH 0 4.60517 3 77 1 TH 1 4.60517 0 206 2 TL 0 3.91202 3 78 1 ACL 0 5.29832 0 207 2 TL 0 3.91202 3 79 1 ACL 1 5.29832 3 208 2 TL 0 3.91202 3 80 1 ACH 0 5.99146 2 209 2 ACL 0 5.29832 2 81 1 ACL 0 5.29832 0 210 2 ACH 0 5.99146 3 82 1 P 0 0.00000 0 211 2 TL 1 3.91202 1 83 1 TH 0 4.60517 0 212 2 ACH 0 5.99146 3 84 1 ACH 0 5.99146 1 213 2 P 0 0.00000 2 85 1 TL 0 3.91202 0 214 2 P 0 0.00000 0 86 1 TH 0 4.60517 3 215 2 TL 0 3.91202 0 87 1 ACH 0 5.99146 0 216 2 TH 0 4.60517 4 88 1 P 0 0.00000 0 217 2 ACH 0 5.99146 2 89 1 ACH 0 5.99146 1 218 2 P 0 0.00000 0 90 1 TL 0 3.91202 0 219 2 TH 0 4.60517 2 91 1 ACL 0 5.29832 1 220 2 TL 0 3.91202 1 92 1 TH 0 4.60517 0 221 2 ACH 0 5.99146 2 93 1 ACL 0 5.29832 1 222 2 TL 0 3.91202 4 94 1 TL 1 3.91202 1 223 2 TH 0 4.60517 2 95 1 TL 1 3.91202 3 224 2 TH 1 4.60517 1 96 1 P 0 0.00000 0 225 2 ACH 0 5.99146 1 97 1 TH 0 4.60517 0 226 2 ACL 0 5.29832 4 98 1 ACL 0 5.29832 0 227 2 P 1 0.00000 3 99 1 P 1 0.00000 0 228 2 ACL 0 5.29832 2 100 1 ACH 0 5.99146 2 229 2 TL 0 3.91202 0 101 1 TH 0 4.60517 0 230 2 ACL 0 5.29832 1 102 1 TL 0 3.91202 0 231 2 ACH 0 5.99146 3 103 1 ACL 0 5.29832 1 232 2 ACL 0 5.29832 3 104 1 TL 0 3.91202 0 233 2 P 0 0.00000 1 105 1 P 0 0.00000 1 234 2 ACL 0 5.29832 4 106 1 ACL 0 5.29832 0 235 2 ACH 0 5.99146 1 107 1 TH 1 4.60517 2 236 2 TH 0 4.60517 1 108 1 P 0 0.00000 0 237 2 ACL 0 5.29832 0 109 1 ACH 0 5.99146 1 238 2 ACL 1 5.29832 4 110 1 TH 1 4.60517 1 239 2 ACL 0 5.29832 3 111 1 TL 0 3.91202 0 240 2 P 0 0.00000 0 112 1 ACH 0 5.99146 0 241 2 P 1 0.00000 3 113 1 TL 0 3.91202 0 242 2 TL 0 3.91202 2 114 1 ACH 0 5.99146 1 243 2 P 0 0.00000 0 115 1 P 0 0.00000 0 244 2 TH 0 4.60517 3 116 1 ACL 0 5.29832 0 245 2 TL 1 3.91202 4 117 1 P 0 0.00000 0 246 2 ACH 1 5.99146 1 118 1 ACH 0 5.99146 3 247 2 P 0 0.00000 1 119 1 TH 0 4.60517 3 248 2 TH 0 4.60517 4 120 1 ACL 0 5.29832 2 249 2 TL 0 3.91202 0 121 1 TH 0 4.60517 2 250 2 TL 0 3.91202 3 122 1 TH 0 4.60517 1 251 2 ACH 0 5.99146 3 123 1 TL 0 3.91202 0 252 2 TH 0 4.60517 3 124 1 ACH 0 5.99146 0 253 2 ACH 0 5.99146 1 125 1 ACL 0 5.29832 0 254 2 TH 0 4.60517 3 126 1 TH 1 4.60517 0 255 2 ACL 0 5.29832 0 127 1 ACL 0 5.29832 0 256 2 TL 0 3.91202 3 128 1 ACH 0 5.99146 1 257 2 P 0 0.00000 1 129 1 ACL 0 5.29832 3 258 2 ACH 0 5.99146 3 ;

Operating System and Release Information

Product FamilyProductSystemSAS Release
ReportedFixed*
SAS SystemSAS/STATAlln/a

***** For software releases that are not yet generally available, the Fixed Release is the software release in which the problem is planned to be fixed.