The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.
The aim of the R package survstan is to provide a toolkit for fitting survival models using Stan.
The R package survstan can be used to fit right-censored survival data under independent censoring. The implemented models allow the fitting of survival data in the presence/absence of covariates. All inferential procedures are currently based on the maximum likelihood (ML) approach.
You can install the released version of survstan from CRAN with:
You can install the development version of survstan from GitHub with:
Let \((t_{i}, \delta_{i})\) be the observed survival time and its corresponding failure indicator, \(i=1, \cdots, n\), and \(\boldsymbol{\theta}\) be a \(k \times 1\) vector of parameters. Then, the likelihood function for right-censored survival data under independent censoring can be expressed as:
\[ L(\boldsymbol{\theta}) = \prod_{i=1}^{n}f(t_{i}|\boldsymbol{\theta})^{\delta_{i}}S(t_{i}|\boldsymbol{\theta})^{1-\delta_{i}}. \]
The maximum likelihood estimate (MLE) of \(\boldsymbol{\theta}\) is obtained by directly maximization of \(\log(L(\boldsymbol{\theta}))\) using the rstan::optimizing() function. The function rstan::optimizing() further provides the hessian matrix of \(\log(L(\boldsymbol{\theta}))\), needed to obtain the observed Fisher information matrix, which is given by:
\[ \mathscr{I}(\hat{\boldsymbol{\theta}}) = -\frac{\partial^2}{\partial \boldsymbol{\theta}\boldsymbol{\theta}'} \log L(\boldsymbol{\theta})\mid_{\boldsymbol{\theta}=\hat{\boldsymbol{\theta}}}, \]
Inferences on \(\boldsymbol{\theta}\) are then based on the asymptotic properties of the MLE, \(\hat{\boldsymbol{\theta}}\), that state that:
\[ \hat{\boldsymbol{\theta}} \asymp N_{k}(\boldsymbol{\theta}, \mathscr{I}^{-1}(\hat{\boldsymbol{\theta}})). \]
Some of the most popular baseline survival distributions are implemented in the R package survstan. Such distributions include:
The parametrizations adopted in the package survstan are presented next.
If \(T \sim \mbox{Exp}(\lambda)\), then
\[ f(t|\lambda) = \lambda\exp\left\{-\lambda t\right\}I_{[0, \infty)}(t), \] where \(\lambda>0\) is the rate parameter.
The survival and hazard functions in this case are given by:
\[ S(t|\lambda) = \exp\left\{-\lambda t\right\} \] and \[ h(t|\lambda) = \lambda. \]
If \(T \sim \mbox{Weibull}(\alpha, \gamma)\), then
\[ f(t|\alpha, \gamma) = \frac{\alpha}{\gamma^{\alpha}}t^{\alpha-1}\exp\left\{-\left(\frac{t}{\gamma}\right)^{\alpha}\right\}I_{[0, \infty)}(t), \] where \(\alpha>0\) and \(\gamma>0\) are the shape and scale parameters, respectively.
The survival and hazard functions in this case are given by:
\[ S(t|\alpha, \gamma) = \exp\left\{-\left(\frac{t}{\gamma}\right)^{\alpha}\right\} \] and \[ h(t|\alpha, \gamma) = \frac{\alpha}{\gamma^{\alpha}}t^{\alpha-1}. \]
If \(T \sim \mbox{LN}(\mu, \sigma)\), then
\[ f(t|\mu, \sigma) = \frac{1}{\sqrt{2\pi}t\sigma}\exp\left\{-\frac{1}{2}\left(\frac{log(t)-\mu}{\sigma}\right)^2\right\}I_{[0, \infty)}(t), \] where \(-\infty < \mu < \infty\) and \(\sigma>0\) are the mean and standard deviation in the log scale of \(T\).
The survival and hazard functions in this case are given by:
\[S(t|\mu, \sigma) = \Phi\left(\frac{-log(t)+\mu}{\sigma}\right)\] and \[h(t|\mu, \sigma) = \frac{f(t|\mu, \sigma)}{S(t|\mu, \sigma)},\] where \(\Phi(\cdot)\) is the cumulative distribution function of the standard normal distribution.
If \(T \sim \mbox{LL}(\alpha, \gamma)\), then
\[ f(t|\alpha, \gamma) = \frac{\frac{\alpha}{\gamma}\left(\frac{t}{\gamma}\right)^{\alpha-1}}{\left[1 + \left(\frac{t}{\gamma}\right)^{\alpha}\right]^2}I_{[0, \infty)}(t), ~ \alpha>0, \gamma>0, \]
where \(\alpha>0\) and \(\gamma>0\) are the shape and scale parameters, respectively.
The survival and hazard functions in this case are given by:
\[S(t|\alpha, \gamma) = \frac{1}{1+ \left(\frac{t}{\gamma}\right)^{\alpha}}\] and \[ h(t|\alpha, \gamma) = \frac{\frac{\alpha}{\gamma}\left(\frac{t}{\gamma}\right)^{\alpha-1}}{1 + \left(\frac{t}{\gamma}\right)^{\alpha}}. \]
If \(T \sim \mbox{Gamma}(\alpha, \lambda)\), then
\[f(t|\alpha, \lambda) = \frac{\lambda^{\alpha}}{\Gamma(\alpha)}t^{\alpha-1}\exp\left\{-\lambda t\right\}I_{[0, \infty)}(t),\]
where \(\Gamma(\alpha) = \int_{0}^{\infty}u^{\alpha-1}\exp\{-u\}du\) is the gamma function.
The survival function is given by
\[S(t|\alpha, \lambda) = 1 - \frac{\gamma^{*}(\alpha, \lambda t)}{\Gamma(\alpha)},\] where \(\gamma^{*}(\alpha, \lambda t)\) is the lower incomplete gamma function, which is available only numerically. Finally, the hazard function is expressed as:
\[h(t|\alpha, \lambda) = \frac{f(t|\alpha, \lambda)}{S(t|\alpha, \lambda)}.\]
If \(T \sim \mbox{ggstacy}(\alpha, \gamma, \kappa)\), then
\[f(t|\alpha, \gamma, \kappa) = \frac{\kappa}{\gamma^{\alpha}\Gamma(\alpha/\kappa)}t^{\alpha-1}\exp\left\{-\left(\frac{t}{\gamma}\right)^{\kappa}\right\}I_{[0, \infty)}(t),\] for \(\alpha>0\), \(\gamma>0\) and \(\kappa>0\).
It can be show that the survival function can be expressed as:
\[S(t|\alpha, \gamma, \kappa) = S_{G}(x|\nu, 1),\] where \(x = \displaystyle\left(\frac{t}{\gamma}\right)^\kappa\), and \(F_{G}(\cdot|\nu, 1)\) corresponds to the distribution function of a gamma distribution with shape parameter \(\nu = \alpha/\gamma\) and scale parameter equals to 1.
Finally, the hazard function is expressed as:
\[h(t|\alpha, \gamma, \kappa) = \frac{f(t|\alpha, \gamma, \kappa)}{S(t|\alpha, \gamma, \kappa)}.\]
If \(T \sim \mbox{ggprentice}(\mu, \sigma, \varphi)\), then
\[f(t | \mu, \sigma, \varphi) = \begin{cases} \frac{|\varphi|(\varphi^{-2})^{\varphi^{-2}}}{\sigma t\Gamma(\varphi^{-2})}\exp\{\varphi^{-2}[\varphi w - \exp(\varphi w)]\}I_{[0, \infty)}(t), & \varphi \neq 0 \\ \frac{1}{\sqrt{2\pi}t\sigma}\exp\left\{-\frac{1}{2}\left(\frac{log(t)-\mu}{\sigma}\right)^2\right\}I_{[0, \infty)}(t), & \varphi = 0 \end{cases} \] where \(w = \frac{\log(t) - \mu}{\sigma}\), for \(-\infty < \mu < \infty\), \(\sigma>0\) and \(-\infty < \varphi < \infty\)$.
It can be show that the survival function can be expressed as:
\[ S(t|\mu, \sigma, \varphi) = \begin{cases} S_{G}(x|1/\varphi^2, 1), & \varphi > 0 \\ 1-S_{G}(x|1/\varphi^2, 1), & \varphi < 0 \\ S_{LN}(x|\mu, \sigma), & \varphi = 0 \end{cases} \] where \(x = \frac{1}{\varphi^2}\exp\{\varphi w\}\), \(S_{G}(\cdot|1/\varphi^2, 1)\) is the distribution function of a gamma distribution with shape parameter \(1/\varphi^2\) and scale parameter equals to 1, and \(S_{LN}(x|\mu, \sigma)\) corresponds to the survival function of a lognormal distribution with location parameter \(\mu\) and scale parameter \(\sigma\).
Finally, the hazard function is expressed as:
\[h(t|\alpha, \gamma, \kappa) = \frac{f(t|\alpha, \gamma, \kappa)}{S(t|\alpha, \gamma, \kappa)}.\]
If \(T \sim \mbox{Gamma}(\alpha, \gamma)\), then
\[f(t|\alpha, \lambda) = \alpha\exp\left\{\gamma x-\frac{\alpha}{\gamma}\left(e^{\gamma x} - 1\right)\right\}I_{[0, \infty)}(t).\]
The survival and hazard functions are given, respectively, by
\[S(t|\alpha, \lambda) = \exp\left\{-\frac{\alpha}{\gamma}\left(e^{\gamma x} - 1\right)\right\}.\] and
\[h(t|\alpha, \lambda) = \alpha\exp\{\gamma x}.\]
Let \(T \sim \mbox{rayleigh}(\sigma)\), where \(\sigma>0\) is a scale parameter. Then, the density, survival and hazard functions are respectively given by:
\[f(t|\sigma) = \frac{x}{\sigma^2}\exp\left\{-\frac{x^2}{2\sigma^2}\right\},\] \[S(t|\sigma) = \exp\left\{-\frac{x^2}{2\sigma^2}\right\}\] and
\[h(t|\sigma) = \frac{x}{\sigma^2}.\]
If \(T \sim \mbox{fatigue}(\alpha, \gamma)\), then
\[ f(t|\alpha, \gamma) = \frac{\sqrt{\frac{t}{\gamma}}+\sqrt{\frac{\gamma}{t}}}{2 \alpha t}\phi\left(\sqrt{\frac{t}{\gamma}}+\sqrt{\frac{\gamma}{t}}\right)(t), ~ \alpha>0, \gamma>0, \]
where \(\phi(\cdot)\) is the probability density function of a standard normal distribution, \(\alpha>0\) and \(\gamma>0\) are the shape and scale parameters, respectively.
The survival function in this case is given by:
\[ S(t|\alpha, \gamma) =\Phi\left(\sqrt{\frac{t}{\gamma}}-\sqrt{\frac{\gamma}{t}}\right)(t) \],
where \(\Phi(\cdot)\) is the cumulative distribution function of a standard normal distribution. The hazard function is given by \[h(t|\mu, \sigma) = \frac{f(t|\alpha, \gamma)}{S(t|\alpha, \gamma)}. \]
When covariates are available, it is possible to fit six different regression models with the R package survstan:
The regression survival models implemented in the R package survstan are briefly described in the sequel. Denote by \(\mathbf{x}\) a \(1\times p\) vector of covariates, and let \(\boldsymbol{\beta}\) and \(\boldsymbol{\phi}\) be \(p \times 1\) vectors of regression coefficients, and \(\boldsymbol{\theta}\) a vector of parameters associated with some baseline survival distribution. To ensure identifiability of the implemented regression models, we shall assume that the linear predictors \(\mathbf{x} \boldsymbol{\beta}\) and \(\mathbf{x} \boldsymbol{\phi}\) do not include an intercept term.
Accelerated failure time (AFT) models are defined as
\[ T = \exp\{\mathbf{x} \boldsymbol{\beta}\}\nu, \] where \(\nu\) follows a baseline distribution with survival function \(S_{0}(\cdot|\boldsymbol{\theta})\) so that
\[ f(t|\boldsymbol{\theta}, \boldsymbol{\beta}, \mathbf{x}) = e^{-\mathbf{x} \boldsymbol{\beta}}f_{0}(te^{-\mathbf{x} \boldsymbol{\beta}}|\boldsymbol{\theta}) \] and
\[ S(t|\boldsymbol{\theta}, \boldsymbol{\beta}, \mathbf{x}) = S_{0}(t e^{-\mathbf{x} \boldsymbol{\beta}}|\boldsymbol{\theta}). \]
Proportional hazards (PH) models are defined as
\[ h(t|\boldsymbol{\theta}, \boldsymbol{\beta}, \mathbf{x}) = h_{0}(t|\boldsymbol{\theta})\exp\{\mathbf{x} \boldsymbol{\beta}\}, \] where \(h_{0}(t|\boldsymbol{\theta})\) is a baseline hazard function so that
\[ f(t|\boldsymbol{\theta}, \boldsymbol{\beta}, \mathbf{x}) = h_{0}(t|\boldsymbol{\theta})\exp\left\{\mathbf{x} \boldsymbol{\beta} - H_{0}(t|\boldsymbol{\theta})e^{\mathbf{x} \boldsymbol{\beta}}\right\}, \] and
\[ S(t|\boldsymbol{\theta}, \boldsymbol{\beta}, \mathbf{x}) = \exp\left\{ - H_{0}(t|\boldsymbol{\theta})e^{\mathbf{x} \boldsymbol{\beta}}\right\}. \]
Proportional Odds (PO) models are defined as
\[ R(t|\boldsymbol{\theta}, \boldsymbol{\beta}, \mathbf{x}) = R_{0}(t|\boldsymbol{\theta})\exp\{\mathbf{x} \boldsymbol{\beta}\}, \] where \(\displaystyle R_{0}(t|\boldsymbol{\theta}) = \frac{1-S_{0}(t|\boldsymbol{\theta})}{S_{0}(t|\boldsymbol{\theta})} = \exp\{H_{0}(t|\boldsymbol{\theta})\}-1\) is a baseline odds function so that
\[ f(t|\boldsymbol{\theta}, \boldsymbol{\beta}, \mathbf{x}) = \frac{h_{0}(t|\boldsymbol{\theta})\exp\{\mathbf{x} \boldsymbol{\beta} + H_{0}(t|\boldsymbol{\theta})\}}{[1 + R_{0}(t|\boldsymbol{\theta})e^{\mathbf{x} \boldsymbol{\beta}}]^2}. \]
and
\[ S(t|\boldsymbol{\theta}, \boldsymbol{\beta}, \mathbf{x}) = \frac{1}{1 + R_{0}(t|\boldsymbol{\theta})e^{\mathbf{x} \boldsymbol{\beta}}}. \]
Accelerated hazard (AH) models can be defined as
\[h(t|\boldsymbol{\theta}, \boldsymbol{\beta},\mathbf{x}) = h_{0}\left(t/e^{\mathbf{x}\boldsymbol{\beta}}|\boldsymbol{\theta}\right)\]
so that
\[S(t|\boldsymbol{\theta}, \boldsymbol{\beta},\mathbf{x}) = \exp\left\{- H_{0}\left(t/ e^{\mathbf{x}\boldsymbol{\beta}}|\boldsymbol{\theta}\right)e^{\mathbf{x}\boldsymbol{\beta}} \right\} \] and \[f(t|\boldsymbol{\theta}, \boldsymbol{\beta}, \mathbf{x}) = h_{0}\left(t/e^{\mathbf{x}\boldsymbol{\beta}}|\boldsymbol{\theta}\right)\exp\left\{- H_{0}\left(t/ e^{\mathbf{x}\boldsymbol{\beta}}|\boldsymbol{\theta}\right)e^{\mathbf{x}\boldsymbol{\beta}} \right\}. \]
The survival function of the extended hazard (EH) model is given by:
\[S(t|\boldsymbol{\theta},\boldsymbol{\beta}, \boldsymbol{\phi}) = \exp\left\{-H_{0}(t/e^{\mathbf{x}\boldsymbol{\beta}}|\boldsymbol{\theta})\exp(\mathbf{x}(\boldsymbol{\beta} + \boldsymbol{\phi}))\right\}. \]
The hazard and the probability density functions are then expressed as:
\[h(t|\boldsymbol{\theta},\boldsymbol{\beta}, \boldsymbol{\phi}) = h_{0}(t/e^{\mathbf{x}\boldsymbol{\beta}}|\boldsymbol{\theta})\exp\{\mathbf{x}\boldsymbol{\phi}\} \] and
\[f(t|\boldsymbol{\theta},\boldsymbol{\beta}, \boldsymbol{\phi}) = h_{0}(t/e^{\mathbf{x}\boldsymbol{\beta}}|\boldsymbol{\theta})\exp\{\mathbf{x}\boldsymbol{\beta}\}\exp\left\{-H_{0}(t/e^{\mathbf{x}\boldsymbol{\beta}}|\boldsymbol{\theta})\exp(\mathbf{x}(\boldsymbol{\beta}+ \boldsymbol{\phi}))\right\}, \]
respectively.
The EH model includes the AH, AFT and PH models as particular cases when \(\boldsymbol{\phi} = \mathbf{0}\), \(\boldsymbol{\phi} = -\boldsymbol{\beta}\), and \(\boldsymbol{\beta} = \mathbf{0}\), respectively.
The survival function of the Yang and Prentice (YP) model is given by:
\[S(t|\boldsymbol{\theta},\boldsymbol{\beta}, \boldsymbol{\phi}) = \left[1+\frac{\kappa_{S}}{\kappa_{L}}R_{0}(t|\boldsymbol{\theta})\right]^{-\kappa_{L}}. \]
The hazard and the probability density functions are then expressed as:
\[h(t|\boldsymbol{\theta},\boldsymbol{\beta}, \boldsymbol{\phi}) = \frac{\kappa_{S}h_{0}(t|\boldsymbol{\theta})\exp\{H_{0}(t|\boldsymbol{\theta})\}}{\left[1+\frac{\kappa_{S}}{\kappa_{L}}R_{0}(t|\boldsymbol{\theta})\right]} \] and
\[f(t|\boldsymbol{\theta},\boldsymbol{\beta}, \boldsymbol{\phi}) = \kappa_{S}h_{0}(t|\boldsymbol{\theta})\exp\{H_{0}(t|\boldsymbol{\theta})\}\left[1+\frac{\kappa_{S}}{\kappa_{L}}R_{0}(t|\boldsymbol{\theta})\right]^{-(1+\kappa_{L})}, \]
respectively, where \(\kappa_{S} = \exp\{\mathbf{x}\boldsymbol{\beta}\}\) and \(\kappa_{L} = \exp\{\mathbf{x}\boldsymbol{\phi}\}\).
The YO model includes the PH and PO models as particular cases when \(\boldsymbol{\phi} = \boldsymbol{\beta}\) and \(\boldsymbol{\phi} = \mathbf{0}\), respectively.
These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.