[Question]: Clarification on a warning message during mcmc_sample() step with exp_surv_dist() function

Question

[Question]: Clarification on a warning message during mcmc_sample() step with exp_surv_dist() function

souravmukherjee2906 opened this issue 2 years ago · 7 comments

souravmukherjee2906 commented 2 years ago

What is your question?

Dear psborrow2 developers,

Great job on making the psborrow2 package. It's really a useful one.

I wanted to reach out to you to clarify regarding an issue that is showing during compilation of mcmc_sample() function with using exp_surv_dist() function in the "create_simulation_obj.R".

I am getting the following warning when I use verbose = T in mcmc_sample():

"Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Chain 1 Exception: exponential_lccdf: Inverse scale parameter is 0, but must be positive finite! (in '/tmp/RtmpXWDRCn/model-38a15dbdb0ce.stan', line 39, column 6 to column 53)"

and sometimes the following warning:

"Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Chain 1 Exception: exponential_lccdf: Inverse scale parameter is Inf, but must be positive finite! (in '/tmp/RtmpXWDRCn/model-38a15dbdb0ce.stan', line 39, column 6 to column 53)"

It seems like the "elp" parameter that you defined in the stan model block, which is also the inverse scale parameter "lambda" in an exponential (lambda) survival distribution, is sometimes taking either 0 or Infinity during the MCMC sampling, whereas theoretically lambda should be a positive finite real number in every iteration of the Markov Chain. My question is does this warning message affect MCMC convergence?

Do you think we should put a constraint on "elp", something like "vector<lower=0>[N] elp;" so that the inverse scale parameter elp remains always positive?

Looking forward to hearing from you soon. Thanks again for this wonderful package.

Code of Conduct

I agree to follow this project's Code of Conduct.

Contribution Guidelines

I agree to follow this project's Contribution Guidelines.

Answer 1 · 2023-06-27T08:46:14.000Z

Hi @souravmukherjee2906 Thanks for the report!

Is this warning coming at the very beginning of sampling? This might be caused by poor initial conditions for the linear predictor (lp) and its exponential (elp). If this only occurs once or twice in the warm up phase, I would not be concerned about convergence.

For sure we can look at adding this constraint and see if it helps.
Another approach would be to pass init values to cmdstanr. These can be given in mcmc_sample(...). Have a look at ?cmdstanr::`model-method-sample`

Answer 2 · 2023-06-28T03:44:57.000Z

Hi @gravesti Thank you so much for your prompt reply.

Yes, this warning is sometimes happening in beginning of sampling phase. However, occasionally for some scenarios, some of the post-warmup (sampling) transitions are getting ended with a divergence.

Warning: 213 of 10000 (2.0%) transitions ended with a divergence.

I tried including init values for elp in the mcmc_sample(...), but it's still showing the same warning. I tried something like the following:

simulation_obj <- create_simulation_obj(
  data_matrix_list = my_sim_data_list,
  covariate = my_covariate_list,  
  outcome = exp_surv_dist("time", "cnsr", baseline_prior = normal_prior(0, 10000)),
  borrowing = my_borrowing_list,
  treatment = treatment_details(trt_flag_col = "trt", trt_prior = normal_prior(0, 10000)),
  quiet = FALSE
)

simulation_res <- mcmc_sample(
  x = simulation_obj,
  iter_warmup = 10000, 
  iter_sampling = 10000,
  chains = 1,
  verbose = TRUE,
  init = list(
    elp = rep(2, n) # chain 1
  )
)

I think we can only pass init values to something which has been defined in the parameter block of Stan. In the case of exp_surv_dist(..), I think elp has been defined in the model block and not in the parameter block. Usually, variables defined in the Stan model block are considered local variables in Stan and local variables may not be defined with constraints because there is no well-defined way to have them be both flexible and easy to validate.

Could you please let me know if there is any way to modify exp_surv_dist(...) to have elp in the parameter block, or fix the above divergence issue?

Thanks again for your help. I really appreciate it.

Answer 3 · 2023-06-28T15:24:09.000Z

@souravmukherjee2906 You're right of course, we can't set init on elp.

Today I tried to set the inits of the components of the linear predictor, which seemed to significantly improve divergences.

mcmc_sample(object,
  iter_warmup = 10000,
  iter_sampling = 2000, 
  chains = 1, 
  verbose = TRUE,
  init = list(list(beta = c(0,0,0,0), beta_trt = 0, alpha = c(0,0), tau = 100))
)

Also I have tried to move the the definitions for lp and elp to transformed parameters and add the constraint to elp but it doesn't seem to make a difference to those warnings which occur at the beginning of warmup about inverse scale == 0. I'm not sure how elp = exp(lp) = 0 as this implies lp = -Inf, which should not happen.
However I don't think these messages are necessarily related to the divergence.

eg I defined lp and elp like this:

transformed parameters {
  vector[N] lp;
  vector<lower=0>[N] elp;
  lp = X * beta + Z * alpha + trt * beta_trt;
  elp = exp(lp) ;
}

If you want to try modifying the model, you can get the stan code from the psborrow object, modify it and then compile it separately and then sample. If you find something that works better, we can integrate it :)

simulation_obj@model_string # copy this to a file and edit
new_model <- cmdstanr::cmdstan_model("edited_model.stan")
results <- new_model$sample(data = psborrow2::prepare_stan_data_inputs(simulation_obj), 
  init = list(list(beta_trt = 0,  ...))
 )
results$summary()

Answer 4 · 2023-06-28T15:47:42.000Z

@souravmukherjee2906 To reduce divergences we also need to consider priors. You could try reducing the variance of the normal priors and a less extreme prior for borrowing, eg gamma(0.1, 0.04)

Answer 5 · 2023-06-30T00:43:52.000Z

Hi @gravesti , Thank you so much for your reply.

I am getting the warning even after specifying the initial values and defining an extra transformed_param_stan_code block in the exp_surv_dist() function where I have used your suggested transformed parameters block in the transformed_param_stan_code. I also reduced the variances in the normal prior and used less extreme prior for borrowing.

I agree that it's okay to avoid the warning on inverse scale parameter if it's showing in the warmup phase, but the problem is I am still getting a divergence message warning apart from that warning, for some chains in the sampling phase. Even if the divergence is occurring for few instances only, that's enough to make the posterior estimates unreliable. Because of this reason, I am getting weird Type 1 Error values.

For example, when drift HR = 1, I am getting highest Type 1 Error (~ 0.06) for "No Borrowing" compared to all other borrowing scenarios (where all Type 1 errors are at 0.04).

I have few more questions:

Does psborrow2 have problem in handling too many covariates? or, types of covariates being used in the data generation? e.g. I have 8 covariates (7 categorical covariates and 1 continuous 'age' covariate)
In the overview_script_simulation.R, Type 1 error for No Borrowing is being calculated using two different lists of data matrices, corresponding to two different scenario of drift HR. Because of this, Type I error values for No borrowing are different in two different drift HR scenario, whereas theoretically it should be same. Could you please let me know if there is any way to fix it? One simple way I can think of is to forcefully use either one of the two list of data matrices to calculate Type 1 error for both cases of drift HR, but is there any way to have same Type 1 error even for both drift HR cases?
What's the fundamental difference between the way of calculating power and type I error with using create_analysis_obj() and create_simulation_obj() in mcmc_sample()?
How to extract R hat from an output of mcmc_sample() where a simulation_obj (created with create_simulation_obj()) is used as an input instead of an analysis_obj?
overview_script_simulation.R only considers the final analysis ("time" and "cnsr" at final) with only one look. What if I have an interim look before the final analysis. Let's say, I also have time-to-event and censoring indicator values at interim with "timeia" and "cnsria". How to calculate the power and type 1 error in the interim analysis by integrating alpha spending function at interim? Also, if the interim is unsuccessful and we don't decide to stop at interim, how to calculate the additional power and type 1 error at the final analysis? Any insights / thoughts on this would be really helpful for me.

Thank you so much for your help.

Answer 6 · 2023-06-30T12:24:13.000Z

Hi @souravmukherjee2906

We might need some input from @mattsecrest on simulations, but I'll try to help.

To help more with divergences, it would be good if you can share a reproducible example.

I'm not aware of a particular problem with large number of covariates. If you have an example simulation, we can investigate further.
As you describe we are getting two estimates from simulations for the same theoretical quantity. If you don't want this, you could break the simulation into two parts:

drift = 1 with No/Full/DBD
drift = 0.6 with Full/BDB.

mcmc_sample() does not calculate power or type I error at all. For analysis objects it merely wraps $sample() from cmdstanr and for simulation objects it calls mcmc_sample on each of the analysis objects in the simulation and summarises the coverage, bias, MSE. I do not believe power or type I error can be calculated from a single model fit.
For simulation objects, there is an argument keep_cmd_stan_models, however this will use a lot of memory to save all objects. Alternatively, you can simply sample the simulation iteration that you are interested in.

my_simulation_object@guide
this_sim_iter <- my_simulation_object@analysis_obj_list[[1]][[1]] # first list is the row from guide, second level is the data
mcmc_sample(this_sim_iter)

I'm afraid I don't have much insight on interim analyses. As is I don't think this is possible, although we have some of the pieces. You could set up the simulation with "timeia" and "cnsria" and then set the posterior_quantiles in mcmc_sample which changes the quantiles which determine if the null hypothesis is rejected to correspond to your alpha spending. With this I think you can at least calculate the power/t1e at interim.

I'm not sure yet how we could calculate values for the overall procedure.

Answer 7 · 2023-07-13T02:50:07.000Z

Hi @gravesti, Thank you so much for your detailed response. It was really helpful.

Is there any way to calculate the number of events (cnsr = 0 & ext = 1) being borrowed from the external data in BDB-conservative and BDB-aggressive methods for both the cases of drift HR = "No drift HR" and drift HR = "Moderate drift HR" in your https://github.com/Genentech/psborrow2/blob/main/vignettes/overview_script_simulation.R ?