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.
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 GaussKronrod 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.
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
Using integration by parts, we can get
This means we can compute I_2 := \int_{0}^t x^{\gamma} \exp(\alpha\beta_1 x)\,dx and then get I via
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) submodeltwophistep_new.stan (2.8 KB)
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.
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.
Phew, that is a relief. Those integrate_1d
problems were frustrating.
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.
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 cmdstan2.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://mcstan.org/rpackages/", getOption("repos")))
install.packages("rstan", repos = c("https://mcstan.org/rpackages/", getOption("repos")))
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 1e6) 
error
(double
, defaults to 1e6)
I have BH
1.75.00 (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(
"submodeltwopsistep.stan"
)
cmd_prefit_psi_step < cmdstan_model(
"submodeltwopsistep.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
)
submodeltwopsistep.stan (2.7 KB)
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…)
Ooo, so the integrate_1d code in standev/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)
Downgrading BH
to 1.72.03
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.750.