Monotone splines

Ramsay, James O. “Monotone regression splines in action.” Statistical science 3.4 (1988): 425-441.

\[h_0(t;\gamma;k) = \sum_{j=1}^{m} M_j(t;k)\gamma_j\]

\[I_j(t;k) = \int_{0}^{t} M_j(u;k) \rm{d}u\]

\[ H_0(t;\gamma;k) = \sum_{j=1}^{m} I_j(t;k)\gamma_j\]

m-Splines

Random linear combinations

i-Splines

Random linear combinations

Stan model

data {
    int<lower=0> N_uncensored;                                   
    int<lower=0> N_censored;                                        
    int<lower=1> m;                                                 
    int<lower=1> NC;                                                
    matrix[N_censored,NC] X_censored;                               
    matrix[N_uncensored,NC] X_uncensored;                                                
    matrix[N_uncensored,m] m_spline_basis_evals_uncensored;                  
    matrix[N_uncensored,m] i_spline_basis_evals_uncensored;   
    matrix[N_censored,m] i_spline_basis_evals_censored;
}
parameters {
    simplex[m] gammas;       
    vector[NC] betas;                                            
    real intercept;   
}
model {
    betas ~ normal(0,1);
    intercept   ~ normal(0,5);
    target += -(i_spline_basis_evals_censored*gammas) .* exp(X_censored*betas + intercept);
    target += -(i_spline_basis_evals_uncensored*gammas) .* exp(X_uncensored*betas + intercept);
      target +=  log(m_spline_basis_evals_uncensored*gammas) + X_uncensored*betas + intercept;
}

Inference

Survival curves

Brier score

Gerds, Thomas A., and Martin Schumacher. “Consistent estimation of the expected Brier score in general survival models with right‐censored event times.” Biometrical Journal 48.6 (2006): 1029-1040. APA

pec: Prediction Error Curves for Risk Prediction Models in Survival Analysis

\[ \rm{MSE}(t,\hat{\pi}_n, S) = \int_{\mathbb{R}^d} \int_0^\infty \left( \mathbf{1}_{\{s>t\}} - \hat{\pi}_n(t\vert \mathbf{x})\right)^2 F(ds\vert \mathbf{x}) H(d\mathbf{x}) \]

\[ \rm{WRSS}(t, \hat{\pi}_n, \hat{G}_n) = \frac{1}{n} \sum_{i=1}^n \left(\mathbf{1}_{\{T_i > t\}} -\hat{\pi}_n(t \vert \mathbf{x}_i)\right)^2 \omega(t, \hat{G}_n, \mathbf{x}_i, T_i, \Delta_i) \]

\[ \omega(t, \hat{G}_n, \mathbf{x}_i, T_i, \Delta_i) = \frac{\mathbf{1}_{\{T_i \leq t\}} \Delta_i}{\hat{G}_n(T_i \vert \mathbf{x}_i)} + \frac{\mathbf{1}_{\{T_i > t\}}}{\hat{G}_n(t \vert \mathbf{x}_i)} \]

All estimators are robust against misspecification of the survival model!

(Conditioned) Posterior predictive checks

How to generate survival time samples (m-Splines)

Reminder: Inverse transform sampling

\[U\sim \rm{Unif}[0,1] \Rightarrow F_{X}^{-1}(U)\sim X\]

Algorithm

  1. Let \(U_{ij}\sim \rm{Exp}\left(e^{\beta' \mathbf{x}_i}\right)\) for \(i=1,2,\dots,n\) and \(j=1,2,\dots\) be an i.i.d. sequence
  2. Let \(T^\star\) be the largest survival / event time in the data.
  3. For each \(i\in 1\dots n\)
    • Let \(U_{ij(i)}\) be the first sample in the i.i.d. sequence \(\{U_{ij}\}_{j=1,2,\dots}\) for which \(U_{ij} \leq H_0(T^\star;\gamma;k)\).
    • Find \(t\) such that \(U_{ij(i)} = H_0(t;\gamma;k)\) (e.g. using R’s uniroot; the loss function is convex, since \(H_0\) is increasing)
    • Return \(\widehat{T}_i = t\)

Comparing distributions

Mean

Median

Standard deviation

Min

Max

Non-parametric Kaplan-Meier estimator

\[ \widehat{S}(t) = \prod_{t_i \leq t} \left( 1- \frac{d_i}{n_i} \right) \]

  • \(t_i\): Time at which at least one event happened
  • \(d_i\): Number of events (e.g. deaths) that happened at time \(t_i\)
  • \(n_i\): Number of individuals known to survive at time \(t_i\) (i.e. have not yet had an event or been censored)

Below we use \(\widehat{S}(70)\) as a statistics for posterior predictive checking. Note that,

  • We only vary the uncensored cases (i.e. we replace them with posterior predictive samples)
  • We keep the censored cases fixed over all replications (our generative model does not contain sufficient structure to explain the censoring process)

Outlook

Survival models in rstanarm

  • More model comparison using Leave-One-Out-Cross-Validation (full loo support)
  • Please test!

Example

library(devtools)
dev_mode(on=T)
library(rstanarm)
fit_stansurv <- stan_surv(Surv(time, event)~metastized, data=df, prior=normal(0,2), seed=42)
dev_mode(on=F)
## # A tibble: 7 x 3
##   term            estimate std.error
##   <chr>              <dbl>     <dbl>
## 1 metastized        0.0906    0.368 
## 2 m-splines-coef1   0.0331    0.0330
## 3 m-splines-coef2   0.129     0.0908
## 4 m-splines-coef3   0.491     0.216 
## 5 m-splines-coef4   0.171     0.165 
## 6 m-splines-coef5   0.202     0.179 
## 7 m-splines-coef6   0.206     0.201
## 
## Model Info:
## 
##  function:        stan_surv
##  baseline hazard: M-splines on hazard scale
##  formula:         Surv(time, event) ~ metastized
##  algorithm:       sampling
##  priors:          see help('prior_summary')
##  sample:          4000 (posterior sample size)
##  observations:    44
##  events:          26 (59.1%)
##  right censored:  18 (40.9%)
##  delayed entry:   no
## 
## Estimates:
##                   mean      sd        2.5%      25%       50%    
## metastized         0.0985    0.3681   -0.6252   -0.1527    0.0906
## m-splines-coef1    0.0453    0.0440    0.0010    0.0140    0.0331
## m-splines-coef2    0.1460    0.0996    0.0094    0.0732    0.1285
## m-splines-coef3    0.5191    0.2283    0.1557    0.3558    0.4909
## m-splines-coef4    0.2233    0.1952    0.0075    0.0742    0.1709
## m-splines-coef5    0.2497    0.2052    0.0111    0.0941    0.2016
## m-splines-coef6    0.2872    0.2701    0.0081    0.0935    0.2059
## log-posterior   -195.6212    2.1610 -200.8249 -196.8072 -195.2351
##                   75%       97.5%  
## metastized         0.3484    0.8257
## m-splines-coef1    0.0623    0.1644
## m-splines-coef2    0.1992    0.3831
## m-splines-coef3    0.6516    1.0562
## m-splines-coef4    0.3118    0.7289
## m-splines-coef5    0.3520    0.7524
## m-splines-coef6    0.3988    1.0087
## log-posterior   -194.0379 -192.5600
## 
## Diagnostics:
##                 mcse   Rhat   n_eff
## metastized      0.0078 0.9998 2210 
## m-splines-coef1 0.0008 0.9999 3356 
## m-splines-coef2 0.0021 1.0004 2267 
## m-splines-coef3 0.0045 1.0000 2583 
## m-splines-coef4 0.0036 1.0011 2878 
## m-splines-coef5 0.0038 1.0000 2869 
## m-splines-coef6 0.0044 1.0005 3712 
## log-posterior   0.0634 1.0052 1161 
## 
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).