Estimating Partial Proportional-Odds Models

From a Partial Proportional-Odds Model to a Likelihood

A partial proportional-odds model lets some predictor effects vary across the cumulative cutpoints while keeping other effects common. Estimation therefore has two connected tasks: maximize the likelihood for the observed categories, and ensure that the resulting cumulative probabilities define valid category probabilities. An optimizer reporting convergence addresses neither task by itself.

Let \(Y_i\) take one of the ordered values \(1,\ldots,K\), for independent observations \(i=1,\ldots,n\). Write \(x_i\) for predictors with common slopes and \(z_i\) for predictors whose slopes may vary by cutpoint. Using the cumulative-logit sign convention, define

\operatorname{logit}\{P(Y_i\leq j\mid x_i,z_i)\} =\eta_{ij} =\alpha_j-x_i^\mathsf{T}\beta-z_i^\mathsf{T}\gamma_j, \qquad j=1,\ldots,K-1.

Here \(\operatorname{logit}(u)=\log\{u/(1-u)\}\), and \(G(t)=1/(1+\exp(-t))\) is its inverse. The vectors \(\gamma_j\) are cutpoint-specific; a predictor omitted from \(z_i\) has a common effect through \(\beta\). The model's cumulative probabilities are \(C_{ij}=G(\eta_{ij})\).

Definition: For an observation with covariates \(x_i,z_i\), the model-implied category probabilities are \(p_{i1}=C_{i1}\), \(p_{ik}=C_{ik}-C_{i,k-1}\) for \(2\leq k\leq K-1\), and \(p_{iK}=1-C_{i,K-1}\). The likelihood is \(L(\theta)=\prod_{i=1}^n p_{i,Y_i}\), where \(\theta\) collects all model parameters, and the log-likelihood is \(\ell(\theta)=\sum_{i=1}^n\log p_{i,Y_i}\).

These expressions require more than a formula for each cumulative logit: adjacent cumulative probabilities must be ordered. With common slopes alone, ordered thresholds ensure this because the common predictor term cancels between cutpoints. With relaxed slopes, the full predictors \(\eta_{ij}\) must be checked as well.

Validity Constraints and Ordered Thresholds

A conventional parameter constraint is \(\alpha_1<\alpha_2<\cdots<\alpha_{K-1}\). One way to enforce it in an optimizer is to set \(\alpha_1=a_1\) and \(\alpha_j=\alpha_{j-1}+\exp(\delta_j)\) for \(j=2,\ldots,K-1\). Every finite value of the unconstrained \(\delta_j\) then produces strictly ordered thresholds. This convention is useful, but does not by itself guarantee valid probabilities when slopes vary across cutpoints.

Theorem: At a fixed covariate value, all category probabilities are strictly positive if and only if \(\eta_{i1}<\eta_{i2}<\cdots<\eta_{i,K-1}\). They are nonnegative if and only if the sequence is nondecreasing.

Proof. The logistic function \(G\) is strictly increasing, since \(G'(t)=G(t)\{1-G(t)\}>0\) for every finite \(t\). Thus \(p_{i1}=G(\eta_{i1})>0\) and \(p_{iK}=1-G(\eta_{i,K-1})>0\) for finite predictors. For each intermediate category,

p_{ik}=G(\eta_{ik})-G(\eta_{i,k-1}).

Strict increase of \(G\) means this difference is positive exactly when \(\eta_{ik}>\eta_{i,k-1}\), and nonnegative exactly when \(\eta_{ik}\geq\eta_{i,k-1}\). Applying this to every adjacent pair proves both claims. Also, the category probabilities sum to one: the differences telescope to \(C_{i,K-1}\), which, together with \(C_{i1}\) and the intermediate differences, gives \(C_{i,K-1}\) in total; adding \(1-C_{i,K-1}\) yields one. \(\square\)

For estimation, it is common to impose ordered thresholds and also constrain the full predictors to be ordered at the covariate values where fitted probabilities are required. With a finite data set, these can be inequality constraints for each observed covariate pattern. Validity at new prediction covariates is not implied by checks only at the training observations; check the intended prediction range too. Depending on the design, full-predictor ordering at a particular covariate value can hold even if thresholds are not ordered. Threshold order is a conventional structural constraint, not a substitute for checking the cumulative predictors.

Worked Example: Evaluating One Observation's Likelihood

Suppose \(K=3\), with one shared predictor \(x\), one relaxed predictor \(z\), and parameters \(\alpha_1=-1\), \(\alpha_2=1\), \(\beta=0.5\), \(\gamma_1=0.2\), and \(\gamma_2=-0.3\). For an observation with \(x=z=1\), calculate the category probabilities and the log-likelihood contribution if \(Y=2\).

The cumulative predictors are \(\eta_1=-1-(1)(0.5)-(1)(0.2)=-1.7\) and \(\eta_2=1-(1)(0.5)-(1)(-0.3)=0.8\). They are ordered, so the category probabilities are positive. Applying \(G\) gives \(C_1=G(-1.7)\approx0.1545\) and \(C_2=G(0.8)\approx0.6900\).

Therefore \(p_1\approx0.1545\), \(p_2\approx0.6900-0.1545=0.5355\), and \(p_3\approx1-0.6900=0.3100\). These sum to one, up to rounding. For the observed category \(Y=2\), the contribution is \(\log(p_2)\), or \(\log(0.5355)\approx-0.625\). In an actual fit, the likelihood uses the unrounded probabilities.

Worked Example: Ordered Thresholds Do Not Guarantee Valid Probabilities

First consider \(\alpha_1=0,\alpha_2=1,\gamma_1=-2,\gamma_2=2\), with \(z=1\) and no shared-predictor contribution. The thresholds are ordered, but the cumulative predictors are \(\eta_1=0-(-2)=2\) and \(\eta_2=1-2=-1\). Thus \(C_1=G(2)\approx0.8808\) and \(C_2=G(-1)\approx0.2689\). The middle-category probability is \(C_2-C_1\approx-0.6119\), which is invalid. Threshold ordering alone has not prevented a failure.

Conversely, consider \(\alpha_1=1,\alpha_2=0,\gamma_1=0,\gamma_2=-2\), again with \(z=1\) and no shared contribution. The thresholds are reversed, but the full predictors are \(\eta_1=1\) and \(\eta_2=2\). Hence the probabilities are \(G(1)\approx0.7311\), \(G(2)-G(1)\approx0.1497\), and \(1-G(2)\approx0.1192\). All are positive and sum to one at this covariate value. This illustrates why validity is determined by the full cumulative predictors; it does not mean that a fit should abandon the usual ordered-threshold constraint.

Score Equations and What an Optimizer Must Maximize

The objective is the sum of the log probabilities of the observed categories, not a separate binary likelihood for each cutpoint. The cutpoint indicators for one ordinal response are dependent, so treating them as independent binary observations would generally produce a different objective. A useful general score formula follows by differentiating the category probability.

Result: Let \(D_{ij}=\partial\eta_{ij}/\partial\theta\), treated as a column vector, and let \(f(t)=G(t)\{1-G(t)\}\). For observation \(i\), the gradient of its log-likelihood contribution is \(f(\eta_{i1})D_{i1}/p_{i1}\) if \(Y_i=1\); it is \(\{f(\eta_{ik})D_{ik}-f(\eta_{i,k-1})D_{i,k-1}\}/p_{ik}\) if \(Y_i=k\) is an intermediate category; and it is \(-f(\eta_{i,K-1})D_{i,K-1}/p_{iK}\) if \(Y_i=K\). The score is the sum of these contributions over observations.

Derivation. For the first category, \(p_{i1}=G(\eta_{i1})\). The chain rule gives \(\partial p_{i1}/\partial\theta=f(\eta_{i1})D_{i1}\), so \(\partial\log p_{i1}/\partial\theta=f(\eta_{i1})D_{i1}/p_{i1}\). For an intermediate category, differentiating \(p_{ik}=G(\eta_{ik})-G(\eta_{i,k-1})\) gives \(f(\eta_{ik})D_{ik}-f(\eta_{i,k-1})D_{i,k-1}\). Dividing by \(p_{ik}\) gives the stated log-likelihood gradient. Finally, \(p_{iK}=1-G(\eta_{i,K-1})\), whose derivative is \(-f(\eta_{i,K-1})D_{i,K-1}\); division by \(p_{iK}\) gives the last expression. Summing the individual derivatives yields the score of \(\ell(\theta)\). \(\square\)

For the parameterization above, the derivatives are explicit: \(\partial\eta_{ij}/\partial\alpha_r\) is one when \(r=j\) and zero otherwise; \(\partial\eta_{ij}/\partial\beta=-x_i\); and \(\partial\eta_{ij}/\partial\gamma_r=-z_i\) when \(r=j\), zero otherwise. A constrained maximum-likelihood estimate maximizes \(\ell(\theta)\) over parameter values satisfying the chosen threshold and probability-validity constraints. If the optimum is interior, the score is zero there. At an active constraint, the appropriate first-order conditions are the constrained (KKT) conditions, not simply a zero unconstrained score.

Worked Example: A Small Partial Proportional-Odds Fit

Suppose there are three categories and a binary predictor \(z\) whose effect is relaxed at both cutpoints. The observed counts are \((10,20,10)\) at \(z=0\) and \((20,10,20)\) at \(z=1\). With \(\eta_j=\alpha_j-z\gamma_j\), this is a saturated model for the two group-specific category distributions. Find its maximum-likelihood estimates.

For a group with counts \(n_1,n_2,n_3\), the likelihood for its category probabilities is proportional to \(p_1^{n_1}p_2^{n_2}p_3^{n_3}\), subject to \(p_1+p_2+p_3=1\). Maximizing \(\sum_k n_k\log p_k\) gives \(n_k/p_k=\lambda\) for each category. Summing the resulting probabilities shows \(\lambda=n_1+n_2+n_3\), so the maximum occurs at \(p_k=n_k/(n_1+n_2+n_3)\). All six counts here are positive, so these estimates are interior and finite.

At \(z=0\), the estimated probabilities are \((0.25,0.50,0.25)\), with cumulative probabilities \(0.25,0.75\). Their logits are \(\alpha_1=\log(0.25/0.75)=-1.0986\) and \(\alpha_2=\log(0.75/0.25)=1.0986\). At \(z=1\), the probabilities are \((0.40,0.20,0.40)\), with cumulative probabilities \(0.40,0.60\) and logits \(-0.4055,0.4055\). Since \(\eta_j=\alpha_j-\gamma_j\) at \(z=1\), the estimates are \(\gamma_1=-1.0986-(-0.4055)=-0.6931\) and \(\gamma_2=1.0986-0.4055=0.6931\).

The cutpoints are ordered, and the full predictors are ordered in both groups: \((-1.0986,1.0986)\) at \(z=0\), and \((-0.4055,0.4055)\) at \(z=1\). As a likelihood check, the maximized log-likelihood, omitting no parameter-dependent terms from the category likelihood, is \(20\log(0.25)+20\log(0.50)+40\log(0.40)+10\log(0.20)\approx-94.3348\). This direct calculation agrees with evaluating the fitted category probabilities against the observed counts.

Convergence, Finite Estimates, and Probability Checks

Maximum-likelihood fitting is usually numerical. In a partial proportional-odds model, there may be inequality constraints from both ordered thresholds and full-predictor ordering. Moreover, sparse patterns or near-separation can push estimates toward very large magnitudes, even when an algorithm terminates. A finite reported iteration limit or a small change in the objective is not, alone, evidence of a reliable finite maximum.

Check the optimization status.
Confirm that the algorithm terminated successfully and that the objective and parameter values are finite. For an interior solution, inspect the score; for an active constraint, inspect constrained optimality diagnostics.

Check stability and curvature.
Compare estimates and log-likelihoods across sensible starting values. At an interior local maximum, the observed information should support local maximum curvature; near-singularity can indicate weak identification or unreliable standard errors.

Check the fitted probabilities.
For every relevant covariate pattern, calculate all cumulative predictors in order, then calculate the category probabilities. Confirm that they are nonnegative (strictly positive for finite, nonboundary fits) and sum to one, allowing only small numerical rounding error.

If probabilities are formed by subtracting two nearly equal cumulative probabilities, direct computation can lose numerical precision. Use stable log-difference calculations when evaluating the log-likelihood, and distinguish a tiny positive probability from a rounded zero. Do not silently replace a negative fitted category probability with zero: it signals an invalid predictor ordering or a numerical failure that requires investigation.

Common Mistake

A frequent error is to impose \(\alpha_1<\cdots<\alpha_{K-1}\), observe that the optimizer converged, and conclude that every fitted category probability is valid. The reasoning fails because relaxed slopes can reverse the order of \(\eta_{ij}\) at particular covariate values even when the thresholds are ordered. The correction is to check the order of the full cumulative predictors at the observed and intended prediction covariates, then form and inspect the category probabilities. Conversely, reversed thresholds do not by themselves prove invalid probabilities at a specified covariate value when cutpoint-specific slopes are present; the full predictors determine validity there. In standard fitting, retain the conventional ordered-threshold constraint and treat full-predictor ordering as an additional, distinct requirement.

Key takeaway: A trustworthy partial proportional-odds fit maximizes the ordinal category likelihood under its constraints, has finite and stable estimates, and yields valid category probabilities wherever the model will be used. Convergence, threshold ordering, and probability validity are separate checks.

Check Your Understanding

Use the model definitions and checks above to answer these questions.

Write the category probabilities for a response with four ordered categories in terms of its three cumulative probabilities.
Why does ordering the thresholds not guarantee valid probabilities when a predictor has cutpoint-specific slopes?
At a fixed covariate value, what condition on the cumulative predictors is equivalent to strictly positive intermediate-category probabilities?
In a constrained fit with an active ordering constraint, why is a zero unconstrained score not the only appropriate optimality check?
What additional probability check is needed before using a fitted model at covariate values outside the observed data?