The LOGISTIC Procedure

Cameron and Trivedi (1998, p. 68) studied the number of doctor visits from the Australian Health Survey 1977–78. The data set contains a dependent variable, dvisits, which contains the number of doctor visits in the past two weeks (0, 1, or 2, where 2 represents two or more visits) and the following explanatory variables: sex, which indicates whether the patient is female; age, which contains the patient’s age in years divided by 100; income, which contains the patient’s annual income (in units of $10,000); levyplus, which indicates whether the patient has private health insurance; freepoor, which indicates that the patient has free government health insurance due to low income; freerepa, which indicates that the patient has free government health insurance for other reasons; illness, which contains the number of illnesses in the past two weeks; actdays, which contains the number of days the illness caused reduced activity; hscore, which is a questionnaire score; chcond1, which indicates a chronic condition that does not limit activity; and chcond2, which indicates a chronic condition that limits activity.

data docvisit;
   input sex age agesq income levyplus freepoor freerepa
         illness actdays hscore chcond1 chcond2 dvisits;
   if ( dvisits > 2) then dvisits = 2;
   datalines;
1 0.19 0.0361 0.55  1  0  0  1  4  1  0  0  1
1 0.19 0.0361 0.45  1  0  0  1  2  1  0  0  1
0 0.19 0.0361 0.90  0  0  0  3  0  0  0  0  1

   ... more lines ...   

1 0.37 0.1369 0.25  0  0  1  1  0  1  0  0  0
1 0.52 0.2704 0.65  0  0  0  0  0  0  0  0  0
0 0.72 0.5184 0.25  0  0  1  0  0  0  0  0  0
;

Because the response variable dvisits has three levels, the proportional odds model constructs two response functions. There is an intercept parameter for each of the two response functions, , and common slope parameters across the functions. The model can be written as

The following statements fit a proportional odds model to this data:

proc logistic data=docvisit;
   model dvisits = sex age agesq income levyplus
                   freepoor freerepa illness actdays hscore
                   chcond1 chcond2;
run;

Selected results are displayed in Output 60.18.1.

The test of the proportional odds assumption in Output 60.18.1 rejects the null hypothesis that all the slopes are equal across the two response functions. This test is very anticonservative; that is, it tends to reject the null hypothesis even when the proportional odds assumption is reasonable.

The proportional odds assumption for ordinal response models can be relaxed by specifying the UNEQUALSLOPES option in the MODEL statement. A general (nonproportional odds) model has different slope parameters for every logit i:

The following statements fit the general model:

proc logistic data=docvisit;
   model dvisits = sex age agesq income levyplus
                   freepoor freerepa illness actdays hscore
                   chcond1 chcond2 / unequalslopes;
   sex:      test sex_0     =sex_1;
   age:      test age_0     =age_1;
   agesq:    test agesq_0   =agesq_1;
   income:   test income_0  =income_1;
   levyplus: test levyplus_0=levyplus_1;
   freepoor: test freepoor_0=freepoor_1;
   freerepa: test freerepa_0=freerepa_1;
   illness:  test illness_0 =illness_1;
   actdays:  test actdays_0 =actdays_1;
   hscore:   test hscore_0  =hscore_1;
   chcond1:  test chcond1_0 =chcond1_1;
   chcond2:  test chcond2_0 =chcond2_1;
run;

The TEST statements test the proportional odds assumption for each of the covariates in the model. The parameter names are constructed by appending the response level that identifies the response function, as described in the section Parameter Names in the OUTEST= Data Set. Selected results from fitting the general model to the data are displayed in Output 60.18.2.

The preceding general model fits slope parameters, and according to the "Linear Hypotheses Testing Results" table in Output 60.18.2, several variables are unnecessarily allowing nonproportional odds. You can obtain a more parsimonious model by specifying a subset of the parameters to have nonproportional odds. The following statements allow the parameters for the variables in the "Linear Hypotheses Testing Results" table that have p-values less than 0.1 (actdays, agesq, and income) to vary across the response functions:

proc logistic data=docvisit;
   model dvisits= sex age agesq income levyplus freepoor
         freerepa illness actdays hscore chcond1 chcond2
   / unequalslopes=(actdays agesq income);
run;

Selected results from fitting this partial proportional odds model are displayed in Output 60.18.3.

You can also specify the following code to let stepwise selection determine which parameters have unequal slopes:

proc logistic data=docvisit;
   model dvisits= sex age agesq income levyplus freepoor
         freerepa illness actdays hscore chcond1 chcond2
   / equalslopes unequalslopes selection=stepwise details;
run;

This selection process chooses sex, freepoor, illness, and hscore as the proportional odds effects and chooses U_actdays and U_agesq as the unconstrained general effects. For more information about model selection for partial proportional odds models, see Derr (2013).

The partial proportional odds model can be written in the same form as the general model by letting and . So the first q parameters have proportional odds and the remaining parameters do not. The last 12 – q parameters can be rewritten to have a common slope plus an increment from the common slope: , where the new parameters contain the increments from the common slopes. The model in this form makes it obvious that the proportional odds model is a submodel of the partial proportional odds model, and both of these are submodels of the general model. This means that you can use likelihood ratio tests to compare models.

You can use the following statements to compute the likelihood ratio tests from the Likelihood Ratio row of the "Testing Global Null hypothesis: BETA=0" tables in the preceding outputs:

data a;
   label p='Pr>ChiSq';
   format p 8.6;
   input Test $10. ChiSq1 DF1 ChiSq2 DF2;
   ChiSq= ChiSq1-ChiSq2;
   DF= DF1-DF2;
   p=1-probchi(ChiSq,DF);
   keep Test Chisq DF p;
   datalines;
Gen vs PO     761.4797 24    734.2971 12
PPO vs PO     752.5512 15    734.2971 12
Gen vs PPO    761.4797 24    752.5512 15
;

proc print data=a label noobs;
   var Test ChiSq DF p;
run;

Therefore, you reject the proportional odds model in favor of both the general model and the partial proportional odds model, and the partial proportional odds model fits as well as the general model. The likelihood ratio test of the general model versus the proportional odds model is very similar to the score test of the proportional odds assumption in Output 60.18.1 because of the large sample size (Stokes, Davis, and Koch, 2000, p. 249).

Note: The proportional odds model has increasing intercepts, which ensures the increasing nature of the cumulative response functions. However, none of the parameters in the partial proportional odds or general models are constrained. Because of this, sometimes during the optimization process a predicted individual probability can be negative; the optimization continues because it might recover from this situation. Sometimes your final model will predict negative individual probabilities for some of the observations; in this case a message is displayed, and you should check your data for outliers and possibly redefine your model. Other times the model fits your data well, but if you try to score new data you can get negative individual probabilities. This means the model is not appropriate for the data you are trying to score, a message is displayed, and the estimates are set to missing.