Nominal Logistic Regression

Nominal Target Variable

It is easy to generalize the binary logistic regression model to the ordinal logistic regression model. However, we need to change the underlying model and math slightly to extend to nominal target variables.

In binary logistic regression we are calculating the probability that an observation has an event. In ordinal logistic regression we are calculating the probability that an observation has at most that event in an ordered list of outcomes. In nominal (or multinomial) logistic regression we are calculating the probability that an observation has a specific event in an unordered list of events.

We will be using the alligator data set to model the association between various factors and alligator food choices - fish, invertebrate, birds, reptiles, and other. The variables in the data set are the following:

Variable	Description
size	small (< 2.3 meters) or large (> 2.3 meters)
lake	lake captured (George, Hancock, Oklawaha, Trafford)
gender	male or female
count	number of observations with characteristics

Generalized Logit Model

Instead of modeling the typical logit from binary or ordinal logistic regression, we will model the generalized logits. These generalized logits are built off target variables with \(m\) categories. The first logit will summarize the first category compared to the reference category. The second logit will summarize the second category compared to the reference. This continues \(m-1\) times as the last logit compares the \(m-1\) category to the reference category. In essence, with a target variable with \(m\) categories, we are building \(m-1\) logistic regressions. However, unlike the ordinal logistic regression where we have different intercepts with the same slope parameters, in nominal logistic regression we have both different intercepts and different slope parameters for each model.

The logits are also not traditional logits. Instead of the natural log of the odds, they are natural logs of relative risks - ratios of two probabilities that do not sum to 1.

Let’s see how to build these generalized logit models in each of our softwares!

The data set is structured slightly differently than previous data sets. Let’s examine the first 10 observations with the head function and the n = 10 option.

Code

head(gator, n = 10)

            size         food    lake gender count
1  <= 2.3 meters         Fish Hancock   Male     7
2  <= 2.3 meters Invertebrate Hancock   Male     1
3  <= 2.3 meters        Other Hancock   Male     5
4   > 2.3 meters         Fish Hancock   Male     4
5   > 2.3 meters         Bird Hancock   Male     1
6   > 2.3 meters        Other Hancock   Male     2
7  <= 2.3 meters         Fish Hancock Female    16
8  <= 2.3 meters Invertebrate Hancock Female     3
9  <= 2.3 meters      Reptile Hancock Female     2
10 <= 2.3 meters         Bird Hancock Female     2

Notice the count variable. In the first observation we can see that we have seven alligators that are small, males captured on Lake Hancock and eat fish. The same applies for each row in our data set. Even with only 56 rows in our data set, we have 219 observations based on the counts.

To account for the multiple counts in a single row of our data set we use the weight = option in the multinom function with the count variable. The formula structure for the multinom function is the same as with previous functions in logistic regression. The target variable food needs to be turned into a factor with the factor function. We set the reference level specified with the relevel function. Here, our reference level is fish because it has the largest proportion of our observations. Similar to before, the summary function will display our needed results.

Code

gator$food <- factor(gator$food)
gator$food <- relevel(gator$food, ref = "Fish")

library(nnet)

glogit.model <- nnet::multinom(food ~ size + lake + gender, weight = count, data = gator)

# weights:  35 (24 variable)
initial  value 352.466903 
iter  10 value 271.274792
iter  20 value 268.935538
final  value 268.932740 
converged

Code

summary(glogit.model)

Call:
nnet::multinom(formula = food ~ size + lake + gender, data = gator, 
    weights = count)

Coefficients:
             (Intercept) size> 2.3 meters lakeHancock lakeOklawaha lakeTrafford
Bird          -2.4321397        0.7300740   0.5754699  -0.55020075     1.237216
Invertebrate   0.1690702       -1.3361658  -1.7805555   0.91304120     1.155722
Other         -1.4309095       -0.2905697   0.7667093   0.02603021     1.557820
Reptile       -3.4161432        0.5571846   1.1296426   2.53024945     3.061087
             genderMale
Bird         -0.6064035
Invertebrate -0.4629388
Other        -0.2524299
Reptile      -0.6276217

Std. Errors:
             (Intercept) size> 2.3 meters lakeHancock lakeOklawaha lakeTrafford
Bird           0.7706720        0.6522657   0.7952303    1.2098680    0.8661052
Invertebrate   0.3787475        0.4111827   0.6232075    0.4761068    0.4927795
Other          0.5381162        0.4599317   0.5685673    0.7777958    0.6256868
Reptile        1.0851582        0.6466092   1.1928075    1.1221413    1.1297557
             genderMale
Bird          0.6888385
Invertebrate  0.3955162
Other         0.4663546
Reptile       0.6852750

Residual Deviance: 537.8655 
AIC: 585.8655

Code

zvalues <- summary(glogit.model)$coefficients / summary(glogit.model)$standard.errors
pnorm(abs(zvalues), lower.tail=FALSE)*2

             (Intercept) size> 2.3 meters lakeHancock lakeOklawaha lakeTrafford
Bird         0.001600209      0.263016749  0.46927947   0.64928120  0.153153204
Invertebrate 0.655313448      0.001155806  0.00427554   0.05514608  0.019011127
Other        0.007834761      0.527539079  0.17749972   0.97330247  0.012782355
Reptile      0.001643580      0.388851410  0.34361579   0.02414336  0.006738219
             genderMale
Bird          0.3786819
Invertebrate  0.2418129
Other         0.5883123
Reptile       0.3597360

In the output above we see similar things to the output we have seen previously. R summarizes the parameter estimates and their standard errors. The main difference in the parameter estimates is the multiple intercepts and slope coefficients for each variable. Since we have five levels of our target variable we have four logistic regression models and therefore four intercepts and four slope parameters for each variable.

The data set is structured slightly differently than previous data sets. Let’s examine the first 10 observations with the head function and the n = 10 option.

Code

gator.head(n = 10)

            size          food     lake  gender  count
0  <= 2.3 meters          Fish  Hancock    Male      7
1  <= 2.3 meters  Invertebrate  Hancock    Male      1
2  <= 2.3 meters         Other  Hancock    Male      5
3   > 2.3 meters          Fish  Hancock    Male      4
4   > 2.3 meters          Bird  Hancock    Male      1
5   > 2.3 meters         Other  Hancock    Male      2
6  <= 2.3 meters          Fish  Hancock  Female     16
7  <= 2.3 meters  Invertebrate  Hancock  Female      3
8  <= 2.3 meters       Reptile  Hancock  Female      2
9  <= 2.3 meters          Bird  Hancock  Female      2

We need to expand out dataset to run the analysis in Python. To do this we just use the repeat and index attributes of the gator dataset. Let’s see the new dataset!

Code

gator = gator.loc[gator.index.repeat(gator['count'])]

gator = gator.reset_index(drop = True)
        
gator.head(n = 10)

            size          food     lake gender  count
0  <= 2.3 meters          Fish  Hancock   Male      7
1  <= 2.3 meters          Fish  Hancock   Male      7
2  <= 2.3 meters          Fish  Hancock   Male      7
3  <= 2.3 meters          Fish  Hancock   Male      7
4  <= 2.3 meters          Fish  Hancock   Male      7
5  <= 2.3 meters          Fish  Hancock   Male      7
6  <= 2.3 meters          Fish  Hancock   Male      7
7  <= 2.3 meters  Invertebrate  Hancock   Male      1
8  <= 2.3 meters         Other  Hancock   Male      5
9  <= 2.3 meters         Other  Hancock   Male      5

Now the dataset is in workable format lengthwise. However, the MNLogit function we need to use for nominal logistic regression doesn’t have a from_formula attribute. Therefore, we need to get all of our predictor variables into a workable format for modeling. Since we have categorical variables, we need to dummy code them. The from_formula framework in earlier sections did that for us before. To do that ourselves we use the get_dummies function on the three categorical variables - size, lake, and gender. We drop the target variable food from this function as well as the original count variable. This set of dummy encoded variables now serves as our X matrix. The target variable food is our y. These are the inputs into the MNLogit function. We use the typical summary function to view the results.

Code

import pandas as pd
import statsmodels.api as sma

y = gator['food']
X = pd.get_dummies(data = gator, columns = ['size', 'lake', 'gender'], drop_first = True).drop(columns = ['food', 'count'])

glogit_model = sma.MNLogit(y, X).fit()

Optimization terminated successfully.
         Current function value: 1.320942
         Iterations 7

Code

glogit_model.summary()

MNLogit Regression Results
Dep. Variable:	food	No. Observations:	219
Model:	MNLogit	Df Residuals:	199
Method:	MLE	Df Model:	16
Date:	Thu, 17 Aug 2023	Pseudo R-squ.:	0.04267
Time:	09:13:22	Log-Likelihood:	-289.29
converged:	True	LL-Null:	-302.18
Covariance Type:	nonrobust	LLR p-value:	0.05705

food=Fish	coef	std err	z	P>\|z\|	[0.025	0.975]
size_> 2.3 meters	0.1289	0.520	0.248	0.804	-0.890	1.147
lake_Hancock	1.2721	0.515	2.470	0.014	0.262	2.282
lake_Oklawaha	1.6785	1.087	1.544	0.123	-0.452	3.809
lake_Trafford	-0.2345	0.699	-0.335	0.737	-1.605	1.136
gender_Male	1.6597	0.529	3.137	0.002	0.623	2.697
food=Invertebrate	coef	std err	z	P>\|z\|	[0.025	0.975]
size_> 2.3 meters	-1.4893	0.607	-2.452	0.014	-2.680	-0.299
lake_Hancock	-0.2610	0.695	-0.375	0.707	-1.623	1.101
lake_Oklawaha	2.9920	1.077	2.778	0.005	0.881	5.103
lake_Trafford	1.3399	0.674	1.989	0.047	0.020	2.660
gender_Male	1.2814	0.564	2.274	0.023	0.177	2.386
food=Other	coef	std err	z	P>\|z\|	[0.025	0.975]
size_> 2.3 meters	-0.4827	0.585	-0.825	0.409	-1.629	0.664
lake_Hancock	0.9376	0.550	1.705	0.088	-0.140	2.015
lake_Oklawaha	0.9792	1.205	0.813	0.416	-1.383	3.341
lake_Trafford	0.6902	0.718	0.962	0.336	-0.716	2.097
gender_Male	0.7256	0.584	1.243	0.214	-0.419	1.870
food=Reptile	coef	std err	z	P>\|z\|	[0.025	0.975]
size_> 2.3 meters	-0.0520	0.604	-0.086	0.931	-1.236	1.132
lake_Hancock	-0.5174	0.748	-0.692	0.489	-1.983	0.949
lake_Oklawaha	1.9040	1.123	1.696	0.090	-0.297	4.105
lake_Trafford	0.6295	0.739	0.852	0.394	-0.818	2.077
gender_Male	0.1429	0.639	0.223	0.823	-1.110	1.396

In the output above we see similar things to the output we have seen previously. Python summarizes the parameter estimates and their standard errors. The main difference in the parameter estimates is the multiple intercepts and slope coefficients for each variable. Since we have five levels of our target variable we have four logistic regression models and therefore four intercepts and four slope parameters for each variable.

The data set is structured slightly differently than previous data sets. Let’s examine the first 10 observations with PROC PRINT and the obs = 10 option.

Code

proc print data=logistic.gator (obs = 10);
run;

To account for the multiple counts in a single row of our data set we use the FREQ statement in PROC LOGISTIC with the count variable. The CLASS statement still defines our categorical variables for lake, size, and gender. The target variable food has a reference level specified with the ref = option. Here, our reference level is fish because it has the largest proportion of our observations. The link = glogit option tells PROC LOGISTIC to use the generalized logit model.

Code

proc logistic data=Logistic.Gator plot(only)=oddsratio(range=clip);
    freq count;
    class lake(param=ref ref='George') size(param=ref ref='<= 2.3 meters') gender(param=ref ref='Male');
    model food(ref='Fish') = lake size gender / link=glogit clodds=pl;
    title 'Model on Alligator Food Choice';
run;
quit;

In the output above we see similar things to the output we have seen previously. PROC LOGISTIC summarizes the model used, the number of observations read and used in the model, and builds a frequency table for our target variable. It then summarizes the categorical variables created in the CLASS statement followed by telling us that the model converged. It also provides us with the same global tests, type 3 analysis of effects, parameter estimates, and odds ratios we have discussed in previous sections. The main difference in the parameter estimates is the multiple intercepts and slope coefficients for each variable. Since we have five levels of our target variable we have four logistic regression models and therefore four intercepts and four slope parameters for each variable.

Interpretation

Similar to binary and ordinal logistic regression models, we exponentiate the coefficients in our nominal logistic regression model to make them interpretable. However, the interpretation changes since these are not odds ratios anymore, but relative risk ratios. Let’s use the coefficient for the birds model and the size variable. It would be incorrect to say that the probability of eating birds is 2.076 times more likely for large alligators compared to small alligators. The correct interpretation would be that the predicted relative probability of eating birds rather than fish is 2.076 times more likely in large alligators compared to small alligators. Notice how our interpretation is in comparison to the reference level. Sometimes these are called conditional interpretations.

Let’s see how to get these from each of our softwares!

As with previous models, we just need to use the exp function to exponentiate the coefficients from our nominal logistic regression model to get the relative risk ratios.

Code

100*(exp(coef(glogit.model))-1)

             (Intercept) size> 2.3 meters lakeHancock lakeOklawaha lakeTrafford
Bird           -91.21513        107.52341    77.79659   -42.316600     244.6005
Invertebrate    18.42033        -73.71484   -83.14555   149.188934     217.6316
Other          -76.09086        -25.21626   115.26708     2.637196     374.8458
Reptile        -96.71612         74.57506   209.45502  1155.663805    2035.0755
             genderMale
Bird          -45.46914
Invertebrate  -37.05689
Other         -22.30894
Reptile       -46.61400

As with previous models, we just need to use the exp function to exponentiate the coefficients from our nominal logistic regression model to get the relative risk ratios.

Code

import numpy as np

100*(np.exp(glogit_model.params)-1)

                            0            1           2           3
size_> 2.3 meters   13.752398   -77.446008  -38.289298   -5.068399
lake_Hancock       256.831218   -22.969462  155.374032  -40.393110
lake_Oklawaha      435.740934  1892.512701  166.237876  571.276040
lake_Trafford      -20.904901   281.867313   99.419636   87.676014
gender_Male        425.787487   260.177574  106.600266   15.358571

By default, SAS produces the odds ratio table that shows the exponentiated coefficients. Even those these are not actually considered odds ratios in the traditional sense, the calculation in the odds ratio table is still correct.

Code

ods html select CLOddsPL;
proc logistic data=Logistic.Gator plot(only)=oddsratio(range=clip);
    freq count;
    class lake(param=ref ref='George') size(param=ref ref='<= 2.3 meters') gender(param=ref ref='Male');
    model food(ref='Fish') = lake size gender / link=glogit clodds=pl;
    title 'Model on Alligator Food Choice';
run;
quit;

Predictions & Diagnostics

Nominal logistic regression has a lot of similarities to binary logistic regression:

Multicollinearity still exists
Non-convergence problems still exists
AIC and BIC metrics can still be calculated
Generalized \(R^2\) remains the same

There are some inherent differences though between binary and nominal logistic regression. A lot of the diagnostics cannot be calculated for nominal logistic regression. ROC curves and residuals cannot typically be calculated because there is actually more than one logistic regression occurring.

Predicted probabilities are actually for each category though. Let’s see how we can get predicted probabilities from each of our softwares!

In R we have the same process of getting predicted probabilities in nominal logistic regression as we did for binary and ordinal logistic regression. The predict function will score any data set. The type = probs option is used to get individual predicted probabilities for each category.

Code

pred_probs <- predict(glogit.model, newdata = gator, type = "probs")
head(pred_probs)

       Fish       Bird Invertebrate     Other    Reptile
1 0.6006304 0.05115737   0.07545645 0.2401706 0.03258518
2 0.6006304 0.05115737   0.07545645 0.2401706 0.03258518
3 0.6006304 0.05115737   0.07545645 0.2401706 0.03258518
4 0.6236286 0.11022853   0.02059329 0.1864858 0.05906375
5 0.6236286 0.11022853   0.02059329 0.1864858 0.05906375
6 0.6236286 0.11022853   0.02059329 0.1864858 0.05906375

The five variables are the individual predicted probabilities for each of the categories of our target variable. Unlike binary logistic regression where we typically did not use a cut-off of 0.5 to determine which category we predict, in nominal logistic regression we typically pick the category with the highest probability as our predicted category.

Accuracy is typically not as high in logistic regressions with more than two categories because there are more chances for incorrect predictions.

In R we have the same process of getting predicted probabilities in nominal logistic regression as we did for binary and ordinal logistic regression. The predict function will score any data set.

Code

glogit_model.predict(X)

            0         1         2         3         4
0    0.035088  0.658310  0.097350  0.185125  0.024127
1    0.035088  0.658310  0.097350  0.185125  0.024127
2    0.035088  0.658310  0.097350  0.185125  0.024127
3    0.035088  0.658310  0.097350  0.185125  0.024127
4    0.035088  0.658310  0.097350  0.185125  0.024127
..        ...       ...       ...       ...       ...
214  0.254486  0.289484  0.057397  0.157045  0.241588
215  0.254486  0.289484  0.057397  0.157045  0.241588
216  0.254486  0.289484  0.057397  0.157045  0.241588
217  0.254486  0.289484  0.057397  0.157045  0.241588
218  0.254486  0.289484  0.057397  0.157045  0.241588

[219 rows x 5 columns]

Accuracy is typically not as high in logistic regressions with more than two categories because there are more chances for incorrect predictions.

In PROC LOGISTIC we have the same process of getting predicted probabilities in nominal logistic regression as we did for binary and ordinal logistic regression. The OUTPUT statement will score the training set, while the SCORE statement will score the validation data set. The predprobs = I option is used to get individual predicted probabilities for each category.

Code

ods html select none;
proc logistic data=Logistic.Gator plot(only)=oddsratio(range=clip);
    freq count;
    class lake(param=ref ref='George') size(param=ref ref='<= 2.3 meters') gender(param=ref ref='Male');
    model food(ref='Fish') = lake size gender / link=glogit clodds=pl;
    output out=pred predprobs=I;
    title 'Model on Alligator Food Choice';
run;
quit;
ods html select all;

proc print data=pred (obs = 5);
run;

proc freq data=pred;
    weight count;
    tables _FROM_*_INTO_;
    title 'Crosstabulation of Observed by Predicted Responses';
run;

The IP variables are the individual predicted probabilities for each of the categories of our target variable. Unlike binary logistic regression where we typically did not use a cut-off of 0.5 to determine which category we predict, in nominal logistic regression we typically pick the category with the highest probability as our predicted category. The FROM and INTO variables are the actual target variable category and the predicted category for the target respectively.

To see how the predicted categories compare with the actual target categories we can use PROC FREQ with a TABLES statement on the FROM and INTO variables.

Accuracy is typically not as high in logistic regressions with more than two categories because there are more chances for incorrect predictions.