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