Delta-Method Intervals For Predicted Risk

Why Differentiate a Risk Prediction?

A model-based cumulative-incidence prediction combines information from the hazard for the event of interest with survival from all event types. Its uncertainty therefore depends not just on the estimated hazard for that event, but also on how every model parameter changes the predicted risk. The delta method provides a systematic first-order way to carry that parameter uncertainty through to a predicted risk.

The central calculation is a gradient: the vector of derivatives of the predicted cumulative incidence with respect to the estimated model parameters. Multiplying this gradient by the parameters’ covariance matrix gives an approximate variance for the prediction. The result is useful when the model estimate and its covariance matrix are available, but its reliability depends on the model, the covariance estimate, and the quality of the normal approximation.

Prediction and Its Derivative

Consider a specified covariate profile and a prediction horizon $t$. Let $\lambda_k(u;\theta)$ be the cause-specific hazard for cause $k$ at time $u$, where $\theta$ is the vector of model parameters. Suppose there are $K$ competing causes, and define the total cumulative hazard through time $u$ as $H(u;\theta)=\int_0^u\sum_{j=1}^K\lambda_j(v;\theta)\,dv$. The event-free survival probability is $S(u;\theta)=\exp\{-H(u;\theta)\}$.

Definition: The model-based cumulative incidence for cause $k$ by time $t$, at the specified covariate profile, is $$ F_k(t;\theta)=\int_0^t S(u;\theta)\lambda_k(u;\theta)\,du. $$ If $\hat\theta$ is the fitted parameter vector, the predicted risk is $\hat F_k(t)=F_k(t;\hat\theta)$.

This formula makes explicit why differentiating only the cause-$k$ hazard is not enough. A change in a parameter can alter $\lambda_k$, the total hazard, or both. Because the total hazard determines $S(u;\theta)$, it can change the predicted incidence even when it belongs to a competing cause.

Assume the hazards are differentiable in $\theta$, the integrals are finite, and differentiation can be interchanged with integration. For parameter component $\theta_r$, differentiating $S(u;\theta)=\exp\{-H(u;\theta)\}$ gives $\partial S/\partial\theta_r=-S\,\partial H/\partial\theta_r$. Therefore,

\frac{\partial F_k(t;\theta)}{\partial\theta_r} = \int_0^t S(u;\theta) \left[ \frac{\partial\lambda_k(u;\theta)}{\partial\theta_r} - \lambda_k(u;\theta)\frac{\partial H(u;\theta)}{\partial\theta_r} \right]du, \qquad \frac{\partial H(u;\theta)}{\partial\theta_r} = \int_0^u\sum_{j=1}^K \frac{\partial\lambda_j(v;\theta)}{\partial\theta_r}\,dv.

The first term inside the integral is the direct change in the cause-$k$ hazard. The second is the change in survival up to time $u$, which can reflect changes in any cause-specific hazard. This derivative is evaluated at the fitted parameter values when constructing the interval.

Delta-Method Standard Error and Interval

Let $\hat\theta$ estimate $\theta$, with estimated covariance matrix $\widehat V_\theta$. If $\hat\theta$ is approximately multivariate normal around $\theta$, and $F_k(t;\theta)$ is sufficiently smooth near the fitted value, a first-order Taylor expansion gives the delta-method variance.

Theorem: Let $g$ be the gradient vector of $F_k(t;\theta)$, evaluated at $\hat\theta$. The delta-method estimate of the prediction variance and standard error is $$ \widehat{\operatorname{Var}}\{\hat F_k(t)\} = g^\mathsf{T}\widehat V_\theta g, \qquad SE\{\hat F_k(t)\} = \sqrt{g^\mathsf{T}\widehat V_\theta g}. $$ An approximate two-sided $100(1-\alpha)\%$ Wald interval on the risk scale is $$ \hat F_k(t)\ \pm\ z_{1-\alpha/2}\,SE\{\hat F_k(t)\}, $$ where $z_{1-\alpha/2}$ is the corresponding standard normal quantile.

In this calculation, the covariance matrix must correspond to the parameter vector used in the prediction. If a model includes estimated baseline-hazard quantities as well as regression coefficients, whether their uncertainty is included depends on how $\widehat V_\theta$ was constructed. Treating estimated quantities as fixed yields uncertainty conditional on those quantities; it does not automatically account for their estimation.

Worked Example: A Constant-Hazard Two-Cause Prediction

Suppose, hypothetically, that a fixed profile has constant cause-specific hazards $\lambda_1=0.12$ per year and $\lambda_2=0.08$ per year. Predict the cumulative incidence of cause 1 by $t=2$ years. Assume the estimated hazards have standard errors 0.020 and 0.015 per year, with correlation 0.25. Use a 95% Wald interval.

With total hazard $\Lambda=\lambda_1+\lambda_2$, event-free survival is $\exp(-\Lambda u)$, so the cause-1 risk is $$ F_1(t)=\int_0^t e^{-\Lambda u}\lambda_1\,du =\frac{\lambda_1}{\Lambda}(1-e^{-\Lambda t}). $$

Here $\Lambda=0.20$ per year, and $$ \hat F_1(2)=\frac{0.12}{0.20}(1-e^{-0.20(2)}) =0.6(1-e^{-0.4}) =0.1978. $$ Thus the predicted cause-1 risk by 2 years is about 19.8%.

To calculate the gradient, write $q(\Lambda)=(1-e^{-\Lambda t})/\Lambda$, so $F_1=\lambda_1q(\Lambda)$. Its derivative is $$ q'(\Lambda)=\frac{\Lambda t e^{-\Lambda t}-(1-e^{-\Lambda t})}{\Lambda^2}. $$ At $\Lambda=0.20$ and $t=2$, $q=1.6484$ years and $q'=-1.5388$ years squared. Since $\Lambda$ changes with either hazard, $$ \frac{\partial F_1}{\partial\lambda_1}=q+\lambda_1q'=1.4637\text{ years}, \qquad \frac{\partial F_1}{\partial\lambda_2}=\lambda_1q'=-0.1847\text{ years}. $$ The negative derivative for $\lambda_2$ has a useful interpretation: increasing the competing-cause hazard reduces the chance of remaining event-free long enough to experience cause 1.

The covariance between the hazards is $0.25(0.020)(0.015)=0.000075$ per year squared. Therefore, $$ \widehat{\operatorname{Var}}(\hat F_1) =(1.4637)^2(0.020)^2+(-0.1847)^2(0.015)^2 +2(1.4637)(-0.1847)(0.000075) =0.0008242. $$ The standard error is $\sqrt{0.0008242}=0.02871$. The 95% interval is $$ 0.1978\pm1.96(0.02871)=(0.1415,\ 0.2541). $$ Under the stated model and covariance assumptions, this Wald interval gives an approximate range of 14.2% to 25.4% for the cause-1 risk by 2 years.

This calculation assumes the constant-hazard model is an appropriate description over the 2-year horizon, the hazard estimates and their covariance are on the stated scale, and their sampling distribution is adequately approximated by a bivariate normal distribution. The covariance term is included because the two hazard estimates are correlated.

Using the Chain Rule for Regression Parameters

In regression models, predictions are usually expressed through model coefficients rather than directly through hazard values at the target profile. The delta-method gradient must then be taken with respect to those coefficients. A convenient way to do this is to first differentiate with respect to the hazards, then apply the chain rule from hazards to coefficients.

Worked Example: A Log-Linear Hazard Model at a Fixed Profile

Suppose, hypothetically, that two constant cause-specific hazards follow $\lambda_j(x)=\exp(\alpha_j+\beta_jx)$, where $x$ is a continuous covariate. At the target profile $x=1.5$, let the fitted hazards be $0.12$ and $0.08$ per year. Use a 2-year horizon and the fitted coefficient order $(\alpha_1,\beta_1,\alpha_2,\beta_2)$. Assume the coefficient covariance matrix is diagonal, with standard errors $0.10, 0.08, 0.12, 0.10$, respectively.

The predicted risk is the same constant-hazard expression as in the previous calculation, so $\hat F_1(2)=0.1978$, $\partial F_1/\partial\lambda_1=1.4637$ years, and $\partial F_1/\partial\lambda_2=-0.1847$ years. Since $\partial\lambda_j/\partial\alpha_j=\lambda_j$ and $\partial\lambda_j/\partial\beta_j=x\lambda_j$, the coefficient gradient is $$ g= \left( \lambda_1\frac{\partial F_1}{\partial\lambda_1}, x\lambda_1\frac{\partial F_1}{\partial\lambda_1}, \lambda_2\frac{\partial F_1}{\partial\lambda_2}, x\lambda_2\frac{\partial F_1}{\partial\lambda_2} \right) =(0.1756,\ 0.2635,\ -0.01477,\ -0.02215). $$

With the stated diagonal covariance matrix, the estimated variance is $$ (0.1756)^2(0.10)^2+(0.2635)^2(0.08)^2 +(-0.01477)^2(0.12)^2+(-0.02215)^2(0.10)^2 \approx0.0007609. $$ Thus the standard error is about $0.02759$, and the 95% raw-scale Wald interval is $0.1978\pm1.96(0.02759)=(0.1437,\ 0.2519)$. This calculation assumes that the supplied coefficient covariance matrix is appropriate; if coefficient estimates are correlated, their off-diagonal covariance terms must also be included.

Intervals on a Transformed Scale

The raw-scale Wald interval is simple, but it can extend below 0 or above 1, even though a cumulative incidence is a probability. A common alternative is to transform the risk, calculate an interval on the transformed scale, and transform the endpoints back. For a risk strictly between 0 and 1, define $\eta=\operatorname{logit}(F)=\log\{F/(1-F)\}$. The delta-method gradient for $\eta$ is the risk gradient divided by $F(1-F)$.

Worked Example: A Logit-Scale Interval

Use the cause-1 prediction from the constant-hazard example, but suppose the covariance matrix is for the log hazards $(\log\lambda_1,\log\lambda_2)$. Assume their standard errors are 0.20 and 0.15, with correlation 0.30. This is a different covariance specification from the one used for the raw-scale interval in that example.

Because $\partial F_1/\partial\log\lambda_j=\lambda_j(\partial F_1/\partial\lambda_j)$, the gradient for risk with respect to the log hazards is $(0.1756,-0.01477)$. The covariance between the log hazards is $0.30(0.20)(0.15)=0.009$. Propagating this covariance gives an approximate risk standard error of $0.03453$. For $\eta=\operatorname{logit}(F_1)$, the standard error is approximately $$ SE(\hat\eta)=\frac{0.03453}{0.1978(1-0.1978)}=0.2175. $$

Here $\operatorname{logit}(0.1978)\approx-1.400$. A 95% interval on the logit scale is $-1.400\pm1.96(0.2175)=(-1.826,-0.974)$. Applying the inverse-logit to both endpoints gives an approximate risk interval of $(0.1388,0.2742)$. Its endpoints lie between 0 and 1, but it remains an approximation and depends on the same model and covariance assumptions.

Diagnostics, Comparisons, and Limitations

A delta-method interval is most credible when the prediction is a smooth function of well-estimated parameters and the fitted parameter distribution is reasonably close to normal. Several checks help identify problems:

Check the gradient. Compare analytic derivatives with small numerical perturbations of each parameter. A disagreement can reveal a sign, scale, or chain-rule error.
Check the covariance calculation. Confirm the parameter order and scale match the gradient, include relevant covariance terms, and verify that the resulting variance is nonnegative apart from negligible rounding error.
Check local sensitivity. A large standard error or a gradient that changes sharply under small parameter perturbations signals that a first-order approximation may be fragile.
Check the interval endpoints. Raw-scale endpoints outside the probability range indicate that the untransformed Wald approximation is unsuitable as a reported risk interval.

The raw-scale Wald interval is easy to calculate and interpret, but its symmetry can be a poor description of uncertainty for risks near 0 or 1. A logit-scale interval respects the probability range after back-transformation and can be asymmetric, but it still relies on a normal approximation after transformation and is not defined at exactly 0 or 1. Bootstrap intervals are an alternative way to approximate prediction uncertainty by resampling and refitting; they may better reflect nonlinearity, but their reliability depends on the resampling scheme and the stability of the fitted model.

The delta method is a first-order approximation, not a guarantee of accurate coverage in a finite sample. It can be unreliable when parameters are weakly identified, estimates are near constraints, the prediction is highly nonlinear over the plausible parameter range, or the predicted risk is close to a boundary. It also does not repair model misspecification. An interval quantifies uncertainty under the chosen model and covariance estimate; it does not include uncertainty about whether that model is appropriate.

A Common Mistake

Common pitfall: Differentiating the cause-$k$ hazard while treating event-free survival as fixed. This omits the contribution of $-\lambda_k(u;\theta)\,\partial H(u;\theta)/\partial\theta_r$ from the gradient. For a competing-cause parameter, the direct derivative of $\lambda_k$ may be zero even though the predicted cause-$k$ risk changes through survival. Differentiate the full cumulative-incidence expression, including the total hazard inside survival, before applying the covariance formula.

Another practical error is to report $\hat F_k\pm1.96SE$ as though it were an exact interval. It is an approximate confidence interval: under repeated sampling and the stated model and regularity conditions, the procedure is intended to cover the true model-based risk at approximately its nominal rate. It is not a probability statement that the fixed risk lies in this particular realized interval.

Check Your Understanding

Answer the questions, then compare your reasoning with the solutions below.

In a competing-risks prediction, why can a parameter for a competing cause affect the predicted incidence of the event of interest?
If the fitted parameter covariance is $\widehat V_\theta$ and the prediction gradient is $g$, what quantity gives the delta-method variance?
For a log-linear hazard $\lambda(x)=\exp(\alpha+\beta x)$, what are the derivatives of the hazard with respect to $\alpha$ and $\beta$?
What does it suggest if a raw-scale Wald interval for a predicted risk extends above 1?

Solutions. (1) The competing-cause hazard contributes to the total hazard and therefore changes event-free survival; the cause-specific cumulative incidence integrates its event hazard against that survival. (2) The variance estimate is $g^\mathsf{T}\widehat V_\theta g$, and the standard error is its square root. (3) The derivatives are $\partial\lambda/\partial\alpha=\lambda$ and $\partial\lambda/\partial\beta=x\lambda$; prediction derivatives then follow by the chain rule. (4) The symmetric raw-scale normal approximation is not suitable as a probability interval in this case; a transformed interval may be preferable, while recognizing that it too is approximate.