Multinomial logit models allow us to model membership in a group based on known variables. For example, operating system preference of a university’s students could be classified as “Windows”, “Mac”, or “Linux”. Perhaps we would like to better understand why students choose one OS versus another. We might want to build a statistical model that allows us to predict the probability of selecting an OS based on information such as sex, major, financial aid and so on. Multinomial logit modeling allows us to propose and fit such models.

It’s important to note that multinomial logit models are best suited for *nominal categories*. These are categories with no natural ordering. The example presented above has no natural ordering. Despite the opinions of some computer users, Linux is not *measurably* “higher” or “lower” than Windows, nor is Mac “higher” or “lower” than Linux. They are distinct groups. Examples of ordered categories include sizes, pain scales, and age groups. These are categories in which one level might be considered greater than or less than an adjacent category. A common approach to modeling ordered categories is the ordered logit model, also known as the proportional odds model. For more information, see our related post, Fitting and Interpreting a Proportional Odds Model

To illustrate multinomial logit models, we’ll use the R statistical programming environment and data on alligator food choice from the text, *An Introduction to Categorical Data Analysis* by Alan Agresti (1996, p. 207). The data contain information on 59 alligators sampled from a lake in Florida. It has the length of the alligator in meters and the primary food type found in the alligator’s stomach. The food type was classified into three categories: Fish, Invertebrate and Other. Of interest is whether or not the length of an alligator is associated with the primary food type. Does knowing the length of an alligator give us some indication about its primary food type? If so, how does length affect the choice of food type?

Let’s read in the data and examine it. (The code below should work for you if have an active internet connection.) Notice we make the food column a factor and set the reference (or baseline) level to “Other” by placing it first in the levels and labels arguments to `factor`

. We’ll explain why we do this in a bit.

```
gators <- read.csv('https://static.lib.virginia.edu/statlab/materials/data/table_8-1.csv')
gators$food <- factor(gators$food,
levels = c("O", "F", "I"),
labels = c("Other", "Fish", "Invertebrates"))
```

Looking at the counts and proportions, we see there is certainly variability between the food groups.

```
# counts
table(gators$food)
```

```
##
## Other Fish Invertebrates
## 8 31 20
```

```
# proportions
prop.table(table(gators$food))
```

```
##
## Other Fish Invertebrates
## 0.1355932 0.5254237 0.3389831
```

We can also explore the data visually using the ggplot2 package. First we make a boxplot of length by food choice.

```
library(ggplot2)
ggplot(gators) +
aes(x = food, y = length) +
geom_boxplot()
```

Next we make a one-dimensional scatter plot of length by food choice. Notice how we use `position_jitter`

to randomly scatter the points in a horizontal direction so they’re not overplotted.

```
ggplot(gators) +
aes(y = length, x = food) +
geom_point(position = position_jitter(width = 0.1, height = 0))
```

In both plots it seems smaller alligators are more likely to eat primarily Invertebrates while bigger alligators are more likely to eat Fish. It’s hard to infer much about Other since we only have 8 observations which appear uniformly scattered over the range of data.

So we can see there appears to be *some* relationship between length and food choice, but how can we quantify it? How does primary food choice change as length changes? This is where multinomial logit modeling comes in.

In R, we can use the `nnet`

package that comes installed with base R. It has the `multinom`

function which fits multinomial logit models via neural networks. Don’t worry, you don’t need to know anything about neural networks to use the function. In fact it works much like the workhorse modeling functions, `lm`

and `glm`

.

To fit our model we specify food be modeled as a function of length using `food ~ length`

. We specify `data = gators`

so the function knows where to find the food and length variables. Finally we specify `Hess = TRUE`

to save something called the Hessian matrix which helps with the calculation of standard errors, and set `trace = FALSE`

to suppress the printing of convergence progress.

```
library(nnet)
m <- multinom(food ~ length, data = gators,
Hess = TRUE, trace = FALSE)
```

After fitting the model we can use the `summary`

function to see coefficient estimates along with standard errors and test statistics. To see the Wald test statistics, we need to also set `Wald.ratios = TRUE`

. We explain what those are below.

`summary(m, Wald.ratios = TRUE)`

```
## Call:
## multinom(formula = food ~ length, data = gators, Hess = TRUE,
## trace = FALSE)
##
## Coefficients:
## (Intercept) length
## Fish 1.617952 -0.1101836
## Invertebrates 5.697543 -2.4654695
##
## Std. Errors:
## (Intercept) length
## Fish 1.307291 0.5170838
## Invertebrates 1.793820 0.8996485
##
## Value/SE (Wald statistics):
## (Intercept) length
## Fish 1.237637 -0.2130865
## Invertebrates 3.176206 -2.7404808
##
## Residual Deviance: 98.34124
## AIC: 106.3412
```

The output is a little different from what you may be accustomed to seeing. In the Coefficient section we’re basically seeing two models, one modeling the effect of length on the log odds that an alligator prefers Fish to Other, and the other that an alligator prefers Invertebrates to Other. This is why multinomial logit models are sometimes called *baseline logit models*. They model each category relative to some baseline level. In this case the baseline level is “Other”, which we specified when we set the food variable as a factor above. In general if you have *J* categories, you will have *J-1* baseline logit models.

The next section provides standard errors of the coefficient estimates. These quantify uncertainty about the estimates. For example, the length coefficient for Invertebrates is -2.46 with a standard error of 0.89. The smaller the standard error in relation to the coefficient the more precise the estimate. The section following the standard errors, Wald statistics, presents the ratio of the coefficients to the standard errors. The farther away the ratio is from 0, the more precise the estimate and the more likely it’s either positive or negative. This is no different from the output of other regression models. It’s just presented differently due to the fact we’re summarizing two models.

The `summary`

method for `multinom`

does not report p-values. However that’s no hindrance to making inference about the “significance” of the coefficients. Wald statistics larger than 2 are evidence of effects different from 0. If you must report p-values, you’ll have to calculate them manually. Here’s one way using the standard normal approximation:

```
s <- summary(m, Wald.ratios = TRUE)
p <- 2 * pnorm(abs(s$Wald.ratios), lower.tail = FALSE)
p
```

```
## (Intercept) length
## Fish 0.215850665 0.831259475
## Invertebrates 0.001492149 0.006134937
```

Despite their ubiquity, p-values don’t tell us anything about the size or direction of coefficients. They simply tell us the probability of our sample producing a test statistic (in this case, a Wald ratio) that far away from 0 (or farther) if the statistic really was 0.

Interpreting coefficients in logit models is difficult. The only thing that’s easy is determining whether or not the coefficient is positive or negative. For example, the length coefficient for Invertebrates is negative (-2.46) and has a large Wald statistic (greater than 2). It appears that the longer an alligator is the less likely it is to eat Invertebrates instead of Other. But that’s nothing we couldn’t have figured out from just looking at our exploratory plots.

The coefficient of -2.46 says that each additional meter of length reduces the log odds of choosing Invertebrates over Other by -2.46. Unfortunately humans don’t think in terms of log odds. It can be shown that exponentiating the coefficient produces an odds ratio, which is a little easier to understand. Here we use the `coef`

function to extract the coefficients from the model, use indexing brackets to specify we want the “length” coefficient for the “Invertebrate” vs Other model, and then use the `exp`

function to exponentiate.

`exp(coef(m)['Invertebrates','length'])`

`## [1] 0.08496894`

The odds ratio is about 0.08. One way to interpret this is to imagine two groups of alligators, one that has gators 2 meters long (*x*) and another that has gators 3 meters long (*x + 1*). For alligators that are 3 meters long, the estimated odds that they’re food choice is Invertebrate rather than Other is (1 – 0.08) 92% less than the odds for alligators that are 2 meters in length. This interpretation is the same no matter the value of x. Obviously you should choose a value in the range of your data.

Both the coefficient and odds ratio are just point estimates. Given a different sample, we would likely get different values. To get a sense of the range of values we might get with repeated samples and models, we can calculate 95% confidence intervals using the `confint`

function.

`confint(m)`

```
## , , Fish
##
## 2.5 % 97.5 %
## (Intercept) -0.9442917 4.1801963
## length -1.1236492 0.9032821
##
## , , Invertebrates
##
## 2.5 % 97.5 %
## (Intercept) 2.181720 9.2133667
## length -4.228748 -0.7021908
```

The confidence interval for the length coefficient in the Invertebrates model is about (-4.23, -0.70). We can exponentiate that to get a confidence interval for the odds ratio. Below we use index brackets to pull out the confidence interval from the 3-dimensional array that `confint(m)`

returns. The first dimension is rows, so we specify “length”. The second dimension is columns, but we leave it blank to indicate we want all columns. The third dimension is the model, and we specify we want the “Invertebrate” model. Finally we exponentiate.

`exp(confint(m)['length',,'Invertebrates'])`

```
## 2.5 % 97.5 %
## 0.01457062 0.49549858
```

This suggests that if we repeatedly sampled and modeled alligators from this lake, 95% of the time we would get odds ratios ranging in value from about 0.01 to 0.50.

What if we wanted to estimate the odds ratio that a gator prefers Fish to Invertebrate? We could refit the model with Invertebrate as the baseline level and use the same process just described. Alternatively we can subtract the Invertebrates model coefficients from the Fish coefficients as follows:

`exp((coef(m)[1,] - coef(m)[2,])["length"])`

```
## length
## 10.54114
```

For alligators of length *x + 1* meters, the estimated odds that the food choice is Fish rather than Invertebrates is about 10 times the estimated odds for alligators of length *x*.

Expressing these relationships in odds ratios is a little better than using log odds, but still not ideal for many people. Probability is something that is more natural to work with. Fortunately we can use multinomial logit models to estimate probabilities. For example, what’s the probability an alligator that is 2.5 meters long prefers Fish? The `predict`

function allows us to easily estimate this. We just need to use the `newdata`

argument to provide the model input as a data frame and set the `type='probs'`

argument to ensure we get probabilities.

`predict(m, newdata = data.frame(length = 2.5), type = 'probs')`

```
## Other Fish Invertebrates
## 0.1832844 0.7017184 0.1149973
```

Our model estimates that a gator 2.5 meters long has 70% probability of preferring Fish. Notice the `predict`

function also returns estimated probabilities for the two other categories as well. Like the coefficients and odds ratios, the probabilities are just point estimates. Once again confidence intervals help us understand the uncertainty in our estimated probabilities.

A convenient function for obtaining confidence intervals for probabilities is the `ggeffect`

function from the `ggeffects`

package. Here we use it to get estimated probabilities and confidence intervals for lengths ranging from 1 to 4 in steps of 0.5. We do that by specifying `terms = "length[1:4,by=0.5]"`

.

```
library(ggeffects)
ggeffect(m, terms = "length[1:4,by=0.5]")
```

```
##
## # Predicted probabilities of food
## # x = length
##
## # Response Level = Other
##
## x | Predicted | SE | 95% CI
## --------------------------------------
## 1.00 | 0.03 | 0.91 | [0.01, 0.17]
## 1.50 | 0.08 | 0.59 | [0.03, 0.22]
## 2.00 | 0.14 | 0.44 | [0.06, 0.28]
## 2.50 | 0.18 | 0.40 | [0.09, 0.33]
## 3.00 | 0.21 | 0.50 | [0.09, 0.41]
## 4.00 | 0.23 | 0.91 | [0.05, 0.65]
##
## # Response Level = Fish
##
## x | Predicted | SE | 95% CI
## --------------------------------------
## 1.00 | 0.15 | 0.69 | [0.04, 0.40]
## 1.50 | 0.34 | 0.38 | [0.19, 0.52]
## 2.00 | 0.56 | 0.32 | [0.41, 0.71]
## 2.50 | 0.70 | 0.38 | [0.53, 0.83]
## 3.00 | 0.75 | 0.47 | [0.55, 0.88]
## 4.00 | 0.76 | 0.90 | [0.35, 0.95]
##
## # Response Level = Invertebrates
##
## x | Predicted | SE | 95% CI
## --------------------------------------
## 1.00 | 0.82 | 0.69 | [0.54, 0.95]
## 1.50 | 0.58 | 0.38 | [0.40, 0.75]
## 2.00 | 0.30 | 0.36 | [0.17, 0.47]
## 2.50 | 0.11 | 0.66 | [0.03, 0.32]
## 3.00 | 0.04 | 1.03 | [0.01, 0.23]
## 4.00 | 0.00 | 1.80 | [0.00, 0.11]
```

`## Standard errors are on link-scale (untransformed).`

Looking in the table, under the section “Response Level = Fish”, we see the 95% confidence interval on the estimated probability of 0.7 is (0.53, 0.83). We can also save this result and use it to visualize all the estimated probabilities and their associated confidence bands as an *effect plot*. First we save the result as `eff`

. Next we use ggplot and specify we want length (x) on the x-axis, predicted probabilities (predicted) on the y-axis, and all colors and fill aesthetics mapped to food group (response.level). Then we draw a line for the predicted values and a ribbon for the lower and upper confidence limits (conf.low and conf.high).

```
eff <- ggeffect(m, terms = "length[1:4,by=0.5]")
ggplot(eff) +
aes(x = x, y = predicted, fill = response.level, color = response.level) +
geom_line() +
geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = 1/3) +
labs(x = 'Length of Alligator', y = 'Predicted Probability') +
ylim(c(0, 1))
```

Now we clearly see that bigger alligators prefer Fish, especially when they are larger than 2 meters. Likewise smaller alligators seem to prefer Invertebrates, probably because they may not be big enough to easily consume fish. We can also see there is a great deal of uncertainty in predicted probabilities of food category beyond 3 meters of length. As expected, the predicted probability of the Other food group as a function of length is inconclusive since we only had 8 gators in that group who widely varied in their lengths.

Finally we sometimes want to test whether a predictor variable is adding any value to a model. In this case we only have one predictor, length, and we’ve shown through Wald statistics and confidence intervals that it seems to be useful in understanding an alligator’s primary food choice. However if we wanted to formally test whether a model with length is better than the null model with just an intercept, we can carry out a Likelihood Ratio Test using the `Anova`

function in the `car`

package.

```
library(car)
Anova(m)
```

```
## Analysis of Deviance Table (Type II tests)
##
## Response: food
## LR Chisq Df Pr(>Chisq)
## length 16.801 2 0.0002248 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

The result is a test of whether or not the length coefficient in both baseline logit models is plausibly 0. The null is that a model with just an intercept is equally as good as a model with an intercept and length. (An intercept-only model is the same as simply using the observed proportions in each food category as a basis for making predictions about an alligator’s food choice.) The small p-value (p = 0.0002248) provides evidence against the null and leads us to keep length in our model.

References

- Agresti, A. (1996).
*An Introduction to Categorical Data Analysis*. Wiley. - Bilder, C. and Loughin, T. (2015).
*Analysis of Categorical Data with R*. CRC Press. - Friendly, M. and Meyer, D. (2016).
*Discrete Data Analysis with R*. CRC Press.

For questions or clarifications regarding this article, contact the UVA Library StatLab: statlab@virginia.edu

View the entire collection of UVA Library StatLab articles.

Clay FordStatistical Research Consultant

University of Virginia Library

September 1, 2020