Multinomial logistic regression - multiple dependent variable group comparisons in the same model

[vpalzes]

I am trying to run a multinomial logistic regression model in SAS using PROC LOGISTIC and would like to know if it is possible to produce multiple dependent variable group comparisons in the same single model.

I am predicting the odds that an individual is in an alcohol use group (see groups below) with a few predictor variables (e.g., age, gender, race/ethnicity, and whether they have asthma). I am primarily interested in calculating the odds ratio of being in the abstinent group (no alcohol use) in people with asthma versus people without asthma.

Alcohol use group ~ Age + Gender + Race + Asthma

Alcohol use group

0 = abstinent

1 = moderate

2 = excessive

3 = very excessive

4 = hazardous

I use the following code to estimate odds ratios with the abstinent group as my reference group (1 vs 0, 2 vs 0, 3 vs 0, 4 vs 0). However, I want to also make these group comparisons:

1) 0 vs 1

2) 0 vs 2/3/4

3) 1 vs 2

4) 1 vs 3

5) 1 vs 4

proc logistic data=&mydata;
    class female(ref='0')

             agegrp(ref='1')

             racegrp6(ref='1')

             asthma(ref='0')
             / param=ref;
    model alcgrp(ref='0') = female agegrp racegrp6 asthma / link=glogit expb cl;
    estimate "log asthma_1" asthma 1 /exp cl category="1";
    estimate "log asthma_1" asthma 1 /exp cl category="2";

    estimate "log asthma_1" asthma 1 /exp cl category="3";

    estimate "log asthma_1" asthma 1 /exp cl category="4";
run;

I basically want to change my reference group for each comparison. I don't know how to do this without running multiple models where each one has a different reference group (code in bold above). I also want to group people with excessive drinking or higher into one group for one comparison, so ideally I would not be running two different sets of models (one with a 3-level dependent variable, and one with a 4- or 5-level dependent variable).

These are my primary questions:

Q1) Is there a way to calculate the inverse of an odds ratio and its corresponding 95% confidence interval and p-value using the parameters estimated from a single model in SAS? I could calculate these by hand, but it would be useless for conducting the significance z-test. For example, my model gives me the odds ratio comparing 1 vs 0, but how do I test the inverse of 0 vs 1, without running another model?

Q2) How do I create a contrast for multinomial logistic regression models in SAS where I would be "adding the probability" of being in multiple groups? For example, if I have a 5-level dependent variable, how do I use the parameters and collapse across levels to generate the test I want (0 vs 2/3/4)?

[StatDave]

Computing a similar, and somewhat more intuitive and useful statistic, is discussed and illustrated in this note. See the second section that discusses comparing probabilities of two response levels.

[vpalzes]

Thanks for your help.

I want to also combine the probabilities of multiple response-level groups for one of my comparisons (e.g., group 0 vs groups 2-4). It seems like the NLMeans macro cannot handle such a contrast for ratios.

I'm also afraid that I'm getting much different OR estimates with this method versus the reference-coding method that I have used previously.

This method:

OR comparing 0 vs 1 = 1.65 (1.52-1.78)

OR comparing 1 vs 2 = 1.08 (1.01-1.14)

OR comparing 1 vs 3 = 1.22 (1.15-1.29)

OR comparing 1 vs 4 = 1.50 (1.33-1.66)

Previous reference-coding method with multiple models (changing the reference group):

OR comparing 0 vs 1 = 1.25 (1.20-1.28)

OR comparing 1 vs 2 = 0.99 (0.93-1.06)

OR comparing 1 vs 3 = 0.99 (0.93-1.05)

OR comparing 1 vs 4 = 1.09 (0.98-1.20)

Here is my code based on the tutorial here:

proc logistic data=&mydata;
    class female agegrp racegrp6 asthma / param=glm;
    model alcgrp = female agegrp racegrp6 asthma / link=glogit;
    lsmeans asthma / e ilink;
    ods output coef=coeffs;
    output out=LogOR xbeta=xb stdxbeta=s;
    store out=logmod;
run;

data cont;
    length label $40;
    infile datalines missover;
    input label k1-k10;
    set=1;
    datalines;
P(dx,0)/P(dx,1)    0  0  0  0  0    1 -1  0  0  0
P(dx,1)/P(dx,2)    0  0  0  0  0    0  1 -1  0  0
P(dx,1)/P(dx,3)    0  0  0  0  0    0  1  0 -1  0
P(dx,1)/P(dx,4)    0  0  0  0  0    0  1  0  0 -1
;

%NLMeans(instore=logmod, coef=coeffs, link=glogit, options=ratio, contrasts=cont, title=Alcohol Group ORs for Patients with Asthma);

NOTE: These contrast lines can only have one 1 and one -1. For example, I cannot do:

P(dx,0)/P(dx,2-4)     0 0 0 0 0     1 0 -0.334 -0.333 -0.333

Otherwise, I recieve this error: 

ERROR: OPTIONS=RATIO requires all contrasts to contain a single 1

for the numerator, a single -1 for the denominator, and

zeros otherwise.

ERROR: The above message is for estimate set 1.

[StatDave]

If you want to contrast the 0 and 2-4 levels of response, then you really just want a binary response model. The easiest way would be to create a dichotomized version of your response and fit a binary logistic model to it.  However, if you are determined to do it in the context of the multinomial model, then you can use the NLEstimate macro to estimate a suitable expression for the desired function. Using the schools and instruction styles example in the note I referred to, suppose you want to estimate the relative risk that contrasts the self and team styles vs. the class style in school 1. That is a dichotomization of the response levels like what you want. Writing this in terms of the response probabilities, suppose you want to estimate the relative risk (p_self+p_team)/p_class in school 1. Writing this in terms of the parameters of the multinomial model, the function to estimate is

exp(intercept_self+school1_self) + exp(intercept_team+school1_team)

which can then be used in the FDATA= data set used by NLEstimate:

      [P(self,S1)+P(team,S1)]/P(class,S1) | exp(b_p1+b_p3)+exp(b_p2+b_p4)

You mentioned getting different OR estimates… The NLMeans macro with options=ratio will estimate relative risks as stated in the note I referred to, not odds ratios. If you want to estimate odds ratios, you will need to estimate suitable expressions in the NLEstimate macro.