Stan Model of the Week: PK Calculation of IV and Oral Dosing
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
We’re going to be starting a Stan “model of the P” (for some time period P) column, so I thought I’d kick things off with one of my own. I’ve been following the Wingvoet blog, the author of which is identified only by the Blogger handle Wingfeet; a couple of days ago this lovely post came out:
Wingfeet’s post implemented an answer to question 6 from chapter 6 of problem from Rowland and Tozer’s 2010 book, Clinical Pharmacokinetics and Pharmacodynamics, Fourth edition, Lippincott, Williams & Wilkins.
So in the grand tradition of using this blog to procrastinate, I thought I’d take a break from C++ coding and paper writing and translate Wingfeet’s model to Stan. First, let me repeat Wingfeet’s model in JAGS, exactly as posted:
tau <- 1/pow(sigma,2)
sigma ~ dunif(0,100)
# IV part
kIV ~dnorm(.4,1)
cIV ~ dlnorm(1,.01)
for (i in 1:nIV) {
predIV[i] <- c0*exp(-k*timeIV[i]) +cIV*exp(-kIV*timeIV[i])
concIV[i] ~ dnorm(predIV[i],tau)
}
c0 <- doseIV/V
V ~ dlnorm(2,.01)
k <- CL/V
CL ~ dlnorm(1,.01)
AUCIV <- doseIV/CL+cIV/kIV
# oral part
for (i in 1:nO) {
predO[i] <- c0star*(exp(-k*timeO[i])-exp(-ka*timeO[i]))
concO[i] ~ dnorm(predO[i],tau)
}
c0star <- doseO*(ka/(ka-k))*F/V
AUCO <- c0star/k
F ~ dunif(0,1)
ka ~dnorm(.4,1)
ta0_5 <- log(2) /ka
t0_5 <- log(2)/k
I still find JAGS models pretty hard to follow on their own, but Wingfeet provided the data, and the call to R2Jags, so I could work backwards from there. Here's what I came up with as a translation to Stan.
data {
int nIV;
int nOral;
real doseIV;
real doseOral;
vector[nIV] timeIV;
vector[nIV] concIV;
vector[nOral] timeOral;
vector[nOral] concOral;
}
parameters {
real CL;
real V;
real F;
real ka;
real kIV;
real cIV;
real sigma;
}
transformed parameters {
real k;
real c0;
real AUCIV;
real c0star;
real AUCOral;
real ta0_5;
real t0_5;
vector[nIV] predIV;
vector[nOral] predOral;
k <- CL / V;
c0 <- doseIV / V;
AUCIV <- doseIV / CL + cIV / kIV;
c0star <- doseOral * (ka / (ka - k)) * F / V;
AUCOral <- c0star / k;
ta0_5 <- log(2) / ka;
t0_5 <- log(2) / k;
predIV <- c0 * exp(-k * timeIV) + cIV * exp(-kIV * timeIV);
predOral <- c0star * (exp(-k * timeOral) - exp(-ka * timeOral));
}
model {
// IV component
kIV ~ normal(0.4, 1);
cIV ~ normal(1,10);
V ~ normal(2,10);
CL ~ normal(1,10);
concIV ~ normal(predIV, sigma);
// oral component
ka ~ normal(0.4,1);
concOral ~ normal(predOral, sigma);
}
Do I ever appreciate the emacs mode for Stan! The Stan model is a direct translation of Wingfeet's JAGS model, with exactly the same priors (Stan uses a sd parameterization of the normal, whereas JAGS follows bugs in using inverse variance). I moved all the deterministic nodes in the JAGS model to the transformed parameters block in the Stan model so we could inspect them in the output. I also vectorized the prediction calculations, which seems much cleaner than using loops, but opinions may vary on this point.
It was critical to put the lower-bounds on all the concentration and volume parameters in Stan, as well as upper and lower bounds on the noise term (which had a broad uniform distribution in the JAGS model).
Here's the data in Stan's dump format, which I dropped in a file named data.R:
nIV <- 12 nOral <- 11 doseIV <- 500 doseOral <- 500 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) timeOral <- c(0.5, 1, 1.5, 2, 4, 6, 10, 16, 24, 32, 48) concOral <- c(2.4, 3.8, 4.2, 4.6, 8.1, 5.8, 5.1, 4.1, 3, 2.3, 1.3)
Here's the code to run it in RStan, with number of warmup iterations and sampling iterations matching Wingfeet's call to the R2Jags package.
library('rstan');
source('data.R');
fit <- stan('pk_iv_oral.stan',
data=c("nIV","nOral","doseIV","doseOral",
"timeIV","concIV","timeOral","concOral"),
chains=4, warmup=5000, iter=14000);
It took about 10s to compile the model and 8s to run 56K iterations in the release version of RStan 2.2.0 (on my several years old Macbook Pro with 2.3GHz Intel Core i7). This uses the default diffuse initializations drawn uniformly on (-2,2) on the unconstrained scale corresponding to roughly (0.1, 10) on the positive-constrained scale
The model converged almost instantly and mixes very well, so this number of iterations is rather an overkill for Stan.
Inference for Stan model: pk_iv_oral.
4 chains, each with iter=14000; warmup=5000; thin=1;
post-warmup draws per chain=9000, total post-warmup draws=36000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
CL 2.41 0.00 0.19 2.07 2.28 2.40 2.53 2.81 16552 1
V 53.42 0.02 2.31 48.73 51.95 53.45 54.93 57.90 15682 1
F 0.86 0.00 0.05 0.75 0.82 0.86 0.89 0.96 14672 1
ka 0.62 0.00 0.10 0.46 0.56 0.62 0.68 0.84 15872 1
kIV 2.28 0.00 0.49 1.36 1.94 2.26 2.60 3.30 13662 1
cIV 10.99 0.02 2.53 6.47 9.27 10.84 12.53 16.48 12323 1
sigma 0.59 0.00 0.12 0.41 0.50 0.57 0.65 0.87 11027 1
k 0.05 0.00 0.00 0.04 0.04 0.04 0.05 0.05 15355 1
c0 9.38 0.00 0.41 8.64 9.10 9.35 9.62 10.26 15374 1
AUCIV 213.50 0.13 16.30 182.75 202.70 213.04 223.79 246.72 16616 1
c0star 8.69 0.01 0.68 7.49 8.23 8.64 9.08 10.18 14116 1
AUCOral 192.83 0.10 14.23 164.46 183.62 192.85 201.99 220.96 22358 1
ta0_5 1.14 0.00 0.17 0.83 1.02 1.12 1.24 1.51 15611 1
t0_5 15.46 0.01 1.49 12.63 14.48 15.43 16.41 18.50 15824 1
predIV[1] 14.33 0.00 0.52 13.25 14.00 14.34 14.68 15.32 18211 1
predIV[2] 12.63 0.00 0.33 11.98 12.41 12.63 12.84 13.28 36000 1
predIV[3] 11.46 0.00 0.34 10.82 11.23 11.45 11.68 12.17 21663 1
predIV[4] 9.17 0.00 0.33 8.60 8.94 9.14 9.36 9.90 13319 1
predIV[5] 8.72 0.00 0.31 8.17 8.51 8.69 8.90 9.40 14229 1
predIV[6] 7.83 0.00 0.27 7.31 7.64 7.82 8.00 8.40 16367 1
predIV[7] 7.15 0.00 0.24 6.69 6.98 7.14 7.30 7.65 17242 1
predIV[8] 5.96 0.00 0.23 5.52 5.82 5.96 6.11 6.41 17905 1
predIV[9] 4.55 0.00 0.25 4.04 4.39 4.55 4.71 5.03 17120 1
predIV[10] 3.17 0.00 0.27 2.63 3.00 3.18 3.35 3.69 16452 1
predIV[11] 2.22 0.00 0.26 1.70 2.05 2.22 2.39 2.72 16205 1
predIV[12] 1.09 0.00 0.20 0.71 0.95 1.08 1.21 1.49 16090 1
predOral[1] 2.11 0.00 0.21 1.73 1.97 2.10 2.24 2.55 21720 1
predOral[2] 3.60 0.00 0.29 3.04 3.41 3.59 3.78 4.20 24051 1
predOral[3] 4.65 0.00 0.31 4.04 4.45 4.64 4.84 5.27 28282 1
predOral[4] 5.37 0.00 0.30 4.77 5.18 5.37 5.56 5.96 32197 1
predOral[5] 6.46 0.00 0.28 5.91 6.28 6.46 6.64 7.01 26222 1
predOral[6] 6.37 0.00 0.31 5.74 6.17 6.37 6.57 6.98 19048 1
predOral[7] 5.49 0.00 0.31 4.89 5.29 5.49 5.69 6.10 18489 1
predOral[8] 4.21 0.00 0.26 3.68 4.04 4.21 4.37 4.70 21074 1
predOral[9] 2.93 0.00 0.23 2.45 2.79 2.94 3.09 3.38 21675 1
predOral[10] 2.05 0.00 0.22 1.60 1.91 2.05 2.19 2.47 20102 1
predOral[11] 1.00 0.00 0.17 0.68 0.89 1.00 1.11 1.34 18362 1
lp__ -9.22 0.03 2.31 -14.78 -10.46 -8.82 -7.52 -5.94 7951 1
Samples were drawn using NUTS(diag_e) at Mon Mar 10 14:49:24 2014.
Here's the fit that Wingfeet got from JAGS:
Inference for Bugs model at "C:/...BEnL/model1b6027171a7a.txt", fit using jags,
4 chains, each with 14000 iterations (first 5000 discarded), n.thin = 2
n.sims = 18000 iterations saved
mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat n.eff
AUCIV 219.959 16.975 189.743 208.636 218.846 230.210 256.836 1.002 2800
AUCO 200.468 15.102 171.908 190.523 199.987 209.928 231.722 1.001 7500
CL 2.346 0.182 1.995 2.226 2.344 2.461 2.712 1.002 2900
F 0.876 0.052 0.773 0.841 0.876 0.911 0.976 1.002 3200
V 56.463 2.752 51.610 54.583 56.280 58.112 62.465 1.002 3100
c0 8.876 0.426 8.005 8.604 8.884 9.160 9.688 1.002 3100
c0star 8.314 0.615 7.119 7.911 8.304 8.704 9.574 1.002 2000
cIV 11.229 2.485 7.029 9.511 11.003 12.654 16.804 1.003 1100
k 0.042 0.004 0.034 0.039 0.042 0.044 0.050 1.002 2200
kIV 2.064 0.510 1.132 1.704 2.041 2.396 3.126 1.004 800
ka 0.659 0.105 0.489 0.589 0.646 0.715 0.903 1.002 3200
predIV[1] 14.343 0.532 13.240 14.009 14.355 14.690 15.378 1.002 2500
predIV[2] 12.637 0.331 11.990 12.422 12.632 12.851 13.298 1.001 12000
predIV[3] 11.435 0.357 10.755 11.196 11.427 11.667 12.147 1.002 1600
predIV[4] 8.923 0.326 8.333 8.709 8.902 9.116 9.638 1.002 2900
predIV[5] 8.410 0.298 7.845 8.213 8.401 8.594 9.021 1.001 14000
predIV[6] 7.522 0.286 6.943 7.340 7.525 7.713 8.064 1.001 5900
predIV[7] 6.910 0.255 6.385 6.749 6.915 7.077 7.399 1.001 6200
predIV[8] 5.848 0.222 5.406 5.704 5.849 5.993 6.285 1.001 7800
predIV[9] 4.557 0.231 4.098 4.409 4.555 4.707 5.012 1.001 4500
predIV[10] 3.270 0.252 2.781 3.107 3.267 3.433 3.773 1.002 3100
predIV[11] 2.349 0.251 1.867 2.183 2.343 2.509 2.865 1.002 2700
predIV[12] 1.217 0.207 0.839 1.076 1.204 1.343 1.662 1.002 2500
predO[1] 2.138 0.214 1.750 1.995 2.124 2.263 2.599 1.001 8700
predO[2] 3.626 0.293 3.070 3.432 3.616 3.807 4.234 1.001 12000
predO[3] 4.653 0.306 4.050 4.452 4.652 4.847 5.264 1.001 16000
predO[4] 5.351 0.294 4.754 5.161 5.354 5.541 5.930 1.001 15000
predO[5] 6.375 0.277 5.801 6.202 6.383 6.562 6.894 1.002 3200
predO[6] 6.274 0.308 5.629 6.079 6.286 6.485 6.842 1.002 2800
predO[7] 5.455 0.292 4.852 5.270 5.461 5.648 6.018 1.002 3900
predO[8] 4.262 0.245 3.764 4.106 4.265 4.423 4.736 1.001 8600
predO[9] 3.057 0.227 2.605 2.910 3.058 3.207 3.504 1.001 8200
predO[10] 2.195 0.218 1.773 2.050 2.192 2.335 2.638 1.001 5200
predO[11] 1.135 0.180 0.805 1.013 1.127 1.247 1.519 1.002 3500
t0_5 16.798 1.713 13.830 15.655 16.652 17.770 20.617 1.002 2200
ta0_5 1.077 0.164 0.768 0.969 1.072 1.177 1.419 1.002 3200
deviance 38.082 5.042 30.832 34.357 37.207 40.918 50.060 1.003 1200
It's evident Stan did a way better job of mixing given the effective sample sizes. Wingfeet reports that the "correct" answers (presumably given in the book) were the following, here reported with the 95% posterior intervals produced by Stan and JAGs (we've been meaning to reduce RStan's posterior intervals to 90% to match CmdStan's, but it hasn't happened yet):
Parameter Correct Stan 95% JAGS 95% --------- ----------- ------------ ------------ t0_5 16.7 hr (12.6, 18.5) (13.8, 20.6) AUCiv 217 mg-hr/L ( 183, 247) ( 190, 256) AUCoral 191 mg-hr/L ( 164, 221) ( 172, 232) CL 2.3 L/hr ( 2.1, 2.8) ( 2.0, 2.7) V 55 L (48.8, 57.9) (51.6, 62.5) F 0.88 (0.75, 0.96) (0.77, 0.98)
I'm not sure what's up with the rather different posterior intervals from JAGS and Stan, as both samplers had enough effective samples to estimate these reasonably well (for instance, JAGS had 2200 effective samples and Stan 15,800 for t0_5). The width of the intervals is similar, so I don't think it was an issue of one system or the other not getting into the tails of the distributions. Of course, these are fairly wide posteriors (as noted by Wingfeet) and the "correct" answers are within the 95% posterior intervals for both systems.
The only "suspicious" parameter is F, which was pinned very near its upper bound of 1.0 in Stan for many samples and has a very wide and skewed posterior; here's a histogram:
The post Stan Model of the Week: PK Calculation of IV and Oral Dosing appeared first on Statistical Modeling, Causal Inference, and Social Science.
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.