Stan simple model for pharmacokinetics

Hello everyone. I am new to Stan and R.
I am trying to run a model for a 1 compartment model in pharmacokinetics.
I found an example on the web and adapted the model .
The model code is saved as stan_model_code in the working directory

data {
        int<lower=0> nIV;
        real<lower=0> doseIV;
        vector<lower=0>[nIV] timeIV;
        vector<lower=0>[nIV] concIV;
      }
      parameters {
        real<lower=0> CL;
        real<lower=0> V;
        real<lower=0> kIV; 
        real<lower=0> cIV; 
        real<lower=0,upper=100> sigma;
      }
      transformed parameters {
        real k;
        real c0;
        real AUCIV;
        real t0_5;
        vector[nIV] predIV;
        k = CL / V;
        c0 = doseIV / V;
        AUCIV = doseIV / CL + cIV / kIV;
        t0_5 = log(2) / k;
        predIV = c0 * exp(-k * timeIV) + cIV * exp(-kIV * timeIV);
      }
      
      model {
        kIV ~ normal(0.4, 1);
        cIV ~ lognormal(3,10);
        V ~ lognormal(2,10);
        CL ~ lognormal(1,10);
        concIV ~ normal(predIV, sigma);
        }

} 

The data file below is named data.R and saved in the working folder

nIV <- 12
doseIV <- 1000
timeIV <- 
  c(0.33, 0.5, 0.67, 1.5, 2, 4, 6, 10, 16, 24, 32, 48)
concIV <- 
  c(14.7, 12.6, 11, 9, 8.2, 7.9, 6.6, 6.2, 4.6, 3.2, 2.3, 1.2)

The main R script below

library(rstan)
source('data.R');

fit <- stan('stan_model_code.stan', data=c("nIV","doseIV", 
            "timeIV","concIV"),
            chains=4, warmup=5000, iter=14000);

The errors I get:
Warning messages:
1: There were 31809 divergent transitions after warmup. See


to find out why this is a problem and how to eliminate them.
2: Examine the pairs() plot to diagnose sampling problems

Any ideas what is wrong?
I adapted the code from:


I eliminated the oral dosing stuff from the example because i don’t need it.
I also ran the full example from the link and got the same errors.
Thank you for any help and also a big thanks to the Stan developers. I was looking for weeks for this kind of software, now I just need to learn to use it
1 Like

Welcome @dan_r. These models are a bit outside my area but a few things:

The divergent transition primer is here:

and contains lots of good information.

Here Bounding posterior output of pharmacokinetic model - #30 by dbarneche some discussion about the same type of model.

Those two posts might get you going in the right direction. I suspect a slightly larger simulated data set will help troubleshoot the model. We generally recommend starting with simulated data so you know the true parameters.

3 Likes

Thank you for replying.
I think my problem is, I don’t understand how to chose priors with a lognormal distribution.
Let’s say that I find in a study that a certain drug has Mean Volume of distribution (V) of 22 L, a Mean Elimination Rate Constant (k) of 0.421, and a Mean Clearance (CL) of 6.6 L/h.
How do I estimate that on a lognormal distribution for my priors? Given of course that my priors are k, V and CL.
In the example above V ~ lognormal(2,10) , and CL ~ lognormal (1,10). How does that relate to a mean V of 22 L and a CL of 6.6 L/h?

1 Like

If I recall correctly lognormal priors can be a bit tricky to define. @Max_Mantei has a post here Understanding shifted lognormal parameters and priors about that and they would know more about setting this kind of priors.

1 Like

Also this post Setting informative priors from pilot study for lognormal model

1 Like

If you can’t find priors directly from domain expertise, the current best way are prior predictive checks - discussed within the recent workflow preprint and the Visualisation paper.

In the most simple approach I can just look what my priors do on the parameters. My understanding of PK is relatively shallow, so I hope I am not messing up, but say I take the drug you describe as at least a bit similar to the one I am studying, so - to stay on the safe side - I can assume V to be a priori within two orders of magnitude of 22, i.e. roughly between 0.2 and 2000. If you are good at math, you can solve (or have your computer solve) for the parameters of the lognormal that would match 95% or 99% prior interval to those values. Or you can just do a bit of trial and error (using R here to draw 30 values from a lognormal distribution with given params ):

set.seed(522364)
center <- 20
rlnorm(30, meanlog = log(center), sdlog = 5)
#  [1] 3.523340e+01 5.241080e+00 4.183498e+04 1.943986e+00 5.580033e+02 9.400129e+06 5.330368e+00 1.387772e+06 5.595159e+02 2.086984e+00
# [11] 4.520441e+02 8.767097e+01 6.771238e+01 3.838290e+02 3.606970e-01 1.013651e+01 1.483323e+05 1.497356e+04 6.265917e-01 5.003035e+00
# [21] 1.370033e+04 9.536687e-01 1.896423e+02 1.322343e-02 1.429113e-02 1.831241e+02 7.630693e+01 8.471494e+00 2.835572e-01 1.170597e+06

OK so I see a bunch of values on the order of 1e6 and there is also some at 1e-2, so my sd is probably too large, let’s try a bit smaller:

rlnorm(30, meanlog = log(center), sdlog = 4)
#  [1] 9.228254e-01 5.863253e+01 1.388731e+00 2.389572e-03 1.147771e+05 3.814239e+00 1.206882e+02 5.894352e+01 3.165645e+00 9.371014e-04
# [11] 3.983885e+00 2.498825e+01 3.555730e-01 4.597215e-01 6.055531e-02 4.430602e+01 1.066487e-01 4.154448e+05 5.076735e+01 5.629688e-01
# [21] 1.001600e+03 1.065909e+01 6.832729e+01 3.306314e+00 1.321551e+03 4.618202e-01 1.773967e+02 1.534803e+03 1.972578e+02 1.018798e+01

Now most of the values are within the range AND the whole range is represented, I would say that’s good for a start.

You could then do a next step and use the params drawn form priors to compute the concentration curve and tweak the priors until many (but not all) curves look somewhat sensible.

Does that make sense?

3 Likes

Thank you for the explanation. It is making sense now. Also I found this article very useful for beginners like me https://www.generable.com/blog/2018/07/optimizing-sampling-and-choosing-priors/.
With my priors set correctly now, I was able to get a very good model in r2jags.
I need to translate it into stan to see if it works

1 Like