### 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).

_{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).

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

Dear Kees Duineveld,

ReplyDeleteI just found your Blog on R2Jags approaches to PK analysis. What I'm trying to do is individual parameter estimation of a subject with 2 measurements and use an established PopPK model. The Bayesian approach seems suitable to use prior information of the published parameters to obtain individual posterior estimates, am I right? Unfortunately, the model is a bit complicated, i.e. a three-compartment model of a infusion. I tried to adapt your Wiekvoet code, but received an error message:

Error in jags.model(model.file, data = data, inits = init.values, n.chains = n.chains, :

Error in node V2

Invalid parent values

Can you image what is wrong?

Many thanks and best wishes

Andreas Meid

R code:

library(R2jags)

individual <- data.frame(

time=c(0,216,240),

conc=c(0,0.00473, 0.00409))

datain <- as.list(individual)

datain$n <- nrow(ch10sp4)

#datain$Rinf <- 1*100

model1 <- function() {

# concentration formula

tau <- 1/pow(sigma,2)

sigma ~ dunif(0,100)

# infusion time = 0

for (i in 1:n) {

pred[i] <- ifelse(time[i] <= 0,

(A/alpha * (1 - exp(-alpha*time[i]))) + (B/beta * (1 - exp(-beta*time[i]))) + (C/gamma * (1 - exp(-gamma*time[i]))),

(A/alpha * (1 - exp(-alpha*0)) * exp(-alpha*(time[i]-0))) +

(B/beta * (1 - exp(-beta*0)) * exp(-beta*(time[i]-0))) +

(C/gamma * (1 - exp(-gamma*0)) * exp(-gamma*(time[i]-0))))

conc[i] ~ dnorm(pred[i],tau)

}

# micro constants

V <- V1

k <- CL/V1

k12 <- Q2/V1

k21 <- Q2/V2

k13 <- Q3/V1

k31 <- Q3/V3

# macro constants

a0 <- k*k21*k31

a1 <- k*k31 + k21*k31 + k21*k13 + k*k21 + k31*k12

a2 <- k + k12 + k13 + k21 + k31

p <- a1 - (a2^2)/3

q <- ((2* a2^3)/27 ) - (a1*(a2^3)) + a0

r1 <- sqrt(-((p^3)/27))

r2 <- 2 * (r1^(1/3))

phi <- (acos(-(q)/(2*r1))) / 3

alpha <- -( (cos(phi) * r2) - a2/3 )

beta <- -( (cos(phi + (2*3.141593/3) ) * r2) - a2/3 )

gamma <- -( (cos(phi + (4*3.141593/3) ) * r2) - a2/3 )

A <- 1/V * ((k21-alpha)/(alpha-beta)) * ((k31-alpha)/(alpha-gamma))

B <- 1/V * ((k21-beta)/(beta-alpha)) * ((k31-beta)/(beta-gamma))

C <- 1/V * ((k21-gamma)/(gamma-beta)) * ((k31-gamma)/(gamma-alpha))

# prior information from published model

CL~dlnorm(log(44.1),log(44.1*0.06))

V1~dlnorm(log(8.9),log(8.9*0.09))

Q2~dlnorm(log(6.1),log(6.1*0.07))

V2~dlnorm(log(7.3),log(7.3*0.09))

Q3~dlnorm(log(14.4),log(14.4*0.12))

V3~dlnorm(log(388),log(388*0.11))

}

parameters <- c('CL','V1','Q2','V2','Q3','V3')

inits <- function()

list(CL=rlnorm(1,log(44),.1),V1=rlnorm(1,log(9),.1),

Q2=rlnorm(1,log(6),.1), V2=rlnorm(1,log(7),.1),

Q3=rlnorm(1,log(14),.1), V3=rlnorm(1,log(388),.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)