Parameter Uncertainty in Landmark State Predictions

Why Predictions Need Uncertainty Measures

A predicted state-occupation probability is an estimate, not a known quantity. Even when a model gives a precise-looking prediction such as a 30% probability of occupying a particular state at a future horizon, the transition probabilities used to calculate it have been estimated from data. Different samples could produce different transition estimates and therefore different predictions.

At a landmark, a prediction is formed by combining the state occupied at that landmark with one or more transition matrices. The central task here is to quantify how estimation uncertainty in those transition parameters affects the resulting state probabilities. The landmark time, prediction horizon, and conditional interpretation follow the earlier material on landmarking and dynamic prediction. We focus on the uncertainty calculation rather than on constructing the prediction itself.

Definition: Let \(\theta\) be a vector of parameters that determine the transition probabilities used for a landmark prediction, and let \(\hat{\theta}\) be an estimate of that vector. Write the predicted state-occupation probability vector at horizon \(u\) as \(\hat{\boldsymbol p}(u)=g(\hat{\theta})\), where \(g\) includes the transition-matrix calculations for the chosen landmark and horizon. Parameter uncertainty is the sampling uncertainty in \(\hat{\boldsymbol p}(u)\) induced by estimating \(\theta\).

Covariance Propagation by the Delta Method

Suppose the estimated parameter vector has estimated covariance matrix \(\widehat V\). If \(g\) is differentiable near the true parameter value, its local sensitivity is described by the Jacobian matrix \(D\). The row for a particular predicted state gives the derivatives of that state probability with respect to the parameters. A first-order Taylor expansion gives \(\,g(\hat{\theta})\approx g(\theta)+D(\hat{\theta}-\theta)\). Consequently, the approximate covariance matrix of the predicted probabilities is obtained by carrying the parameter covariance through that linear transformation.

\widehat{\operatorname{Cov}}\{\hat{\boldsymbol p}(u)\} \approx D\,\widehat V\,D^{\mathsf T}, \qquad D=\left.\frac{\partial g(\theta)}{\partial\theta^{\mathsf T}}\right|_{\theta=\hat{\theta}}.

The diagonal entries of this matrix estimate the variances of the predicted probabilities; their square roots are standard errors. Its off-diagonal entries describe how errors in different state probabilities tend to move together. For one predicted component \(g_j\), the same calculation is \(\widehat{\operatorname{Var}}\{g_j(\hat{\theta})\}\approx \nabla g_j(\hat{\theta})^{\mathsf T}\widehat V\nabla g_j(\hat{\theta})\). This expression includes covariance between parameter estimates, which must not be discarded unless independence is justified.

Assumptions for the usual approximation: The parameter estimator is approximately unbiased and jointly normal in large samples, with a consistent covariance estimate; the prediction function is differentiable in a neighborhood of the fitted parameter; and the fitted model and its conditioning set are appropriate for the stated prediction. The approximation is first-order: it can be poor with small samples, weakly identified parameters, probabilities near a boundary, or a prediction function that changes sharply over the plausible parameter range.

For a single predicted probability \(\hat p_j\), a familiar approximate Wald interval is \(\hat p_j\pm 1.96\,SE(\hat p_j)\) for 95% confidence. It is an interval-procedure statement: over repeated samples under the assumptions, intervals formed this way cover the fixed target about 95% of the time. It is not a 95% probability that the fixed target lies in this particular interval. Because probabilities must lie between zero and one, a Wald interval that crosses a boundary is a warning that this approximation is questionable.

Worked Example: Two-Step Prediction With Two Estimated Transition Probabilities

Worked Example: Propagating Uncertainty Through Two Transitions

Consider a hypothetical three-state process. State 1 can progress to state 2, state 2 can progress to absorbing state 3, and there are two equal prediction intervals. In each interval, the transition matrix is assumed to have the same form: \[ P(a,b)= \begin{pmatrix} 1-a&a&0\\ 0&1-b&b\\ 0&0&1 \end{pmatrix}. \] At the landmark, the subject is known to be in state 1. Suppose the estimates are \(\hat a=0.20\) and \(\hat b=0.30\), with estimated standard errors 0.04 and 0.05 and estimated correlation 0.25. These are hypothetical estimates and covariance values. The model assumes these transition probabilities apply over both intervals and that the stated covariance adequately describes their joint estimation uncertainty.

The predicted distribution after two intervals is the first row of \(P(a,b)^2\). Enumerating the paths gives \(p_1=(1-a)^2\), \(p_2=a(1-b)+(1-a)a=a(2-a-b)\), and \(p_3=ab\). At the fitted values, \[ \hat{\boldsymbol p}= \big((0.80)^2,\;0.20(2-0.20-0.30),\;0.20(0.30)\big) =(0.640,\;0.300,\;0.060). \] The values sum to one, as required for mutually exclusive and exhaustive states.

The parameter covariance matrix has variances \(0.04^2=0.0016\) and \(0.05^2=0.0025\), and covariance \(0.25(0.04)(0.05)=0.0005\): \[ \widehat V= \begin{pmatrix} 0.0016&0.0005\\ 0.0005&0.0025 \end{pmatrix}. \] The derivatives, with columns ordered as \(a,b\), are \[ D= \begin{pmatrix} -2(1-a)&0\\ 2-2a-b&-a\\ b&a \end{pmatrix}_{(a,b)=(0.20,0.30)} = \begin{pmatrix} -1.6&0\\ 1.3&-0.2\\ 0.3&0.2 \end{pmatrix}. \]

For example, the variance estimate for state 2 is \[ (1.3,-0.2)\widehat V(1.3,-0.2)^{\mathsf T} =1.69(0.0016)+0.04(0.0025)+2(1.3)(-0.2)(0.0005) =0.002544. \] Thus \(SE(\hat p_2)=\sqrt{0.002544}=0.0504\). The full propagated covariance matrix is \[ D\widehat V D^{\mathsf T}= \begin{pmatrix} 0.004096&-0.003168&-0.000928\\ -0.003168&0.002544&0.000624\\ -0.000928&0.000624&0.000304 \end{pmatrix}. \] The estimated standard errors for states 1, 2, and 3 are respectively \(0.0640\), \(0.0504\), and \(0.0174\). For state 3, the approximate 95% Wald interval is \(0.060\pm1.96\sqrt{0.000304}=[0.0258,0.0942]\), rounded to four decimal places. Under the model and large-sample approximation, this gives a plausible range for the state 3 occupation probability two intervals after the landmark. It is not an interval for an individual’s eventual absorption time.

Why the Predicted Probabilities Are Dependent

State probabilities are not separate unrestricted quantities. If the possible states are mutually exclusive and exhaustive, their probabilities add to one at every horizon. Therefore, a change in one component generally requires a compensating change in one or more others. In the two-step example, the derivatives of the three probabilities with respect to \(a\) sum to \(-1.6+1.3+0.3=0\), and the derivatives with respect to \(b\) sum to \(0-0.2+0.2=0\). The propagated covariance matrix accordingly has rows that sum to zero.

This dependence matters when reporting joint uncertainty or uncertainty in a contrast. For example, the variance of \(\hat p_1+\hat p_2\) is not the sum of their variances: it also includes twice their covariance. In fact, because \(p_1+p_2=1-p_3\), its variance equals the variance of \(p_3\), not the sum of the two marginal variances. Reporting only separate standard errors hides this structure.

Derivation: If the state probabilities sum to one for every parameter value, then \(\boldsymbol 1^{\mathsf T}g(\theta)=1\). Differentiating with respect to \(\theta\) gives \(\boldsymbol 1^{\mathsf T}D=0\). Thus \(\boldsymbol 1^{\mathsf T}D\widehat V D^{\mathsf T}=0\): the rows of the propagated covariance matrix sum to zero. This is a useful check on calculations, not evidence that the individual probabilities have no uncertainty.

Worked Example: A One-Interval Multinomial Estimate

Worked Example: Uncertainty From Observed State Counts

In a hypothetical sample of 100 independent individuals, all known to occupy a specified landmark state, the estimated one-interval probabilities of occupying states A, B, and C are based on counts 70, 20, and 10. Assume the individuals are sampled from the same conditional population and the outcome has exactly these three categories. Then \(\hat{\boldsymbol p}=(0.70,0.20,0.10)\).

For multinomial sampling, the estimated covariance matrix of the category proportions is \(\{\operatorname{diag}(\hat{\boldsymbol p})-\hat{\boldsymbol p}\hat{\boldsymbol p}^{\mathsf T}\}/n\). Here it is \[ \frac{1}{100} \begin{pmatrix} 0.21&-0.14&-0.07\\ -0.14&0.16&-0.02\\ -0.07&-0.02&0.09 \end{pmatrix} = \begin{pmatrix} 0.0021&-0.0014&-0.0007\\ -0.0014&0.0016&-0.0002\\ -0.0007&-0.0002&0.0009 \end{pmatrix}. \] The standard error for state C is \(\sqrt{0.0009}=0.0300\). A normal-approximation 95% interval is \(0.10\pm1.96(0.03)=[0.0412,0.1588]\). It describes uncertainty in the conditional state C probability for this one-interval target under the sampling assumptions. The negative covariance between A and B, for instance, reflects the fact that a larger estimated share in one category leaves less share for the others.

This example is a direct estimate of a single transition row, so the multinomial covariance is available immediately. In a multi-step prediction, the estimated transition probabilities may enter repeatedly and may be correlated across rows or time intervals. The general Jacobian calculation is needed in that setting; substituting independent standard errors for the full covariance matrix can give the wrong answer.

Worked Example: Including Uncertainty in Landmark Composition

Worked Example: An Estimated Starting-State Distribution

Suppose a hypothetical target population at the landmark has two states. The estimated fraction in state 1 is \(\hat q=0.60\), based on 100 independent individuals, and state 2 has fraction \(0.40\). For one prediction interval, assume the transition probabilities to state 1 are 0.80 from state 1 and 0.30 from state 2, treated here as known. The predicted probability of state 1 is \[ \hat p_1=0.60(0.80)+0.40(0.30)=0.60. \]

The uncertainty in \(\hat q\) is approximately \(\hat q(1-\hat q)/100=0.60(0.40)/100=0.0024\), so \(SE(\hat q)=0.0490\). Since \(p_1(q)=0.30+0.50q\), its derivative with respect to \(q\) is 0.50. Propagation gives \[ \widehat{\operatorname{Var}}(\hat p_1)=0.50^2(0.0024)=0.000600, \qquad SE(\hat p_1)=\sqrt{0.000600}=0.0245. \] This standard error reflects estimated landmark composition only; the transition probabilities were stipulated as known. If transition estimates were also uncertain, their covariance and any covariance with \(\hat q\) would need to be included in the parameter vector and propagated jointly. The result also clarifies the target: a prediction conditional on a known current state does not have this composition component, whereas a prediction averaged over an estimated landmark population does.

Diagnostics, Alternatives, and Limitations

Before interpreting a propagated standard error, check that the fitted transition rows are valid probability distributions, that the predicted components sum to one, and that the covariance matrix is symmetric and positive semidefinite up to rounding. For a full state-probability vector, the covariance rows should sum approximately to zero when the probabilities are constrained to sum to one. A failure can reveal a derivative, indexing, or covariance-entry error. Also check whether an approximate interval falls outside \([0,1]\), and whether the estimated standard error is large relative to the probability itself.

The delta method is efficient and makes clear which parameters drive uncertainty, but it relies on a local linear approximation and an adequate covariance estimate. It can be unreliable near boundaries or with sparse transitions. The bootstrap provides another way to form intervals by resampling and recomputing the prediction; that approach is especially useful when the prediction function is complicated, although its validity depends on an appropriate resampling scheme and adequate data. A Bayesian analysis instead propagates a posterior distribution through the prediction function and interprets resulting intervals using posterior probability. These approaches answer related uncertainty questions, but their assumptions and interpretations are not interchangeable.

A common mistake is to calculate a standard error for each transition probability, treat them all as independent, and then combine those standard errors while ignoring row-sum constraints and covariance. This can overstate or understate prediction uncertainty. Use the joint covariance matrix of the estimated parameters, including covariance between estimated transition probabilities, and propagate it through the complete prediction function. A second pitfall is to report uncertainty conditional on a fixed landmark state distribution when the intended prediction averages over a distribution that was itself estimated. State clearly which sources of uncertainty the calculation includes.

Key takeaway: For a differentiable landmark prediction \(g(\theta)\), approximate its covariance by \(D\widehat V D^{\mathsf T}\). Preserve parameter covariances, check the probability constraints, and state whether uncertainty in the landmark composition is included.

Check Your Understanding

Answer each question, then compare with the brief solutions below.

If a predicted probability is \(g(\hat\theta)\), what information is needed in addition to its fitted value to obtain a delta-method standard error?
Why should the covariance between two estimated transition parameters be retained when calculating prediction uncertainty?
For probabilities of three exhaustive states, what constraint should their propagated covariance matrix satisfy?
In the two-state composition example, what would the predicted state 1 probability be if \(q=0.40\), with transition probabilities still 0.80 and 0.30?
Name one warning sign that a normal-based Wald interval may be unreliable for a state probability.

Solutions. (1) The gradient of the prediction with respect to the estimated parameters and their joint covariance matrix; the variance estimate is the gradient times covariance times its transpose. (2) Parameter estimates can move together, and the resulting cross-products contribute to the variance of the prediction. (3) Its rows should sum to zero, up to rounding, because the predicted probabilities sum to one. (4) \(0.40(0.80)+0.60(0.30)=0.50\). (5) An interval extending below zero or above one, or a probability near a boundary with a strongly nonlinear prediction function, signals that the normal approximation may be poor.