Modeling of semi-competing risks by means of first passage times of a stochastic process

Abstract

In semi-competing risks one considers a terminal event, such as death of a person, and a non-terminal event, such as disease recurrence. We present a model where the time to the terminal event is the first passage time to a fixed level c in a stochastic process, while the time to the non-terminal event is represented by the first passage time of the same process to a stochastic threshold S, assumed to be independent of the stochastic process. In order to be explicit, we let the stochastic process be a gamma process, but other processes with independent increments may alternatively be used. For semi-competing risks this appears to be a new modeling approach, being an alternative to traditional approaches based on illness-death models and copula models. In this paper we consider a fully parametric approach. The likelihood function is derived and statistical inference in the model is illustrated on both simulated and real data.

Appendix

1.1 Joint density of $(T_d,T_c)$ for $d<c$

We will first calculate the joint density of $(T_d,D(T_d))$. This is done by integration of the joint density of $(T_d,D(T_d^-),D(T_d))$, which is given in Kahle et al. (2016, Theorem 2.37). The result is, in our notation,

$$\begin{aligned} \tilde{f}(t,z)= & {} \frac{e^{-z}v'(t)}{\Gamma (v(t))} \int _0^d \frac{x^{v(t)-1}}{z-x} dx \\= & {} \frac{e^{-z}v'(t)}{\Gamma (v(t))} \cdot \frac{d^{v(t)}\left( v(t)d \; _2F_1(1,v(t)+1;v(t)+2;d/z) + v(t)z+z \right) }{v(t)(v(t)+1)z^2} \end{aligned}$$

From this we get the joint density of $(T_d,T_c)$ for $d < c$,

$$\begin{aligned} f(t_1,t_2&;&v(t_1),v(t_2),d,c)dt_1dt_2 = P(t_1 \le T_d \le t_1 + dt_1, t_2 \le T_c \le t_2 + dt_2) \\= & {} \int _{z=d}^c P(t_1 \le T_d \le t_1 + dt_1, z \le D(T_d) \le z +dz, t_2 \le T_c \le t_2 + dt_2) \\= & {} \left[ \int _{z=d}^c \tilde{f}(t,z) P(t_2 \le T_c \le t_2 + dt_2|D(T_d)=z) dz \right] dt_1\\= & {} \left[ \int _{z=d}^c \tilde{f}(t,z) f(t_2 -t_1; v(t_2) - v(t_1),c -z)dz \right] dt_1 dt_2 \end{aligned}$$

There is, however, also a possibility that the process crosses both d and c at the same time, giving a tie between $T_d$ and $T_c$. In this case, the relevant density is

$$\begin{aligned} f(t_1,t_1&;&v(t_1),v(t_1),d,c)dt_1 = P(t_1 \le T_d \le t_1 + dt_1, t_1 \le T_c \le t_1 + dt_1) \\= & {} P(t_1 \le T_d \le t_1 + dt_1, D(T_d)>c) \\= & {} \left[ \int _c^\infty \tilde{f}(t_1,z)dz \right] dt_1 \end{aligned}$$

In our computer calculations we have used the simplifying approximation of the joint density of $(T_d,T_c)$, which assumes that $D(T_d)=d$ and $D(T_c)=c$. In this case we have, for $d<c$,

$$\begin{aligned} f(t_1,t_2&;&v(t_1),v(t_2),d,c)dt_1dt_2 = P(t_1 \le T_d \le t_1 + dt_1, t_2 \le T_c \le t_2 + dt_2) \\= & {} P(t_1 \le T_d \le t_1 + dt_1)P(t_2 \le T_c \le t_2 + dt_2|T_d = t_1)\\= & {} f(t_1; v(t_1), d)dt_1f(t_2 -t_1; v(t_2) - v(t_1),c -d)dt_1dt_2 \end{aligned}$$

1.2 Identifiability of the model

We first prove identifiability of the parameters $c,\alpha ,\beta$ from the distribution of X. Note that we have $X=T_c$. Thus, from (4) we have,

$$\begin{aligned} P(X>t) = P(T_c > t) = \gamma (v(t),c)/\Gamma (v(t)) \end{aligned}$$

(11)

where $\gamma (a,c) = \int _0^c z^{a-1}e^{-z}dz$. We first show as a digression that if c is unknown, then the function v(t) is not nonparametrically identifiable. This follows since in the right hand expression of (11) we may for any given $c>0$ solve for v(t) for each fixed t. To see this, note from (11) that $P(X>t)$ equals $P(W>c)$ where $W \sim$ gamma(v(t), 1). Since the gamma distribution is stochastically increasing in the shape parameter, here v(t), we may always adjust the v(t) to get a given value for $P(W>c)$.

Thus suppose instead that (10) holds. Now we use a result from Temme (1975) to see that as $a \rightarrow \infty$ we have

$$\begin{aligned} \gamma (a,c)/\Gamma (a) \sim \frac{c^a e^{-c}}{\Gamma (1+a)} \sim (2\pi a)^{-1/2} e^{a-c} \left( \frac{c}{a}\right) ^{a} . \end{aligned}$$

Here the last expression is obtained by using Stirling’s formula. Letting $a=\alpha t^\beta$ and taking the logarithm we get for large t,

$$\begin{aligned} \log P(X>t)\sim & {} -(1/2) \log 2\pi -(1/2) \log \alpha -(1/2) \beta \log t + \alpha t^\beta - c + \alpha t^\beta \log c \nonumber \\- & {} \alpha t^\beta \log \alpha - \alpha t^\beta \beta \log t \end{aligned}$$

(12)

Suppose now there is another combination of $c, \alpha , \beta$, denoted $c^*,\alpha ^*,\beta ^*$, for which the same $P(X>t)$ is obtained for all t. Then letting $t \rightarrow \infty$, the dominant term in (12) is $\alpha t^\beta \log t$ which hence must be equal for the two parametrizations, implying $\beta =\beta ^*$ and hence also $\alpha =\alpha ^*$. Finally, this clearly implies $c=c^*$ and we are done.

For identifiability of the full threshold model, it remains to show that the distribution of S conditional on $S<c$ is identifiable when c and the parameters of the process D(t) are given. We are in fact able to show that this distribution is nonparametrically identifiable for any given v(t) and c, which we for simplicity will assume to be strictly increasing and continuous, with $v(0)=0$, $v(\infty )=\infty$. Suppose first that $c=\infty$, so that $T_S$ is always observed. We will show that the distribution of $T_S$ uniquely determines the distribution of S. Now

$$\begin{aligned} P(T_S>t) = P(S>D(t)) = E [P(S>D(t)|D(t))] = E[ \bar{F}_S(W) ] \end{aligned}$$

(13)

where $W \sim$ gamma(v(t), 1) and $\bar{F}_S(s)=P(S>s)$. Since this is to hold for all $t>0$, from the fact that the family of $W \sim$ gamma($\theta ,1)$ is a complete family of distributions, it follows that $\bar{F}_S$ is uniquely given and hence that the distribution of $T_S$ uniquely determines the distribution of S Casella and Berger (2002, Chapter 6.2).

Next, for a given $c < \infty$ we need to show that the (observable) distribution $P(T_S>t|S<c)$ uniquely determines the distribution $P(S>s|S<c)$. This follows directly from the above argument which had $c=\infty$ by considering only distributions for S with support in (0, c).

Note finally that (13) is a result of interest in itself if the distribution of the threshold S is given and one wants the distribution of $T_S$. Paroissin and Salami (2014) consider the cases where S is, respectively, exponentially and gamma distributed.

1.3 Nonparametric estimation of crude quantities

Consider competing risks with latent variables X and Z. Suppose that n units are observed, either until (independent) right censoring or until time $T = \min (X,Z)$, whatever comes first. Let $t_{1}< \cdots < t_{k}$ be the sorted event times, i.e., observations of T. Let further $\hat{S}(\cdot )$ be the Kaplan–Meier estimator of the survival function of T. Then the so-called Aalen–Johansen estimator of the sub-distribution functions are (Borgan 1998):

$$\begin{aligned} \hat{F}^{*}_{X}(t)= & {} \sum _{i; t_i \le t} \hat{S}(t_i)\frac{\delta _{iX}}{n_i},\nonumber \\ \hat{F}^{*}_{Z}(t)= & {} \sum _{i; t_i \le t} \hat{S}(t_i)\frac{\delta _{iZ}}{n_i}. \end{aligned}$$

(14)

Here $n_i$ is the number at risk at time $t_{i}$ while $\delta _{iX} = 1$ ($\delta _{iZ} = 1$) if the observation at time $t_i$ is an X (Z). The natural estimates of the conditional sub-distribution functions $\tilde{F}_X(t)$ and $\tilde{F}_Z(t)$ are hence

$$\begin{aligned} \hat{\tilde{F}}_X(t) = \frac{\hat{F}_X^{*}(t)}{\hat{F}_X^{*}(\infty )} \quad \text {and} \quad \hat{\tilde{F}}_Z(t) = \frac{\hat{F}_Z^{*}(t)}{\hat{F}_Z^{*}(\infty )}. \end{aligned}$$

(15)

With the same notation we have nonparametric estimates of the cumulative cause-specific hazard functions for the two risks given by

$$\begin{aligned} \hat{\Lambda }^{*}_{X}(t)= & {} \sum _{i; t_i \le t}\frac{\delta _{iX}}{n_i},\nonumber \\ \hat{\Lambda }^{*}_{Z}(t)= & {} \sum _{i; t_i \le t}\frac{\delta _{iZ}}{n_i}. \end{aligned}$$

(16)