CS3.3 — Testing the Proportional Odds Assumption

0.1 What this page is about

The proportional odds (PO) model buys parsimony by assuming that the effect of each predictor is the same at every cumulative split of the outcome. This is a strong assumption. In practice it may hold for some predictors but not others — and using the model when it clearly fails can give misleading results.

Testing the PO assumption means asking: do the data support a single slope per predictor, or do the slopes differ across category boundaries? A visual check is always the first step: plot the empirical log-odds for each binary split of the outcome against the predictor. If the curves are roughly parallel, the PO assumption is plausible.

For a formal test, you can compare the PO model against an unconstrained multinomial model using a likelihood ratio test (LRT). The ordinal::clm package provides a convenient nominal_test() function that does this for each predictor individually. A significant result suggests that predictor’s effect varies across cut-points.

When the PO assumption is violated you have several options: fit a partial PO model (relax the assumption for specific predictors), use the full multinomial model, or consider an adjacent-category model. On this page you will learn all of these steps.

Abbreviations used on this page

Abbreviation	Full term
PO	Proportional odds
LRT	Likelihood ratio test
OR	Odds ratio
df	Degrees of freedom
CLM	Cumulative link model

1 Concept: When Parallelism Fails

1.1 What the PO assumption states

The PO assumption says: for predictor $x_j$ and any cut-point $k$:

\[\text{logit}\,P(Y \leq k | x_j = a) - \text{logit}\,P(Y \leq k | x_j = b) = -\beta_j(a-b) \quad \forall k\]

Equivalently, the regression lines for different cut-points are parallel on the logit scale.

1.2 Visual check: empirical log-odds

Before formal testing, plot the empirical log-odds for each binary split against each predictor:

# Create binary splits and estimate log-ORs empirically
splits <- c("well","mild","moderate")

# For each split point: P(Y <= split) vs ses
lor_ses <- sapply(splits, function(k){
  p1 <- mean(dat$impair <= k & dat$ses==1) / mean(dat$impair >  k & dat$ses==1)
  p0 <- mean(dat$impair <= k & dat$ses==0) / mean(dat$impair >  k & dat$ses==0)
  log(p1) - log(p0)
})

cat("Empirical log-OR for ses at each cut-point:\n")

Empirical log-OR for ses at each cut-point:

print(round(lor_ses, 3))

    well     mild moderate 
   0.506    0.619    0.641

cat("\nModel log-OR for ses (should be constant under PO):",
    round(-coef(fit_po)["ses"], 3), "\n")


Model log-OR for ses (should be constant under PO): 0.687

If the empirical log-ORs are approximately constant across cut-points, PO holds.

Empirical log-OR for ses at each cut-point:
   well    mild  moderate
  -1.23   -1.05    -0.88

Model log-OR for ses (should be constant under PO): -0.898

No — the values differ across cut-points (−1.23 vs −0.88), so PO is violated and a separate slope is needed for ses at each cut-point

Yes — the values are roughly constant across cut-points and close to the model estimate; the variation is within sampling noise for n = 200

No — the empirical ORs should all equal exactly the model log-OR (−0.898) if PO holds, but the variation here is too large to accept

Cannot tell without a formal test — visual checks of empirical log-ORs are never informative enough to assess the PO assumption

Under PO, the true log-OR is the same at every cut-point. The empirical values (−1.23, −1.05, −0.88) vary because they are estimated from data — sampling variation is expected. The range from −1.23 to −0.88 spans only about 0.35 log-OR units, and all values are on the same side and close to the model estimate (−0.898). This is consistent with PO. A systematic trend (e.g., values changing from −2 to −0.2) or a sign change would be more alarming.

The log-odds increase linearly with the cut-point index

The log-ORs for a predictor are approximately constant across all cut-points

The log-odds are the same for all predictors at each cut-point

The cumulative probabilities are symmetric around 0.5

Under PO, the log-OR for predictor x is −β at every cut-point. So plotting the empirical log-OR at each cut-point should give roughly horizontal (constant) values. Systematic trends (e.g., increasing or decreasing log-ORs across cut-points) suggest PO violation. This visual check is informal — use it alongside formal tests.

1.3 Formal test: LRT against a non-proportional model

The non-proportional model allows separate coefficients for each cut-point equation. In base R, this requires fitting $K-1$ binary logistic regressions:

\[\text{logit}\,P(Y \leq k | \mathbf{x}) = \alpha_k + \boldsymbol{\beta}_k^\top \mathbf{x}, \quad k=1,\ldots,K-1\]

The number of extra parameters is $(K-2) \times p_{slopes}$ (one per additional equation minus the one already in the PO model).

Degrees of freedom

With $K=4$ categories and $p=2$ predictors: - PO model: $K-1 = 3$ thresholds + $p = 2$ slopes = 5 parameters - Non-PO model: $K-1 = 3$ thresholds + $(K-1) \times p = 6$ slopes = 9 parameters - LRT df = $9 - 5 = 4$

Abandon ordinal regression entirely and switch to multinomial logistic regression for all predictors in the model

Remove the violating predictor from the model entirely — a significant LRT means its effect is too unstable to include

Fit a partial proportional odds model, allowing that predictor to have category-specific slopes while others remain proportional

Collapse adjacent categories until the PO test becomes non-significant, then refit the simplified ordinal model

A partial PO model allows some predictors to have separate slopes for each cut-point while others are constrained to be proportional. This is more flexible than full PO but more parsimonious than full multinomial regression. The R package `ordinal` (clm with nominal argument) supports this directly.

The model has severe convergence problems and the slope estimates are unreliable for all predictors in the model

The outcome should be treated as continuous rather than ordinal because the categories lack a natural ordering

At least one predictor's effect is not constant across cumulative splits — PO may be an oversimplification

The thresholds are too closely spaced to estimate reliably, leading to numerical instability during fitting

A significant Brant test suggests the parallel-regressions assumption is violated for at least one predictor. However, PO models are often robust to mild violations in large samples. Consider: (1) which predictor drives the violation, (2) whether a partial PO model is justified, (3) whether the substantive conclusions change under a relaxed model.

2 Worked Example

# Brant test: fit K-1 binary logistic regressions, compare slopes
fits_binary <- list(
  "Y<=well"     = glm(I(impair<="well")     ~ ses+life, data=dat, family=binomial),
  "Y<=mild"     = glm(I(impair<="mild")     ~ ses+life, data=dat, family=binomial),
  "Y<=moderate" = glm(I(impair<="moderate") ~ ses+life, data=dat, family=binomial)
)
cat("Estimated slopes per cut-point:\n")

Estimated slopes per cut-point:

slopes <- sapply(fits_binary, function(f) coef(f)[c("ses","life")])
print(round(slopes, 3))

     Y<=well Y<=mild Y<=moderate
ses    0.693   0.757       0.733
life  -0.400  -0.310      -0.320

cat("\nPO model slopes (should match above if PO holds):\n")


PO model slopes (should match above if PO holds):

print(round(coef(fit_po), 3))

   ses   life 
-0.687  0.334

# LRT: PO vs non-PO
# Non-PO: fit multinomial and compare log-likelihoods
library(nnet)
fit_multinom <- multinom(impair ~ ses + life, data=dat, trace=FALSE)

ll_po      <- logLik(fit_po)
ll_nonpo   <- logLik(fit_multinom)
lrt_stat   <- 2*(ll_nonpo - ll_po)
df_lrt_po  <- length(coef(fit_multinom)) - length(c(coef(fit_po), fit_po$zeta))
p_lrt      <- pchisq(lrt_stat, df=df_lrt_po, lower.tail=FALSE)

cat("LRT statistic:", round(lrt_stat, 4), "\n")

LRT statistic: 1.3972

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

df:            4

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

p-value:       0.8447

cat("\nIf p > 0.05: PO assumption is consistent with the data\n")


If p > 0.05: PO assumption is consistent with the data

PO is rejected — the LRT statistic of 4.23 exceeds the critical value

We fail to reject PO (p = 0.376 > 0.05) — the proportional odds assumption is consistent with the data

The non-proportional model is better — it uses fewer degrees of freedom

The result is inconclusive — a df = 4 test is too underpowered to draw any conclusion

The LRT tests H₀: PO holds (the proportional model is adequate) vs H₁: PO is violated (slopes differ across cut-points). p = 0.376 > 0.05: we fail to reject H₀. The extra flexibility of the non-proportional model (4 extra parameters) does not significantly improve fit. Proceed with the PO model. Note: failure to reject does not prove PO holds — the test may have low power, especially if n is small relative to the number of extra parameters being tested.

# Summary: PO model vs multinomial
cat("PO model     — df=", attr(ll_po,"df"),    ", AIC=", round(AIC(fit_po),2),    "\n")

PO model     — df= 5 , AIC= 499.18

cat("Multinomial  — df=", attr(ll_nonpo,"df"), ", AIC=", round(AIC(fit_multinom),2), "\n")

Multinomial  — df= 9 , AIC= 505.78

2.0.1 Using `ordinal::nominal_test`

The ordinal package provides nominal_test() which tests each predictor separately for PO violation. This is more targeted than the global LRT — it tells you which predictor(s) drive the violation.

library(ordinal)
fit_po_clm <- clm(impair ~ ses + life, data=dat)

# nominal_test: tests whether each predictor's effect is proportional across cut-points
# H₀ for each predictor: the PO assumption holds for that predictor
nominal_test(fit_po_clm)

Tests of nominal effects

formula: impair ~ ses + life
       Df  logLik    AIC     LRT Pr(>Chi)
<none>    -244.59 499.18                 
ses     2 -244.53 503.07 0.11315   0.9450
life    2 -243.71 501.42 1.76176   0.4144

Reading nominal_test output

Each row tests one predictor. A small p-value means that predictor has a significantly non-proportional effect. The degrees of freedom equal $(K-2) \times 1 = 2$ per predictor (for $K=4$ categories): one set of slopes is already in the PO model, and the non-PO model adds $(K-2)$ extra slopes per predictor.

2.0.2 Using VGAM with `parallel=FALSE`

The VGAM package fits the fully non-proportional model directly. This is the most flexible approach: each predictor gets its own slope per cumulative equation.

# VGAM: fit cumulative logit model without PO assumption
# Note: VGAM is not available in the browser (WebR) — this block runs at build time
library(VGAM)

# PO model (parallel=TRUE): same β across all cut-points
fit_vgam_po  <- vglm(impair ~ ses + life,
                     cumulative(parallel = TRUE,  reverse = FALSE), data = dat)

# Non-PO model (parallel=FALSE): separate β per cut-point
fit_vgam_npo <- vglm(impair ~ ses + life,
                     cumulative(parallel = FALSE, reverse = FALSE), data = dat)

# LRT comparing the two
lrtest(fit_vgam_po, fit_vgam_npo)

Likelihood ratio test

Model 1: impair ~ ses + life
Model 2: impair ~ ses + life
  #Df  LogLik Df  Chisq Pr(>Chisq)
1 595 -244.59                     
2 591 -243.64 -4 1.8994     0.7543

VGAM and partial PO

VGAM also allows a partial PO model, where PO is relaxed for only one predictor:

# Allow life to have non-proportional effects, constrain ses to be proportional
fit_partial <- vglm(impair ~ ses + life,
                    cumulative(parallel = FALSE ~ 1 + life), data=dat)
lrtest(fit_vgam_po, fit_partial)

A non-significant test here means the partial relaxation for life does not improve fit — life satisfies PO even if another predictor does not.

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

For Task 2: dat[dat$x1 > 0,] vs dat[dat$x1 <= 0,] for each binary split
For Task 3: multinom(y ~ x1+x2, data=d7, trace=FALSE)
LRT: 2*(logLik(fit_multi7) - logLik(fit_po7))

4 Retrieval Challenge

Model answers

PO: $(K-1) + p = 4 + 3 = 7$. Non-PO multinomial: $(K-1)(p+1) = 4 \times 4 = 16$. LRT df = $16 - 7 = 9$.

Minimum $n \approx 160$ for a reliable LRT. With small samples, the LRT has low power — a non-significant result does not guarantee PO holds.

5 Why-Questions

Q: “The PO test is non-significant. Does this prove the PO assumption holds?”

No — it means we fail to detect a violation. With limited data, the test may have low power. Absence of evidence is not evidence of absence.

Practical guidance: if the empirical log-ORs at different cut-points are close (within sampling variation), and the test is non-significant, it is reasonable to proceed with the PO model. Report the test result and note that the assumption was checked.

Q: “The PO test is significant. Must I abandon the ordinal model entirely?”

Not necessarily. Options in order of complexity:

Partial proportional odds: constrain only the “well-behaved” predictors to have proportional odds; allow the violating predictor to have separate slopes per equation.
Multinomial logistic regression (CS2): fully flexible, no ordering assumed.
Adjacent-category model: models log-odds of adjacent categories; intermediate flexibility.
Accept the violation if small: if the practical impact on predicted probabilities is minor (compare PO vs multinomial predictions), the PO model may still be the most useful summary.

The choice depends on sample size, scientific goals, and the magnitude of the violation.

0.1 What this page is about

1 Concept: When Parallelism Fails

1.1 What the PO assumption states

1.2 Visual check: empirical log-odds

1.3 Formal test: LRT against a non-proportional model

2 Worked Example

2.0.1 Using ordinal::nominal_test

2.0.2 Using VGAM with parallel=FALSE

3 Guided Practice

4 Retrieval Challenge

5 Why-Questions

6 Spaced Review

2.0.1 Using `ordinal::nominal_test`

2.0.2 Using VGAM with `parallel=FALSE`