PK calculations for infusion at constant rate

In this third PK posting I move to chapter 10, study problem 4 of Rowland and Tozer (Clinical pharmacokinetics and pharmacodynamics, 4th edition). In this problem one subject gets a 24 hours continuous dose. In many respects the Jags calculation is not very different from before, but the relation between parameters in the model and PK parameters is a bit different.

Data

Data is from infusion at constant rate, (1.5 µm/hr, but the calculations come out on the expected values with 1.5 µg/hr, and no molecular weight is given in the question) for 24 hours. Drug concentration (µg/L) in plasma is measured for 28 hours.
ch10sp4 <- data.frame(
    time=c(0,1,2,4,6,7.5,10.5,14,18,24,25,26,27,28),
    conc=c(0,4.2,14.5,21,23,19.8,22,20,18,21,18,11.6,7.1,4.1))

Model

The model is straightforward, but does suffer a bit from sensitivity of init values. The construction (time<24 and time>=24 is just a quick way to ensure that the up and down doing parts of the model sit at the correct spots in time.
library(R2jags)
datain <- as.list(ch10sp4)
datain$n <- nrow(ch10sp4)
datain$Rinf <- 1.5*1000
model1 <- function() {
  for (i in 1:n) {
    pred[i] <- Css*((1-exp(-k*time[i]))*(time[i]<24)
                  +exp(-k*(time[i]-24))*(time[i]>=24))
    conc[i] ~ dnorm(pred[i],tau)
  }
  tau <- 1/pow(sigma,2)
  sigma ~ dunif(0,100)
  Css <- Rinf/CL
  k <- CL/V
  V ~ dlnorm(2,.01)
  CL ~ dlnorm(1,.01)
  t0_5 <- log(2)/k
  Cssalt <- 3000/CL
  h1 <- Cssalt*(1-exp(-k))
  h2 <- Cssalt*(1-exp(-2*k))
  CLpermin <- CL/60
}

parameters <- c('k','CL','V','t0_5','Cssalt','h1','h2','Css','CLpermin')
inits <- function() 
  list(V=rlnorm(1,log(100),.1),CL=rlnorm(1,log(70),.1),sigma=rlnorm(1,1,.1))
jagsfit <- jags(datain, model=model1, 
    inits=inits, 
    parameters=parameters,progress.bar="gui",
    n.chains=4,n.iter=14000,n.burnin=5000,n.thin=2)
jagsfit

Inference for Bugs model at "C:/Users/.../modelf7824253314.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

CL        68.021   3.252  62.069  65.902  67.845  69.932  75.065 1.001  6600

CLpermin   1.134   0.054   1.034   1.098   1.131   1.166   1.251 1.001  6600

Css       22.102   1.046  19.983  21.449  22.109  22.761  24.167 1.001  6600

Cssalt    44.204   2.092  39.965  42.899  44.218  45.522  48.333 1.001  6600

V        170.798  23.922 127.164 155.118 169.417 184.662 222.523 1.001 18000

h1        14.679   1.698  11.572  13.575  14.586  15.685  18.367 1.001 18000

h2        24.419   2.307  20.014  22.913  24.360  25.849  29.249 1.001 18000

k          0.406   0.058   0.306   0.368   0.401   0.438   0.536 1.001 18000

t0_5       1.742   0.243   1.293   1.583   1.730   1.886   2.267 1.001 18000

deviance  66.267   3.024  62.809  64.041  65.477  67.692  74.131 1.002  2000

For each parameter, n.eff is a crude measure of effective sample size,

and Rhat is the potential scale reduction factor (at convergence, Rhat=1).

DIC info (using the rule, pD = var(deviance)/2)

pD = 4.6 and DIC = 70.8

DIC is an estimate of expected predictive error (lower deviance is better).

Expected values: CL=1.2 L/min, t_1/2=1.6 hr, V=164 L. The values are a bit off, it was given in the text that Css is 21 mg/L while the model finds 22.
Three additional parameters; Cssalt, h1, h2 correspond to part (c) of the question; what concentrations if the dosage is doubled. Expected answers: Cssalt=42, h1=14.8, h2=24.

Plot of the model

Last week I made a similar plot in classical R graphics, so this week using ggplot2. What I personally notice from a modelling perspective is that the narrowness of the band near time=0 which misses the point at time 1 h. The whole pattern of points there almost suggests a small delay and steeper increase (lower t1/2) would make for a better fit. Interestingly, the same can be seen in the decreasing part of the curve. From a data perspective, the wide variability around the plateau concentration may be important. Perhaps that particular variation is the cause of the steepness which I mentioned before. 

part (b) Is the curve in the up and down going part equally steep?

This requires a different model, since right now by definition they are equally shaped. The parameter k12, difference between k1 and k2, is the parameter of interest.
datain2 <- as.list(ch10sp4)
datain2$n <- nrow(ch10sp4)
model2 <- function() {
  for (i in 1:n) {
    pred[i] <- Css*((1-exp(-k1*time[i]))*(time[i]<24)
          +exp(-k2*(time[i]-24))*(time[i]>=24))
    conc[i] ~ dnorm(pred[i],tau)
  }
  tau <- 1/pow(sigma,2)
  sigma ~ dunif(0,100)
  k1 ~ dlnorm(0,0.001)
  k2 ~ dlnorm(0,0.001)
  Css ~ dnorm(21,.01)
  k12 <- k1-k2
}

parameters2 <- c('k1','k2','k12','Css')
jagsfit2 <- jags(datain2, model=model2,  
    parameters=parameters2,progress.bar="gui",
    n.chains=4,n.iter=14000,n.burnin=5000,n.thin=2)

Inference for Bugs model at "C:/Users/...modelf783e02e64.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
Css       21.275   1.131 19.089 20.559 21.258 21.982 23.578 1.001  6800
k1         0.529   0.138  0.333  0.446  0.512  0.588  0.814 1.002  8300
k12        0.194   0.163 -0.074  0.099  0.182  0.272  0.520 1.015  3800
k2         0.334   0.068  0.215  0.291  0.329  0.374  0.480 1.001  8600
deviance  64.761   3.663 60.215 62.073 63.929 66.528 73.979 1.002  3900

For each parameter, n.eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor (at convergence, Rhat=1).

DIC info (using the rule, pD = var(deviance)/2)
pD = 6.7 and DIC = 71.5

DIC is an estimate of expected predictive error (lower deviance is better).

It seems the difference is is just not large enough to be convincing, which is what was supposed to be found. Css has now decreased compared to the first model, closer to the value which was given.

Previous posts in this series:

Simple Pharmacokinetics with Jags

PK calculation of IV and oral dosing in JAGS Translated in Stan by Bob Carpenter here: Stan Model of the Week: PK Calculation of IV and Oral Dosing