Initial time in ode_integrate?

[elisa.b...@universitadipavia.it]

Hi! 

I've tried to implement a very simple PK single subject model with one compartment and i.v. bolus dose using the ODE solver integrate_ode. 

The model is:

log(V) ~normal(log(muv), sigma2k)

log(k) ~ normal(log(muk), sigma2k)

sigma2e^-1 ~ gamma(a,b) 

for (i in 1:ntimes){

z[i] ~ normal(c[i], sigma2e)

}

where c[i]<- q[i]/V  (V is the apparent volume of distribution) and q[i] comes from the resolution of this ODE:

dq/dt=-k*q (k is the rate elimination constant and q represents the amount in the compartment)

The model implemented in stan (one_comp_ode.stan) and the R code (one_comp_ode.R) are in the attachments.

Problems:

 - Why does Stan ask me to set a t0 less than 0?

I try to set t0 less than 0, then I've eliminated the first measure and the first time instance and reset t0 equal to 0 but the algorithm seems to reject always the proposal

(and it is extremely slow). 

I would ask if you can tell me where is the problem in my code or if there is another way to do the same thing. 

Thanks. 

Elisa

[Bob Carpenter]

We should make this clearer in the manual, because
you're not the only one who's raised this question.

* t0 is the initial time
* y0 is the initial state
* ts is the array of requested solution times

There is a requirement that ts[i] > 0 (strictly) for all i.

If there's measurement error on the initial state y0,
you need to make y0 a parameter and then model the measurement
error directly.

Otherwise, the model looks OK to me --- you have all the necessary
constraints on parameters. You could make sigmae
and/or theta local variables in the model, unless you want to save
sigmae in addition to taue. I don't see any reason to save theta since
it's just an array containing k.

- Bob

> On Dec 10, 2014, at 5:57 AM, elisa.b...@universitadipavia.it wrote:
>
> Hi!
> I've tried to implement a very simple PK single subject model with one compartment and i.v. bolus dose using the ODE solver integrate_ode.
>
> The model is:
> log(V) ~normal(log(muv), sigma2k)
> log(k) ~ normal(log(muk), sigma2k)
> sigma2e^-1 ~ gamma(a,b)

> z[i] ~ normal(c[i], sigma2e)
>

> c[i] comes from this ODE:
>
> dq(t)/dt=-k*q/t) (k is the rate elimination constant and q represents the amount in the compartment)
>
> c(t)=q(t)/V (V is the apparent volume of distribution)

>
> The model implemented in stan (one_comp_ode.stan) and the R code (one_comp_ode.R) are in the attachments.
>
> Problems:
> - Why does Stan ask me to set a t0 less than 0?
> I try to set t0 less than 0, then I've eliminated the first measure and the first time instance and reset t0 equal to 0 but the algorithm seems to reject always the proposal
> (and it is extremely slow).
>
> I would ask if you can tell me where is the problem in my code or if there is another way to do the same thing.
> Thanks.
>
> Elisa
>
>
>
>

> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
> <one_comp_ode.stan><one_comp_ode.R>

[Elisa Borella]

Even if I translate the time vector of a certain constant, i.e. t=t+10, Stan requires me to set a t0 minor then 10. So, is it right that t0 can't be fixed at the first element of the time vector? 

As regarding the initial condition y0, in my case I want to fix the amount in the compartment equal to the dose amount to simulate a bolus injection. Neverthless when I run the script, the algorithm seems not to move from the first iteration also after a long time. 

This is the message that appears:

Iteration:    1 / 10000 [  0%]  (Warmup)

Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

stan::prob::normal_log(N4stan5agrad3varE): Location parameter is 1.#INF:0, but must be finite!

If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

I've attached the scripts with the changes I've made. 

2014-12-10 20:26 GMT+01:00 Bob Carpenter <ca...@alias-i.com>:

We should make this clearer in the manual, because
you're not the only one who's raised this question.

  * t0 is the initial time
  * y0 is the initial state
  * ts is the array of requested solution times

There is a requirement that ts[i] > 0 (strictly) for all i.

If there's measurement error on the initial state y0,
you need to make y0 a parameter and then model the measurement
error directly.

Otherwise, the model looks OK to me --- you have all the necessary
constraints on parameters.  You could make sigmae
and/or theta local variables in the model, unless you want to save
sigmae in addition to taue.  I don't see any reason to save theta since
it's just an array containing k.

- Bob

> On Dec 10, 2014, at 5:57 AM, elisa.borella02@universitadipavia.it wrote:
>
> Hi!
> I've tried to implement a very simple PK single subject model with one compartment and i.v. bolus dose using the ODE solver integrate_ode.
>
> The model is:
> log(V) ~normal(log(muv), sigma2k)
> log(k) ~ normal(log(muk), sigma2k)
> sigma2e^-1 ~ gamma(a,b)
> z[i] ~ normal(c[i], sigma2e)
>
> c[i] comes from this ODE:
>
> dq(t)/dt=-k*q/t) (k is the rate elimination constant and q represents the amount in the compartment)
>
> c(t)=q(t)/V (V is the apparent volume of distribution)
>
> The model implemented in stan (one_comp_ode.stan) and the R code (one_comp_ode.R) are in the attachments.
>
> Problems:
>  - Why does Stan ask me to set a t0 less than 0?
> I try to set t0 less than 0, then I've eliminated the first measure and the first time instance and reset t0 equal to 0 but the algorithm seems to reject always the proposal
> (and it is extremely slow).
>
> I would ask if you can tell me where is the problem in my code or if there is another way to do the same thing.
> Thanks.
>
> Elisa
>
>
>
>
> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.

> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.

> For more options, visit https://groups.google.com/d/optout.
> <one_comp_ode.stan><one_comp_ode.R>

--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/GsjoA9NL1z0/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

For more options, visit https://groups.google.com/d/optout.

--

Elisa Borella

[Bob Carpenter]

> On Dec 11, 2014, at 4:21 AM, Elisa Borella <elisa.b...@universitadipavia.it> wrote:
>
> Even if I translate the time vector of a certain constant, i.e. t=t+10, Stan requires me to set a t0 minor then 10. So, is it right that t0 can't be fixed at the first element of the time vector?

Correct. All elements of the time vector must be strictly
greater than t0.

> As regarding the initial condition y0, in my case I want to fix the amount in the compartment equal to the dose amount to simulate a bolus injection.

That's what we've done, as well.

That means that y0 is the value for t0. There's nothing
to solve, which is why it's not in the list of times.

> Neverthless when I run the script, the algorithm seems not to move from the first iteration also after a long time.
> This is the message that appears:
>
> Iteration: 1 / 10000 [ 0%] (Warmup)
>
> Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
> stan::prob::normal_log(N4stan5agrad3varE): Location parameter is 1.#INF:0, but must be finite!
> If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
> but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

This message is saying that you have a normal(...) somewhere that
is getting an infinite mean. Sorry we don't have line-by-line
error messages yet --- it's on our to-do list.

It's hard to say more without seeing the model. You might want to
try running in diagnostic mode just to make sure your init is
well formed.

Then, try running with a small initial step size, like 0.001 to see
if that helps you move. It controls the discretization of the
Hamiltonian simulation (itself solving an ODE within an iteration),
and making it smaller lets it follow the curvature of the posterior
more closely.

A problem that can arise is that the parameters will wander into
values that cause the ODE to become stiff. Stan's going to grind to
a halt when that happens, at least until we add a stiff solver, which
is one of our highest priorities at the moment and should happen in
the next six to nine months if not sooner. The smaller initial step
size can help (as can a higher target acceptance rate for adaptation,
which will keep the step size smaller).

Also, the more informative you can make the priors, the better --- it
will help prevent the parameters from wandering into stiff territory.

- Bob

[Elisa Borella]

1. It's hard to say more without seeing the model. -> The model is described in the first mail where I provided also the .stan and .R codes. 

2.the more informative you can make the priors -> Setting the priors more strictly (not only the lower bound at 0 but also the upper bound), the algorithm runs, but despite this the parameters estimated are different from those used to simulate the data (so the real parameters). I implemented the model also in another software for Bayesian analysis (using always a ODE solver) and in that case the parameters are estimated correctly. Couldn't be the integrate_ode function the responsible for this? If I don't use the ODE, but the analytical solution (c(t)=D/V * exp(-k*T) ) the parameters are estimated correctly! 

Here is a comparison between the estimated parameters in both the cases and the real ones:

Parameters estimated with ODE solver:

Inference for Stan model: one_comp_ode.

1 chains, each with iter=10000; warmup=5000; thin=1; 

post-warmup draws per chain=5000, total post-warmup draws=5000.

        mean se_mean    sd   2.5%    25%    50%    75%  97.5% n_eff Rhat

V       8.53    0.27 12.14   0.75   2.55   4.97   9.91  38.42  2030    1

k[1]    1.67    0.00  0.03   1.61   1.65   1.67   1.69   1.74  2022    1

taue    0.05    0.00  0.02   0.01   0.03   0.04   0.06   0.10  2333    1

lp__  -14.13    0.03  1.45 -17.69 -14.90 -13.80 -13.07 -12.29  1828    1

Parameters estimated with analytical solution:

Inference for Stan model: one_comp.

3 chains, each with iter=14000; warmup=5000; thin=1; 

post-warmup draws per chain=9000, total post-warmup draws=27000.

       mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat

V      5.02    0.00 0.17   4.68   4.91   5.01   5.12   5.38  9037    1

k      2.07    0.00 0.19   1.71   1.94   2.06   2.18   2.50  8804    1

taue   0.08    0.00 0.04   0.02   0.05   0.07   0.10   0.18  7806    1

lp__ -16.89    0.02 1.47 -20.66 -17.56 -16.50 -15.81 -15.22  6205    1

Real values: 

V 5.022

k 2.063

taue 0.07909

3. try running with a small initial step size -> If I would set the integration step where should I put the value? In an argument of the integrate_ode function? Which one? (I didn't found nothing about this on the manuals)

Thanks for the help. If you need the scripts, just ask it. 

--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/GsjoA9NL1z0/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

For more options, visit https://groups.google.com/d/optout.

--

Elisa Borella

[Bob Carpenter]

> On Dec 11, 2014, at 2:03 PM, Elisa Borella <elisa.b...@universitadipavia.it> wrote:
>
> 1. It's hard to say more without seeing the model. -> The model is described in the first mail where I provided also the .stan and .R codes.

Thanks. I forgot you'd posted it already.

I think the problem is in your definition of the ODE system.

Your code is this:

dydt[1] <- -theta[1] * y[1];

But your math is this:

> dq(t)/dt=-k*q/t) (k is the rate elimination constant and q represents the amount in the compartment)

So it looks like you're missing the division by t in the system definition.

> 2.the more informative you can make the priors -> Setting the priors more strictly (not only the lower bound at 0 but also the upper bound), the algorithm runs, but despite this the parameters estimated are different from those used to simulate the data (so the real parameters). I implemented the model also in another software for Bayesian analysis (using always a ODE solver) and in that case the parameters are estimated correctly. Couldn't be the integrate_ode function the responsible for this? If I don't use the ODE, but the analytical solution (c(t)=D/V * exp(-k*T) ) the parameters are estimated correctly!

My guess is that it's in the coding of the ODE system, not in
the ODE solver.

As we discussed earlier, you also shouldn't have t = 0 in your
array of solution times if t0=0 is your initial time. If you
want to model noise in the initial state, you need a separate
statement in Stan.

> Here is a comparison between the estimated parameters in both the cases and the real ones:
>

...

> 3. try running with a small initial step size -> If I would set the integration step where should I put the value? In an argument of the integrate_ode function? Which one? (I didn't found nothing about this on the manuals)

I meant the initial step size for the Hamiltonian integrator we use
for sampling on the outside, not for the ode_integrate() function on the
inside. The step size determines how long the discretization interval
is for the Hamiltonian simulation over the parameters. How to set it
depends on whether you're using RStan, PyStan, CmdStan, etc. In RStan, it's
another argument to stan() or sampling():

control = list(stepsize=0.001)

You may also find it helpful to up the adaptation acceptance target to keep
the step size low at convergence, maybe something like:

control = list(stepsize=0.001, adapt_delta=0.95)

This is all independent of the step size in the Runge-Kutta integrator used
for integrate_ode(). That one can't be controlled yet, but we'll make it
configurable in the future.

- Bob

[Elisa Borella]

> dq(t)/dt=-k*q/t) (k is the rate elimination constant and q represents the amount in the compartment)

Sorry, I typed wrong; I meant dq(t)/dt=-k*q(t). Now I can try to transate time vector, estimate also the initial condition and set the step. I'll give you a feedback as soon as I have the results. 

2014-12-12 23:23 GMT+01:00 Bob Carpenter <ca...@alias-i.com>:

- Bob

--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/GsjoA9NL1z0/unsubscribe.

To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

For more options, visit https://groups.google.com/d/optout.

--

Elisa Borella