Covariance Of Estimated State-Occupation Probabilities

Why Joint Covariance Matters

At a fixed horizon, the estimated probabilities of occupying the possible states form a vector, not a collection of unrelated estimates. If a person is more likely to be in one state, the probabilities assigned to other states must adjust. Consequently, the estimates are generally dependent, and their covariances matter whenever we compare states or calculate uncertainty for a combination of probabilities.

This tutorial develops a plug-in covariance estimate for the Aalen–Johansen state-occupation vector. The calculation follows the same product updates that produce the estimates, while carrying each subject’s influence through those updates. It applies to a finite-state model observed on independent subjects, under the conditions that support the Aalen–Johansen estimator and its usual large-sample approximation. In particular, transition and censoring records must be observed consistently, and censoring must not selectively conceal future state transitions after conditioning on the observed history used by the model. The details of censoring assumptions belong to the introductory material on Aalen–Johansen estimation and its dedicated treatment of independent censoring.

The state probabilities are obtained from the Aalen–Johansen transition-probability matrix. With an initial row distribution \(\hat{\boldsymbol p}(0)\), the estimate at time \(t\) is \(\hat{\boldsymbol p}(0)\hat{\mathbf P}(t)\). We use row vectors for probabilities and transition matrices; when writing a covariance matrix, probability vectors inside outer products are treated as columns. The estimates sum to one, so their covariance matrix cannot have full rank when all \(K\) states are included.

Propagating Influence Through Product Updates

Suppose there are \(n\) independent subjects, and consider distinct transition event times \(u\) up to \(t\). Let \(Y_i(u)\) be the number at risk in state \(i\) immediately before \(u\), and let \(dN_{ij}(u)\) be the number of observed transitions from state \(i\) to state \(j\) at that time. For \(i\ne j\), the estimated cumulative transition-intensity increment is \(d\hat A_{ij}(u)=dN_{ij}(u)/Y_i(u)\); the diagonal increment is \(d\hat A_{ii}(u)=-\sum_{j\ne i}d\hat A_{ij}(u)\). With distinct event times, the product update is \(\mathbf I+d\hat{\mathbf A}(u)\).

For subject \(r\), let \(Y_{ir}(u)\) indicate that the subject is at risk in state \(i\) just before \(u\), and let \(dN_{ij,r}(u)\) indicate their transition from \(i\) to \(j\) at \(u\). A subject-level contribution to the estimated intensity increment is

Set the diagonal contribution to \(d\hat A_{ii,r}(u)=-\sum_{j\ne i}d\hat A_{ij,r}(u)\). These contributions sum to zero over subjects for every increment: \(\sum_r dN_{ij,r}=dN_{ij}\) and \(\sum_r Y_{ir}=Y_i\). They represent each subject’s centered contribution to the estimated intensity increment.

Write \(\boldsymbol\phi_r(u)\) for subject \(r\)’s row-vector contribution to the state-probability estimate immediately after time \(u\). If the initial distribution is fixed, set \(\boldsymbol\phi_r(0)=\boldsymbol 0\). If it is estimated from the same subjects, include its contribution: for a subject beginning in state \(i\), the initial contribution is \((\boldsymbol e_i-\hat{\boldsymbol p}(0))/n\), where \(\boldsymbol e_i\) is the row vector with a one in position \(i\) and zeros elsewhere. Propagate the contributions by

The first term carries the subject’s existing probability contribution through the estimated update. The second adds the subject’s contribution to that update, weighted by the estimated probability of being in each origin state just before the event time. Continue through all event times up to \(t\).

Why This Formula Works

The justification is a first-order expansion of the product updates. If the estimated state vector immediately before an event is perturbed by a small row vector \(\boldsymbol\phi_r(u-)\), and the estimated intensity increment by \(d\hat{\mathbf A}_r(u)\), then the updated vector is approximately

Multiplying and retaining the baseline and first-order terms gives \(\hat{\boldsymbol p}(u-)\{\mathbf I+d\hat{\mathbf A}(u)\}+\boldsymbol\phi_r(u-)\{\mathbf I+d\hat{\mathbf A}(u)\}+\hat{\boldsymbol p}(u-)d\hat{\mathbf A}_r(u)\). The first term is the fitted update; the remaining terms are exactly the stated recursion. The omitted product \(\boldsymbol\phi_r(u-)d\hat{\mathbf A}_r(u)\) is second order in the perturbations and is discarded in the first-order approximation.

Under the usual large-sample conditions, the estimator’s error is approximated by the sum of its subject-level contributions. Independent subjects then give the covariance estimate as the sum of their contribution outer products. This is the same first-order logic used for delta-method uncertainty, expressed here in a form that retains all state-to-state covariance terms.

The Unit-Sum Constraint and Linear Combinations

At every time, the state probabilities sum to one. Therefore, if \(\boldsymbol 1\) is the column vector of \(K\) ones, then \(\hat{\mathbf V}(t)\boldsymbol 1=\boldsymbol 0\) in exact arithmetic. This also follows from the contribution formula: each subject’s contribution sums to zero, so every outer product contributes zero when multiplied by \(\boldsymbol 1\). Thus, the covariance matrix is singular, with at least one zero eigenvalue. This is a structural feature, not a failure of the estimator.

For any fixed column vector \(\boldsymbol c\), the variance of the estimated linear combination \(\boldsymbol c^{\mathsf T}\hat{\boldsymbol p}(t)^{\mathsf T}\) is estimated by \(\boldsymbol c^{\mathsf T}\hat{\mathbf V}(t)\boldsymbol c\). For example, to compare the probabilities of states 1 and 2, take \(\boldsymbol c=(1,-1,0)^{\mathsf T}\). The variance then includes both marginal variances and their covariance; omitting the covariance gives the wrong uncertainty for the difference.

Diagnostics, Alternatives, and Limitations

Several checks help detect implementation errors. The estimated covariance matrix should be symmetric, its diagonal entries should be nonnegative, and its row sums should be approximately zero. It should also be positive semidefinite, apart from small numerical errors. A materially nonzero row sum may indicate that contributions were not centered, a diagonal increment was formed incorrectly, or the state vector and matrix multiplication conventions were mixed. Tiny negative eigenvalues can arise from rounding; substantial negative eigenvalues require investigation.

A common mistake is to calculate a separate standard error for each state and then treat the estimates as independent. That discards the off-diagonal covariances and produces incorrect uncertainty for contrasts and sums. The correction is to retain the whole covariance matrix and calculate uncertainty for a linear combination using \(\boldsymbol c^{\mathsf T}\hat{\mathbf V}\boldsymbol c\). Another mistake is to interpret singularity as a defect. With all states included, the unit-sum constraint necessarily makes one direction have zero variance.

The covariance estimate relies on the large-sample, first-order approximation and on valid estimation of the transition increments. Sparse risk sets, few observed transitions, heavy censoring, or a state with almost no observations can make the approximation unstable. The distinct-event-time recursion also needs modification when transitions are tied: the grouped update and its covariance contribution must respect the actual joint increment structure. For a small sample or an especially complicated estimator, resampling can provide an alternative uncertainty assessment; bootstrap confidence intervals are treated separately.

Unlike reporting only marginal standard errors, the joint matrix supports direct inference about differences and other prespecified combinations. Unlike a covariance calculation based on a simple endpoint multinomial, it accommodates the estimated transition history in the Aalen–Johansen product. Its limitations are equally important: it does not repair model misspecification, informative censoring, inadequate state definitions, or a weak large-sample approximation.