Stabilising survival model integrand for use in integrate_1d

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

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

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:


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

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

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

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

3 Likes

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.

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

3 Likes

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 :/.

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

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

Yes, I am following the build instructions from the PR

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.

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)

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?

@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.

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

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

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

Downgrading BH to 1.72.0-3 solves the rstan (2.26) issue. Enormously helpful, thanks.

Probably worth making sure those PRs get merged before a binary gets out / goes to CRAN? At the moment, install.packages('rstan') into an R session with an empty library will install BH 1.75-0.

4 Likes