6.5 Multinomial response variable
This model arises when the response variable can take a small number of discrete choices. Otherwise, if the number of possibilities is very high, it will be then more suitable to use the count data model (which will be discussed later on). It should be noted that the response values are assumed to be not ordered.
As we did previously, we will use again the logit function to model the conditional probability of each choice with respect to the regressors values. Assuming the response can take \(C\) choices, the probability of the individual \(i\) to take the choice \(c\) is given by:
\[\begin{equation} Pr\big(y_i=c\big|x_i)=\frac{exp\big(x_i\beta_c\big)}{\sum\limits^C_{c=1}exp\big(x_i\beta_c\big)} \tag{6.59} \end{equation}\]
This function called Multinomial logistic function, or softmax function.
Notice that \(\beta\) now is a \((K\times C)\) matrix, such that the values of each column are the coefficients of the corresponding choice.
Following the notation of (Murphy 2012), the log likelihood function is given by:
\[\begin{align} LL(\beta)&=log\prod\limits^n_{i=1}\prod\limits^C_{c=1}u_{ic}^{y_{ic}} \\ &=\sum\limits^{n}_{i=1}\sum\limits^{C}_{c=1}y_{ic}log(u_{ic}) \\ &=\sum\limits^{n}_{i=1}\Bigg[\bigg(\sum\limits^{C}_{c=1}y_{ic}x_i\beta_c\bigg)-log\bigg(\sum\limits^{C}_{c=1}exp(x_i\beta_c\bigg)\Bigg] \tag{6.60} \end{align}\]
Where \(u_{ic}=Pr(y_i=c|x_i, \beta)\), \(y_{ic}=I_{(y_i=c)}\).
Before computing the gradient we should first set the coefficients of the last choice \(\beta_C=0\), and stuck all the remaining columns of that matrix in one column vector denoted \(\mathrm w\). The gradient that will be used in optimization then will be:
\[\begin{equation} \mathrm g(\mathrm w)=\sum\limits^{n}_{i=1}(\mathrm u_i-\mathrm y_i)\otimes x_i \tag{6.61} \end{equation}\]
Where \(\mathrm u_i=\bigg(\begin{matrix}Pr(y_i=1|x_i,\beta)&..&Pr(y_i=C-1|x_i,\beta)\end{matrix}\bigg)^t\), \(\mathrm y_i=\bigg(\begin{matrix}I(y_i=1)&..&I(y_i=C-1)\end{matrix}\bigg)^t\)
For instance, if we have 3 choices and 2 features the gradient will be:
\[\begin{equation} \mathrm g(\mathrm w)=\sum\limits^{n}_{i=1}\begin{pmatrix}(u_{i1}-y_{i1})x_{i1}\\ (u_{i1}-y_{i1})x_{i2}\\ (u_{i2}-y_{i2})x_{i1}\\ (u_{i2}-y_{i2})x_{i2} \end{pmatrix} \tag{6.62} \end{equation}\]
Example 6.3 We will use Python package scikit-learn library to simulate 500 data points using multinomial distribution and assuming having 3 significant regressors.
In Python:
from sklearn.datasets import make_classification
import pandas as pd
import numpy as np
# simulate 500 data points
X, y = make_classification(n_samples=500, n_features=3, n_redundant=0, n_clusters_per_class=1, n_classes=3, random_state=123)
# store values in a dataframe
mn_data=pd.DataFrame(np.around(X,3), columns=[f'x{i}' for i in range(1,4)])
mn_data["y"]=np.around(y)
# display the dependent variable
mn_data.y.value_counts()
[out] y
[out] 1 168
[out] 2 166
[out] 0 166
[out] Name: count, dtype: int64We have three balanced classes. To fit a multinomial model to this data, we will again use the stasmodels package.
import statsmodels.formula.api as smf
# using only the three significant regressors
mn_model_py=smf.mnlogit(f"y~{'+'.join(mn_data.columns[:-1])}", mn_data).fit()[out] Optimization terminated successfully.
[out] Current function value: 0.334163
[out] Iterations 8[out] <class 'statsmodels.iolib.summary.Summary'>
[out] """
[out] MNLogit Regression Results
[out] ==============================================================================
[out] Dep. Variable: y No. Observations: 500
[out] Model: MNLogit Df Residuals: 492
[out] Method: MLE Df Model: 6
[out] Date: Thu, 17 Sept 2026 Pseudo R-squ.: 0.6958
[out] Time: 11:13:45 Log-Likelihood: -167.08
[out] converged: True LL-Null: -549.30
[out] Covariance Type: nonrobust LLR p-value: 7.435e-162
[out] ==============================================================================
[out] y=1 coef std err z P>|z| [0.025 0.975]
[out] ------------------------------------------------------------------------------
[out] Intercept -1.0891 0.294 -3.698 0.000 -1.666 -0.512
[out] x1 0.0713 0.222 0.321 0.748 -0.364 0.506
[out] x2 -3.2838 0.386 -8.518 0.000 -4.039 -2.528
[out] x3 -2.0509 0.346 -5.922 0.000 -2.730 -1.372
[out] ------------------------------------------------------------------------------
[out] y=2 coef std err z P>|z| [0.025 0.975]
[out] ------------------------------------------------------------------------------
[out] Intercept -2.0399 0.361 -5.643 0.000 -2.748 -1.331
[out] x1 0.1051 0.203 0.518 0.604 -0.292 0.503
[out] x2 1.0214 0.288 3.545 0.000 0.457 1.586
[out] x3 -3.6829 0.367 -10.037 0.000 -4.402 -2.964
[out] ==============================================================================
[out] """Every coefficient of the above output is the relative log odds of the associated choice over the base choice (where it is here \(y=0\)) if the associated regressor increases by one (if it is a continuous regressor). To understand the interpretation, it is better to convert the log odds into odds (probabilities) by exponentiating those coefficients. For instance, let us take the coefficient \(1.0214\) of \(x_2\) for the choice \(y=2\), and after exponentiation, we get:
[out] 2.7772058125599464This means that if \(x_2\) increases by one, the probability of being with the choice \(y=2\) is \(2.7\) times of that of the base choice \(y=0\). However, if the regressor \(x_2\) is a specific category from categorical regressor, this category then will be compared to the base category of that regressor.
In R:
In R, There exist bunch of packages that handle multinomial models, but we will use the nnet. Since converting the response variable from numeric to categorical type in R is not automatic, we should do it by hand:
suppressPackageStartupMessages(library(nnet))
# Move the data to R space
mn_dataR <-py$mn_data
# convert the response variable to categorical type
mn_dataR$y <- factor(mn_dataR$y)
#fit the model
mn_dataR_model <- multinom(y~., mn_dataR)[out] # weights: 15 (8 variable)
[out] initial value 549.306144
[out] iter 10 value 206.751573
[out] iter 20 value 167.083897
[out] iter 30 value 167.081527
[out] iter 30 value 167.081527
[out] iter 30 value 167.081527
[out] final value 167.081527
[out] converged[out] Call:
[out] multinom(formula = y ~ ., data = mn_dataR)
[out]
[out] Coefficients:
[out] (Intercept) x1 x2 x3
[out] 1 -1.089116 0.07129945 -3.283815 -2.050855
[out] 2 -2.039883 0.10510994 1.021445 -3.682891
[out]
[out] Std. Errors:
[out] (Intercept) x1 x2 x3
[out] 1 0.2944956 0.2218684 0.3855269 0.3463240
[out] 2 0.3614999 0.2028569 0.2880969 0.3669349
[out]
[out] Residual Deviance: 334.1631
[out] AIC: 350.1631If we want to change the base reference to \(1\) for instance instead of \(0\), we can use the following script:
As we see, we get the same results as in python.