Fitting a model to data always succeeds in some technical sense — the algorithm converges and produces coefficients. But how well does the model actually fit? And if you have two competing models, which one is better?
Goodness of fit refers to how closely the model’s predicted probabilities match the observed outcomes. The primary measure is deviance (−2 × log-likelihood): a perfect model has zero deviance; a model that ignores all predictors (the null model) has a high deviance. The difference in deviance between two nested models is the likelihood ratio test (LRT) statistic.
For comparing non-nested models, or penalising unnecessary complexity, the Akaike Information Criterion (AIC) is used. AIC rewards fit but penalises the number of parameters — a model with lower AIC is preferred.
On this page you will learn how to read deviance output, compare models with LRT and AIC, and think critically about what “good fit” means for a multinomial model.
NoteAbbreviations used on this page
Abbreviation
Full term
AIC
Akaike Information Criterion
BIC
Bayesian Information Criterion
LRT
Likelihood ratio test
df
Degrees of freedom
OR
Odds ratio
MLE
Maximum likelihood estimator
1 Concept: Assessing Model Fit
1.1 Deviance and the null model
The deviance of a model is \(D = -2\ell(\hat{\boldsymbol{\beta}})\). Smaller deviance = better fit.
The null deviance is the deviance of the intercept-only model. The difference:
is the likelihood ratio statistic for the overall model against the null. Under \(H_0\) (no predictors matter): \(G^2 \overset{a}{\sim} \chi^2_{df}\) where \(df\) = number of non-intercept parameters.
The model is overfit — a large G² is a sign of overfitting and the model should be simplified
The predictors collectively explain a significant amount of the variation in the outcome
The model has more parameters than can be reliably estimated given the available sample size
The Pearson chi-square assumption is violated because the expected counts are too large
G² = 2[ℓ(model) − ℓ(null)] is the LRT for overall model fit. Under H₀ (no predictors matter), G² ~ χ²_{df} where df = number of non-intercept parameters. A large G² (small p-value) means at least one predictor improves fit significantly over guessing the most common category.
1.2 Pearson chi-square statistic
For grouped data (cells defined by covariate patterns with enough observations each):
where \(O_{ik}\) = observed count, \(E_{ik} = n_i \hat{\pi}_{ik}\) = expected count. Under a well-fitting model: \(X^2 \overset{a}{\sim} \chi^2_{(c-p)(K-1)}\) where \(c\) = number of covariate patterns.
WarningSmall expected counts
The \(\chi^2\) approximation requires \(E_{ik} \geq 5\) in each cell. With continuous predictors or small \(n\), this is rarely met for individual data. In that case, use the LRT instead, or the Hosmer-Lemeshow approach (grouping observations into \(g\) bins).
Continuous predictors violate the multinomial distribution, which requires only discrete outcomes
Each observation forms its own covariate pattern, making expected cell counts much less than 5
The χ² approximation breaks down unless all groups have balanced numbers of observations
Continuous predictors produce non-identifiable deviance because the likelihood becomes unbounded
The Pearson χ² test requires E_ik ≥ 5 in each cell. With continuous predictors, each observation has a unique covariate pattern, giving c = n covariate patterns each with expected count near 1. The χ² approximation breaks down completely. Use LRT instead, or group observations (Hosmer-Lemeshow approach).
1.3 Hosmer-Lemeshow test for individual data
When predictors are continuous, each observation has a unique covariate pattern and the Pearson test breaks down. The Hosmer-Lemeshow approach groups observations into \(g\) bins (typically 10) by their predicted probability and computes a chi-square statistic comparing observed vs expected counts within each bin.
where \(O_{bk}\) = observed count of category \(k\) in bin \(b\), and \(E_{bk}\) = sum of predicted probabilities for category \(k\) across observations in bin \(b\).
Warning
The Hosmer-Lemeshow test has low power in small samples. A non-significant result does not guarantee good fit — it may simply reflect low sample size.
Pearson X² test — it works for any multinomial data
Deviance test comparing to the saturated model
Hosmer-Lemeshow test — group observations into bins, then compare observed vs expected counts
No GOF test is possible — multinomial models cannot be evaluated for fit
Pearson X² and deviance-based GOF tests require grouped data with large expected cell counts (≥5 per cell). With a continuous covariate, each observation forms its own covariate pattern, giving n=500 cells each with expected count near 1 — the chi-squared approximation breaks down. The Hosmer-Lemeshow approach pools observations into bins to create adequate expected counts.
1.4 Grouped vs individual data format
The same multinomial model can be fitted to two different data formats:
Format
Structure
R syntax
Individual
One row per observation; outcome is a factor
multinom(y ~ x1 + x2, data=dat)
Grouped
One row per covariate pattern; outcome is a count matrix
multinom(cbind(y1,y2,y3,y4) ~ x1 + x2, data=dat)
ML estimates are identical either way. The difference matters for goodness-of-fit tests: Pearson \(X^2\) and deviance are valid GOF statistics only for grouped data with sufficient expected counts (≥5 per cell). For individual data, use Hosmer-Lemeshow or LRT instead.
We can test whether two categories share the same predictor effects by fitting a constrained model where their coefficients are equal. This is an LRT with \(df\) = number of equality constraints.
Whether one category can be collapsed into the reference category
Whether the predictor effect is homogeneous across those two categories
Whether the model is identifiable with the given sample size
Whether the proportional odds assumption holds for ordinal outcomes
Constraining β_k = β_j for two categories k and j tests H₀: the predictor has the same log-OR effect on both. If the LRT is non-significant, the simpler constrained model is supported — the two categories respond identically to the predictor. The df = number of equality constraints imposed.
2 Worked Example
# Overall model fitcat("Null deviance: ", round(deviance(fit0), 4), "\n")
The model fits the data well — p > 0.05 confirms good fit because non-significant means adequacy of the fitted multinomial model
The predictors are clearly useless — the model is no better than the null since p = 0.16 far exceeds the α = 0.05 threshold
We fail to reject H₀ (p = 0.16), but with n = 40 the test has low power — a non-significant result does not rule out meaningful effects
The model should be rejected — a reduction of 9.24 in deviance for 6 additional parameters is too small to justify the added complexity
G² = 9.24 on df = 6 gives p = 0.16: we cannot reject H₀ that all non-intercept parameters are zero at α = 0.05. However, with only n = 40 observations the test has very limited power. The deviance decreased by 9.24 points (from 93.45 to 84.21), which is substantively meaningful even if not quite statistically significant. This is a classic case where the study is underpowered to detect real effects.
# Information criteriacat("AIC:", round(AIC(fit), 4), "\n")
AIC: 93.6416
cat("BIC:", round(BIC(fit), 4), "\n")
BIC: 108.8415
# BIC penalises complexity more strongly: n*log(n) vs 2*k
# Test whether mild and moderate share the same coefficients# (i.e., can we merge them?)# Fit a model where mild and moderate have the same ses and life slopes:# This is done by creating a new merged outcomedat2 <- datlevels(dat2$impair)[levels(dat2$impair) %in%c("mild","moderate")] <-"mild_mod"fit_constrained <-multinom(impair ~ ses + life, data=dat2, trace=FALSE)# The constrained model has fewer parameterscat("Full model df: ", length(coef(fit)), "\n")
Full model df: 9
cat("Constrained model df: ", length(coef(fit_constrained)), "\n")
If \(p > 0.05\): the data are consistent with mild and moderate having the same covariate effects — they may be combinable, simplifying the model. If \(p < 0.05\): the categories respond differently to the predictors and should remain separate.
# Hosmer-Lemeshow test for multinomial individual data# (manual implementation — no external package needed)pred_probs <-predict(fit, type ="probs") # n × K matrix of predicted probsK_cats <-ncol(pred_probs)n_obs <-nrow(pred_probs)g_bins <-5# number of bins (small here because n=40 is tiny)# Group by overall predicted risk: use the fitted probability of the most likely categorymax_prob <-apply(pred_probs, 1, max)bin_breaks <-quantile(max_prob, probs =seq(0, 1, length.out = g_bins +1))hl_groups <-cut(max_prob, breaks = bin_breaks, include.lowest =TRUE, labels =FALSE)# For each bin and each category: observed vs expectedhl_stat <-0cat("Bin-level comparison (observed vs expected per category):\n")
Bin-level comparison (observed vs expected per category):
cat("\nNote: n =", n_obs, "— this sample is very small; interpret cautiously.\n")
Note: n = 40 — this sample is very small; interpret cautiously.
NoteWhy the HL test has low power here
With \(n = 40\) and \(K = 4\) categories, there are only ~8 observations per bin. Expected cell counts are far below 5 in many cells, so the chi-squared approximation is unreliable. The Hosmer-Lemeshow test is most trustworthy with \(n \geq 200\) and at least 5 expected observations per cell.
# Grouped data format: one row per covariate pattern, outcome = count matrix# This mirrors the lung cancer example from the course (therapy × sex)lung_grouped <-data.frame(therapy =c("sequential", "sequential", "alternative", "alternative"),sex =c("male", "female", "male", "female"),progressive =c(28L, 4L, 41L, 12L),no_change =c(45L, 12L, 44L, 7L),partial =c(29L, 5L, 20L, 3L),complete =c(26L, 2L, 20L, 1L))# Fit multinomial with count matrix on the left-hand sidefit_grouped <-multinom(cbind(progressive, no_change, partial, complete) ~ therapy + sex,data = lung_grouped, trace =FALSE)summary(fit_grouped)
cat("\nNote: with only", n_patterns, "covariate patterns,\n","the model is nearly saturated — df is very small.\n")
Note: with only 4 covariate patterns,
the model is nearly saturated — df is very small.
ImportantGrouped data: when does GOF testing work?
The chi-squared approximation for deviance and Pearson \(X^2\) requires each expected cell count \(E_{ik} \geq 5\). With only 4 covariate patterns and 4 response categories (16 cells total), you need at least 80 total observations. The lung cancer data has 304 total observations — borderline acceptable. Check: colSums(lung_grouped[,3:6]) to see if any category has very few total observations.
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.
For Task 3: relevel the merged outcome, fit the constrained model, LRT
For Task 4: predict(fit4_full, type="probs") gives the probability matrix; group by apply(...,1,max) and use cut(..., quantile(...))
4 Retrieval Challenge
NoteModel answers
Q1: Model B is preferred (lower AIC). The difference of 7 is substantial — commonly, \(\Delta\)AIC > 4 is considered meaningful.
Q2: BIC’s penalty is \(\log(n)\) per parameter. When \(n > 8\), \(\log(n) > 2\), so BIC penalises complexity more than AIC and selects simpler models. With large \(n\) (e.g., \(n=1000\), \(\log(n)\approx 6.9\)), BIC is much more conservative.
5 Why-Questions
CautionQ: “A non-significant goodness-of-fit test means the model fits well. Right?”
Not exactly. A non-significant GOF test means we fail to detect lack of fit — but with small \(n\) (like \(n=40\)), the test has low power. A poorly fitting model can produce a non-significant test simply because we don’t have enough data to detect the misspecification.
Conversely, with very large \(n\), even trivial deviations from the model will produce significant GOF statistics. In those cases, GOF testing is less useful; focus on predictive accuracy and practical significance instead.
CautionQ: “I’m comparing AIC across 8 models. Do I just pick the lowest AIC?”
Approximately, yes — but: 1. AIC differences < 2 are not meaningful (models are essentially equivalent) 2. AIC assumes models are estimated on the same data — don’t compare across datasets or after removing missing data 3. AIC optimises predictive accuracy; BIC approximates the Bayesian posterior model probability under equal priors — choose based on your goal 4. If you are doing inference (not prediction), the LRT on a specific pre-specified hypothesis is more principled than AIC selection
6 Spaced Review
NoteFurther reading
Agresti, A. (1996). An Introduction to Categorical Data Analysis. Wiley. — Chapter 6, Section 1.