[This article was first published on Statistical Modeling, Causal Inference, and Social Science » R, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
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: