**Statistical Modeling, Causal Inference, and Social Science » R**, and kindly contributed to R-bloggers)

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 { intnIV; 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.

**leave a comment**for the author, please follow the link and comment on his blog:

**Statistical Modeling, Causal Inference, and Social Science » R**.

R-bloggers.com offers

**daily e-mail updates**about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...