Time-to-event, with only fraction of events observed - All post-warmup steps divergent

Modeling
1

Hello,

I have a dataset where each entry contains:
(Time to event (NA if event not observed), Time observed), where Time to event <= Time observed.

I expect only a fraction of the total population to ever observe the event, so I’m trying to build a model with the following likelihood:

For observed events (ie, time to event is not N/A):
L = p(event observed during lifetime) * p(Time to event), as:
L = Bernoulli(1 | theta) * Weibull(Time to event | alpha, sigma)

For non-observed events:
L = p(event not observed during lifetime) + p(event observed during lifetime) * p(Time observed < Time to event) as
L = Bernoulli (0 | theta) + Bernoulli(1 | Theta) * Weibull(Imputed time to event | alpha, sigma) T [Time observed, Inf]

I have tried encoding this as the attached stan-model, but when I run the model, all post-warmup steps are divergent, which in my experience points to a serious flaw in the model.

Is there something obviously wrong with the model? I’m not sure trying to impute time to event from the truncated distribution in the missing observation case is the ideal way of expressing p(time observed < time to event). Any ideas how to do this directly using the CDF?

data {
  int<lower=0> N_observed;
  vector<lower=0>[N_observed] y_observed; // time to event
  int<lower=0> N_not_observed;
  vector[N_not_observed] y_not_observed; // time observed without event
}

parameters {
  real<lower=0> alpha;
  real<lower=0> sigma;
  real<lower=0, upper=1> theta;
  vector<lower=0>[N_not_observed] not_observed_t;
}

model {
  theta ~ beta(1, 1);
  alpha ~ exponential(1);
  sigma ~ exponential(1);
  for(i in 1:N_observed) {
    target += bernoulli_lpmf(1 | theta);
    target += weibull_lpdf(y_observed[i] | alpha, sigma);
  }
  
  for(i in 1:N_not_observed) {
    not_observed_t[i] ~ weibull(alpha, sigma) T[y_not_observed[i],];
    target += log_sum_exp(
      bernoulli_lpmf(0 | theta),
      bernoulli_lpmf(1 | theta)
      + weibull_lpdf(not_observed_t[i] | alpha, sigma) 
      - weibull_lccdf(y_not_observed[i] | alpha, sigma));
  }
}
2

Hello and welcome,

This type of model is outside my area. But a few things to help folks out

  1. Can you supply a snippet of data or some simulated data for the model? This will help folks run it and troubleshoot.
  2. Can you post what version of Stan you using along with OS, python or R or something else?

There is a time to event model in rstanarm here:

and a linger read here with a version in brms (uses Stan and you can dump out the Stan code):
https://www.rensvandeschoot.com/tutorials/discrete-time-survival/

3

Hi Ara,

Thanks, and thank you for the reply! I have updated the model to the following (integrating out the censored non-observed events, as suggested in 4.3 Censored data | Stan User’s Guide). It now looks as follows:

  int<lower=0> N_observed;
  int<lower=0> N_not_observed;
  vector<lower=0> [N_observed] t_observed;
  vector<lower=0>[N_not_observed] t_not_observed;
}

parameters {
  real<lower=0> alpha;
  real<lower=0> sigma;
  real<lower=0, upper=1> theta;
}

model {
  alpha ~ exponential(3);
  sigma ~ exponential(3);
  theta ~ beta(1,1);

  for(i in 1:N_observed) {
    target += bernoulli_lpmf(1 | theta);
    target += weibull_lpdf(t_observed[i] | alpha, sigma);
  }
  
  for(i in 1:N_not_observed) {
    target += log_sum_exp(
      bernoulli_lpmf(0 | theta),
      bernoulli_lpmf(1 | theta) * weibull_lccdf(t_not_observed[i] | alpha, sigma)
    );
  }
}

Simulated data:

n_samples <- 1000
sim_alpha <- 0.3
sim_sigma <- 0.7
sim_theta <- 0.2
y = tibble(
  obs_time = runif(n_samples, 0, 900),
  event_time = rweibull(n_samples, sim_alpha, sim_sigma),
  converted = rbernoulli(n_samples, sim_theta)) %>%
  mutate(event_time = if_else(converted, event_time, as.double(NA)))

t_observed <- y %>% filter(event_time <= obs_time) %>% pull(event_time)
t_not_observed <- y %>% filter(is.na(event_time) | event_time > obs_time) %>% pull(obs_time)

fit <- stan("mix2.stan",
            control = list(adapt_delta = 0.995),
            data = list(N_observed = length(t_observed),
                        t_observed = t_observed,
                        N_not_observed = length(t_not_observed),
                        t_not_observed = t_not_observed))

With adapt_delta = 0.9, it still gets the occasional divergence, but also ESS is extremely small (n_eff about 2-3). Looking at the chains for all three parameters in shinystan, all three parameters are non-stationary.

Admittedly, I’m not very knowledgeable about survival/cure models (will check out your links!), but most seem to assume the subject always dies (either by uncured disease, or later, by natural causes). My use case is for user conversion for e-commerce (some clients never end up converting), that’s why I’m trying to “mix” with the bernoulli distribution.

My system:

> sessionInfo()
R version 4.0.2 (2020-06-22)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Catalina 10.15.7

Matrix products: default
BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] shinystan_2.5.0      shiny_1.5.0          rstan_2.21.2         StanHeaders_2.21.0-6 lubridate_1.7.9      forcats_0.5.0        stringr_1.4.0        dplyr_1.0.2         
 [9] purrr_0.3.4          readr_1.3.1          tidyr_1.1.2          tibble_3.0.3         ggplot2_3.3.2        tidyverse_1.3.0     

loaded via a namespace (and not attached):
 [1] matrixStats_0.57.0 fs_1.5.0           xts_0.12-0         threejs_0.3.3      httr_1.4.2         tools_4.0.2        backports_1.1.9    DT_0.15            R6_2.5.0          
[10] DBI_1.1.0          colorspace_2.0-0   withr_2.2.0        tidyselect_1.1.0   gridExtra_2.3      prettyunits_1.1.1  processx_3.4.5     curl_4.3           compiler_4.0.2    
[19] cli_2.0.2          rvest_0.3.6        shinyjs_2.0.0      xml2_1.3.2         labeling_0.4.2     colourpicker_1.1.0 scales_1.1.1       dygraphs_1.1.1.6   ggridges_0.5.2    
[28] callr_3.5.1        digest_0.6.27      rmarkdown_2.6      base64enc_0.1-3    pkgconfig_2.0.3    htmltools_0.5.1.1  dbplyr_1.4.4       fastmap_1.0.1      htmlwidgets_1.5.1 
[37] rlang_0.4.10       readxl_1.3.1       rstudioapi_0.13    farver_2.0.3       generics_0.0.2     zoo_1.8-8          jsonlite_1.7.2     gtools_3.8.2       crosstalk_1.1.0.1 
[46] inline_0.3.16      magrittr_2.0.1     loo_2.3.1          bayesplot_1.7.2    Rcpp_1.0.5         munsell_0.5.0      fansi_0.4.1        lifecycle_0.2.0    stringi_1.5.3     
[55] yaml_2.2.1         pkgbuild_1.1.0     plyr_1.8.6         grid_4.0.2         blob_1.2.1         parallel_4.0.2     promises_1.1.1     crayon_1.3.4       miniUI_0.1.1.1    
[64] lattice_0.20-41    haven_2.3.1        hms_0.5.3          knitr_1.30         ps_1.5.0           pillar_1.4.6       igraph_1.2.5       markdown_1.1       reshape2_1.4.4    
[73] codetools_0.2-16   stats4_4.0.2       reprex_0.3.0       glue_1.4.2         evaluate_0.14      V8_3.2.0           RcppParallel_5.0.2 modelr_0.1.8       vctrs_0.3.4       
[82] httpuv_1.5.4       cellranger_1.1.0   gtable_0.3.0       assertthat_0.2.1   xfun_0.20          mime_0.9           xtable_1.8-4       broom_0.7.0        later_1.1.0.1     
[91] rsconnect_0.8.16   shinythemes_1.1.2  ellipsis_0.3.1
2 Likes
4

Thanks for taking the time to write that all up. That will helps folks troubleshoot.