PKPD modeling

[Linas Mockus]

Hi,

Can somebody explain which is the right way to model non negative concentrations? In one case I had:
functions {
   real [] pbpk (...) {
...
   }
}
transformed parameters {
   real y_hat [T,N];
   y_hat = integrate_ode_bdf (pbpk, ...);
}

In another case I enforced nonnegative constraint on concentration y_hat:
transformed parameters {
   real <lower=0> y_hat [T,N];
   y_hat = integrate_ode_bdf (pbpk, ...);
}

While in the first case convergence was good but the posterior of y_hat included negative values. In the second case I am getting divergent__ all over the place and some chains didn't converge at all.
The full stan model for the second case is attached. Specifically I am concerned about urine for which upper bound of 95% credible interval is higher than the amount of drug administered.

Thanks for any advice,
Linas

[Michael Betancourt]

Your best option is to rewrite the diff eq in terms of log y

using the relationship

(d logy )/dt = 1/y dy/dt(y) = e^{- log y } dy/dt(exp(log y));

--
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.
To post to this group, send email to stan-...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
<pbpkauto.stan>

[Linas Mockus]

Michael,

Thanks a lot. Why
   ly ~ normal (log (to_vector (y_hat [ts_idx, idx])), sigma);

doesn't take care of it since only positive y_hat are considered.

Linas

[Michael Betancourt]

Because you still need to transform the ODE to be unconstrained.

[Linas Mockus]

Thanks Micael, I thought that I have to bite the bullet... Are thare any libraries/codes to do this? I am sure I will mess up when doing manually.

Thank you,

Linas

[Sebastian Weber]

Hi!

This should be easy to do. Maybe the ascii math is not so obvious, but what you need to do is to integrate log(y_n) for the nth state of the ODE instead of y_n. To do that take you current ODE system

dy_n/dt = f(y_1, ... y_n, ...)

then divide the n-th ODE rhs by y_n, so you get

1/y_n dy_n/dt = dly_n/dt =  f(y_1, ... y_n, ...)/y_n

This will make your states y_n to be log-states ly_n=log(y_n). All you have to do now is to transform from the log(y_n) to y_n just before the ODE rhs in the functor which you wrote in Stan. That's it.

However, there is a catch here - calculating things on the log scale means that you cannot represent a 0. You can only represent values close to zero which is something to consider when initializing. So you will have to call 10^-5 or something like that to be zero; which number here makes sense depends on your problem specifics. Do not make this number crazy small as this will make the system terribly hard to solve (going from nothing to something is a huge jump on the log-scale).

Sebastian

[lmockus]

Hi,

 

For some reason I don't think that we need to go to log scale. If ODEs are correct and initial concentrations are non negative the the posterior of concentration (and urine) is non negative - assuming integration time points are tight enough. I think it was my mistake to contaminate posterior of concentration with sigma (where y stands for concentration and l for log). How do you think?

generated quantities {
   real ly_pred [T];
   real lurine_pred;

 

   for (t in 1 : T)
      ly_pred [t] = normal_rng (log (y_hat [t, idx]), sigma);
   lurine_pred = normal_rng (log (y_hat [T, idx_urine]), sigma_urine);
}

 

Linas

[Sebastian Weber]

Hi!

Of course, if your residual error model allows for negative values, then you can get these.

Sebastian

[Michael Betancourt]

But so will numerical errors, both from the approximate ODE solvers

and floating point arithmetic.  The point is that if you are trying to

solve an ODE system whose states are positive but are very close

to zero then you are much better off working in log space.

[lmockus]

Michael and Sebastian,

 

Thanks for your suggestions. I already rewrote the model (and messed up somewhere - hope to figure out myself) but my concern is not with negative concentration, rather it is with elimination via urine. The model I am working with is detailed enough to capture elimination via kidney. The error model for urine is assumed to be lognormal.  Actually data supports this. While I am getting posterior of urine which makes physical sense (urine < dose administered) when I contaminate the posterior of urine with measurement error (my understanding it is called posterior predictive simulation) I am getting upper 95% confidence bound way above administered dose. Is this as expected or I should modify error model somewhat? The model is attached, see generated quantities block.

 

Thanks again for very good suggestions,

Linas

 

 

[Sebastian Weber]

This can happen, yes. BUT: Then I would not trust at all either your model or the measurement device - either one or both are not doing a good thing here.