CS2.2 — Model Fitting & Interpretation

0.1 What this page is about

Multinomial logistic regression extends binary logistic regression to outcomes with three or more unordered categories — for example, treatment choice (surgery / radiation / medication), or level of psychiatric impairment (well / mild / moderate / impaired).

The core idea is straightforward: instead of one equation comparing “event vs no event”, the model fits one equation per non-reference category, each modelling the log-odds of being in that category compared to the reference. For K categories you get K−1 equations, and each equation has its own set of coefficients.

This flexibility is what makes the model powerful: the effect of a predictor such as socioeconomic status on the odds of being mildly impaired vs well does not have to equal its effect on the odds of being severely impaired vs well. The coefficients can — and often do — differ across equations.

On this page you will see how the model is estimated (maximum likelihood), how to read the coefficient matrix, and how to compute and interpret odds ratios per category comparison.

1 Concept: Estimation by Maximum Likelihood

1.1 The multinomial log-likelihood

For $n$ independent observations, the log-likelihood of the multinomial model is:

\[\ell(\boldsymbol{\beta}) = \sum_{i=1}^{n} \sum_{k=1}^{K} y_{ik} \log \pi_{ik}(\boldsymbol{\beta})\]

where $y_{ik} = \mathbf{1}[Y_i = k]$ is the indicator that observation $i$ is in category $k$, and $\pi_{ik}$ is the softmax probability.

There is no closed-form solution — $\ell$ is maximised iteratively using Newton-Raphson or quasi-Newton algorithms. nnet::multinom uses BFGS.

Because the multinomial log-likelihood is not differentiable

Because the design matrix is always rank-deficient in multinomial models

Because the log-likelihood is nonlinear in the parameters and must be maximised iteratively

Because the outcome categories are mutually exclusive

The multinomial log-likelihood involves log(π_ik(β)) where π_ik is a nonlinear (softmax) function of β. Unlike ordinary linear regression, there is no analytic solution. Algorithms like Newton-Raphson or BFGS (used by nnet::multinom) iteratively update β until convergence.

1.2 The information matrix and standard errors

At convergence, the MLE $\hat{\boldsymbol{\beta}}$ has approximate covariance:

\[\widehat{\text{Var}}(\hat{\boldsymbol{\beta}}) = \left[-\frac{\partial^2 \ell}{\partial \boldsymbol{\beta} \partial \boldsymbol{\beta}^\top}\bigg|_{\hat{\boldsymbol{\beta}}}\right]^{-1} = \mathcal{I}(\hat{\boldsymbol{\beta}})^{-1}\]

Standard errors are square roots of diagonal elements. This is the observed information matrix.

The covariance matrix of the predictor variables, measuring how they relate to each other

The matrix of first derivatives of the log-likelihood, known as the score vector

The negative matrix of second derivatives of the log-likelihood, evaluated at the MLE

The matrix of squared residuals from the fitted model, used to estimate unexplained variance

The observed information matrix I(β̂) = −∂²ℓ/∂β∂β'|_{β̂} measures the curvature of the log-likelihood at the MLE. Its inverse approximates the covariance matrix of β̂. Square roots of the diagonal give standard errors. Larger curvature (sharper peak) = more information = smaller SEs.

1.3 Interpreting coefficients

For predictor $x_j$ and category $k$ (vs reference $r$), the coefficient $\hat{\beta}_{kj}$ is:

The estimated log-odds ratio: $\log \frac{\pi_k / \pi_r}{\pi_k^* / \pi_r^*}$ when $x_j$ increases by 1 (all else constant)
Exponentiated: $\exp(\hat{\beta}_{kj})$ is the odds ratio for category $k$ vs $r$ per unit increase in $x_j$

What the OR is NOT

The OR in a multinomial model is not marginal — it conditions on all other predictors. It also does not directly give the change in probability (that depends on the current probability level via the softmax).

The predicted probability of being in category k

The probability ratio between category k and the reference for predictor j

The odds ratio for category k vs reference per one-unit increase in predictor j

The relative risk of category k vs reference

β_kj is a log-odds ratio, so exp(β_kj) is the odds ratio for category k vs the reference, holding other predictors fixed. An OR > 1 means predictor j is associated with higher odds of being in category k rather than the reference. Note: this is not a relative risk or a probability ratio.

2 Worked Example

summary(fit)

Call:
multinom(formula = impair ~ ses + life, data = dat, trace = FALSE)

Coefficients:
         (Intercept)        ses      life
mild      -0.7728136 -0.5573249 0.1453964
moderate   0.2431942 -2.0353372 0.5123313
impaired   0.8872316 -2.7933927 0.6554945

Std. Errors:
         (Intercept)       ses      life
mild       0.6612680 0.7178526 0.1339639
moderate   0.4768376 0.5294046 0.1026607
impaired   0.4430075 0.5112686 0.1005892

Residual Deviance: 532.6075 
AIC: 550.6075

# Full coefficient table with ORs and 95% Wald CIs
betas <- coef(fit)
ses_m <- summary(fit)$standard.errors

cat("=== Coefficients (log-OR scale) ===\n")

=== Coefficients (log-OR scale) ===

round(betas, 4)

         (Intercept)     ses   life
mild         -0.7728 -0.5573 0.1454
moderate      0.2432 -2.0353 0.5123
impaired      0.8872 -2.7934 0.6555

cat("\n=== Odds Ratios ===\n")


=== Odds Ratios ===

round(exp(betas), 4)

         (Intercept)    ses   life
mild          0.4617 0.5727 1.1565
moderate      1.2753 0.1306 1.6692
impaired      2.4284 0.0612 1.9261

cat("\n=== 95% Wald CIs (OR scale) ===\n")


=== 95% Wald CIs (OR scale) ===

ci_lo <- exp(betas - 1.96*ses_m)
ci_hi <- exp(betas + 1.96*ses_m)
for(cat_name in rownames(betas)){
  cat(cat_name, ":\n")
  cat("  ses:  OR =", round(exp(betas[cat_name,"ses"]),3),
      " [", round(ci_lo[cat_name,"ses"],3), ",",
      round(ci_hi[cat_name,"ses"],3), "]\n")
  cat("  life: OR =", round(exp(betas[cat_name,"life"]),3),
      " [", round(ci_lo[cat_name,"life"],3), ",",
      round(ci_hi[cat_name,"life"],3), "]\n")
}

mild :
  ses:  OR = 0.573  [ 0.14 , 2.339 ]
  life: OR = 1.156  [ 0.889 , 1.504 ]
moderate :
  ses:  OR = 0.131  [ 0.046 , 0.369 ]
  life: OR = 1.669  [ 1.365 , 2.041 ]
impaired :
  ses:  OR = 0.061  [ 0.022 , 0.167 ]
  life: OR = 1.926  [ 1.581 , 2.346 ]

Coefficients:
          (Intercept)    ses    life
mild           -0.82  -0.52   0.21
moderate       -0.14  -1.19   0.44
impaired        0.48  -1.73   0.59

The model has been mis-specified — ses coefficients should be equal in every row

Higher SES is increasingly protective against more severe impairment: the log-odds advantage of higher SES grows with severity

Higher SES reduces the probability of all impairment categories by the same amount

The ses coefficient in the impaired row (−1.73) should be interpreted relative to the mild row (−0.52), not relative to the reference category

Each row is a separate logistic equation comparing that category to the reference (well). The ses coefficient in each row is the log-OR for that specific comparison. Here: OR for mild vs well = exp(−0.52) ≈ 0.60; OR for impaired vs well = exp(−1.73) ≈ 0.18. The pattern — increasingly negative ses as severity increases — shows that higher SES is most strongly protective against the most severe outcomes. This variation across rows is exactly what the multinomial model allows and what a proportional odds model would constrain to be equal across all comparisons.

# Deviance and AIC
cat("Deviance (= -2*logLik):", round(deviance(fit), 4), "\n")

Deviance (= -2*logLik): 532.6075

cat("AIC:                   ", round(AIC(fit), 4), "\n")

AIC:                    550.6075

cat("Log-likelihood:        ", round(logLik(fit), 4), "\n")

Log-likelihood:         -266.3037

3 Guided Practice

R is loading in the background — the exercises below will be interactive once the status bar disappears. You can read the concept sections above immediately.

Hints

multinom(party ~ age + educ, data=vote_dat, trace=FALSE)
OR for educ in Liberal eq: exp(coef(fit_vote)["Liberal","educ"])
predict(fit_vote, newdata=data.frame(age=40,educ=3), type="probs")

4 Retrieval Challenge

Model answers

Q1: OR = exp(1.8) ≈ 6.05 — treated patients have 6× higher odds of improving vs remaining stable, compared to untreated.

Q2: All three probabilities computed above should sum to 1.

5 Why-Questions

Q: “The coefficient for ‘ses’ differs across the three equations. Doesn’t this complicate interpretation?”

Yes — and this is the price of flexibility. The multinomial model allows the effect of SES to differ depending on which category you’re comparing to the reference.

If you expect the effect to be the same across all comparisons, consider testing this as a constraint (CS2.5). If the constrained model fits as well as the unconstrained, you gain parsimony; if not, the different slopes are a substantive finding.

In ordinal data (CS3), the proportional odds model imposes equal slopes by construction — a special case of the constrained multinomial.

Q: “Why does multinom use BFGS and not Newton-Raphson? Does it matter?”

Both are gradient-based optimisers. Newton-Raphson uses exact second derivatives (Hessian); BFGS approximates the Hessian. For multinomial regression the log-likelihood is concave (unique maximum), so both converge to the same point. BFGS is preferred in practice because computing the full Hessian for high-dimensional $\boldsymbol{\beta}$ is expensive.

The SEs from summary(fit)$standard.errors are based on the inverse Hessian at convergence — regardless of which algorithm was used.

0.1 What this page is about

1 Concept: Estimation by Maximum Likelihood

1.1 The multinomial log-likelihood

1.2 The information matrix and standard errors

1.3 Interpreting coefficients

2 Worked Example

3 Guided Practice

4 Retrieval Challenge

5 Why-Questions

6 Spaced Review