Delta-Method Uncertainty For Aalen-Johansen Estimates

Why Uncertainty Must Follow the Product Integral

An Aalen–Johansen estimate is built by multiplying estimated transition matrices across event times. Its uncertainty therefore depends not just on the uncertainty of each estimated transition intensity, but also on how the product integral combines those estimates. The delta method provides a systematic way to propagate that uncertainty to a selected transition probability, such as the probability of occupying state 3 at a fixed horizon when starting in state 1.

The key object is the covariance matrix of the estimated intensity increments. Ignoring covariance, or treating the final transition probability as though it were a single intensity estimate, can give an incorrect standard error. This tutorial focuses on uncertainty for selected entries of the estimated transition matrix; it does not derive covariance formulas for an entire vector of state-occupation probabilities.

Definition: Let $\widehat{\theta}$ be a vector of estimated transition-intensity increments, with estimated covariance matrix $\widehat{\Sigma}$. Let $\widehat{p}=f(\widehat{\theta})$ be a selected Aalen–Johansen transition probability, expressed as a differentiable function of those increments. The delta-method variance estimate is $$ \widehat{\operatorname{Var}}(\widehat{p}) = \widehat{g}^{\mathsf T}\widehat{\Sigma}\widehat{g}, \qquad \widehat{g} = \left.\frac{\partial f(\theta)}{\partial \theta}\right|_{\theta=\widehat{\theta}}. $$ The square root of this variance estimate is the delta-method standard error.

This approximation relies on the estimated increments being approximately jointly normal in large samples, with a covariance estimate that consistently estimates their sampling covariance. The selected function must be differentiable near the true increment vector. In the usual finite-grid Aalen–Johansen calculation, the transition matrices are formed from estimated intensity increments and multiplied in time order, so their entries—and hence any selected entry of the product—are differentiable functions of those increments.

The underlying transition-intensity estimates must also be valid for the data and risk sets being analyzed. In particular, the risk sets must be defined appropriately, and the assumptions supporting the Aalen–Johansen estimator must be reasonable. The estimated covariance matrix must reflect the joint uncertainty of the increments used in the calculation. For distinct transition types, standard counting-process variance calculations often give zero estimated cross-covariances under their usual conditions, but this is not a rule to impose blindly: if the covariance calculation or data structure implies nonzero covariance, it must be retained.

Differentiating the Product-Integral Recursion

Suppose the event-time updates, in chronological order, are $\widehat{B}_1,\ldots,\widehat{B}_m$. Each update is an identity matrix plus the estimated transition-intensity increment matrix at that time. The estimated transition matrix at the horizon is $\widehat{P}=\widehat{B}_1\widehat{B}_2\cdots\widehat{B}_m$. For a selected entry $\widehat{P}_{ij}$, the delta method requires its derivative with respect to each uncertain increment.

The derivative can be computed by differentiating the same recursion used to calculate the estimate. Write $P_0=I$ and $P_r=P_{r-1}B_r$. For a scalar increment parameter $\theta_q$, the product rule gives

\frac{\partial P_r}{\partial \theta_q} = \frac{\partial P_{r-1}}{\partial \theta_q}B_r + P_{r-1}\frac{\partial B_r}{\partial \theta_q}.

If $\theta_q$ does not occur in $B_r$, its derivative there is the zero matrix. Starting with $\partial P_0/\partial\theta_q=0$, repeated application yields $\partial P_m/\partial\theta_q$; the derivative of the selected probability is the $(i,j)$ entry of that matrix. This is a direct justification of the gradient calculation: each uncertain update contributes both through the product accumulated before it and through all later updates that multiply it.

Once the gradient is available, the covariance calculation is a quadratic form. When the increments are independent, this reduces to summing the squared derivatives times their variances. When they are correlated, the cross-product terms matter. The following examples show both calculations.

Worked Example: A Two-State Probability With Two Updates

Consider a hypothetical two-state model in which state 1 can transition to absorbing state 2. At two ordered event times, the estimated transition-intensity increments are $\widehat{a}_1=0.10$ and $\widehat{a}_2=0.20$. Suppose their estimated variances are $0.0025$ and $0.0064$, respectively, and their estimated covariance is zero. These are assumed to be valid covariance estimates from the intensity-estimation procedure.

At an update with increment $a_r$, the probability of remaining in state 1 is multiplied by $1-a_r$. Thus the estimated probability of transitioning from state 1 to state 2 by the horizon is

\widehat{p}_{12} = 1-(1-\widehat{a}_1)(1-\widehat{a}_2) = 1-(0.90)(0.80) = 0.28.

For $f(a_1,a_2)=1-(1-a_1)(1-a_2)$, the gradient is $\nabla f=(1-a_2,1-a_1)$. At the estimates, it is $(0.80,0.90)$. The estimated variance is therefore

\widehat{\operatorname{Var}}(\widehat{p}_{12}) = (0.80)^2(0.0025)+(0.90)^2(0.0064) = 0.001600+0.005184 = 0.006784.

The standard error is $\sqrt{0.006784}=0.08236$. A nominal 95% Wald interval is $0.28\mathbin{\pm}1.96(0.08236)$, or approximately $(0.1186,0.4414)$. Under the delta-method normal approximation, this interval is a procedure with approximately 95% repeated-sampling coverage under the model and regularity conditions; it is not a statement that there is a 95% probability that the fixed true probability lies in this particular realized interval.

A Multistate Calculation

For a multistate model, define each event-time update using its estimated transition increments. The diagonal entry in a row is one minus the sum of the outgoing increments from that state; off-diagonal entries are the corresponding increments. This construction preserves each row sum at one when the increments are valid. The gradient for a selected probability can be found by differentiating the product of these update matrices, as in the recursion above.

Worked Example: Uncertainty for a Selected Three-State Probability

Consider a hypothetical three-state model with state 3 absorbing. At the first update, state 1 has estimated increments $a=0.20$ to state 2 and $b=0.10$ to state 3. At a later update, state 1 has increment $c=0.10$ to state 3, and state 2 has increment $d=0.30$ to state 3. Suppose the estimated covariance matrix for $(a,b,c,d)$ is diagonal with entries $0.002, 0.001, 0.0015, 0.004$. These hypothetical inputs assume that the intensity estimates and their covariance matrix satisfy the conditions described above.

The two update matrices are

B_1= \begin{pmatrix} 1-a-b & a & b\\ 0 & 1 & 0\\ 0 & 0 & 1 \end{pmatrix}, \qquad B_2= \begin{pmatrix} 1-c & 0 & c\\ 0 & 1-d & d\\ 0 & 0 & 1 \end{pmatrix}.

The probability of being in state 3 at the horizon, starting in state 1, is the $(1,3)$ entry of $B_1B_2$. There are three routes into state 3: remain in state 1 at the first update and transition later, transition to state 2 and then to state 3, or transition directly to state 3 at the first update. Consequently,

p_{13}=(1-a-b)c+ad+b = (0.70)(0.10)+(0.20)(0.30)+0.10 =0.23.

Differentiate the expression to obtain $\nabla p_{13}=(d-c,\,1-c,\,1-a-b,\,a)$. At the estimates, this is $(0.20,0.90,0.70,0.20)$. Using the stated diagonal covariance matrix,

\begin{aligned} \widehat{\operatorname{Var}}(\widehat{p}_{13}) &=(0.20)^2(0.002)+(0.90)^2(0.001) +(0.70)^2(0.0015)+(0.20)^2(0.004)\\ &=0.000080+0.000810+0.000735+0.000160\\ &=0.001785. \end{aligned}

Thus $SE=\sqrt{0.001785}=0.04225$. The nominal 95% Wald interval is $0.23\mathbin{\pm}1.96(0.04225)$, approximately $(0.1472,0.3128)$. The interval describes uncertainty in the estimated probability of being in state 3 at the stated horizon, starting from state 1. It does not describe the probability of a particular individual taking one of these paths.

Why the Full Covariance Matrix Matters

A diagonal covariance matrix is convenient, but it is an assumption about the estimated increments, not a universal property of the product integral. In general, for gradient components $g_r$ and covariance entries $\widehat{\Sigma}_{rs}$, the variance is $\sum_r\sum_s g_rg_s\widehat{\Sigma}_{rs}$. Covariances can either increase or decrease the resulting variance, depending on their signs and the corresponding gradient components.

Worked Example: Including a Nonzero Covariance

In a hypothetical two-update, two-state calculation, let $\widehat{a}_1=0.15$ and $\widehat{a}_2=0.25$. Suppose the estimated variances are $0.0016$ and $0.0036$, and the estimated covariance is $0.0012$. The covariance matrix is valid: its determinant is $0.0016(0.0036)-(0.0012)^2=0.00000432>0$, and its diagonal entries are positive.

The estimated transition probability is $1-(0.85)(0.75)=0.3625$. Its gradient with respect to $(a_1,a_2)$ is $(1-a_2,1-a_1)=(0.75,0.85)$. The variance, including the off-diagonal terms, is

\begin{aligned} \widehat{\operatorname{Var}}(\widehat{p}_{12}) &=(0.75)^2(0.0016)+(0.85)^2(0.0036) +2(0.75)(0.85)(0.0012)\\ &=0.000900+0.002601+0.001530\\ &=0.005031. \end{aligned}

The standard error is $0.07093$, giving a nominal 95% Wald interval of $0.3625\mathbin{\pm}1.96(0.07093)$, or approximately $(0.2235,0.5015)$. If the covariance were wrongly set to zero, the estimated variance would be $0.003501$ and the standard error about $0.05917$, understating uncertainty in this example. The covariance term is part of the calculation, not an optional correction.

Diagnostics, Alternatives, and Limitations

Several checks help reveal when a delta-method calculation may be going wrong:

Check the covariance matrix. It should be symmetric and positive semidefinite, apart from small numerical rounding errors. Implausible entries or a materially negative eigenvalue can signal a coding or covariance-estimation problem.
Check the gradient against the recursion. For a small perturbation $\epsilon$ in one increment, compare the analytic derivative with the finite difference $\{f(\widehat{\theta}+\epsilon e_r)-f(\widehat{\theta})\}/\epsilon$. Use a sufficiently small perturbation that the update remains valid.
Check the interval against the parameter space. A transition probability must lie between zero and one. A Wald interval extending beyond those limits may be a warning that a symmetric normal approximation is poor, especially near a boundary or when the standard error is large relative to the probability.
Check the inputs and time order. The updates must match the estimated increments and be multiplied in chronological order. A mistaken transition direction, omitted increment, or reversed product can change both the estimate and its gradient.

The delta method is attractive because it is transparent, computationally efficient, and directly tied to the estimated intensity covariance. A Greenwood-type calculation is a special-case style of variance calculation for particular survival or transition settings; it does not replace the general gradient-and-covariance calculation when the selected probability depends on several increments. Bootstrap confidence intervals are another alternative, particularly when a normal approximation appears questionable, but they require a suitable resampling scheme and do not automatically resolve model or data problems.

The approximation is local: it uses the slope of $f$ near the estimated increments and a normal approximation for their joint distribution. It can be unreliable with small effective risk sets, highly skewed estimates, probabilities near zero or one, or constraints that make the sampling distribution strongly asymmetric. Transforming an interval to keep it within the probability range can address an impossible endpoint, but does not by itself establish good coverage. When reporting results, state the selected origin and destination states, horizon, point estimate, standard error or interval, and the assumptions supporting the calculation.

Common Mistake

Common pitfall: Treating the transition probability as a single estimated intensity increment, or adding increment variances without accounting for derivatives and covariances. The correction is to write the selected probability as a function of all relevant increments, differentiate that function, and calculate $\widehat{g}^{\mathsf T}\widehat{\Sigma}\widehat{g}$. A product-integral probability is generally not a simple sum of the increment estimates.

Key takeaway: For a selected Aalen–Johansen transition probability, uncertainty is propagated through the product integral by differentiating the product with respect to its estimated intensity increments and applying the full covariance matrix. The resulting standard error and Wald interval are useful when the joint normal and local linear approximations are adequate.

Check Your Understanding

Answer each question, then compare your work with the solutions below.

For a selected transition probability $f(\theta)$, what two ingredients determine its delta-method variance?
In a two-state model with two increments $a_1=0.10$ and $a_2=0.20$, calculate the transition probability by the horizon.
Why can setting an estimated covariance to zero produce a wrong standard error?
A nominal Wald interval for a transition probability is $(-0.04,0.18)$. What is a diagnostic concern, and what does this interval not imply?

Solutions. (1) The gradient of the probability with respect to the increments and their covariance matrix; the variance is $g^{\mathsf T}\Sigma g$. (2) $1-(1-0.10)(1-0.20)=1-(0.90)(0.80)=0.28$. (3) The variance includes cross-products $2g_rg_s\Sigma_{rs}$; omitting nonzero covariance terms can either increase or decrease the variance. (4) The interval extends below zero, so the symmetric normal approximation may be poor or the standard error may be large relative to the estimate. It does not imply that a probability outside $[0,1]$ is possible; it is an approximate interval procedure whose endpoints can fall outside the parameter space.