[This article was first published on Wiekvoet, 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.
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.Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
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).
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).
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 JagsPK 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
To leave a comment for the author, please follow the link and comment on their blog: Wiekvoet.
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.