Stabilising survival model integrand for use in integrate_1d

Modeling
,
1

I’m trying to use integrated_1d to numerically integrate the following

I = \int_{0}^{t} u^{\gamma - 1} \exp\{\alpha \beta_{1} u\} \text{d}u

for t, \gamma > 0 and \alpha, \beta \in \mathbb{R}. Note that this domain means I cannot use the gamma_p function in Stan. For some parameter configurations, integrate_1d does not like this integrand. The only reason I can think of is that in some reasonable configurations u^{\gamma - 1} might overflow (for \gamma < 1 and small u), so I rewrote I as

I = \int_{0}^{t} \exp\{(\gamma - 1) \log(u) + \alpha \beta_{1} u\} \text{d}u

but this is also produces integrate_1d errors, and occur with such frequency that I can’t get Stan to initialise.

suppressMessages(library(cmdstanr))

cmdstan_model <- cmdstan_model( 
  "some-path-to-make-reprex-happy/submodel-two-phi-step.stan"
)
#> Compiling Stan program...
stan_data <- list(
  n_patients = 16L, 
  baseline_measurement = c(
   -0.935695521120025, -2.12761603089136, -0.876325387003469, 0.931350429946891,
    0.422731897638725, -1.30375210258178,  1.85158267797611, -0.289448169455647, 
    0.325903831278964, -1.2757407765375,   1.58794213975781,  1.59307842403321, 
   -0.0584097720952887, 1.80357260097634, -2.30428952138529,  0.655115279462301
 ), 
 log_crude_event_rate = 3.17644992456411
)

res <- cmdstan_model$sample(
  data = stan_data,
  chains = 1,
  iter_warmup = 1,
  iter_sampling = 0
)
#> Running MCMC with 1 chain...
#> 
#> Chain 1 event_time: 3.15863 hazard_gamma: 3.00293 alpha: -1.14303 beta_one: -1.82362 
#> Chain 1  surv_prob_term_two: 2578.11 
#> Chain 1 event_time: 0.181123 hazard_gamma: 3.00293 alpha: -1.14303 beta_one: 1.38322 
#> Chain 1  surv_prob_term_two: 0.00159067 
#> Chain 1 event_time: 0.529865 hazard_gamma: 3.00293 alpha: -1.14303 beta_one: 1.01519 
#> Chain 1  surv_prob_term_two: 0.0314054 
#> Chain 1 event_time: 4.02702 hazard_gamma: 3.00293 alpha: -1.14303 beta_one: -0.887469 
#> Chain 1  surv_prob_term_two: 599.159 
#> Chain 1 event_time: 0.581211 hazard_gamma: 3.00293 alpha: -1.14303 beta_one: -0.830041 
#> Chain 1  surv_prob_term_two: 0.0992707 
#> Chain 1 event_time: 1.74709 hazard_gamma: 3.00293 alpha: -1.14303 beta_one: 0.952865 
#> Chain 1  surv_prob_term_two: 0.459929 
#> Chain 1 event_time: 0.137359 hazard_gamma: 3.00293 alpha: -1.14303 beta_one: -1.10788 
#> Chain 1  surv_prob_term_two: 0.000978182 
#> Chain 1 event_time: 0.216454 hazard_gamma: 3.00293 alpha: -1.14303 beta_one: 0.00135613 
#> Chain 1  surv_prob_term_two: 0.00336121 
#> Chain 1 event_time: 4.41431 hazard_gamma: 3.00293 alpha: -1.14303 beta_one: -1.39034 
#> Chain 1  surv_prob_term_two: 10356.2 
#> Chain 1 event_time: 0.847361 hazard_gamma: 3.00293 alpha: -1.14303 beta_one: 0.840636 
#> Chain 1  surv_prob_term_two: 0.111383 
#> Chain 1 event_time: 1.12224 hazard_gamma: 3.00293 alpha: -1.14303 beta_one: 0.780646 
#> Chain 1  surv_prob_term_two: 0.226627 
#> Chain 1 event_time: 0.628329 hazard_gamma: 3.00293 alpha: -1.14303 beta_one: -1.17817 
#> Chain 1  surv_prob_term_two: 0.157637 
#> Chain 1 event_time: 0.516754 hazard_gamma: 3.00293 alpha: -1.14303 beta_one: -1.66055 
#> Chain 1  surv_prob_term_two: 0.0973701 
#> Chain 1 event_time: 4.18989 hazard_gamma: 3.00293 alpha: -1.14303 beta_one: 0.899539 
#> Chain 1  surv_prob_term_two: 1.48168 
#> Chain 1 event_time: 0.798516 hazard_gamma: 3.00293 alpha: -1.14303 beta_one: -1.63808 
#> Chain 1  surv_prob_term_two: 0.540543 
#> Chain 1 event_time: 0.204784 hazard_gamma: 3.00293 alpha: -1.14303 beta_one: -1.27114 
#> Chain 1  surv_prob_term_two: 0.00356428 
#> Chain 1 event_time: 3.15863 hazard_gamma: 3.00293 alpha: -1.14303 beta_one: -1.82362 
#> Chain 1  surv_prob_term_two: 2578.11 
#> Chain 1 event_time: 0.181123 hazard_gamma: 3.00293 alpha: -1.14303 beta_one: 1.38322 
#> Chain 1 Exception: integrate: error estimate of integral 0.000158824 exceeds the given relative tolerance times norm of integral (in '/var/folders/9q/xrp79_nd2d7dhnp5j6nncymd6zq2rx/T/Rtmp12DnEu/model-116b06733f3da.stan', line 16, column 4 to column 16)
#> Chain 1 Exception: integrate: error estimate of integral 0.000158824 exceeds the given relative tolerance times norm of integral (in '/var/folders/9q/xrp79_nd2d7dhnp5j6nncymd6zq2rx/T/Rtmp12DnEu/model-116b06733f3da.stan', line 16, column 4 to column 16)
#> Chain 1 Exception: integrate: error estimate of integral 0.000158824 exceeds the given relative tolerance times norm of integral (in '/var/folders/9q/xrp79_nd2d7dhnp5j6nncymd6zq2rx/T/Rtmp12DnEu/model-116b06733f3da.stan', line 16, column 4 to column 16)
#> Chain 1 Exception: integrate: error estimate of integral 0.000158824 exceeds the given relative tolerance times norm of integral (in '/var/folders/9q/xrp79_nd2d7dhnp5j6nncymd6zq2rx/T/Rtmp12DnEu/model-116b06733f3da.stan', line 16, column 4 to column 16)
#> Warning: Chain 1 finished unexpectedly!
#> Warning: No chains finished successfully. Unable to retrieve the fit.

Here are some examples I’ve encountered. Initially I thought the asymptote at u = 0 was going to be the issue, but this also happens when the integrand is equal to zero at u = 0.

# some prefab examples:
# Example - decaying behaviour
event_time <- 0.333964
hazard_gamma <- 0.280632
alpha <- 1.71436
beta_one <- 1.641

# Chain 1 Rejecting initial value:
# Chain 1   Error evaluating the log probability at the initial value.
# Chain 1 Exception: integrate: error estimate of integral 0.0568744 exceeds the
# given relative tolerance times norm of integral

f_integrand <- function(x) {
  x^(hazard_gamma - 1) * exp(alpha * beta_one * x)
}

f_alt <- function(x) {
  exp((hazard_gamma - 1) * log(x) + alpha * beta_one * x)
}

integrate(f_integrand, lower = 0, upper = event_time)
#> 3.338441 with absolute error < 0.00011
integrate(f_alt, lower = 0, upper = event_time)
#> 3.338441 with absolute error < 0.00011

curve(f_integrand, from = 0, to = event_time)
curve(f_alt, from = 0, to = event_time, add = T, col = "blue")


# And another - increasing behaviour
event_time <- 0.141916
hazard_gamma <- 4.43969
alpha <- -1.20662 
beta_one <- -0.81199 

#Chain 1 Exception: integrate: error estimate of integral 7.81929e-08 exceeds

f_integrand <- function(x) {
  x^(hazard_gamma - 1) * exp(alpha * beta_one * x)
}

f_alt <- function(x) {
  exp((hazard_gamma - 1) * log(x) + alpha * beta_one * x)
}

integrate(f_integrand, lower = 0, upper = event_time)
#> 4.33822e-05 with absolute error < 1e-14
integrate(f_alt, lower = 0, upper = event_time)
#> 4.33822e-05 with absolute error < 1e-14

curve(f_integrand, from = 0, to = event_time)
curve(f_alt, from = 0, to = event_time, add = T, col = "blue")

Created on 2021-01-26 by the reprex package (v0.3.0)

The Stan code is attached if you want to play around:
submodel-two-phi-step.stan (2.4 KB)

1 Like
2

I’m going to think a bit more about this problem, which is old nemesis of mine. For now, you might want to look at what others have done here in the forum.

1 Like
3

I appreciate any time/brain power you can spend on it. It does make me think that this is one reason @sambrilleman might have opted for Gauss-Kronrod quadrature, done in R, in rstanarm::stan_jm. The difference here is that my event time is uncertain, which means the region of integration varies and hence the evaluation points and weights would change each iteration.

1 Like
4

Right. I’ve given this a bit more thought but I unfortunately still can’t make it reliably work. Here’s what I have up to now: our integral of interest is

I := \int_{0}^t x^{\gamma-1} \exp(\alpha\beta_1x)\,dx.

Using integration by parts, we can get

\int_{0}^t x^{\gamma} \exp(\alpha\beta_1 x)\,dx = \frac{t^\gamma \exp(\alpha\beta_1 t)}{\alpha\beta_1} - \frac{\gamma}{\alpha\beta_1}\int_{0}^t x^{\gamma-1} \exp(\alpha\beta_1 x)\,dx.

This means we can compute I_2 := \int_{0}^t x^{\gamma} \exp(\alpha\beta_1 x)\,dx and then get I via

I = \frac{t^\gamma \exp(\alpha\beta_1 t)- \alpha\beta_1I_2}{\gamma}.

We haven’t solved much of anything, really, but notice that the integrand in I_2 is a tad nicer.
For instance, it evaluates to 0 at x=0.
Here’s a plot for the same configuration you used:

image
image902×651 29.4 KB

Notice how there’s no running off to infinity as x \to 0 and how the whole range of the function on [0, t] is narrower. I take these as signs of greater numerical robustness, but I might be wrong.
Now the problem is: I implemented this in R and got integrals with much lower error (see attached code) but when I tried to implement stuff in Stan I still got the norm convergence problem from integrate_1d.
Maybe @nhuurre knows of a trick that will work here.
Now, @bbbales2, is this not an indication that integrate_1d is brittle? Am I missing something? Like a configuration parameter that makes all of this go away?

robust_survival.r (1.0 KB) submodel-two-phi-step_new.stan (2.8 KB)

1 Like
5

This is very neat and I would think the same.

I see there is a news item for cmdstan 2.26.0 which is

Fixed a bug with integrate_1d tolerances.

so I might upgrade and give it a whirl.

2 Likes
6

I did upgrade and it seemed to make a world of difference. Both programs now finish fine with a minor warning that the integrand gives inf. For the first one it’s obvious what’s going on; for the second, I don’t know where it could be evaluating to inf, but that can be investigated later.

1 Like
7

Phew, that is a relief. Those integrate_1d problems were frustrating.

3 Likes
8

Thanks @maxbiostat and @bbbales2, using Stan 2.26 seems to solve the integrate_1d errors to a sufficient degree.

Unfortunately, for reasons specific to my use case, I need to use rstan::log_prob on the resulting Stan model. Which means I have to use rstan and hence Stan 2.21.3 :(. I guess I can have a crack at installing this experimental branch of rstan by @hsbadr that seems to incorporate 2.26.

9

For posterity - this is less difficult than I imagined. I can now use rstan::log_prob with Stan 2.26. Neat.

3 Likes
10

Never mind, I can “install” StanHeaders_2.26.0.9000 and rstan_2.26.0.9000, but using this version of rstan still gives the relative tolerance errors, but using cmdstanr with cmdstan-2.26.0 gives a small number of quadrature got inf errors. I can’t figure out what is out of sync :/.

11

@rok_cesnovar @hsbadr – any idea why I would get such different behaviour between rstan and cmdstanr using the same Stan version (2.26)?

12

Did you do?

remove.packages(c("StanHeaders", "rstan"))
install.packages("StanHeaders", repos = c("https://mc-stan.org/r-packages/", getOption("repos")))
install.packages("rstan", repos = c("https://mc-stan.org/r-packages/", getOption("repos")))
1 Like
13

Yes, I am following the build instructions from the PR

14

The only difference I can see is the different Boost library. Cmdstanr uses Boost 1.72.0. I am not sure what version Boost headers R package is at.

15

I think the binaries need to be updated.

Also, you may want to tune the control parameters in rstan:

control A named list of parameters to control the sampler’s behavior. It defaults to NULL so all the default values are used. First, the following are adaptation parameters for sampling algorithms. These are parameters used in Stan with similar names here.

  • adapt_engaged ( logical )
  • adapt_gamma ( double , positive, defaults to 0.05)
  • adapt_delta ( double , between 0 and 1, defaults to 0.8)
  • adapt_kappa ( double , positive, defaults to 0.75)
  • adapt_t0 ( double , positive, defaults to 10)
  • adapt_init_buffer ( integer , positive, defaults to 75)
  • adapt_term_buffer ( integer , positive, defaults to 50)
  • adapt_window ( integer , positive, defaults to 25)

In addition, algorithm HMC (called ‘static HMC’ in Stan) and NUTS share the following parameters:

  • stepsize ( double , positive)
  • stepsize_jitter ( double , [0,1])
  • metric ( string , one of “unit_e”, “diag_e”, “dense_e”)

For algorithm NUTS, we can also set:

  • max_treedepth ( integer , positive)

For algorithm HMC, we can also set:

  • int_time ( double , positive)

For test_grad mode, the following parameters can be set:

  • epsilon ( double , defaults to 1e-6)
  • error ( double , defaults to 1e-6)
16

I have BH 1.75.0-0 (that is the R package version)

The only relevant control parameters will be the ones controlling initialisation – I can’t get the model to initialise (in rstan) because I get the integration tolerance errors on each of the 100 initialisation attempts. AFAIK cmdstanr uses the same initialisation strategy?

17

@rok_cesnovar I was going to report a similar issue with cmdstandr, which almost always has undefined values in the first iteration; could be an CSV I/O sync or loops start from 0 vs 1.

18

Here’s a reprex and some revised Stan code. Also I think we’re getting off topic here and I should probably post this information somewhere else.

library(rstan)
library(cmdstanr)

prefit_psi_step <- stan_model(
  "submodel-two-psi-step.stan"
)

cmd_prefit_psi_step <- cmdstan_model(
  "submodel-two-psi-step.stan"
)

psi_step_data <- list(n_patients = 16L, baseline_measurement = c(-0.935695521120025, 
-2.12761603089136, -0.876325387003469, 0.931350429946891, 0.422731897638725, 
-1.30375210258178, 1.85158267797611, -0.289448169455647, 0.325903831278964, 
-1.2757407765375, 1.58794213975781, 1.59307842403321, -0.0584097720952887, 
1.80357260097634, -2.30428952138529, 0.655115279462301), log_crude_event_rate = 3.17644992456411, 
    event_time = c(3.88219330975027, 17.5607555174655, 8.14097264393642, 
    0.874471581905152, 0.875182734259227, 1.39419961264823, 1.13884032148061, 
    43.2464501785995, 1.01706489078387, 19.9898187874785, 1.1633097689907, 
    0.769646147838224, 9.44064603011716, 1.49631432330604, 22.3176589776014, 
    1.04117377075568), long_beta = structure(c(0.931386697031433, 
    0.865137767714317, 0.973505472094455, 0.880036367203131, 
    0.935119573179882, 0.910063465843014, 1.34798084760589, 0.919172808940375, 
    1.01141192543184, 0.890762064216965, 1.05317362916398, 1.22448828904961, 
    0.843895992304378, 0.872039892731229, 0.854977513431319, 
    0.958580645647191, -0.511886090506962, 0.660548718396229, 
    1.21418574954005, 1.57400898285864, 1.48149041414514, 0.176628427838126, 
    1.46367493523126, 0.437751542691176, 1.21862018076306, 0.117778624770988, 
    0.988969750444339, 0.940177413736985, 0.142150803498799, 
    1.96856722070831, 0.366973246331529, 1.02998350811324), .Dim = c(16L, 
    2L)))

res <- sampling(
  prefit_psi_step,
  data = psi_step_data,
  chains = 1,
  warmup = 100,
  iter = 110,
  control = list(adapt_delta = 0.9)
)

res <- cmd_prefit_psi_step$sample(
  data = psi_step_data,
  chains = 1,
  iter_warmup = 100,
  iter_sampling = 110,
  adapt_delta = 0.9
)

submodel-two-psi-step.stan (2.7 KB)

2 Likes
19

Sorry I’m only getting round to replying to this now. And seems like you’ve worked through solutions to your issue anyway.

But fwiw, I think I did look at integrate_1d briefly for stan_surv and/or stan_jm. But it was after a lot of the GK quadrature stuff had already been written. From memory (although the details are a bit vague for me now) the major problem was that integrate_1d had to take the single variable over which we want to integrate and then the other unknowns (i.e. input arguments) in a very structured way, and it was going to be near impossible to conform to that. In stan_jm even the longitudinal submodel itself is a function of t and so I would to have essentially had to put all of the jm.stan work, transformations, etc, inside the function definition that was going to be passed to integrate_1d. (Or something like that…)

2 Likes
20

Ooo, so the integrate_1d code in stan-dev/math hasn’t been updated to work with Boost 1.75 yet. These are the two pulls relevant pulls: here (upgrading boost) and here (removing integrate_1d patch that was needed for 1.72 but not for 1.75)

2 Likes