Score Tests For Adding Multinomial Predictors

Why Test at the Restricted Fit?

Suppose a baseline-category multinomial logit model already includes a set of predictors, and a scientific question asks whether another set should be added. A score test answers this question by measuring how strongly the likelihood is pushing away from the restricted model, while evaluating the score and information at the fit under the null hypothesis. Unlike a likelihood-ratio test, it does not require fitting the expanded model.

The key is to account for the coefficients already in the model. The added coefficients can be statistically entangled with those existing nuisance parameters. Consequently, the information relevant to the test is not generally just the information block for the added coefficients.

Model, Hypothesis, and Score Statistic

Let $Y_i$ be the response for independent observational units $i=1,\ldots,n$, taking one of $K$ categories. Use category $K$ as the baseline, and write $p_{ik}=P(Y_i=k)$. The baseline-category logits have the form $ \log(p_{ik}/p_{iK})=\eta_{ik} $ for $k=1,\ldots,K-1$. Let $\gamma$ collect the coefficients in the restricted model and let $\beta$ collect the $q$ added, category-specific coefficients. The null hypothesis is $H_0:\beta=0$; the alternative allows $\beta$ to vary freely.

Definition: The score test statistic for $H_0:\beta=0$ is calculated from the score for $\beta$ and the information for $\beta$ adjusted for estimation of $\gamma$, both evaluated at the restricted maximum-likelihood estimate $(\tilde\gamma,0)$. Write the Fisher information there in blocks as $$ I= \begin{pmatrix} I_{\gamma\gamma}&I_{\gamma\beta}\\ I_{\beta\gamma}&I_{\beta\beta} \end{pmatrix}. $$ The adjusted information and score statistic are $$ J_{\beta\beta\cdot\gamma} =I_{\beta\beta}-I_{\beta\gamma}I_{\gamma\gamma}^{-1}I_{\gamma\beta}, \qquad Q=U_\beta(\tilde\gamma,0)^\mathsf{T} J_{\beta\beta\cdot\gamma}^{-1} U_\beta(\tilde\gamma,0). $$

Here $U_\beta$ is the derivative of the log-likelihood with respect to $\beta$. The restricted maximum-likelihood estimate satisfies $U_\gamma(\tilde\gamma,0)=0$, but the added-coefficient score $U_\beta(\tilde\gamma,0)$ need not be zero. The matrix $J_{\beta\beta\cdot\gamma}$ is the efficient information for $\beta$ after accounting for the nuisance parameters.

For example, if $z_i$ is an added predictor and $\beta_{kr}$ is its coefficient for nonbaseline category $k$ and predictor column $r$, the score component is $$ U_{\beta_{kr}}=\sum_{i=1}^n z_{ir}\{\mathbf{1}(Y_i=k)-p_{ik}\}, \qquad k=1,\ldots,K-1. $$ For grouped multinomial observations, the indicator contribution is replaced by the observed category count, and its expectation by the group size times $p_{ik}$. These scores and the information blocks are evaluated using the probabilities from the restricted fit.

Worked Example: One Added Predictor and Two Category Coefficients

There are two groups, each with 40 observations, with predictor $x=0$ or $x=1$. The observed counts in categories 1, 2, and baseline category 3 are:

Predictor $x$	Category 1	Category 2	Category 3	Total
0	20	10	10	40
1	10	20	10	40

The restricted model has intercepts only. Its pooled proportions are $30/80=0.375$, $30/80=0.375$, and $20/80=0.25$. Thus the restricted fitted probabilities are $p_1=p_2=0.375$ and $p_3=0.25$. The intercept score is zero because the pooled observed counts equal their fitted totals.

The added slope scores, one for each nonbaseline category, come from the $x=1$ group: $U_\beta=(10-40(0.375),\,20-40(0.375))^\mathsf{T}=(-5,5)^\mathsf{T}$. For a multinomial observation, the covariance matrix for the two nonbaseline indicators is $$ C=\begin{pmatrix} p_1(1-p_1)&-p_1p_2\\ -p_1p_2&p_2(1-p_2) \end{pmatrix} =\begin{pmatrix} 0.234375&-0.140625\\ -0.140625&0.234375 \end{pmatrix}. $$

Because the intercepts are nuisance parameters, the adjusted predictor information is $C$ times the weighted sum of squared deviations of $x$ from its weighted mean. The mean is $0.5$, so that sum is $40(0-0.5)^2+40(1-0.5)^2=20$. Therefore $J_{\beta\beta\cdot\gamma}=20C$, with diagonal entries $4.6875$ and off-diagonal entries $-2.8125$. The vector $(-5,5)^\mathsf{T}$ lies in the direction with information eigenvalue $7.5$; its squared length is $50$. Hence $$ Q=\frac{50}{7.5}=\frac{20}{3}\approx 6.667. $$

This is a joint test of the two added slopes, so the reference distribution has 2 degrees of freedom. Notice that using the unadjusted slope information $40C$ would give $50/15=10/3$, not $20/3$. The restricted intercepts must be accounted for.

Why the Adjusted Information Appears

A finite change in an added coefficient can be partly mimicked by changing existing coefficients. The subtraction in $J_{\beta\beta\cdot\gamma}$ removes the information associated with that overlap. The following block-matrix fact also gives a convenient way to calculate the statistic.

Result: If $I_{\gamma\gamma}$ and $J_{\beta\beta\cdot\gamma}$ are nonsingular, then the $\beta,\beta$ block of $I^{-1}$ is $J_{\beta\beta\cdot\gamma}^{-1}$. Thus the statistic can also be computed as $U_\beta^\mathsf{T}[I^{-1}]_{\beta\beta}U_\beta$, using the full information matrix at the restricted fit.

To verify the identity, solve $I(a^\mathsf{T},b^\mathsf{T})^\mathsf{T}=(0^\mathsf{T},v^\mathsf{T})^\mathsf{T}$. The first block equation gives $a=-I_{\gamma\gamma}^{-1}I_{\gamma\beta}b$. Substitution into the second gives $J_{\beta\beta\cdot\gamma}b=v$, and hence $b=J_{\beta\beta\cdot\gamma}^{-1}v$. This is precisely the lower-right block action of $I^{-1}$, proving the claim. $\square$

In the multinomial logit model, information blocks can be assembled from the category-indicator covariance. For an existing design column $x_r$ and an added column $z_s$, the block entry for categories $k,l<K$ is $\sum_i x_{ir}z_{is}p_{ik}\{\mathbf{1}(k=l)-p_{il}\}$. Replacing $x_r$ by $z_r$ gives added-coefficient information. This formula makes explicit that correlations between existing and added design columns contribute to the adjustment.

Worked Example: Testing Two Added Predictors Jointly

Consider three groups, each containing three observations, with predictor values $(x_1,x_2)=(-1,1),(0,-2),(1,1)$. The observed counts in categories 1, 2, and 3 are, respectively, $(0,1,2),(1,1,1),(2,0,1)$. Test the four category-specific slopes for $x_1$ and $x_2$ jointly, retaining only category intercepts in the restricted model.

Each category has three observations overall, so the restricted fit has $p_1=p_2=p_3=1/3$. The covariance matrix for the two nonbaseline indicators is $C=\frac{1}{9}\begin{pmatrix}2&-1\\-1&2\end{pmatrix}$. For $x_1$, the score vector across categories 1 and 2 is $(2,-1)^\mathsf{T}$: its first component is $(-1)(0-1)+(1)(2-1)=2$, and its second is $(-1)(1-1)+(1)(0-1)=-1$. For $x_2$, the score vector is $(2,-1)^\mathsf{T}$: the corresponding sums are $(-1)(0-1)+(-2)(1-1)+(1)(2-1)=2$ and $(1)(1-1)+(-2)(1-1)+(1)(0-1)=-1$.

Therefore the $x_2$ score is $(2,-1)^\mathsf{T}$. The weighted predictor sums of squares are $6$ for $x_1$ and $18$ for $x_2$, and their weighted cross-product is zero. The information separates by predictor: it is $6C$ for $x_1$ and $18C$ for $x_2$. Inverting these blocks gives $(6C)^{-1}=\begin{pmatrix}1&0.5\\0.5&1\end{pmatrix}$ and $(18C)^{-1}=\begin{pmatrix}1/3&1/6\\1/6&1/3\end{pmatrix}$.

The $x_1$ contribution to the statistic is $ (2,-1)(6C)^{-1}(2,-1)^\mathsf{T}=4+1-2=3$. The $x_2$ contribution is $ (2,-1)(18C)^{-1}(2,-1)^\mathsf{T}=1$. Thus $Q=4$, with 4 degrees of freedom because four coefficients are tested. The explicit score calculation is important: each residual must be multiplied by its predictor value before summation.

Asymptotic Reference Distribution

Theorem: Suppose the observations are independent and follow the stated multinomial logit model; the null parameter $(\gamma_0,0)$ is an interior point; the model is identifiable near that point; the log-likelihood is sufficiently smooth for a second-order Taylor expansion; the information identity and central limit theorem apply; and the per-observation Fisher information converges to a finite positive-definite matrix. Assume also that the restricted maximum-likelihood estimator is consistent and asymptotically normal, and that $q$ is fixed. Then, under $H_0$, $Q\xrightarrow{d}\chi^2_q$.

Here is the justification. Write the information per observation at the null as $i=I/n$, with blocks $i_{\gamma\gamma},i_{\gamma\beta},i_{\beta\beta}$. The central limit theorem and the information identity give $n^{-1/2}(U_\gamma,U_\beta)^\mathsf{T}\xrightarrow{d}N(0,i)$, where the scores on the right are evaluated at the true null parameter. A Taylor expansion of the restricted score equation $U_\gamma(\tilde\gamma,0)=0$ yields $$ \sqrt n(\tilde\gamma-\gamma_0) =i_{\gamma\gamma}^{-1}\frac{U_\gamma(\gamma_0,0)}{\sqrt n}+o_p(1). $$

Expand the added-coefficient score around $(\gamma_0,0)$. Its derivative with respect to $\gamma$, divided by $n$, converges to $-i_{\beta\gamma}$. Substituting the restricted-estimator expansion gives $$ \frac{U_\beta(\tilde\gamma,0)}{\sqrt n} =\frac{U_\beta(\gamma_0,0)}{\sqrt n} -i_{\beta\gamma}i_{\gamma\gamma}^{-1} \frac{U_\gamma(\gamma_0,0)}{\sqrt n}+o_p(1). $$ This converges to a mean-zero normal vector with covariance $j=i_{\beta\beta}-i_{\beta\gamma}i_{\gamma\gamma}^{-1}i_{\gamma\beta}$. The restricted-fit adjusted information divided by $n$ converges in probability to $j$. Therefore $Q$ converges in distribution to $Z^\mathsf{T}j^{-1}Z$, where $Z\sim N(0,j)$. Since $j$ is positive definite, $j^{-1/2}Z\sim N(0,I_q)$, and the sum of the squares of its $q$ components has the $\chi^2_q$ distribution. This proves the result. $\square$

Worked Example: A Joint Test With Unequal Category Probabilities

Suppose the restricted model has intercepts only and two groups of 20 observations, with $x=0$ and $x=1$. Their counts in categories 1, 2, and baseline category 3 are $(12,5,3)$ and $(8,7,5)$. The pooled counts are $(20,12,8)$, so the restricted probabilities are $(0.5,0.3,0.2)$.

For the added $x$ slopes, the score is the residual in the $x=1$ group: $U_\beta=(8-20(0.5),\,7-20(0.3))^\mathsf{T}=(-2,1)^\mathsf{T}$. The indicator covariance is $C=\begin{pmatrix}0.25&-0.15\\-0.15&0.21\end{pmatrix}$. The weighted mean of $x$ is $0.5$, and the weighted sum of squared deviations is $20(0.5)^2+20(0.5)^2=10$. Consequently, $J_{\beta\beta\cdot\gamma}=10C=\begin{pmatrix}2.5&-1.5\\-1.5&2.1\end{pmatrix}$. Its determinant is $2.5(2.1)-1.5^2=3$, so its inverse is $\begin{pmatrix}0.7&0.5\\0.5&5/6\end{pmatrix}$.

The score statistic is $Q=4(0.7)+2(-2)(1)(0.5)+5/6=49/30\approx1.633$. There are two added coefficients, so the reference distribution is $\chi^2_2$. Its upper-tail probability is $\exp(-Q/2)=\exp(-49/60)\approx0.442$. These data therefore provide little evidence against the restricted model at conventional significance levels.

Common Pitfall and Practical Interpretation

A frequent error is to calculate $U_\beta^\mathsf{T}I_{\beta\beta}^{-1}U_\beta$ and call it the score test. That expression treats the nuisance parameters as though they were known. In a fitted restricted model they are estimated, so the correct information is the Schur complement $I_{\beta\beta}-I_{\beta\gamma}I_{\gamma\gamma}^{-1}I_{\gamma\beta}$. The first worked example shows that the unadjusted calculation can change the statistic substantially.

The degrees of freedom equal the number of independently estimable added coefficients being tested, not the number of response categories alone. For $K$ categories and $r$ added predictor columns, testing all category-specific slopes ordinarily tests $q=(K-1)r$ coefficients under the baseline-category parameterization. The score test is an asymptotic test, not an exact finite-sample calculation; sparse counts, separation, or nearly singular information can make its chi-square approximation unreliable.

Key takeaway: Fit the restricted model, evaluate the added-coefficient score there, and standardize it using information adjusted for the nuisance parameters. Under the stated regularity conditions, the resulting quadratic statistic has an asymptotic $\chi^2_q$ null distribution.

Check Your Understanding

Use the definitions and calculations above to answer these questions.

In a model with three response categories and two added predictor columns, how many coefficients are tested if every nonbaseline category receives both slopes?
Why is the added-coefficient score generally evaluated at the restricted fit rather than at an expanded-model estimate?
What component of the statistic changes when nuisance parameters are estimated rather than treated as known?
For the first worked example, explain in words why $40C$ is not the appropriate adjusted information matrix.
State two regularity conditions needed for the chi-square reference distribution to be justified.