Why MNLogit returns `classes_num - 1` params and how get them all?

user_7599215 0
Viewed 526

If I have few classes, i.e. 3. I am expecting to get 3 Generalized Linear Regression coefficients arrays, as in sklearn.linear_model.LogisticRegression but statsmodels.discrete.discrete_model.MNLogit provides classes_num - 1 coefficients (in this case -- 2).

Example:

import statsmodels.api as st
from sklearn.linear_model import LogisticRegression


iris = st.datasets.get_rdataset('iris','datasets')

y = iris.data.Species
x = iris.data.iloc[:, :-1]
mdl = st.MNLogit(y, x)
# mdl_fit = mdl.fit()
mdl_fit = mdl.fit(method='bfgs' , maxiter=1000)
print(mdl_fit.params.shape)  # (4, 2)

model = LogisticRegression(fit_intercept = False, C = 1e9)
mdl = model.fit(x, y)
print(model.coef_.shape)  # (3, 4)

How I am supposed to get regression coefficients for all 3 classes using MNLogit?

scikit-learn statsmodels python
2 Answers

Those coefficients are not computed to force identifiability of the model. In other words, not computing them ensures that the coefficients for the other classes are unique. If you had three sets of coefficients, there are an infinite number of models that give the same predictions but have different values for the coefficients. And this is bad if you want to know standard errors, p-values and so on.

The logit of the missing class is assumed to be zero. Demo:

mm = st.MNLogit(
    np.random.randint(1, 5, size=(100,)),
    np.random.normal(size=(100, 3))
)

res = mm.fit()
xt = np.random.normal(size=(2, 3))
res.predict(xt)

Results in:

array([[0.19918096, 0.34265719, 0.21307297, 0.24508888],
       [0.33974178, 0.21649687, 0.20971884, 0.23404251]])

Now these are the logits, plus the zeros for the first class

logits = np.hstack([np.zeros((xt.shape[0], 1)), xt.dot(res.params)])

array([[ 0.        ,  0.54251673,  0.06742093,  0.20740715],
       [ 0.        , -0.45060978, -0.4824181 , -0.37268309]])

And the predictions via a softmax:

np.exp(logits) / (np.sum(np.exp(logits), axis=1, keepdims=1))

array([[0.19918096, 0.34265719, 0.21307297, 0.24508888],
       [0.33974178, 0.21649687, 0.20971884, 0.23404251]])

which match the predictions from the model.

To reiterate: you cannot find those coefficients. Use a constant logit of zero for the first class. And you cannot find how much features are influential for the first class. That is actually an ill-posed question: the features cannot be influential for the reference class, because the reference class is never predicted directly. What the coefficients tell you is how much the log-odds for a given class, compared the reference class, change as a result of an unit increase for a particular feature.

user_5217760

The predicted probabilities have to add up to 1 over classes. So, we loose one free parameter in the full model in order to impose this constraint.

The predicted probability for the reference class is one minus the sum of probabilities of all other classes.

This is similar to the binary case, where we don't have separate parameters for success and failure because one probability is just one minus the other probability.

user_3337000
Related