Cefamandole: PK after IV with six subjects in JAGS
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
After last week’s post on Theoph I noticed the MEMSS library contained more PK data sets. The Cefamandole data is an IV data set with six subjects. It is somewhat challenging, since there seem to be several elimination rates.
Data
Data are displayed below. Notice that the lines seem somewhat curved, there is clearly more going on that just simple elimination.
Models
The first model fits concentration as function of two eliminations. The problem here is that the eliminations need to be ordered. All faster eliminations should have one prior, all slower another. After much fiddling around I found that I cannot force the ordering at subject level with priors or limits. Hence I chose a different way out. The model gets an extra component, it will predict 0 if the parameters are incorrectly ordered. This indeed gave the desired result in terms of model behavior.
The model itself though, was less to my linking. It goes way too low at the later times. Upon consideration, since the model assumes homoscedastic error it does not matter to the likelihoood if the fit at the lower end is not so good. After all, an error of 1 is nothing compared to errors of 10 or 20 at the upper end of the scale. In addition it can be argued that the error is not homoscedastic. If it were homoscedastic, the underlying data would show way more variation at the lower end of the scale. It would need confirmation from the measuring team, but for now I will assume proportional error.
Model 2: Proportional error
The proportional error will be implemented by using a log normal distribution for concentration. This was pretty simple to code. The log of the expected concentration is used in combination with dlnorm(). The only issue remaining is that expectation 0 for incorrectly ordered parameters has to be shifted a bit, log(0) is a bit of an issue.
This model looks pretty decent. There are still some things to do, similar to last week, translate the parameters c1star and c2star into the parameters which can be interpreted from PK perspective.
Model outputs
Model 1
Inference for Bugs model at “C:/Users/Kees/AppData/Local/Temp/RtmpSe5nAn/modele207c187661.txt”, fit using jags,
4 chains, each with 10000 iterations (first 5000 discarded), n.thin = 10
n.sims = 2000 iterations saved
mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat n.eff
C[1] 91.473 24.508 51.690 76.444 88.708 103.600 144.987 1.016 160
C[2] 324.947 132.978 147.584 253.433 304.913 366.175 626.711 1.007 460
Ke[1] 0.020 0.002 0.016 0.018 0.019 0.021 0.024 1.045 76
Ke[2] 0.141 0.192 0.081 0.116 0.132 0.147 0.202 1.072 180
c1star[1] 58.783 8.381 43.892 53.222 58.514 63.406 76.601 1.047 64
c1star[2] 59.841 18.962 32.887 47.487 56.324 67.606 108.087 1.048 72
c1star[3] 92.337 10.569 70.647 85.720 92.805 99.469 112.338 1.041 74
c1star[4] 89.554 12.223 67.975 81.086 88.831 96.932 116.699 1.027 110
c1star[5] 146.419 21.504 106.395 132.153 145.199 159.149 193.558 1.032 84
c1star[6] 121.228 14.872 87.712 112.706 122.379 130.819 148.106 1.092 42
c2star[1] 426.241 107.777 278.572 355.121 406.845 471.583 703.121 1.038 92
c2star[2] 201.118 35.861 149.478 179.179 197.081 217.084 278.248 1.033 830
c2star[3] 489.294 218.162 268.493 349.208 419.178 537.325 1125.062 1.088 42
c2star[4] 449.402 68.164 340.627 402.410 440.666 487.652 603.535 1.038 75
c2star[5] 402.518 43.391 335.760 373.724 397.703 423.902 498.253 1.017 380
c2star[6] 140.893 46.737 79.548 114.343 133.240 157.967 237.080 1.024 650
ke1[1] 0.020 0.002 0.016 0.018 0.020 0.021 0.025 1.042 74
ke1[2] 0.021 0.004 0.016 0.019 0.020 0.023 0.032 1.060 57
ke1[3] 0.018 0.002 0.014 0.017 0.018 0.019 0.022 1.031 99
ke1[4] 0.020 0.002 0.017 0.019 0.020 0.021 0.025 1.023 120
ke1[5] 0.020 0.002 0.016 0.019 0.020 0.021 0.024 1.036 74
ke1[6] 0.019 0.002 0.015 0.018 0.019 0.020 0.022 1.071 46
ke2[1] 0.167 0.025 0.128 0.150 0.164 0.181 0.225 1.050 64
ke2[2] 0.101 0.026 0.067 0.085 0.096 0.111 0.164 1.050 91
ke2[3] 0.181 0.042 0.122 0.151 0.173 0.200 0.283 1.076 43
ke2[4] 0.143 0.019 0.110 0.130 0.141 0.155 0.181 1.040 74
ke2[5] 0.112 0.017 0.086 0.100 0.110 0.121 0.151 1.027 130
ke2[6] 0.118 0.038 0.062 0.093 0.114 0.137 0.205 1.051 61
deviance 464.708 7.996 451.264 458.919 464.029 469.608 482.206 1.006 550
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 = 31.8 and DIC = 496.6
DIC is an estimate of expected predictive error (lower deviance is better).
Model 2
Code used
library(MEMSS)
library(ggplot2)
library(R2jags)
png(‘Cefamandole.png’)
qplot(y=conc,x=Time,col=Subject,data=Cefamandole) +
scale_y_log10(‘cefamandole (mcg/ml)’)+
geom_line()+
theme(legend.position=’bottom’)
dev.off()
library(R2jags)
datain <- list(
time=Cefamandole$Time,
conc=Cefamandole$conc,
n=nrow(Cefamandole),
subject =c(1:nlevels(Cefamandole$Subject))[Cefamandole$Subject],
nsubject = nlevels(Cefamandole$Subject)
)
model1 <- function() {
tau <- 1/pow(sigma,2)
sigma ~ dunif(0,100)
for (i in 1:n) {
pred[i] <- (ke1[subject[i]]
(c1star[subject[i]]*exp(-ke1[subject[i]]*time[i])+
c2star[subject[i]]*exp(-ke2[subject[i]]*time[i]))
conc[i] ~ dnorm(pred[i],tau)
}
for(i in 1:nsubject) {
ke1[i] ~ dlnorm(ke[1],kemtau[1])
ke2[i] ~ dlnorm(ke[2],kemtau[2])
c1star[i] ~ dlnorm(cm[1],ctau[1])
c2star[i] ~ dlnorm(cm[2],ctau[2])
}
for (i in 1:2) {
kem[i] ~ dunif(-10,10)
kemtau[i] <- 1/pow(kesd[i],2)
kesd[i] ~ dunif(0,10)
Ke[i] <- exp(ke[i])
cm[i] ~ dunif(-10,20)
ctau[i] <- 1/pow(csd[i],2)
csd[i] ~ dunif(0,100)
C[i] <- exp(cm[i])
}
ke <- sort(kem)
}
parameters <- c('Ke','ke1','ke2','c1star','c2star','C')
inits <- function() {
list(ke1 = rnorm(6,0.13,.03),
ke2= rnorm(6,0.0180,.005),
kem= log(c(rnorm(1,0.13,.03),rnorm(1,0.0180,.005))),
kesd =runif(2,.1,.4),
cm = log(rnorm(2,25,5)),
c1star=rnorm(6,25,3),
c2star=rnorm(6,25,3)
)}
jagsfit1 <- jags(datain, model=model1,
inits=inits,
parameters=parameters,progress.bar=”gui”,
n.iter=10000,
n.chains=4,n.thin=10)
jagsfit1
#plot(jagsfit1)
#traceplot(jagsfit1)
Time <- seq(0,max(Cefamandole$Time))
la1 <- sapply(1:nrow(jagsfit1$BUGSoutput$sims.matrix),
function(i) {
C1 <- jagsfit1$BUGSoutput$sims.matrix[i,'C[1]']
C2 <- jagsfit1$BUGSoutput$sims.matrix[i,'C[2]']
k1 <- jagsfit1$BUGSoutput$sims.matrix[i,'Ke[1]']
k2 <- jagsfit1$BUGSoutput$sims.matrix[i,'Ke[2]']
y=C1*exp(-k1*Time)+C2*exp(-k2*Time)
})
res1 <- data.frame(
Time=Time,
Conc=apply(la1,1,mean),
C05=apply(la1,1,FUN=function(x) quantile(x,0.05)),
C95=apply(la1,1,FUN=function(x) quantile(x,0.95)))
png(‘cef1.png’)
ggplot(res1, aes(x=Time)) +
geom_line(aes(y=Conc))+
scale_y_log10(‘cefamandole (mcg/ml)’)+
geom_line(aes(y=conc,x=Time,col=Subject),alpha=1,data=Cefamandole)+
geom_ribbon(aes(ymin=C05, ymax=C95),alpha=.2)+
theme(legend.position=’bottom’)
dev.off()
#########
#
# model 2: proportional error
#
#########
model2 <- function() {
tau <- 1/pow(sigma,2)
sigma ~ dunif(0,100)
for (i in 1:n) {
pred[i] <- log(
(ke1[subject[i]]
(c1star[subject[i]]*exp(-ke1[subject[i]]*time[i])+
c2star[subject[i]]*exp(-ke2[subject[i]]*time[i]))
+.001*(ke1[subject[i]]>ke2[subject[i]]))
conc[i] ~ dlnorm(pred[i],tau)
}
for(i in 1:nsubject) {
ke1[i] ~ dlnorm(ke[1],kemtau[1])
ke2[i] ~ dlnorm(ke[2],kemtau[2])
c1star[i] ~ dlnorm(cm[1],ctau[1])
c2star[i] ~ dlnorm(cm[2],ctau[2])
}
for (i in 1:2) {
kem[i] ~ dunif(-10,10)
kemtau[i] <- 1/pow(kesd[i],2)
kesd[i] ~ dunif(0,10)
Ke[i] <- exp(ke[i])
cm[i] ~ dunif(-10,20)
ctau[i] <- 1/pow(csd[i],2)
csd[i] ~ dunif(0,100)
C[i] <- exp(cm[i])
}
ke <- sort(kem)
}
parameters <- c('Ke','ke1','ke2','c1star','c2star','C')
inits <- function() {
list(ke1 = rnorm(6,0.13,.03),
ke2= rnorm(6,0.0180,.005),
kem= log(c(rnorm(1,0.13,.03),rnorm(1,0.0180,.005))),
kesd =runif(2,.1,.4),
cm = log(rnorm(2,25,5)),
c1star=rnorm(6,25,3),
c2star=rnorm(6,25,3)
)}
jagsfit2 <- jags(datain, model=model2,
inits=inits,
parameters=parameters,progress.bar=”gui”,
n.iter=10000,
n.chains=4,n.thin=10)
jagsfit2
la2 <- sapply(1:nrow(jagsfit2$BUGSoutput$sims.matrix),
function(i) {
C1 <- jagsfit2$BUGSoutput$sims.matrix[i,'C[1]']
C2 <- jagsfit2$BUGSoutput$sims.matrix[i,'C[2]']
k1 <- jagsfit2$BUGSoutput$sims.matrix[i,'Ke[1]']
k2 <- jagsfit2$BUGSoutput$sims.matrix[i,'Ke[2]']
y=C1*exp(-k1*Time)+C2*exp(-k2*Time)
})
res2 <- data.frame(
Time=Time,
Conc=apply(la2,1,mean),
C05=apply(la2,1,FUN=function(x) quantile(x,0.05)),
C95=apply(la2,1,FUN=function(x) quantile(x,0.95)))
png(‘cef2.png’)
ggplot(res2, aes(x=Time)) +
geom_line(aes(y=Conc))+
scale_y_log10(‘cefamandole (mcg/ml)’)+
geom_line(aes(y=conc,x=Time,col=Subject),alpha=1,data=Cefamandole)+
geom_ribbon(aes(ymin=C05, ymax=C95),alpha=.2)+
theme(legend.position=’bottom’)
dev.off()
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.