A comparison of statistical methods to predict the residual lifetime risk

Abstract

Lifetime risk measures the cumulative risk for developing a disease over one’s lifespan. Modeling the lifetime risk must account for left truncation, the competing risk of death, and inference at a fixed age. In addition, statistical methods to predict the lifetime risk should account for covariate-outcome associations that change with age. In this paper, we review and compare statistical methods to predict the lifetime risk. We first consider a generalized linear model for the lifetime risk using pseudo-observations of the Aalen-Johansen estimator at a fixed age, allowing for left truncation. We also consider modeling the subdistribution hazard with Fine-Gray and Royston-Parmar flexible parametric models in left truncated data with time-covariate interactions, and using these models to predict lifetime risk. In simulation studies, we found the pseudo-observation approach had the least bias, particularly in settings with crossing or converging cumulative incidence curves. We illustrate our method by modeling the lifetime risk of atrial fibrillation in the Framingham Heart Study. We provide technical guidance to replicate all analyses in R.

Appendix

Variance of non-parametric estimator of lifetime risk

The variance of the lifetime risk is estimated using Greenwood’s formula, [4, 18]

$$\begin{aligned}&\widehat{\mathrm {var}} \left\{ \hat{F_1}(\tau | T>\tau _0) \right\} \\&\quad = \sum _{\tau _0< t_j \le \tau } \left\{ {\hat{S}}(t_{j-1}) \frac{\varDelta N_1(t_j)}{Y(t_j)} \right\} ^2 \\&\quad \quad \times \left\{ \frac{Y(t_j) - \varDelta N_1(t_j)}{Y(t_j) \varDelta N_1(t_j)} + \sum _{t_l< t_{t_{j}}} \frac{\varDelta N(t_l)}{Y(t_l) (Y(t_l) - \varDelta N(t_l))} \right\} \\&\qquad + 2 \sum _{\tau _0< t_j< \tau } \sum _{ t_k \in {(}t_j, \tau {]}} {\hat{S}}(t_{j-1}) \frac{\varDelta N_1(t_j)}{Y(t_j)} {\hat{S}}(t_{k-1}) \\&\quad \quad \times \frac{\varDelta N_1(t_k)}{Y(t_k)} \left\{ -\frac{1}{Y(t_j)} + \sum _{t_l<t_j} \frac{\varDelta N(t_l)}{Y(t_l) (Y(t_l) - \varDelta N(t_l))} \right\} . \end{aligned}$$

Weighted Breslow estimator of the baseline cumulative subdistribution hazard

The weighted Breslow estimator of the baseline cumulative subdistribution hazard is given by

$$\begin{aligned} {\hat{\varLambda }}_{1,0}(t | T_i>\tau _0)&= \int _{\tau _0}^{\tau } {\hat{\lambda }}_{1,0}(t | T_i>\tau _0) dt \\&= \frac{1}{n} \sum _{i=1}^n \int _{\tau _0}^{\tau } \frac{w_i (u)}{{\hat{S}}^0(\hat{\varvec{\pi }}, u)} dN_i(u) \end{aligned}$$

where ${\hat{S}}^0(\hat{\varvec{\pi }}, u) = \sum _{i=1}^n w_i(u) T_i(u) \exp \{ \hat{\varvec{\pi }}^T Z_i(u) \}$ [11].

Illustration of non-proportional subdistribution hazards with proportional cause-specific hazards

Let the cause-specific hazards take Weibull form, $\alpha _1(t; Z)=\frac{a_1}{b_1^{a_1}}t^{a_1 - 1} e^{\gamma _1 Z}$ and $\alpha _2(t; Z)=\frac{a_2}{b_2^{a_2}}t^{a_2 - 1} e^{\gamma _2 Z}$, where Z is a binary covariate with a proportional effect on both cause-specific hazards. The all-cause cumulative hazard is $A(t; Z)=\frac{t}{b_1}^{a_1} e^{\gamma _1 Z} + \frac{t}{b_2}^{a_2} e^{\gamma _1 Z}$. The proportional effect of Z on the cause-specific hazard can be verified by deriving the cause-specific hazard ratio (CSHR),

$$\begin{aligned} \frac{\alpha _1(t;Z=1)}{\alpha _1(t;Z=0)}&= \frac{\frac{a_1}{b_1^{a_1}}t^{a_1 - 1} e^{\gamma _1}}{\frac{a_1}{b_1^{a_1}}t^{a_1 - 1}}\\&=e^{\gamma _1}. \end{aligned}$$

The subdistribution hazard can be represented as a function of the cause-specific hazard, $\lambda _1(t) = \alpha _1(t)/\big ( 1 + \frac{F_2 (t)}{S(t)} \big )$, where $F_k(t) = \int _0^t S(u) \alpha _k(u) du = 1 - \exp (-\int _0^t \lambda _k (u) du)$. The subdistribution hazard ratio (SHR) is then

$$\begin{aligned} \frac{\lambda _1(t; Z=1)}{\lambda _1(t; Z=0)}&= \frac{ \alpha _1(t;Z=1) / \left( 1 + \frac{F_2 (t;Z=1)}{S(t; Z=1)} \right) }{ \alpha _1(t;Z=0) / \left( 1 + \frac{F_2 (t; Z=0)}{S(t; Z=0)} \right) } \\&= e^{\gamma _1} \cdot \left( \frac{1 + \frac{F_2 (t; Z=0)}{S(t; Z=0)} }{ 1 + \frac{F_2 (t; Z=1)}{S(t; Z=1)}} \right) . \end{aligned}$$

The SHR is not constant, but equal to the CSHR multiplied by a function of time.

Supplemental tables and figures

Table 8 Participant characteristics at index age 55 in the Framingham Heart Study

Table 9 Sensitivity analysis: predicted differences in lifetime risk of atrial fibrillation at age 95, from index age 55 using BIC

Table 10 Multivariable models for incident atrial fibrillation in the Framingham Heart Study, from index age 55 to age 95

Table 11 Multivariable models for incident death without atrial fibrillation in the Framingham Heart Study, from index age 55 to age 95

Fig. 4

True cumulative incidence functions in simulation study

Fig. 5

Nested loop plot showing the simulation study results: relative bias

Fig. 6

Nested loop plot showing the simulation study results: coverage

Fig. 7

Nested loop plot showing the simulation study results: type I error and power of pseudo-observation method