Propagating Baseline Hazard Uncertainty

Why Baseline-Hazard Uncertainty Matters

An absolute-risk prediction is built from estimated cause-specific hazards. Even if a patient’s covariates are fixed and the regression coefficients were known exactly, the estimated baseline hazards would still vary from sample to sample. That variation changes predicted survival and the cumulative incidence for each cause. A prediction interval that treats the baseline estimates as fixed therefore omits a potentially important source of uncertainty.

The focus here is how to propagate uncertainty in estimated baseline hazards through a competing-risk prediction. Earlier material on uncertainty in absolute-risk predictions introduced the broader problem, while the delta-method and bootstrap material covered general interval strategies. The emphasis now is on identifying the baseline-hazard component, retaining dependence among its estimates, and recognizing what a baseline-only interval does—and does not—include.

Definition: Baseline-hazard uncertainty is the sampling uncertainty in estimated baseline hazard quantities, such as cause-specific baseline hazard increments or rates. Propagating it means carrying their estimated covariance through the function that maps the hazards to a predicted cumulative incidence.

Prediction as a Function of Estimated Hazards

For a specified covariate profile, a cause-specific model combines a baseline hazard with a covariate multiplier. For cause \(k\), write the hazard at time \(t\) as \(\lambda_k(t\mid z)=\lambda_{0k}(t)\exp(\beta_k^\mathsf{T}z)\). Here \(\lambda_{0k}(t)\) is the cause-specific baseline hazard, \(z\) is the specified covariate vector, and \(\beta_k\) is the corresponding coefficient vector. A cumulative-incidence prediction depends on the hazards for the target cause and its competing causes: increases in one cause’s hazard can change both the target cause’s event contribution and the probability of remaining event-free.

At a finite set of event times, a fitted model typically supplies estimated baseline cumulative-hazard increments. Denote the collection of baseline quantities used in a prediction by the vector \(\widehat{\boldsymbol{\theta}}_0\), and let \(\widehat{\Sigma}_0\) be its estimated covariance matrix. With regression coefficients temporarily held fixed, write the predicted cumulative incidence as \(g(\widehat{\boldsymbol{\theta}}_0)\). A first-order Taylor expansion around the estimated baseline quantities gives the approximate variance

\widehat{\operatorname{Var}}_0\{g(\widehat{\boldsymbol{\theta}}_0)\} \approx \nabla g(\widehat{\boldsymbol{\theta}}_0)^\mathsf{T} \widehat{\Sigma}_0 \nabla g(\widehat{\boldsymbol{\theta}}_0).

The subscript 0 emphasizes that this is the variance contribution from baseline-hazard estimation alone, conditional on the fitted regression coefficients and the specified profile. The gradient contains the prediction’s sensitivity to each baseline quantity. The off-diagonal entries of \(\widehat{\Sigma}_0\) matter whenever baseline estimates are correlated; they are not optional adjustments.

A One-Interval Derivation

The following calculation gives an explicit sensitivity formula. Suppose there are two causes and, over a single interval of length one, each cause-specific hazard is constant. Let \(\lambda_{01}\) and \(\lambda_{02}\) be the baseline rates, and suppose the specified profile multiplies them by 1.5 and 0.8. The profile-specific rates are \(x=1.5\lambda_{01}\) and \(y=0.8\lambda_{02}\), with total rate \(d=x+y\).

Under these constant-hazard assumptions, the probability of any event during the interval is \(1-\exp(-d)\). Conditional on an event in the interval, the probability that its cause is 1 is \(x/d\). Thus cause 1’s cumulative incidence over the interval is

F_1=\frac{x}{d}\{1-\exp(-d)\}.

To differentiate with respect to the profile-specific rate \(x\), apply the product rule, using \(\partial d/\partial x=1\). Differentiating \(x/d\) gives \(y/d^2\), while differentiating \(1-\exp(-d)\) gives \(\exp(-d)\). Therefore

\frac{\partial F_1}{\partial x} = \frac{y\{1-\exp(-d)\}}{d^2} + \frac{x}{d}\exp(-d).

Similarly, differentiating with respect to \(y\) changes \(d\) but not the numerator \(x\). Hence

\frac{\partial F_1}{\partial y} = x\frac{(d+1)\exp(-d)-1}{d^2}.

The chain rule converts these into sensitivities to the baseline rates: \(\partial F_1/\partial\lambda_{01}=1.5\,\partial F_1/\partial x\) and \(\partial F_1/\partial\lambda_{02}=0.8\,\partial F_1/\partial y\). This makes explicit why a competing cause’s baseline hazard contributes to uncertainty in the target cause’s risk.

Worked Example: A Baseline-Only Interval for One Cause

Consider hypothetical fitted baseline rates \(\widehat{\lambda}_{01}=0.12\) and \(\widehat{\lambda}_{02}=0.08\) per year for the one-year interval, with the profile multipliers above. Suppose their estimated covariance matrix is

\widehat{\Sigma}_0= \begin{pmatrix} 0.0004 & 0.0001\\ 0.0001 & 0.0009 \end{pmatrix}.

The profile-specific rates are \(x=1.5(0.12)=0.18\) per year and \(y=0.8(0.08)=0.064\) per year, so \(d=0.244\) per year. The predicted cause 1 risk is

\widehat{F}_1 = \frac{0.18}{0.244}\{1-\exp(-0.244)\} \approx 0.1597.

For the gradient, substitute \(x=0.18\), \(y=0.064\), and \(d=0.244\) into the derivative formulas:

\begin{aligned} g_1 &= 1.5\left[ \frac{0.064\{1-\exp(-0.244)\}}{0.244^2} + \frac{0.18}{0.244}\exp(-0.244) \right] \approx 1.2161,\\ g_2 &= 0.8(0.18) \frac{(0.244+1)\exp(-0.244)-1}{0.244^2} \approx -0.06129. \end{aligned}

The baseline-only variance includes both marginal variances and their covariance:

\begin{aligned} \widehat{\operatorname{Var}}_0(\widehat{F}_1) &=g_1^2(0.0004)+g_2^2(0.0009) +2g_1g_2(0.0001)\\ &\approx 0.0005800, \qquad SE_0=\sqrt{0.0005800}\approx 0.02408. \end{aligned}

A normal-approximation 95% interval is \(0.1597\pm1.96(0.02408)\), or approximately \((0.1125,\,0.2069)\). In context, this is an approximate interval for the one-year cause 1 risk at the stated profile, accounting for estimated baseline rates while holding the regression multipliers fixed. It is not a full interval for all model uncertainty.

Keeping the Covariance and Defining the Scope

Baseline-hazard estimates are generally not independent. Estimates across causes can be related because they use overlapping risk sets, and increments at different times can also be dependent. The covariance term in the variance formula can increase or decrease the result, depending on the signs of the sensitivities and covariances. In the example, \(g_1\) is positive and \(g_2\) is negative, so the positive covariance contributes a negative term to the approximate variance.

For a prediction using many baseline increments, the same rule applies: collect all relevant estimates in a vector, calculate the gradient of the predicted risk with respect to that vector, and use the full covariance matrix. The gradient can be obtained by differentiating the prediction formula or by numerical perturbations. If numerical differences are used, perturb each component on a scale appropriate to that component and check that smaller perturbations give similar derivatives.

Important distinction: A baseline-only variance treats regression coefficients as fixed. A full-model variance instead propagates uncertainty in both the baseline estimates and the regression coefficients, including their cross-covariances when available. The two quantities answer different questions and should be labelled accordingly.

In a joint calculation, place the baseline quantities and regression coefficients together in a parameter vector \(\widehat{\boldsymbol{\theta}}\), use its joint covariance matrix \(\widehat{\Sigma}\), and differentiate the predicted risk with respect to every component. The resulting variance is \(\nabla g^\mathsf{T}\widehat{\Sigma}\nabla g\). Simply adding baseline-only and coefficient-only variances is justified only when the relevant covariance contributions are zero or deliberately ignored.

Worked Example: Adding Coefficient Uncertainty

Continue the hypothetical calculation above. Suppose a scalar coefficient also affects the prediction, and the fitted model’s derivative of risk with respect to that coefficient is 0.080. For illustration, assume its estimated variance is 0.040 and that its covariance with both baseline estimates is zero. These assumptions make the joint variance the sum of the baseline contribution and the coefficient contribution.

The coefficient contribution is \(0.080^2(0.040)=0.000256\). Adding it to the previously computed baseline contribution gives

\widehat{\operatorname{Var}}(\widehat{F}_1) =0.0005800+0.000256 =0.0008360, \qquad SE=\sqrt{0.0008360}\approx0.02891.

The corresponding normal-approximation interval is \(0.1597\pm1.96(0.02891)\), approximately \((0.1031,\,0.2164)\). Under the stated zero-covariance assumption, this interval is wider because it includes an additional source of sampling variability. If the coefficient and baseline estimates are correlated, their cross-covariance terms must also be included; the simple sum would then be incomplete.

Worked Example: Estimating Baseline Covariance From Paired Replicates

Suppose three hypothetical paired resampling fits give baseline-rate estimates \((0.10,0.07)\), \((0.12,0.08)\), and \((0.14,0.09)\), in inverse years. Their means are \(0.12\) and \(0.08\). The deviations from those means are \((-0.02,-0.01)\), \((0,0)\), and \((0.02,0.01)\).

Using the sample covariance divisor \(3-1=2\), the estimated variances are \(0.0008/2=0.0004\) and \(0.0002/2=0.0001\). Their covariance is \(\{(-0.02)(-0.01)+(0)(0)+(0.02)(0.01)\}/2=0.0002\). Thus, for these replicates, the covariance matrix is

\widehat{\Sigma}_0= \begin{pmatrix} 0.0004 & 0.0002\\ 0.0002 & 0.0001 \end{pmatrix}.

Using the one-interval sensitivities \(g_1=1.2161\) and \(g_2=-0.06129\), the propagated variance is \(1.2161^2(0.0004)+(-0.06129)^2(0.0001)+2(1.2161)(-0.06129)(0.0002)\approx0.0005621\). Thus \(SE_0\approx\sqrt{0.0005621}=0.02371\). This calculation shows why paired estimates matter: calculating two marginal variances alone would discard the covariance information carried by the same resampling fits.

Diagnostics, Alternatives, and Limitations

Before reporting a propagated interval, check that its inputs and scope match the prediction:

Check the parameterization. A covariance matrix for cumulative-hazard increments cannot be used as though it were for hazard rates without the appropriate transformation.
Check the covariance matrix. It should be symmetric and positive semidefinite up to numerical tolerance. A materially negative propagated variance signals a coding, scaling, or covariance-estimation problem.
Check the gradient. Compare analytic derivatives with finite differences where feasible. A sign error can be especially consequential for competing causes, whose sensitivities need not be positive.
Check the approximation. The delta method is a local, first-order approximation. Strong curvature, sparse events, or highly variable baseline estimates can make a symmetric normal interval inaccurate.
Check the interval’s scope and scale. A normal interval may extend outside the possible risk range. That warns that the approximation may be unsuitable; it does not justify interpreting impossible endpoints as valid probabilities.

A baseline-only delta-method interval is useful when the goal is specifically to isolate the effect of baseline estimation, or when other parameters are being treated as fixed by design. A joint delta-method interval is more appropriate when sampling uncertainty in all fitted parameters is the target and a reliable joint covariance estimate is available. A bootstrap can propagate uncertainty through the complete fitting and prediction procedure, including nonlinearities that a first-order approximation may miss, provided each resample refits the relevant model and preserves the sampling structure. These approaches are alternatives with different targets and assumptions, not interchangeable labels for the same interval.

Limitations remain even when the covariance calculation is correct. The interval reflects sampling variation under the fitted model and the assumptions supporting its estimation; it does not automatically account for model misspecification, unmeasured predictors, or changes in the population where predictions will be used. A baseline-only interval can be narrow simply because it excludes coefficient uncertainty. State explicitly which components were varied and which were held fixed.

Common Mistake

A common mistake is to treat the estimated baseline hazards as known once they have been fitted, then propagate only coefficient uncertainty—or to add baseline and coefficient variances without considering their dependence. The correction is to identify the complete parameter vector relevant to the stated inferential target, use its covariance structure, and label a deliberately partial calculation as baseline-only or coefficient-only. For a baseline-only calculation, include covariance among the baseline estimates; for a full-model calculation, include cross-covariances as well.

Key takeaway: Predicted competing-risk incidence is a function of all relevant cause-specific hazards. Baseline uncertainty enters through the gradient of that function and the full covariance matrix of the estimated baseline quantities; the resulting interval is only as comprehensive as the uncertainty components it includes.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions.

Why can uncertainty in a competing cause’s baseline hazard affect the predicted cumulative incidence of the target cause?
In the one-interval example, why is the covariance term \(2g_1g_2\operatorname{Cov}(\widehat{\lambda}_{01},\widehat{\lambda}_{02})\) negative when the covariance is positive?
What does a baseline-only prediction interval hold fixed, and what does it vary?
A propagated delta-method variance is negative. What should be investigated before reporting an interval?

Solutions. (1) The competing hazard changes the total hazard and therefore the probability of remaining event-free and the allocation of events across causes. (2) The cause 1 sensitivity is positive and the cause 2 sensitivity is negative, so their product is negative; multiplying it by a positive covariance gives a negative covariance contribution. (3) It holds the regression coefficients and specified covariate profile fixed, while accounting for sampling uncertainty in the estimated baseline-hazard quantities and their covariance. (4) Check the covariance matrix, parameter ordering and scaling, gradient calculations, and numerical implementation; a valid covariance matrix and gradient cannot produce a negative quadratic-form variance except for negligible numerical error.