This package enables us to simulate survival data with several levels of complexity from different survival models: Cox model (Cox 1972), the Accelerated Failure Time (AFT) model, and the Accelerated Hazard (AH) model (Chen and Wang 2000). To simulate data from a Cox model, we consider the procedure of (Bender, Augustin, and Blettner 2005). This model is the most popular in survival analysis, but other models exist. The package also enables us to simulate survival data from an AFT model based on (Leemis, Shih, and Reynertson 1990). Considering different models for data simulation is interesting because the assumptions associated with these models are different. Indeed, the Cox model is a proportional risk model, not the AFT model. In the AFT model, the variables will have an accelerating or decelerating effect on on individuals’ survival. Despite the survival functions of an AFT model never intersect as in the Cox model. However, having data with intersecting survival curves allows the data to be more complex and makes it more difficult for methods to predict survival. Besides, Cox and AFT models produce data with survival curves not crossing. To have crossing survival curves, we considered two approaches. The first approach consists of modifying the AFT model to have intersecting survival curves. The second approach concerns using an AH model to generate the survival data. The AH model is more flexible than the two models mentioned above. In the AH model, the variables will accelerate or decelerate the instantaneous risk of death. The survival curves of the AH model can therefore cross each other. This package also allows us to simulate survival data from modified AFT and AH models. The generation of survival times carries out from the different models mentioned above. The models’ baseline risk function of the models is assumed to be known and follows a specific probability distribution.
1.2 Survival times generation
The models considered are composed of one function, \(\alpha_0(t)\) the baseline risk and a parameter \(\beta,\) reflecting the variables’ effect on survival times. For data generation, we assume that the baseline risk function \(\alpha_0(t)\) is known and therefore follows a probability distribution. For this initial version of the package, the baseline hazard distributions are Weibull or log-normale. We summarize the characteristics of the distributions (Weibull/Log-normale) in the table below.
| Parameters |
\[\lambda > 0\] \[a > 0\] |
\[\mu \in ]-\infty,+\infty[\] \[\sigma > 0\] |
| Support |
\[\mathbb{R}^+\] |
\[\mathbb{R}^+_{\star}\] |
| Baseline hazard |
\[\alpha_0(t) = \lambda a t^{(a-1)}\] |
\[\alpha_0(t) = \frac{\frac{1}{\sigma\sqrt{2\pi t}} \exp\left[-\frac{(\log t - \mu)^2 }{2 \sigma^2}\right]}{1 - \Phi\left[\frac{\log t - \mu}{\sigma}\right]}\] |
| Cumulative Hazard |
\[H_0(t) = \lambda t^{a}\] |
\[H_0(t) = - \log(1 - \Phi\left[\frac{\log t - \mu}{\sigma}\right])\] |
| Inverse cumulative hazard |
\[H_0^{-1}(u) = \left( \frac{u}{\lambda} \right)^{1/a}\] |
\[H_0^{-1}(u) = \exp(\sigma\Phi^{-1}(1-\exp(-u))+\mu)\] |
| Density function |
\[f(t) = \lambda a t^{(a-1)} \exp(-\lambda t^{a})\] |
\[f(t) = \exp\left[-\frac{(\log t - \mu)^2 }{2 \sigma^2}\right] \frac{1}{\sigma t \sqrt{2\pi }}\] |
| Cumulative function |
\[F(t) = \exp(-\lambda t^{a})\] |
\[F(t) = 1 - \Phi\left[\frac{\log t - \mu}{\sigma}\right]\] |
| Expectation |
\[\mathbb{E}(T) = \Gamma(\frac{1}{a} + 1) \frac{1}{\sqrt[a\,]{\lambda}}\] |
\[\mathbb{E}(T) = \exp (\mu + \frac{\sigma^2}{2})\] |
| Variance |
\[\mathbb{V}(T) = \left[ \Gamma(\frac{2}{a} + 1) - \Gamma^2(\frac{1}{a} + 1)\right] \frac{1}{\sqrt[a\,]{\mu^2} }\] |
\[ \mathbb{V}(T) = (\exp(\sigma^2) -1) \exp(2\mu+\sigma^2)\] |
with \(\Gamma\) is the gamma function and \(\Phi\) is the cumulative distribution function of the standard normal distribution.
The distribution function is deduced from the survival function from the following formula:
\[\begin{equation}
F(t\mid X) = 1 - S(t\mid X).
\end{equation}\]
For data generation, if
\(Y\) is a random variable that follows a probability distribution
\(F,\) then
\(U = F(Y)\) follows a uniform distribution over the interval
\([0,1],\) and
\((1-U)\) also follows a uniform distribution
\(\mathcal{U}[0,1].\) We finally obtain that:
\[\begin{align}
1 - U &= S(t\mid X) \\ %\sim \mathcal{U} [0,1] \\
&= \exp(-H_0(\psi_1(X)t)\psi_2(X)) % \sim \mathcal{U} [0,1].
\end{align}\]
If
\(\alpha_0(t)\) is strictly positive for any
\(t,\) then
\(H_0(t)\) can be inverted and the survival time of each of the considered models (Cox, AFT, and AH) expresses for
\(H_0^{-1}(u)\). The expression of the survival times for each of the models writes in a general way:
\[\begin{equation}
T = \frac{1}{\psi_1(X)} H^{-1}_0 \left( \frac{\log(1-U)}{\psi_2(X)} \right), \text{with}
\end{equation}\]
\[\begin{equation}
(\psi_1(X), \psi_2(X)) = \left\{
\begin{array}{ll}
(1, \exp(\beta^TX)) & \mbox{for the Cox model } \\
(\exp(\beta^TX), \exp(-\beta^TX)) & \mbox{for the AH model} \\\
(\exp(\beta^TX), 1) & \mbox{for the AFT model. }
\end{array}
\right.
\end{equation}\]
Two distributions are proposed for the cumulative risk function \(H_0(t)\) to generate the survival data. If the survival times are distributed according to a Weibull distribution \(\mathcal{W}(a, \lambda),\) the baseline risk is of the form:
\[\begin{equation}
\alpha_0(t) = a\lambda t^{a-1}, \lambda > 0, a > 0.
\end{equation}\]
The cumulative risk function is therefore written as follows:
\[\begin{equation}
H_0(t) = \lambda t^{a}, \lambda > 0, a > 0
\end{equation}\]
and the inverse of this function is expressed as follows:
\[\begin{equation}
H_0^{-1}(u) = \left( \frac{u}{\lambda} \right)^{1/a}.
\end{equation}\]
In a second step, we considered that the survival times followed a log-normal
\(\mathcal{LN}(\mu, \sigma)\) distribution of mean
\(\mu\) and standard deviation
\(\sigma\). The basic risk function is therefore written as:
\[\begin{equation}
\alpha_0(t) = \frac{\frac{\frac{1}{\sigma\sqrt{2\pi t}}} \exp\left[-\frac{(\log t - \mu)^2 }{2 \sigma^2}\right]}{1 - \Phi\left[\frac{\log t - \mu}{\sigma}\right]},
\end{equation}\]
with
\(\Phi(t)\) the distribution function of a centered and reduced normal distribution. The cumulative risk function is written as:
\[\begin{equation}
H_0(t) = - \log\left [1 - \Phi\left(\frac{\log t - \mu}{\sigma}\right)\right]
\end{equation}\]
and therefore the inverse of this function is expressed by:
\[\begin{equation}
H_0^{-1}(u) = \exp(\sigma\Phi^{-1}(1-\exp(-u))+\mu),
\end{equation}\]
with \(\Phi^{-1}(t)\) the inverse of the distribution function of a centered and reduced normal distribution.
