In this post I want to analyze a first order pharmocokinetcs problem: the data of study problem 9, chapter 3 of Rowland and Tozer (Clinical pharmacokinetics and pharmacodynamics, 4th edition) with Jags. It is a surprising simple set of data, but still there is enough to play around with.

Data, model and analysis

The data is simple enough. One subject’s concentration of cocaine (μg/L) in plasma as function of time, after a starting dose of 33 mg cocaine hydrochlorine. The model is simple enough, C=C0*exp(-k*t). A number of additional parameters are to be calculated, half-live (t1/2), Clearance (CL) and distribution volume (V). Additional information needed; molecular weight of cocaine 303 g/mol, cocaine hydrochloride 330 g/mol.
ch3sp9 <- data.frame(
    time=c(0.16,0.5,1,1.5,2,2.5,3),
    conc=c(170,122,74,45,28,17,10))

Bayesian model

This model has as feature that the (non-informative) priors are directly on the parameters which are in PK analysis important; distribution volume and clearance.
library(R2jags)
datain <- as.list(ch3sp9)
datain$n <- nrow(ch3sp9)
datain$dose <- 33*1000*303/340 # in ug cocaine
model1 <- function() {
  for (i in 1:n) {
    pred[i] <- c0*(exp(-k*time[i]))
    conc[i] ~ dnorm(pred[i],tau)
  }
  tau <- 1/pow(sigma,2)
  sigma ~ dunif(0,100)
  c0 <- dose/V
  V ~ dlnorm(2,.01)
  k <- CL/V
  CL ~ dlnorm(1,.01)
  AUC <- dose/CL
  t0_5 <- log(2)/k
}

parameters <- c('k','AUC','CL','V','t0_5')
inits <- function() 
  list(V=rlnorm(1,1,.1),CL=rlnorm(1,1,.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 “…”, 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
AUC      201.756   0.766 200.258 201.369 201.749 202.143 203.233 1.001 10000
CL       145.767   0.553 144.705 145.485 145.769 146.045 146.855 1.001 10000
V        147.568   0.496 146.607 147.320 147.569 147.819 148.515 1.002 18000
k          0.988   0.005   0.978   0.985   0.988   0.991   0.998 1.001 13000
t0_5       0.702   0.004   0.694   0.700   0.702   0.704   0.709 1.001 13000
deviance   6.769   4.166   2.236   3.801   5.682   8.557  17.538 1.001  6300

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 = 8.7 and DIC = 15.4

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

Alternative Baysian model

Looking at residuals

Using a package