A fitted proportional odds (PO) model gives you coefficient estimates β but how confident should you be in them? Could the true effect be zero? How wide is the range of plausible values?
Inference answers these questions through hypothesis tests and confidence intervals (CIs). The same two approaches from CS2.3 apply here. A Wald test computes a z-statistic from the coefficient divided by its standard error (SE) and compares it to a normal distribution. A likelihood ratio test (LRT) compares the fit of the full model against a reduced model with the predictor removed; the difference in log-likelihoods follows a chi-squared distribution.
One practical difference from multinomial inference: in the PO model there is only one coefficient per predictor (not Kβ1), so dropping a predictor from the model costs only 1 degree of freedom (df) in the LRT rather than Kβ1. This makes inference simpler and more powerful.
On this page you will learn how to extract Wald CIs from polr() and clm(), how to run LRTs, and why profile likelihood CIs are preferred over Wald CIs in smaller samples.
NoteAbbreviations used on this page
Abbreviation
Full term
PO
Proportional odds
CI
Confidence interval
SE
Standard error
LRT
Likelihood ratio test
OR
Odds ratio
df
Degrees of freedom
MLE
Maximum likelihood estimator
1 Concept: Inference for the Proportional Odds Model
1.1 Wald tests and CIs
The MLE \(\hat{\boldsymbol{\beta}}\) is asymptotically normal. Wald CI for \(\beta_j\):
summary(polr(...)) shows \(t\)-values but not \(p\)-values, because the degrees of freedom for the \(t\)-distribution are ambiguous. Convert manually: \(p = 2(1 - \Phi(|t|))\) using the normal approximation.
Use anova() on the polr object β it automatically reports p-values
Compute p = 2Β·(1 β Ξ¦(|t|)) treating the t-value as a z-statistic
Divide the coefficient by its standard error and compare to a t-distribution with nβp df
p-values are not valid for ordinal regression β use CIs only
polr omits p-values because the df for the t-distribution are ambiguous in ordinal models (same issue as lmer). The standard workaround is to treat the t-statistic as a z-score under the asymptotic normal approximation: p = 2Β·(1 β Ξ¦(|t|)). This is anti-conservative in small samples β use the LRT for small n.
1.2 Profile likelihood CIs
confint.polr uses profile likelihood by default β more accurate than Wald, especially for small samples. It traces the likelihood surface along each parameter direction.
Wald CI β it is the fastest and most widely reported method and is always recommended for ordinal models
Bootstrap CI β it makes no distributional assumptions and is therefore the most robust option in small samples
Profile likelihood CI β it traces the actual likelihood surface and is more accurate than Wald in small samples
Bayesian credible interval β it uses a flat prior for the regression coefficients and is model-free
confint.polr uses profile likelihood by default. It finds the set of Ξ² values where 2[β(Ξ²Μ) β β(Ξ²Μ|Ξ²_j=b)] β€ 3.84 (the 95% contour). This accounts for asymmetry in the likelihood surface, giving more accurate intervals than Wald β especially for small samples or parameters near boundary values.
1.3 Likelihood Ratio Test for a predictor
To test \(H_0: \beta_j = 0\) (dropping predictor \(j\)):
1 β because the PO model has a single coefficient per predictor shared across all cut-points
Kβ1 β one constraint per cumulative equation because each equation has an independent coefficient
n β p β the residual degrees of freedom, just as in linear regression for any single coefficient test
It depends on the predictor type β categorical predictors with L levels cost Lβ1 degrees of freedom
The PO assumption means Ξ²_j is shared across all Kβ1 cumulative equations. Setting Ξ²_j = 0 imposes exactly one constraint, so df = 1. This is a key efficiency gain over multinomial regression, where dropping a predictor imposes Kβ1 constraints (df = Kβ1). For a K-level categorical predictor in PO, df = Kβ1 (one constraint per dummy variable).
# Compare Wald and profile CIs for sesse_ses <- coef_tbl["ses","Std. Error"]beta_ses <-coef(fit)["ses"]wald_ci <-c(beta_ses -1.96*se_ses, beta_ses +1.96*se_ses)prof_ci <- ci_prof["ses",]cat("=== Comparison for ses ===\n")
=== Comparison for ses ===
cat("Wald CI :", round(wald_ci, 4), "\n")
Wald CI : -1.2102 -0.1639
cat("Profile CI:", round(prof_ci, 4), "\n")
Profile CI: -1.2145 -0.1668
cat("\nOR scale:\n")
OR scale:
cat("Wald CI :", round(exp(wald_ci), 4), "\n")
Wald CI : 0.2981 0.8488
cat("Profile CI:", round(exp(prof_ci), 4), "\n")
Profile CI: 0.2969 0.8464
# The same inferences with ordinal::clm β p-values printed directlylibrary(ordinal)fit_clm <-clm(impair ~ ses + life, data=dat)summary(fit_clm)
formula: impair ~ ses + life
data: dat
link threshold nobs logLik AIC niter max.grad cond.H
logit flexible 200 -244.59 499.18 4(0) 3.48e-08 3.8e+02
Coefficients:
Estimate Std. Error z value Pr(>|z|)
ses -0.68704 0.26690 -2.574 0.01 *
life 0.33407 0.04891 6.830 8.49e-12 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Threshold coefficients:
Estimate Std. Error z value
well|mild -0.5288 0.2822 -1.874
mild|moderate 0.7822 0.2842 2.752
moderate|impaired 2.0388 0.3150 6.473
# Profile likelihood CIs from clm (uses observed information matrix)ci_clm <-confint(fit_clm)cat("Profile CIs from clm (log-OR scale):\n")
Profile CIs from clm (log-OR scale):
print(round(ci_clm, 4))
2.5 % 97.5 %
ses -1.2144 -0.1668
life 0.2401 0.4322
cat("\nOR scale:\n")
OR scale:
print(round(exp(ci_clm), 4))
2.5 % 97.5 %
ses 0.2969 0.8464
life 1.2714 1.5406
NoteWhy clm prints p-values but polr does not
ordinal::clm computes p-values using the normal (z) approximation automatically and prints them in summary(). MASS::polr deliberately omits them to force the user to be explicit about their choice of reference distribution. Both approaches yield the same test statistic β the difference is just in how the output is presented.
The clm p-values use the observed information matrix (ββΒ²β/βΞ²Β² evaluated at Ξ²Μ), while polr uses the expected information matrix (E[ββΒ²β/βΞ²Β²]). In the cumulative logit model these differ slightly, so you may see small numerical differences in standard errors between the two packages.
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.
MASS::polr reports \(t\)-values but omits p-values because the appropriate reference distribution is ambiguous:
Strictly, under the Wald framework, the statistic is asymptotically \(\mathcal{N}(0,1)\) (z-test)
But with small samples or many parameters, \(t\)-distribution with some estimated df may be more accurate
There is no obvious df formula for ordinal regression (unlike linear mixed models where lmerTest provides approximations)
The convention in many R packages is to let the user choose. For a z-test p-value: 2*(1-pnorm(abs(t_value))). For a rigorous test, prefer the LRT.
CautionQ: βProfile likelihood CI computation takes much longer than Wald. Is it always worth it?β
In large samples (\(n > 500\), many groups), Wald and profile CIs converge and the difference is negligible. Wald is preferable for its speed.
For small or medium samples (\(n < 200\)), especially for variance components or boundary-case parameters, profile likelihood provides meaningfully more accurate intervals. The computation is asymptotically \(O(p)\) likelihood maximisations β expensive for high-dimensional \(\boldsymbol{\beta}\), but tractable for the 2β5 parameters typical in ordinal models.
6 Spaced Review
NoteFurther reading
Agresti, A. (2010). Analysis of Ordinal Categorical Data (2nd ed.). Wiley. β Chapter 3.