Uncertainty For Predicted Category Probabilities

Why Predicted Probabilities Need Standard Errors

A fitted multinomial logit model can report a probability for every response category at a chosen covariate value. Those probabilities are estimates: the fitted coefficients would generally change if a different sample were observed. To quantify that parameter uncertainty on the probability scale, we propagate the estimated coefficient covariance through the softmax function.

The calculation has two ingredients: the covariance matrix of the fitted coefficients and the derivative of the predicted probability with respect to those coefficients. The multivariate delta method combines them. The resulting standard error describes uncertainty in the estimated conditional probability at the specified covariate value; it does not describe the variability of an individual’s future categorical outcome.

The Softmax Derivative for a Predicted Category

Consider a baseline-category logit model with $K$ response categories, where category $K$ is the baseline. Let $\beta_j$ be the coefficient vector for category $j$, for $j=1,\ldots,K-1$, and let $z_j$ be its design vector at a fixed covariate value. The linear predictors and probabilities are

\eta_j=z_j^\top\beta_j,\qquad p_j=\frac{\exp(\eta_j)}{1+\sum_{m=1}^{K-1}\exp(\eta_m)} \quad (j<K),\qquad p_K=\frac{1}{1+\sum_{m=1}^{K-1}\exp(\eta_m)}.

Stack the coefficient vectors in the same order as the model’s fitted covariance matrix, writing $\beta=(\beta_1^\top,\ldots,\beta_{K-1}^\top)^\top$. For a nonbaseline category $r$, differentiate $p_r$ first with respect to each linear predictor. The softmax derivative is $\partial p_r/\partial\eta_l=p_r(\mathbf{1}\{r=l\}-p_l)$. The chain rule then gives the gradient with respect to the full coefficient vector.

Definition: At a fixed covariate value, the coefficient-gradient block for the predicted probability $p_r$, $r<K$, is $$ g_{r,l}=\frac{\partial p_r}{\partial\beta_l} =p_r\bigl(\mathbf{1}\{r=l\}-p_l\bigr)z_l^\top,\qquad l=1,\ldots,K-1. $$ The full gradient $g_r=\nabla_\beta p_r$ is the column vector formed by stacking the transposes of these row blocks. For the baseline probability, $$ g_{K,l}=\frac{\partial p_K}{\partial\beta_l}=-p_Kp_lz_l^\top,\qquad l=1,\ldots,K-1. $$

The design vector may differ across categories when the model uses category-specific predictors. If the same covariate design is used for each nonbaseline category, the blocks share that design vector, but the probability factors still differ by category. These derivatives are evaluated at the fitted coefficients when calculating estimated standard errors.

Delta-Method Result and Plug-In Standard Error

The following result states the regularity conditions needed for both the distributional approximation and the usual estimated standard error. The continuity condition on the gradient is relevant to consistency of the plug-in variance; differentiability at the true coefficient alone is enough for the distributional delta-method conclusion.

Theorem: Suppose $\beta_0\in\mathbb{R}^q$, $h:\mathbb{R}^q\to\mathbb{R}$ is differentiable at $\beta_0$, and $$ \sqrt n(\hat\beta-\beta_0)\xrightarrow{d}N(0,\Sigma), $$ where $\Sigma$ is a finite positive semidefinite matrix. Then $$ \sqrt n\{h(\hat\beta)-h(\beta_0)\}\xrightarrow{d} N\bigl(0,\nabla h(\beta_0)^\top\Sigma\nabla h(\beta_0)\bigr). $$ If, in addition, $\nabla h$ is continuous at $\beta_0$, and an estimated coefficient covariance matrix satisfies $n\hat V\xrightarrow{p}\Sigma$, then $$ \hat g^\top\hat V\hat g $$ is a consistent estimator of the first-order variance of $h(\hat\beta)$, where $\hat g=\nabla h(\hat\beta)$. In particular, $\sqrt{\hat g^\top\hat V\hat g}$ estimates the standard error of $h(\hat\beta)$.

Proof. Differentiability at $\beta_0$ gives a remainder $r(u)$ such that

h(\beta_0+u)-h(\beta_0) =\nabla h(\beta_0)^\top u+r(u), \qquad \frac{r(u)}{\lVert u\rVert}\longrightarrow 0 \quad\text{as }u\longrightarrow 0.

The assumed convergence in distribution implies $\sqrt n(\hat\beta-\beta_0)=O_p(1)$, hence $\hat\beta-\beta_0=O_p(n^{-1/2})$ and $\hat\beta\xrightarrow{p}\beta_0$. Therefore $\sqrt n\,r(\hat\beta-\beta_0)=o_p(1)$: the remainder divided by $\lVert\hat\beta-\beta_0\rVert$ converges in probability to zero, while $\sqrt n\lVert\hat\beta-\beta_0\rVert=O_p(1)$. Multiplying the expansion by $\sqrt n$ and applying the assumed limit and Slutsky’s theorem proves the distributional statement.

For the plug-in claim, continuity of $\nabla h$ at $\beta_0$ and $\hat\beta\xrightarrow{p}\beta_0$ imply $\hat g\xrightarrow{p}g_0=\nabla h(\beta_0)$. Consequently,

n\hat g^\top\hat V\hat g =\hat g^\top(n\hat V)\hat g \xrightarrow{p}g_0^\top\Sigma g_0,

by continuity of matrix multiplication. This is the stated first-order variance limit. $\square$

For a fitted category probability, take $h(\beta)=p_r(\beta)$ and substitute its softmax gradient $g_r$. If $\hat V$ is the estimated covariance matrix of $\hat\beta$, the estimated variance and standard error are

\widehat{\operatorname{Var}}(\hat p_r)=\hat g_r^\top\hat V\hat g_r, \qquad \widehat{\operatorname{SE}}(\hat p_r) =\sqrt{\hat g_r^\top\hat V\hat g_r}.

For a vector of probabilities, stack their gradients as the rows of a Jacobian matrix $J$. The estimated covariance matrix of the probability vector is $\hat J\hat V\hat J^\top$. Its diagonal entries are the individual estimated variances; its off-diagonal entries retain the dependence between category probabilities.

Worked Probability-Scale Calculations

Worked Example: Three Categories with Two Intercepts

Suppose there are three response categories, category 3 is the baseline, and the model has an intercept for each nonbaseline category. At the prediction point, let $\hat\eta_1=0$ and $\hat\eta_2=\log 2$. The softmax denominator is $1+1+2=4$, so $\hat p_1=1/4$, $\hat p_2=1/2$, and $\hat p_3=1/4$.

For category 1, the gradient with respect to the two intercepts is $\hat g_1=(\hat p_1(1-\hat p_1),-\hat p_1\hat p_2)^\top=(3/16,-1/8)^\top$. Suppose the fitted coefficient covariance matrix, in this same coefficient order, is $\hat V=\left(\begin{smallmatrix}0.04&0\\0&0.09\end{smallmatrix}\right)$. Then

\widehat{\operatorname{Var}}(\hat p_1) =\left(\frac{3}{16}\right)^2(0.04) +\left(-\frac18\right)^2(0.09) =0.00140625+0.00140625=0.0028125,

so the estimated standard error is $\sqrt{0.0028125}\approx0.0530$. A 95% probability-scale Wald interval uses the standard-normal critical value $1.96$: $0.25\pm1.96(0.0530)\approx0.25\pm0.1039$, or approximately $(0.1461,0.3539)$. This interval is an approximation based on the fitted probability’s sampling distribution.

Worked Example: Checking the Probability-Scale Interval Arithmetic

Suppose a three-category fitted model gives probabilities $\hat p_1=4/9$, $\hat p_2=1/3$, and $\hat p_3=2/9$, with category 3 as baseline. The corresponding two linear predictors are $\hat\eta_1=\log 2$ and $\hat\eta_2=\log(3/2)$, because the category-to-baseline probability ratios are $2$ and $3/2$.

Let the estimated covariance matrix of the two intercept estimates be $0.0444I_2$, where $I_2$ is the two-dimensional identity matrix. The probability gradient is $\hat g_1=(20/81,-4/27)^\top$. Its squared length is $(20/81)^2+(4/27)^2=(400+144)/6561=544/6561$. Thus

\widehat{\operatorname{Var}}(\hat p_1) =0.0444\left(\frac{544}{6561}\right) =\frac{24.1536}{6561} \approx0.0036814,\qquad \widehat{\operatorname{SE}}(\hat p_1)\approx0.0607.

Using the displayed rounded standard error, the margin is $1.96(0.0607)=0.118972$. Since $4/9=0.444444\ldots$, the endpoints are $0.444444\ldots-0.118972=0.325472\ldots$ and $0.444444\ldots+0.118972=0.563416\ldots$. The 95% Wald interval is therefore approximately $(0.3255,0.5634)$. Keeping the unrounded standard error through the calculation gives the same four-decimal endpoints.

Intervals on the Logit Scale

A probability-scale Wald interval is simple to calculate, but its endpoints are not constrained to lie between zero and one. A common alternative is to form an interval for the logit of the category probability and then transform the endpoints back. For a probability $p_r$, define $h(\beta)=\log\{p_r(\beta)/(1-p_r(\beta))\}$. The chain rule gives

\nabla h(\beta)=\frac{\nabla p_r(\beta)}{p_r(\beta)\{1-p_r(\beta)\}}, \qquad \widehat{\operatorname{SE}}\{h(\hat\beta)\} =\frac{\widehat{\operatorname{SE}}(\hat p_r)} {\hat p_r(1-\hat p_r)}.

This is the log-odds of category $r$ versus all the other categories combined, not generally one of the model’s individual baseline-category logits. A normal interval for $h(\hat\beta)$ is transformed with the inverse-logit function $s(t)=1/(1+\exp(-t))$. Because $s(t)$ lies strictly between zero and one for every finite $t$, the resulting interval respects the probability bounds.

Worked Example: Transforming a Logit-Scale Interval

Continue with an estimated probability $\hat p_1=4/9$, whose estimated probability-scale standard error is $0.0607$. The estimated logit is $\log\{(4/9)/(5/9)\}=\log(4/5)\approx-0.2231436$. Since $\hat p_1(1-\hat p_1)=(4/9)(5/9)=20/81$, the transformed-scale standard error is approximately $0.0607/(20/81)=0.245835$, which rounds to $0.2458$.

Using the displayed rounded values, the 95% logit interval has margin $1.96(0.2458)=0.481768$. Its endpoints are $-0.2231-0.481768=-0.704868$ and $-0.2231+0.481768=0.258668$, or $(-0.7049,0.2587)$ to four decimal places. Applying the inverse logit gives approximately $(0.3308,0.5643)$. Rounding both scales consistently matters: the interval is obtained by transforming the two endpoints, not by transforming the center and adding a probability-scale margin.

Assumptions, Interpretation, and Common Pitfalls

These intervals are large-sample approximations. Their usual interpretation is that, under the stated model and regularity conditions, the procedure has approximately 95% coverage over repeated samples when the sample size is sufficiently large. The covariance matrix must correspond to the fitted coefficient vector in the same ordering used to build the gradient. If a different covariance estimator is used, such as a robust estimator, the same propagation formula applies provided its asymptotic covariance conditions are appropriate.

A predicted probability is a function of estimated coefficients at a specified covariate value. It is not the outcome of a future categorical draw. For example, if a future response has probability $0.4$ of category 1, the outcome itself is still categorical; the delta-method standard error instead quantifies how much the estimated value $0.4$ would vary across repeated data sets because the coefficients are estimated.

Common mistake: treating category probabilities as independent. They sum to one, so their estimation errors are dependent. For joint uncertainty, use $\hat J\hat V\hat J^\top$, including its off-diagonal entries; do not estimate each category’s uncertainty as if the other probabilities were unrelated.
Common mistake: using the wrong gradient dimension or coefficient order. Each gradient block must match the corresponding coefficient block in $\hat V$. Check the ordering before multiplying; otherwise the quadratic form can be numerically valid but refer to the wrong quantity.
Common mistake: assuming every Wald interval is a valid probability interval. The ordinary interval $\hat p_r\pm1.96\,\widehat{\operatorname{SE}}(\hat p_r)$ can extend below zero or above one. A logit-scale interval avoids that particular defect, though it remains a large-sample approximation and can behave poorly when estimates are very near a boundary.
Common mistake: confusing coefficient standard errors with probability standard errors. The probability-scale standard error requires the softmax derivative and the full covariance matrix. A coefficient’s standard error alone does not give uncertainty for a nonlinear predicted probability.

Key takeaway: At a fixed covariate value, differentiate the softmax probability with respect to every fitted coefficient, then calculate $\hat g_r^\top\hat V\hat g_r$. The square root is the delta-method standard error; an interval must also be chosen and interpreted with its approximation and probability-bound limitations in mind.

Check Your Understanding

Work through these questions using the definitions and calculations above.

In a baseline-category model, why does the derivative of a nonbaseline probability with respect to a different category’s linear predictor have a negative sign?
What information is needed, in addition to a fitted probability, to calculate its delta-method standard error?
For a vector of predicted probabilities, what do the off-diagonal entries of $\hat J\hat V\hat J^\top$ represent?
Explain why differentiability at $\beta_0$ is enough for the delta-method distributional result, while consistency of the plug-in variance estimate was stated with the additional gradient-continuity condition.
Why is a future categorical outcome’s variability not the same as the delta-method standard error of its predicted category probability?