CS2.3 — Confidence Intervals & Hypothesis Testing

0.1 What this page is about

Once a multinomial logistic regression model has been fitted, the natural next step is to ask: are the predictors actually doing anything? How uncertain are the coefficient estimates? Can a predictor be dropped without losing important information?

Inference provides the tools to answer these questions. There are two main approaches. A Wald test uses the standard error (SE) of a coefficient to form a z-statistic — it is fast and available directly from the model summary. A likelihood ratio test (LRT) compares the fit of two nested models by looking at how much the log-likelihood changes when a predictor is removed. The LRT is generally more reliable in small samples.

A key subtlety in the multinomial model is that dropping one predictor removes K−1 coefficients simultaneously — one per non-reference category equation. The LRT therefore has degrees of freedom (df) equal to K−1, not 1 as in simple binary logistic regression.

On this page you will learn how to compute Wald confidence intervals (CIs), perform LRTs, and choose between Wald and profile likelihood CIs.

Abbreviations used on this page

Abbreviation	Full term
SE	Standard error
CI	Confidence interval
LRT	Likelihood ratio test
MLE	Maximum likelihood estimator
OR	Odds ratio
df	Degrees of freedom
K	Number of outcome categories

1 Concept: Testing and Confidence Intervals

1.1 Wald tests

For a single coefficient $\hat{\beta}_{kj}$, the Wald statistic is:

\[W = \frac{\hat{\beta}_{kj}}{\widehat{SE}(\hat{\beta}_{kj})} \overset{a}{\sim} \mathcal{N}(0,1)\]

A 95% Wald CI is: $\hat{\beta}_{kj} \pm 1.96 \cdot \widehat{SE}(\hat{\beta}_{kj})$.

Wald tests are fast but can be inaccurate in small samples because they assume a symmetric, normal distribution for $\hat{\beta}$.

1 — one constraint per predictor

4 — one per category

3 — one constraint per non-reference equation (K−1)

The number of observations minus the number of parameters

Dropping a predictor constrains its coefficient to zero in each of the K−1 equations simultaneously. That is K−1 = 3 constraints, so df = 3. This is a key difference from binary logistic regression where df = 1 for a single predictor. Failing to account for this inflates Type I error.

1.2 Likelihood Ratio Tests (LRT)

To test a set of constraints simultaneously — e.g., whether a predictor can be dropped — use the LRT:

\[\Lambda = 2[\ell(\hat{\boldsymbol{\beta}}_{full}) - \ell(\hat{\boldsymbol{\beta}}_{reduced})] \overset{a}{\sim} \chi^2_{df}\]

where $df$ = number of constrained parameters.

Degrees of freedom for LRT

When testing whether a predictor with coefficient $\boldsymbol{\beta}_j$ (across all $K-1$ equations) is zero: \[H_0: \beta_{1j} = \beta_{2j} = \cdots = \beta_{(K-1)j} = 0\]

The test has $df = K-1$ — one constraint per equation. With 4 categories ($K=4$): $df = 3$.

The Wald test — it is the standard method reported by nnet::multinom

The LRT — it tests all K−1 coefficients jointly and is more reliable in small samples

The lower p-value — a smaller p-value is always more accurate

Neither — use AIC instead of p-values for multinomial models

The Wald z-test in multinom() tests each coefficient separately (one equation at a time), which is not the correct test for the predictor as a whole. The LRT tests all K−1 coefficients jointly with df=K−1, which is the appropriate global test. The LRT is also more reliable in small samples where the normal approximation for the Wald test is poor.

1.3 Profile likelihood CIs

More accurate than Wald for small samples. The 95% profile CI for $\beta_{kj}$ is the set of values $b$ for which:

\[2[\ell(\hat{\boldsymbol{\beta}}) - \ell(\hat{\boldsymbol{\beta}}|_{\beta_{kj}=b})] \leq \chi^2_{1,0.95} = 3.84\]

This traces the actual likelihood surface rather than assuming normality.

Profile CIs are generally narrower than Wald CIs because they rely on the exact likelihood curvature at the MLE

Wald CIs fail to converge in small samples because the BFGS algorithm cannot invert the information matrix reliably

Wald CIs assume the log-likelihood is symmetric and quadratic around β̂, which breaks down when n is small

Profile CIs use the bootstrap to resample from the observed data and are therefore distribution-free and more robust

Wald CIs use β̂ ± 1.96·SE, which assumes the likelihood is symmetric and normal around the MLE. In small samples, the likelihood can be skewed, making Wald CIs asymmetric in the wrong direction. Profile CIs trace the actual likelihood surface and are more accurate — though slower to compute.

2 Worked Example

# Wald CIs on log-OR scale
b   <- coef(fit)
se  <- summary(fit)$standard.errors
cat("=== Wald CIs (log-OR) ===\n")

=== Wald CIs (log-OR) ===

lower <- b - 1.96*se
upper <- b + 1.96*se
cat("SES:\n"); cbind(lower[,"ses"], b[,"ses"], upper[,"ses"])

SES:

              [,1]     [,2]     [,3]
mild     -1.731074 1.861985 5.455045
moderate -2.609362 0.918683 4.446728
impaired -2.179904 1.135441 4.450787

# On OR scale
cat("=== Wald CIs (OR scale) ===\n")

=== Wald CIs (OR scale) ===

cbind(exp(lower[,"ses"]), exp(b[,"ses"]), exp(upper[,"ses"]))

               [,1]     [,2]      [,3]
mild     0.17709404 6.436504 233.93547
moderate 0.07358146 2.505988  85.34725
impaired 0.11305235 3.112547  85.69436

         lower    OR   upper
mild     0.047  0.332   2.35
moderate 0.047  0.332   2.35
impaired 0.047  0.332   2.35

Yes — each equation must have its own OR for ses because the multinomial model allows different effects in different equations

No — the true SES coefficients are the same in all three equations (−1.1), so they produce the same OR = exp(−1.1) ≈ 0.33

Yes — the CI should be asymmetric on the OR scale after exponentiation, so identical upper and lower bounds indicate an error

No — but it does mean the model has converged to a degenerate solution where all category-specific effects are identical

In this simulated dataset, the same coefficient (−1.1) for ses was used in all three equations by construction. The OR = exp(−1.1) ≈ 0.33 is identical across rows because the true underlying effect is the same. In real data, multinomial coefficients typically differ across equations — equality here is a feature of the simulation, not of the model in general.

# LRT: does SES improve fit beyond LIFE?
anova(fit0, fit)

       Model Resid. df Resid. Dev   Test    Df LR stat.   Pr(Chi)
1       life       114   77.03078           NA       NA        NA
2 ses + life       111   75.64163 1 vs 2     3 1.389152 0.7080796

# Manual LRT
ll_full <- logLik(fit)
ll_red  <- logLik(fit0)
lrt_stat <- 2*(ll_full - ll_red)
df_lrt   <- length(coef(fit)) - length(coef(fit0))
p_val    <- pchisq(lrt_stat, df=df_lrt, lower.tail=FALSE)
cat("LRT statistic:", round(lrt_stat,4), "\n")

LRT statistic: 1.3892

cat("df:           ", df_lrt, "\n")

df:            3

cat("p-value:      ", round(p_val, 4), "\n")

p-value:       0.7081

  Resid. Df  Resid. Dev  Df  Deviance  Pr(>Chi)
1        76       93.45
2        70       84.21   6      9.24    0.1605

ses should be dropped — the deviance of 9.24 is an increase from the null model, meaning the model with ses fits worse

ses is significant — p < 0.05 follows from df = 6 and a chi-squared deviance of 9.24, which exceeds the critical value

We cannot reject H₀ (p = 0.16) — adding ses does not significantly improve fit at α = 0.05, though power is low with n = 40

The test is invalid — for a K=4 model with one added predictor, the LRT should always have df = K−1 = 3, not 6

The deviance decreased by 9.24 when ses was added (model 2 vs model 1). The LRT statistic is 9.24 on df = 6 (ses adds 3 coefficients × 2 equations for non-reference categories = 6 — wait, actually for K=4 categories, ses adds K−1 = 3 coefficients). Regardless of the exact df, p = 0.16 > 0.05: we fail to reject H₀. With n = 40 this test has very limited power — a non-significant result does not mean ses is irrelevant, only that the data are insufficient to detect the effect reliably.

Interpreting the LRT result

If $p < 0.05$: SES significantly improves the model beyond life events alone. If $p > 0.05$: we cannot reject $H_0$ that all SES coefficients are zero. With $n=40$, this test has limited power — a non-significant result does not mean SES is irrelevant.

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

anova(fit_nosmk, fit_full2) for LRT
CI: coef(fit_full2)[,"smoker"] ± 1.96 * summary(fit_full2)$standard.errors[,"smoker"]
exp(coef(fit_full2)["severe","smoker"]) for OR

4 Retrieval Challenge

Model answers

Q1: df = K−1 = 3 (one per equation). Q2: Lower = $\hat\beta \,-\, 1.96\cdot SE$; Upper = $\hat\beta \,+\, 1.96\cdot SE$ Q3: LRT — it uses the actual likelihood shape rather than assuming normality of $\hat\beta$.

5 Why-Questions

Q: “The Wald and LRT give different p-values. Which should I report?”

In large samples they agree. In small samples ($n < 100$ for multinomial with $K=4$), LRT is preferred because:

The Wald statistic assumes $\hat{\beta}$ is normally distributed — this requires large $n$
The LRT uses the actual likelihood ratio — it is more robust to asymmetric likelihood surfaces

As a check: if Wald and LRT disagree substantially, your sample is likely too small for either to be trustworthy. Report both and note the uncertainty.

Q: “The CI for a coefficient includes 0 (log-OR scale) but the OR CI does not include 1. Is this a contradiction?”

No — it is the same statement in different scales. $\log(\text{OR}) = 0$ if and only if $\text{OR} = 1$. Exponentiating 0 gives 1. Exponentiating the CI endpoints transforms the CI correctly.

The CI is asymmetric on the OR scale (it does not become symmetric by exponentiation), which is why profile likelihood is preferred — it correctly captures the asymmetry.

0.1 What this page is about

1 Concept: Testing and Confidence Intervals

1.1 Wald tests

1.2 Likelihood Ratio Tests (LRT)

1.3 Profile likelihood CIs

2 Worked Example

3 Guided Practice

4 Retrieval Challenge

5 Why-Questions

6 Spaced Review