Deriving the Score Test for Proportional Odds

Why nuisance adjustment matters

The proportional-odds assumption says that a predictor has the same slope in every cumulative logit. Testing that assumption is not as simple as measuring the slope differences while holding everything else fixed: the thresholds and common slopes are estimated from the same data. Their uncertainty affects how much evidence the data provide about nonparallel slopes. The score test accounts for that uncertainty while requiring only a fit of the proportional-odds model.

Write \(J\) for the number of ordered response categories and \(K=J-1\) for the number of cumulative cutpoints. For observation \(i\), let \(Y_i\in\{1,\ldots,J\}\) be the response and \(x_i\in\mathbb{R}^p\) the predictor vector. We use the cumulative-logit convention \(\operatorname{logit}P(Y_i\leq j\mid x_i)=\alpha_j-x_i^\mathsf{T}\beta_j\). The proportional-odds null sets all the slope vectors equal.

Definition: The proportional-odds null hypothesis is \(H_0:\beta_1=\cdots=\beta_K\). The alternative allows a separate slope vector at each cutpoint. Its number of slope-equality restrictions is \(r=(K-1)p=(J-2)p\).

For a convenient parameterization of the alternative, take the last cutpoint's slope as the reference: set \(\beta_K=\beta\) and \(\beta_j=\beta+\delta_j\) for \(j=1,\ldots,K-1\). The nuisance parameter is \(\eta=(\alpha_1,\ldots,\alpha_K,\beta)\), and the tested parameter is \(\delta=(\delta_1^\mathsf{T},\ldots,\delta_{K-1}^\mathsf{T})^\mathsf{T}\). Thus \(H_0:\delta=0\). A different full-rank choice of slope contrasts gives the same test.

Category probabilities and scores

Let \(G(t)=1/(1+e^{-t})\), and define \(q_{ij}=P(Y_i\leq j\mid x_i)=G(\alpha_j-x_i^\mathsf{T}\beta_j)\), with \(q_{i0}=0\) and \(q_{iJ}=1\). The probability of category \(j\) is \(\pi_{ij}=q_{ij}-q_{i,j-1}\). Under the null, all \(\beta_j=\beta\); ordered thresholds \(\alpha_1<\cdots<\alpha_K\) give positive category probabilities.

For any parameter vector \(\theta\), write \(t_{ij}=\partial q_{ij}/\partial\theta\), with \(t_{i0}=t_{iJ}=0\). If \(Y_i=y\), its log-likelihood contribution is \(\log\pi_{iy}\), so its score is the derivative of the category log-probability.

Result: At a parameter value with positive category probabilities, the score contribution from observation \(i\), whose observed category is \(y\), is \(\displaystyle s_i(\theta)=\frac{t_{iy}-t_{i,y-1}}{\pi_{iy}}\). At the proportional-odds null, put \(h_{ij}=q_{ij}(1-q_{ij})\). The derivatives of \(q_{ij}\) with respect to \(\alpha_l\), \(\beta\), and \(\delta_l\) are, respectively, \(h_{ij}\mathbf{1}(j=l)\), \(-h_{ij}x_i\), and \(-h_{ij}x_i\mathbf{1}(j=l)\) for \(l<K\).

These derivatives follow from \(G'(t)=G(t)\{1-G(t)\}\) and the chain rule. The indicator in the \(\delta_l\) derivative appears because that deviation changes only the slope at cutpoint \(l\). Applying the displayed score formula to the resulting cumulative-probability derivatives gives each observation's threshold, common-slope, and slope-deviation scores.

Worked Example: One observation's score contributions

Consider three response categories, one predictor, and the null parameter values \(\alpha_1=0\), \(\alpha_2=\log 3\), and \(\beta=0\). Let the observed predictor value be \(x=2\). Then \(q_1=G(0)=0.5\), \(q_2=G(\log 3)=0.75\), and the category probabilities are \((\pi_1,\pi_2,\pi_3)=(0.5,0.25,0.25)\). Also, \(h_1=0.25\) and \(h_2=0.75(0.25)=0.1875\).

For an observation in category 2, the cumulative-probability derivatives with respect to \((\alpha_1,\alpha_2)\) are \((h_1,h_2)\) at the two cutpoints. Its threshold score is therefore \((-h_1/\pi_2,h_2/\pi_2)=(-1,0.75)\): the first cumulative probability is subtracted from the second in \(\pi_2\).

The common-slope derivative of \(\pi_2\) is \((-h_2x)-(-h_1x)=(h_1-h_2)x=0.125\). Dividing by \(\pi_2=0.25\) gives a common-slope score of \(0.5\). There is one deviation, \(\delta_1\), because there are two cutpoints. It changes \(q_1\) but not \(q_2\), so the derivative of \(\pi_2\) with respect to \(\delta_1\) is \(h_1x=0.5\), and its score is \(2\). This is one record's score, not a test statistic: the latter uses all records and adjusts for nuisance estimation.

Efficient information and the score statistic

Let \(\ell(\eta,\delta)\) be the full log-likelihood. At the restricted maximum-likelihood estimate \((\hat\eta_0,0)\), the nuisance score \(U_\eta=\partial\ell/\partial\eta\) is zero, while the score for the tested deviations is \(U_\delta=\partial\ell/\partial\delta\). Partition the Fisher information, or its regular observed-information counterpart, into blocks for \(\eta\) and \(\delta\). The information relevant to testing \(\delta=0\) after nuisance adjustment is the Schur complement \(I_{\delta\cdot\eta}=I_{\delta\delta}-I_{\delta\eta}I_{\eta\eta}^{-1}I_{\eta\delta}\).

Definition: The score statistic for proportional odds is \(\displaystyle S=U_\delta^\mathsf{T}I_{\delta\cdot\eta}^{-1}U_\delta\), where the score and information are evaluated under the null, ordinarily at its restricted maximum-likelihood estimate. The inverse requires nonsingular nuisance and efficient information matrices.

The adjustment can be understood directly from the score and information. For one observation with category probabilities \(\pi_{ij}\), its score has conditional mean zero because \(\sum_{y=1}^J\pi_{iy}s_{iy}=\sum_{y=1}^J(\partial\pi_{iy}/\partial\theta)=\partial(1)/\partial\theta=0\). Its expected information is consequently \(\sum_y\pi_{iy}s_{iy}s_{iy}^\mathsf{T}\). Summing over independent observations gives the information matrix used in the test. This formula also makes clear that the information averages over possible categories; it is not just the outer product for the observed category.

Derivation: The adjusted information is the information left for slope deviations after their linear association with the nuisance scores is removed. Indeed, for a vector \(a\), the residualized score is \(U_\delta-I_{\delta\eta}I_{\eta\eta}^{-1}U_\eta\). Its covariance is \[ I_{\delta\delta}-I_{\delta\eta}I_{\eta\eta}^{-1}I_{\eta\delta}=I_{\delta\cdot\eta}. \] This follows by expanding its covariance and using \(\operatorname{Cov}(U_\eta)=I_{\eta\eta}\), \(\operatorname{Cov}(U_\delta,U_\eta)=I_{\delta\eta}\), and \(\operatorname{Cov}(U_\eta,U_\delta)=I_{\eta\delta}\). Thus nuisance information is projected out rather than treated as known.

Worked Example: A three-category score-test calculation

Use three categories, one predictor, and the null values \(\alpha_1=0\), \(\alpha_2=\log 3\), \(\beta=0\). Consider four design points and responses: two records with \(x=0,Y=1\), one with \(x=-3,Y=2\), and one with \(x=1,Y=3\). At these null values the category probabilities are \((0.5,0.25,0.25)\). For categories 1, 2, and 3, respectively, the score contributions for \((\alpha_1,\alpha_2,\beta,\delta_1)\) are \[ (0.5,0,-0.5x,-0.5x),\quad (-1,0.75,0.25x,x),\quad (0,-0.75,0.75x,0). \]

Summing the observed scores gives \(U_{\alpha_1}=2(0.5)-1=0\), \(U_{\alpha_2}=0.75-0.75=0\), \(U_\beta=0.25(-3)+0.75(1)=0\), and \(U_{\delta_1}=-3\). The nuisance score is zero at this point; in this illustration, take the point to be the restricted fit.

For one record at predictor value \(x\), the expected information entries for the two thresholds are \(I_{\alpha\alpha}=\begin{pmatrix}0.375&-0.1875\\-0.1875&0.28125\end{pmatrix}\). The threshold-common-slope entries are \((-0.1875x,-0.09375x)^\mathsf{T}\), the threshold-deviation entries are \((-0.375x,0.1875x)^\mathsf{T}\), and the remaining entries are \(I_{\beta\beta}=0.28125x^2\), \(I_{\beta\delta}=0.1875x^2\), and \(I_{\delta\delta}=0.375x^2\). These values are obtained by summing \(\pi_y s_y s_y^\mathsf{T}\) over the three categories.

Across the four records, \(n=4\), \(\sum_i x_i=-2\), and \(\sum_i x_i^2=10\). In the order \(\eta=(\alpha_1,\alpha_2,\beta)\), this gives \[ I_{\eta\eta}=\frac1{16}\begin{pmatrix}24&-12&6\\-12&18&3\\6&3&45\end{pmatrix},\quad I_{\eta\delta}=\frac1{16}\begin{pmatrix}12\\-6\\30\end{pmatrix},\quad I_{\delta\delta}=\frac{60}{16}. \] To compute the nuisance adjustment, solve \(\begin{pmatrix}24&-12&6\\-12&18&3\\6&3&45\end{pmatrix}z=(12,-6,30)^\mathsf{T}\). Direct substitution gives \(z=(1/6,-1/3,2/3)^\mathsf{T}\). Hence \(I_{\delta\eta}I_{\eta\eta}^{-1}I_{\eta\delta}=(12,-6,30)z/16=24/16=1.5\), and \(I_{\delta\cdot\eta}=60/16-24/16=36/16=2.25\). The statistic is \(S=(-3)^2/2.25=4\), compared with a chi-square distribution on one degree of freedom. Its upper-tail probability is approximately \(0.0455\). In an applied analysis, the nuisance fit and information must be calculated from the data; this example shows the adjustment arithmetic.

Why the reference distribution is chi-square

The chi-square calibration is asymptotic. It requires more than the algebraic definition of the statistic: if the null fit is not identifiable or the probabilities approach zero in a way that breaks regularity, the usual reference distribution need not apply.

Theorem: Suppose observations are independent, the number of categories and predictors is fixed, and the true parameter satisfies \(H_0\) at an interior point with strictly ordered thresholds and positive category probabilities. Assume the model is identifiable, the Fisher information per observation converges to a positive-definite matrix, the score satisfies a central limit theorem (for example, through bounded covariates and a Lindeberg condition), and the restricted maximum-likelihood estimator is consistent and asymptotically normal. Then, under \(H_0\), \(S\xrightarrow{d}\chi^2_r\), where \(r=(K-1)p\).

To see the role of nuisance estimation, let \((\eta_0,0)\) be the true parameter and write \(I\) for the limiting information per observation at that parameter. The score central limit theorem gives \(n^{-1/2}(U_\eta,U_\delta)\xrightarrow{d}N(0,I)\). A first-order expansion of the nuisance score equation \(U_\eta(\hat\eta_0,0)=0\), together with consistency and asymptotic normality, yields \[ \sqrt n(\hat\eta_0-\eta_0)=I_{\eta\eta}^{-1}\frac{U_\eta(\eta_0,0)}{\sqrt n}+o_p(1). \]

Expand the tested score at the restricted estimate. Its derivative with respect to \(\eta\) is asymptotically \(-nI_{\delta\eta}\), so \[ \frac{U_\delta(\hat\eta_0,0)}{\sqrt n} = \frac{U_\delta(\eta_0,0)}{\sqrt n} -I_{\delta\eta}I_{\eta\eta}^{-1}\frac{U_\eta(\eta_0,0)}{\sqrt n} +o_p(1). \] The limiting covariance of this adjusted score is \(I_{\delta\cdot\eta}\), by the covariance calculation above. The plug-in efficient information divided by \(n\) converges in probability to \(I_{\delta\cdot\eta}\). Therefore the quadratic form converges to the sum of squares of \(r\) independent standard normal variables, which is \(\chi^2_r\). \(\square\)

Worked Example: Counting restrictions with several predictors

Suppose the ordinal response has four categories, so there are \(K=3\) cumulative cutpoints, and there are two predictors. The alternative has three slope vectors, each of length two. Under the null they are equal, leaving \(K-1=2\) independent slope deviations, each with two components. The test therefore has \(r=2(2)=4\) degrees of freedom when the theorem's regularity conditions hold.

For a score-contribution check, let \(x=(1,-1)^\mathsf{T}\), \(\beta=0\), and thresholds \((0,\log 3,\log 7)\). Then the cumulative probabilities are \((0.5,0.75,0.875)\), and the category probabilities are \((0.5,0.25,0.125,0.125)\). For an observation in category 2, the threshold score is \((-1,0.75,0)\). The common-slope score is \(0.25x=(0.25,-0.25)^\mathsf{T}\). There are two deviations, \(\delta_1,\delta_2\): their score contributions are \(x=(1,-1)^\mathsf{T}\) and \(-0.75x=(-0.75,0.75)^\mathsf{T}\), respectively. For five categories instead, \(K=4\), and the same two predictors give \(r=(4-1)2=6\) restrictions.

Common mistake and practical interpretation

Ignoring nuisance estimation. Using \(U_\delta^\mathsf{T}I_{\delta\delta}^{-1}U_\delta\) treats thresholds and the common slope as known. The correction is to use \(I_{\delta\cdot\eta}\), which removes the information shared with nuisance parameters.
Using the wrong degrees of freedom. The test concerns slope equality, not the number of thresholds alone. With \(p\) predictors and \(K\) cutpoints, the restriction count is \((K-1)p\).
Evaluating at an unrestricted fit. The score test is evaluated under the null, at the restricted fit. An unrestricted estimate is not a substitute for that fit.
Assuming the asymptotic calibration is automatic. Sparse categories, near-zero fitted probabilities, or weak identification can undermine the regularity conditions behind the chi-square approximation.

The score test is useful when fitting the proportional-odds model is straightforward but fitting the larger alternative is difficult or unnecessary. A large statistic indicates that, after accounting for uncertainty in thresholds and common slopes, the data show evidence against equal slopes across cutpoints. It does not by itself identify which predictor or cutpoint is responsible.

Check Your Understanding

Work through these questions using the parameterization and assumptions above.

For \(J\) response categories and \(p\) predictors, how many slope-equality restrictions does the proportional-odds null impose?
Why does the expected information for one observation sum over all possible categories, rather than use only its observed category?
What does the Schur-complement adjustment remove from the slope-deviation information?
For five categories and three predictors, what are the score test's asymptotic degrees of freedom?
Name two regularity conditions needed for the chi-square reference distribution, and explain why a zero fitted category probability may be problematic.