Statistical Specifications

IMPORTANT

Please note that this document is currently a work-in-progress and does not contain complete information for this package yet.

Survival Model Specification

This package can only be used to fit proportional hazards models of the form:

$$\log(h_i(t \mid \theta, \psi_i)) = \log(h_0(t \mid \theta)) + X_i \beta + \sum_j G_j(t \mid \psi_i)$$

Where:

• $h_0(.)$ is a parametric baseline hazard function
• $t$ is the event time
• $\theta$ is a vector of parameters that parameterise the baseline hazard
• $\psi_i$ is an arbitrary vector of parameters for subject $i$ specified by the longitudinal model
• $G_j(.)$ is a link function that maps $\psi_i$ to a contribution to the log-hazard function where $j$ simply indexes the given link function
• $X_i$ is the subjects covariate design matrix
• $\beta$ is the corresponding coefficients to scale the design matrix covariates contribution to the log-hazard function

The following sections outline the available distributions to users which can be selected for the baseline hazard $h_0(.)$. Please note that some of these distributions do not have the proportional-hazards property meaning that the resulting survival model corresponding to the hazard $h_i()$ will not be of the same parametric family as the baseline distribution with the hazard $h_0(.)$.

Exponential Distribution

$$h(t \mid \lambda) = \lambda$$

Where: - $\lambda > 0$ is the rate parameter

Weibull Distribution (Proportional Hazard Parameterisation)

$$h(t \mid \lambda, \gamma) = \lambda \gamma t^{\gamma - 1 };$$

Where: - $\lambda > 0$ is the rate parameter - $\gamma > 0$ is the shape parameter Note that with $\gamma = 1$ we obtain the exponential distribution as a special case.

Log-Logistic Distribution

$$h(t \mid a, b) = \frac {(b/a)(t/a)^{(b-1)}} {1 + (t/a)^b}$$

Where:

• $a > 0$ is the scale parameter
• $b > 0$ is the shape parameter

Gamma Distribution

$$h(t \mid k, \theta) = \frac{ t^{k - 1} e^{-t/\theta} } { \theta ^ k ( \Gamma(k) - \gamma(k, t/\theta) ) }$$

Where:

• $k > 0 $ is the shape parameter
• $\theta > 0$ is the scale parameter
• $\Gamma(.)$ is the complete gamma function
• $\gamma(., .)$ is the lower incomplete gamma function

Covariate Selection with Regularized Horseshoe Prior

In the survival model the covariate coefficients $\beta$ typically just have a normal prior distribution, using prior_normal() or prior_normal_vector(). However, if the user wishes to perform covariate selection then they can instead use the regularized horseshoe prior via prior_horseshoe(). Here we explain how the prior is defined and how the hyperparameters can be set. Note that we use the same parametrization and notation as in the brms package (Bürkner 2017).

The horseshoe prior is a global-local shrinkage prior (Carvalho, Polson, and Scott 2010). The global parameter shrinks all coefficients towards zero, while coefficient-specific local parameters allow some coefficients to escape the shrinkage. The regularized horseshoe prior (Piironen and Vehtari 2017) additionally introduces a slab component that regularizes coefficients that are far from zero. This is useful for covariate selection because all candidate covariates can be included in one model fit, with irrelevant coefficients strongly shrunk towards zero.

For $p$ covariate coefficients $\beta_{j}, j = 1, \dotsc, p$, the regularized horseshoe prior applied with prior_horseshoe(), and with positive hyperparameters $\nu_{\lambda}, \nu_{\tau}, \nu_c, s_{\tau}, s_c$, can be written as

$$\begin{align*} z_j &\sim \operatorname{Normal}(0, 1), \\ \lambda_j &\sim t^+_{\nu_{\lambda}}(0, 1), \\ \tau &\sim t^+_{\nu_{\tau}}(0, s_{\tau}), \\ \xi &\sim \text{Inv-Gamma}\!\left(\frac{\nu_c}{2}, \frac{\nu_c}{2}\right), \\ c^2 &= s_c^2 \xi, \\ \tilde{\lambda}_j &= \left( \frac{c^2\lambda_{j}^2} {c^2 + \tau^2\lambda_{j}^2} \right)^{1/2}, \\ \beta_j &= z_j\tau\tilde{\lambda}_j, \qquad j = 1, \dotsc, p . \end{align*}$$

Here $t^+_{\nu}(0, s)$ denotes a half-Student-$t$ distribution with degrees of freedom $\nu$ and scale $s$. The $p$ local shrinkage parameters $\lambda_j$ control the amount of shrinkage for each coefficient, while the global shrinkage parameter $\tau$ controls the overall amount of shrinkage across all coefficients. The slab scale $s_c$ and slab degrees of freedom $\nu_c$ determine how strongly large coefficients are regularized.

Equivalently, conditional on $\lambda_j$, $\tau$, and $c$, this implies

$$\beta_j \mid \lambda_j, \tau, c \sim \operatorname{Normal}\!\left( 0,\, \tau^2 \frac{c^2\lambda_j^2}{c^2 + \tau^2\lambda_j^2} \right).$$

The arguments of prior_horseshoe() correspond to these hyperparameters:

• df: local degrees of freedom, $\nu_{\lambda}$
• df_global: global degrees of freedom, $\nu_{\tau}$
• df_slab: slab degrees of freedom, $\nu_c$
• scale_global: global scale, $s_{\tau}$
• scale_slab: slab scale, $s_c$

As a special case, setting $\nu_{\lambda} = 1$ and $\nu_{\tau} = 1$ gives half-Cauchy priors for the local and global shrinkage parameters. The slab part has the interpretation

$$c^2 \sim \text{Inv-Gamma}\!\left( \frac{\nu_c}{2}, \frac{\nu_c s_c^2}{2} \right),$$

which corresponds to a Student-$t$ slab with $\nu_c$ degrees of freedom and scale $s_c$ for coefficients that escape the main shrinkage. Smaller values of $s_{\tau}$ imply stronger prior shrinkage towards zero, while larger values allow more coefficients to remain non-zero. The slab scale should be chosen large enough to permit clinically or scientifically plausible covariate effects, but not so large that unrealistically large effects are left essentially unregularized.

A useful summary of the amount of shrinkage applied to coefficient $j$ is the shrinkage factor

$$\kappa_j = \frac{1}{1 + \tau^2\lambda_j^2 / c^2}.$$

Values of $\kappa_j$ close to 1 indicate strong shrinkage towards zero, while values close to 0 indicate little shrinkage. Regularized horseshoe priors have been applied in PK/PD covariate selection by Pourzanjani and Davis (2026), where the main practical advantage is the ability to assess all candidate covariates in a single Bayesian model fit.

Longitudinal Model Specification

Random-Slope Model

$$y_{ij} = \mu_{l(i)} + s_i t_{ij} + \epsilon_{ij}$$

where:

• $y_{ij}$ is the tumour size for subject $i$ at timepoint $j$
• $\mu_{l(i)}$ is the intercept for subject $i$
• $s_i$ is the random slope for subject $i$ with $s_i \sim N(\mu_{sk(i)}, \sigma_{sk(i)})$
• $\mu_{sk(i)}$ is the mean for the random slope within treatment arm $k(i)$
• $\sigma_{sk(i)}$ is the variance term for the slope in treatment arm $k(i)$
• $\epsilon_{ij}$ is the error term with $\epsilon_{ij} \sim N(0, \sigma)$
• $k(i)$ is the treatment arm index for subject $i$
• $l(i)$ is the study index for subject $i$

Derivative of the SLD Trajectory Link / Growth Parameter Link

$$G(t_{ij} \mid \mu_0, s_i) = s_i$$

Accessible via linkDSLD() & linkGrowth()

Identity Link

$$\begin{align*} G(t_{ij} \mid \mu_0, s_i) = \mu_0 + s_i t_{ij} \end{align*}$$

Accessible via linkIdentity()

Stein-Fojo Model

$$\begin{align*} y_{ij} &\sim N(SLD_{ij},\ SLD_{ij}^2 \sigma^2) \\ \\ SLD_{ij} &= \begin{cases} b_i[e^{-s_it_{ij}} + e^{g_i t_{ij}} - 1] & \text{if } t_{ij}\geq 0 \\ b_i e^{g_i t_{ij}} & \text{if } t_{ij}\lt 0 \end{cases}\\\\ b_i &\sim \text{LogNormal}(\mu_{bl(i)}, \omega_{b l(i)}) \\ s_i &\sim \text{LogNormal}(\mu_{sk(i)}, \omega_{s k(i)}) \\ g_i &\sim \text{LogNormal}(\mu_{gk(i)}, \omega_{g k(i)}) \\ \end{align*}$$

Where:

• $i$ is the subject index
• $j$ is the visit index
• $y_{ij}$ is the observed tumour measurements
• $SLD_{ij}$ is the expected sum of longest diameter for subject $i$ at time point $j$
• $t_{ij}$ is the time since first treatment for subject $i$ at visit $j$
• $b_i$ is the subject baseline SLD measurement
• $s_i$ is the subject kinetics shrinkage parameter
• $g_i$ is the subject kinetics tumour growth parameter
• $\phi_i$ is the subject proportion of cells affected by the treatment
• $k(i)$ is the treatment arm index for subject $i$
• $l(i)$ is the study index for subject $i$
• $\mu_{\theta k(i)}$ is the population mean for parameter $\theta$ in group $k(i)$
• $\omega_{\theta k(i)}$ is the population variance for parameter $\theta$.

If using the non-centred parameterisation then the following alternative formulation is used:

$$\begin{align*} b_i &= exp(\mu_{bl(i)} + \omega_{b l(i)} * \eta_{b i}) \\ s_i &= exp(\mu_{sk(i)} + \omega_{s k(i)} * \eta_{s i}) \\ g_i &= exp(\mu_{gk(i)} + \omega_{g k(i)} * \eta_{g i}) \\ \\ \eta_{b i} &\sim N(0, 1)\\ \eta_{s i} &\sim N(0, 1) \\ \eta_{g i} &\sim N(0, 1) \\ \end{align*}$$

Where:

• $\eta_{\theta i}$ is a random effects offset on parameter $\theta$ for subject $i$

If using the unscaled variance parameterisation then the following alternative formulation is used: $$y_{ij} \sim N(SLD_{ij},\sigma^2)$$

Derivative of the SLD Trajectory Link

$$G(t_{ij} \mid b_i, s_i, g_i) = \begin{cases} b_i(g_i e^{g_i t_{ij}} -s_ie^{-s_it_{ij}} ) & \text{if } t_{ij}\geq 0 \\ b_i g_i e^{g_i t_{ij}} & \text{if } t_{ij}\lt 0 \end{cases}$$

Accessible via linkDSLD()

Time to Growth Link

$$G(t_{ij} \mid b_i, s_i, g_i) = \max \left( \frac{ \text{log}(s_i) - \text{log}(g_i) }{ s_i + g_i }, 0 \right)$$

Accessible via linkTTG()

Identity Link

$$\begin{align*} G(t_{ij} \mid b_i, s_i, g_i) &= SLD_{ij} \end{align*}$$

Accessible via linkIdentity()

Growth Parameter Link

$$\begin{align*} G(t_{ij} \mid b_i, s_i, g_i) &= log(g_i) \end{align*}$$

Accessible via linkGrowth()

Shrinkage Parameter Link

$$\begin{align*} G(t_{ij} \mid b_i, s_i, g_i) &= log(s_i) \end{align*}$$

Accessible via linkShrinkage()

Population Quantities

Note that when generating population quantities for the $b$, $s$ and $g$ parameters the median of the distribution is used. This is because this has the same interpretation as using a non-centred parameterisation and setting the “random effects” term to be 0. For the $\phi$ parameter the mean of the distribution is used.

Generalized Stein-Fojo (GSF) Model

$$\begin{align*} y_{ij} &\sim N(SLD_{ij},\ SLD_{ij}^2 \sigma^2) \\ \\ SLD_{ij} &= \begin{cases} b_i[\phi_i e^{-s_it_{ij}} + (1-\phi_i)e^{g_i t_{ij}}] & \text{if } t_{ij}\geq 0 \\ b_i e^{g_i t_{ij}} & \text{if } t_{ij}\lt 0 \end{cases}\\ \\ b_i &\sim \text{LogNormal}(\mu_{bl(i)}, \omega_{b l(i)}) \\ s_i &\sim \text{LogNormal}(\mu_{sk(i)}, \omega_{s k(i)}) \\ g_i &\sim \text{LogNormal}(\mu_{gk(i)}, \omega_{g k(i)}) \\ \phi_i &\sim \text{LogitNormal}(\mu_{\phi k(i)}, \omega_{\phi k(i)}) \end{align*}$$

Where:

• $i$ is the subject index
• $j$ is the visit index
• $y_{ij}$ is the observed tumour measurements
• $SLD_{ij}$ is the expected sum of longest diameter for subject $i$ at time point $j$
• $t_{ij}$ is the time since first treatment for subject $i$ at visit $j$
• $b_i$ is the subject baseline SLD measurement
• $s_i$ is the subject kinetics shrinkage parameter
• $g_i$ is the subject kinetics tumour growth parameter
• $\phi_i$ is the subject proportion of cells affected by the treatment
• $k(i)$ is the treatment arm index for subject $i$
• $l(i)$ is the study index for subject $i$
• $\mu_{\theta k(i)}$ is the population mean for parameter $\theta$ in group $k(i)$
• $\omega_{\theta k(i)}$ is the population variance for parameter $\theta$ in group $k(i)$.

If using the non-centred parameterisation then the following alternative formulation is used:

$$\begin{align*} b_i &= exp(\mu_{bl(i)} + \omega_{b l(i)} * \eta_{b i}) \\ s_i &= exp(\mu_{sk(i)} + \omega_{s k(i)} * \eta_{s i}) \\ g_i &= exp(\mu_{gk(i)} + \omega_{g k(i)} * \eta_{g i}) \\ \phi_i &= \text{logistic}(\mu_{\phi k(i)} + \omega_{\phi k(i)} * \eta_{\phi i}) \\ \\ \eta_{b i} &\sim N(0, 1)\\ \eta_{s i} &\sim N(0, 1) \\ \eta_{g i} &\sim N(0, 1) \\ \eta_{\phi i} &\sim N(0, 1) \\ \end{align*}$$

Where:

• $\eta_{\theta i}$ is a random effects offset on parameter $\theta$ for subject $i$

If using the unscaled variance parameterisation then the following alternative formulation is used: $$y_{ij} \sim N(SLD_{ij}, \sigma^2)$$

Derivative of the SLD Trajectory Link

$$G(t_{ij} \mid b_i, s_i, g_i, \phi_i) = \begin{cases} b_i(-s_i\phi_ie^{-s_it_{ij}} + (1 - \phi_i)g_i e^{g_i t_{ij}}) & \text{if } t_{ij}\geq 0 \\ b_i g_i e^{g_i t_{ij}} & \text{if } t_{ij}\lt 0 \end{cases}$$

Accessible via linkDSLD()

Time to Growth Link

$$G(t_{ij} \mid b_i, s_i, g_i, \phi_i) = \max \left( \frac{ \log(s_i) + \log(\phi_i) - \log(g_i) - \log(1-\phi_i) }{ s_i + g_i }, 0 \right)$$

Accessible via linkTTG()

Identity Link

$$\begin{align*} G(t_{ij} \mid b_i, s_i, g_i, \phi_i) &= SLD_{ij} \end{align*}$$

Accessible via linkIdentity()

Growth Parameter Link

$$\begin{align*} G(t_{ij} \mid b_i, s_i, g_i, \phi_i) &= log(g_i) \end{align*}$$

Accessible via linkGrowth()

Shrinkage Parameter Link

$$\begin{align*} G(t_{ij} \mid b_i, s_i, g_i, \phi_i) &= log(s_i) \end{align*}$$

Accessible via linkShrinkage()

Population Quantities

Note that when generating population quantities for the $b$, $s$, $g$ and $\phi$ parameters the median of the distribution is used. This is because this has the same interpretation as using a non-centred parameterisation and setting the “random effects” term to be 0.

Claret-Bruno Model

$$\begin{align*} y_{ij} &\sim N(SLD_{ij},\ SLD_{ij}^2 \sigma^2) \\ \\ SLD_{ij} &= \begin{cases} b_i e^{g_i t_{ij}} & \text{if } t_{ij} < 0, \\ b_i \cdot \exp\left(g_i t_{ij} - \frac{p_i}{c_i} \left(1 - e^{-c_i t_{ij} }\right)\right) & \text{if } t_{ij} \geq 0. \end{cases}\\ \\ b_i &\sim \text{LogNormal}(\mu_{bl(i)}, \omega_{b l(i)}) \\ g_i &\sim \text{LogNormal}(\mu_{gk(i)}, \omega_{g k(i)}) \\ c_i &\sim \text{LogNormal}(\mu_{ck(i)}, \omega_{c k(i)}) \\ p_i &\sim \text{LogNormal}(\mu_{pk(i)}, \omega_{p k(i)}) \\ \end{align*}$$

Where:

• $i$ is the subject index
• $j$ is the visit index
• $y_{ij}$ is the observed tumour measurements
• $SLD_{ij}$ is the expected sum of longest diameter for subject $i$ at time point $j$
• $t_{ij}$ is the time since first treatment for subject $i$ at visit $j$
• $b_i$ is the subject baseline SLD value.
• $g_i$ is the subject tumour growth rate.
• $c_i$ is the subject treatment resistance rate.
• $p_i$ is the subject treatment growth inhibition response.
• $k(i)$ is the treatment arm index for subject $i$
• $l(i)$ is the study index for subject $i$
• $\mu_{\theta k(i)}$ is the population mean for parameter $\theta$ in group $k(i)$
• $\omega_{\theta k(i)}$ is the population variance for parameter $\theta$ in group $k(i)$

If using the non-centred parameterisation then the following alternative formulation is used:

$$\begin{align*} b_i &= exp(\mu_{b l(i)} + \omega_{b l(i)} * \eta_{b i}) \\ g_i &= exp(\mu_{g k(i)} + \omega_{g k(i)} * \eta_{g i}) \\ c_i &= exp(\mu_{c k(i)} + \omega_{c k(i)} * \eta_{c i}) \\ p_i &= exp(\mu_{p k(i)} + \omega_{p k(i)} * \eta_{p i}) \\ \\ \eta_{b i} &\sim N(0, 1)\\ \eta_{g i} &\sim N(0, 1) \\ \eta_{c i} &\sim N(0, 1) \\ \eta_{p i} &\sim N(0, 1) \\ \end{align*}$$

If using the unscaled variance parameterisation then the following alternative formulation is used: $$y_{ij} \sim N(SLD_{ij}, \sigma^2)$$

Derivative of the SLD Trajectory Link

$$G(t_{ij} \mid b_i, g_i, c_i, p_i) = \begin{cases} b_i g_i e^{g_i t_{ij}} & \text{if } t_{ij}\lt 0 \\ \left( g_i - p_i e^{-c_i t_{ij}} \right) SLD_{ij} & \text{if } t_{ij}\geq 0 \end{cases}$$

Accessible via linkDSLD()

Time to Growth Link

$$G(t_{ij} \mid b_i, g_i, c_i, p_i) = \max \left( \frac{ln(\frac{p_i}{g_i})}{c_i}, 0 \right)$$

Accessible via linkTTG()

Identity Link

$$\begin{align*} G(t_{ij} \mid b_i, g_i, c_i, p_i) &= SLD_{ij} \end{align*}$$

Accessible via linkIdentity()

Growth Parameter Link

$$\begin{align*} G(t_{ij} \mid b_i, g_i, c_i, p_i) &= log(g_i) \end{align*}$$

Accessible via linkGrowth()

Population Quantities

When generating population quantities for the $b$, $g$, $c$ and $p$ parameters the median of the distribution is used. This is because this has the same interpretation as using a non-centred parameterisation and setting the “random effects” term to be 0.

Post-Processing

Brier Score

The Brier Score is used to measure a models predictive performance. In the case of Survival Models it is the weighted squared difference between whether a subject died within a given time interval and the models estimated probability of them dying within said interval. Within jmpost the following formula (as described in Blanche et al. (2015)) has been implemented:

$$\hat{BS}(t) = \frac{1}{n}\sum_{i=1}^{n} \hat{W}_i(t) \Big(D_i(t) - \pi_i(t) \Big)^2$$

$$\hat{W}_i(t) = \frac{\mathbb{1}_{(T_i \gt t)}}{\hat{G}(t)} + \frac{\mathbb{1}_{(T_i \le t)} \cdot \Delta_i}{\hat{G}(T_i)}$$

Where:

• $T_i$ is the observed event/censor time for subject $i$
• $\Delta_i$ is the event indicator which is 1 if $T_i$ is an event and 0 if $T_i$ is censored
• $\mathbb{1}_{(.)}$ is the indicator which is 1 if ${(.)}$ is true else 0
• $D_i(t)$ is the event indicator function for subject $i$ e.g. $D_i(t) = \mathbb{1}_{(T_i \lt t,\ \delta_i = 1)}$
• $\pi_i(t)$ is a model predicted probability of subject $i$ dying before time $t$
• $\hat{G}(u)$ is the Kaplan-Meier estimator of survival function of the censoring time at $u$

Note that by default $\hat{G}(T_i)$ is estimated by $\hat{G}(T_i-)$ and that in the case of ties event times are always considered to have occurred before censored times (this is in contrast to survival::survfit which regards censored times as coming before event times when estimating the censoring distribution). Both of these default options can be changed if required.

References

Blanche P, Proust-Lima C, Loubère L, Berr C, Dartigues J-F, Jacqmin-Gadda H (2015). “Quantifying and Comparing Dynamic Predictive Accuracy of Joint Models for Longitudinal Marker and Time-to-Event in Presence of Censoring and Competing Risks.” Biometrics, 71(1), 102–113. https://doi.org/https://doi.org/10.1111/biom.12232.

Bürkner P-C (2017). “brms: An R Package for Bayesian Multilevel Models Using Stan.” Journal of Statistical Software, 80(1), 1–28. https://doi.org/10.18637/jss.v080.i01.

Carvalho CM, Polson NG, Scott JG (2010). “The Horseshoe Estimator for Sparse Signals.” Biometrika, 97(2), 465–480. https://doi.org/10.1093/biomet/asq017.

Piironen J, Vehtari A (2017). “Sparsity Information and Regularization in the Horseshoe and Other Shrinkage Priors.” Electronic Journal of Statistics, 11(2), 5018–5051. https://doi.org/10.1214/17-EJS1337SI.

Pourzanjani A, Davis C (2026). “The Regularized Horseshoe for Covariate Selection Improves Convenience and Predictive Performance in Population PK/PD Models.” CPT: Pharmacometrics & Systems Pharmacology, 15(2), e70198. https://doi.org/10.1002/psp4.70198.