Exercise 2.2: Model Calibration
Nan-Hung Hsieh
2019/12/09 (update: 2019-12-08)
1 TASK
After the uncertainty analysis, we can calibrate the model parameters by Markov chain Monte Carlo technique
Use the parameter distributions that we test in uncertainty analysis and conduct MCMC simulation to do model
2 SOLUTION
2.1 Packages
library(simuloR)
library(tidyverse)
library(httk)
library(corrplot)
library(bayesplot)
library(rstan)
theme_set(theme_bw())
color_scheme_set("mix-blue-red")2.2 Single chain modeling
2.2.1 Input simulation
MCMC ("MCMC.default.out","", # name of output file
"", # name of data file
2000,0, # iterations, print predictions flag
1,2000, # printing frequency, iters to print
10101010); # random seed (default)
Level {
# Prior
Distrib (Fgutabs, Uniform, 0.8, 1);
Distrib (kgutabs, LogUniform, 0.218, 21.8);
Distrib (kelim, Uniform, 0.05307864, 1.326966);
Distrib (Vdist, Uniform, 0.15603994, 3.900998);
# Likelihood
Likelihood(Conc, LogNormal, Prediction(Conc) , 1.1);
Simulation { # 1
# Constant
IngDose = 319.992; # ingested dose (mg)
BW = 79.6; # body weight (kg)
# Data
Print (Conc, 0.01 0.25 0.57 1.12 2.02 3.82 5.10 7.03 9.05 12.12 24.37 );
Data (Conc, 0.74 2.84 6.57 10.50 9.66 8.58 8.36 7.47 6.89 5.94 3.28 );
}
}
End.set.seed(1)
ffpk.mcmc.1 <- mcsim(model = "models/ffpk.model", input = "inputs/ffpk.mcmc.in")## * Creating executable program, pleas wait...## * Created executable program 'models/mcsim.ffpk.model'.## Executing...## * Create 'MCMC.check.out' from the last iteration.2.2.2 Posterior checks
par(mfrow = c(2,2))
plot(ffpk.mcmc.1$Fgutabs.1., type = "l")
plot(ffpk.mcmc.1$kgutabs.1., type = "l")
plot(ffpk.mcmc.1$kelim.1., type = "l")
plot(ffpk.mcmc.1$Vdist.1., type = "l")
parms <- httk::parameterize_1comp(chem.name = "theophylline")## Warning in available_rblood2plasma(chem.cas = chem.cas, chem.name =
## chem.name, : Human in vivo Rblood2plasma returned.Fgutabs <- parms$Fgutabs * parms$hepatic.bioavailability
kgutabs <- parms$kgutabs
kelim <- parms$kelim
Vdist <- parms$Vdisti <- 1001:2000
ffpk.mcmc.1$Fgutabs.1.[i] %>% density() %>% plot(main = "Fgutabs")
plot.rng <- seq(par("usr")[1], par("usr")[2], length.out = 1000)
prior.dens <- do.call("dunif", c(list(x=plot.rng), c(0.8,1)))
lines(plot.rng, prior.dens, lty="dashed") # 0.9141961
abline(v=Fgutabs, col="red")
mean(ffpk.mcmc.1$Fgutabs.1.[i])## [1] 0.9119495sd(ffpk.mcmc.1$Fgutabs.1.[i])## [1] 0.05815058ffpk.mcmc.1$kgutabs.1.[i] %>% density() %>% plot(main = "kgutabs")
plot.rng <- seq(par("usr")[1], par("usr")[2], length.out = 1000)
prior.dens <- do.call("dunif", c(list(x=plot.rng), c(0.218,21.8)))
lines(plot.rng, prior.dens, lty="dashed")
abline(v=kgutabs, col="red") # ka = 2.18 
mean(ffpk.mcmc.1$kgutabs.1.[i])## [1] 4.219778sd(ffpk.mcmc.1$kgutabs.1.[i])## [1] 0.5011811ffpk.mcmc.1$kelim.1.[i] %>% density() %>% plot(main = "kelim")
plot.rng <- seq(par("usr")[1], par("usr")[2], length.out = 1000)
prior.dens <- do.call("dunif", c(list(x=plot.rng), c(0.05307864, 1.326966)))
lines(plot.rng, prior.dens, lty="dashed")
abline(v=kelim, col="red") # ke = 0.2653932
mean(ffpk.mcmc.1$kelim.1.[i])## [1] 0.05440021sd(ffpk.mcmc.1$kelim.1.[i])## [1] 0.001261038ffpk.mcmc.1$Vdist.1.[i] %>% density() %>% plot(main = "Vdist")
plot.rng <- seq(par("usr")[1], par("usr")[2], length.out = 1000)
prior.dens <- do.call("dunif", c(list(x=plot.rng), c(0.15603994, 3.900998)))
lines(plot.rng, prior.dens, lty="dashed")
abline(v=Vdist, col="red") # 0.78
mean(ffpk.mcmc.1$Vdist.1.[i])## [1] 0.3577756sd(ffpk.mcmc.1$Vdist.1.[i])## [1] 0.02641781names(ffpk.mcmc.1)## [1] "iter" "Fgutabs.1." "kgutabs.1." "kelim.1." "Vdist.1."
## [6] "LnPrior" "LnData" "LnPosterior"ffpk.mcmc.1[i,c("Fgutabs.1.","kgutabs.1.","kelim.1.","Vdist.1.")] %>% cor()## Fgutabs.1. kgutabs.1. kelim.1. Vdist.1.
## Fgutabs.1. 1.00000000 0.07109448 -0.1136462 0.8783650
## kgutabs.1. 0.07109448 1.00000000 -0.1465097 0.3240788
## kelim.1. -0.11364618 -0.14650965 1.0000000 -0.2469482
## Vdist.1. 0.87836496 0.32407883 -0.2469482 1.0000000ffpk.mcmc.1[i,c("Fgutabs.1.","kgutabs.1.","kelim.1.","Vdist.1.")] %>% cor() %>% corrplot(method = "number")
2.2.3 Posterior predective checks
Single draws
check.out <- read.delim("MCMC.check.out")
Observed <- check.out$Data
Expected <- check.out$Prediction
plot(x=Observed, y=Expected, xlim=c(0,10), ylim=c(0,10))
abline(0,1)
plot(Expected, Observed/Expected, pch='x')
abline(1,0)
Multiple draws
ffpk <- function(IngDose, Fgutabs, kgutabs, kelim, Vdist, BW, t){
A <- (IngDose * Fgutabs * kgutabs)/(Vdist * BW * (kgutabs - kelim))
Conc <- A * exp(-kelim * t) - A * exp(-kgutabs * t)
return(Conc)
}Fgutabs <- ffpk.mcmc.1$Fgutabs.1.
kgutabs <- ffpk.mcmc.1$kgutabs.1.
kelim <- ffpk.mcmc.1$kelim.1.
Vdist <- ffpk.mcmc.1$Vdist.1.
df <- data.frame(Fgutabs, kgutabs, kelim, Vdist)S1 <- subset(Theoph, Subject==1)
t1 <- S1$Time
Dose1 <- mean(S1$Dose)
BW1 <- mean(S1$Wt)Conc.matrix <- matrix(nrow = length(Fgutabs), ncol = length(t1))
for(i in 1:11){
Conc.matrix[,i] <- ffpk(IngDose = Dose1*BW1,
Fgutabs = df$Fgutabs, kgutabs = df$kgutabs,
kelim = df$kelim,
Vdist = df$Vdist,
BW = BW1,
t = t1[i])
}tail(Conc.matrix)## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1996,] 0 6.718464 9.042450 9.454967 9.060106 8.211251 7.655808
## [1997,] 0 6.380599 8.917774 9.539501 9.184477 8.325652 7.762471
## [1998,] 0 6.476178 9.053583 9.689727 9.337618 8.480128 7.916897
## [1999,] 0 6.420394 8.974875 9.603881 9.252128 8.397402 7.836290
## [2000,] 0 6.468327 9.041879 9.675581 9.321202 8.460095 7.894794
## [2001,] 0 6.408414 8.958128 9.585961 9.234864 8.381733 7.821668
## [,8] [,9] [,10] [,11]
## [1996,] 6.888519 6.167681 5.213932 2.667180
## [1997,] 6.984492 6.253611 5.286575 2.704340
## [1998,] 7.137573 6.403943 5.430754 2.813263
## [1999,] 7.060314 6.330321 5.362777 2.766621
## [2000,] 7.113025 6.377582 5.402814 2.787276
## [2001,] 7.047140 6.318509 5.352771 2.761459plot(t1, Observed, pch = "")
for(i in 1001:2000){
lines(t1, Conc.matrix[i,], col="grey")
}
points(t1, Observed)
2.2.4 Use setpoints analysis
SetPoints ("", "MCMC.default.out", 0, Fgutabs, kgutabs, kelim, Vdist);
Simulation {
# Constant
IngDose = 319.992; # ingested dose (mg)
BW = 79.6; # body weight (kg)
Print (Conc, 0.01 0.25 0.57 1.12 2.02 3.82 5.10 7.03 9.05 12.12 24.37 );
}
END.ffpk.setpts.out <- mcsim("models/ffpk.model", input = "inputs/ffpk.setpts.in")## Warning in readLines(input): incomplete final line found on 'inputs/
## ffpk.setpts.in'## * Creating executable program, pleas wait...## * Created executable program 'models/mcsim.ffpk.model'.## Executing...tail(ffpk.setpts.out)## Iter Fgutabs kgutabs kelim Vdist Conc_1.1 Conc_1.2 Conc_1.3
## 1996 1995 0.901786 4.52603 0.0547194 0.362596 0.442296 6.71846 9.04245
## 1997 1996 0.912991 4.03587 0.0547194 0.362596 0.400270 6.38060 8.91777
## 1998 1997 0.926529 4.03587 0.0536925 0.362596 0.406207 6.47618 9.05358
## 1999 1998 0.918593 4.03587 0.0540290 0.362596 0.402727 6.42039 8.97487
## 2000 1999 0.925451 4.03587 0.0540290 0.362596 0.405734 6.46833 9.04188
## 2001 2000 0.916879 4.03587 0.0540290 0.362596 0.401976 6.40841 8.95813
## Conc_1.4 Conc_1.5 Conc_1.6 Conc_1.7 Conc_1.8 Conc_1.9 Conc_1.10
## 1996 9.45497 9.06011 8.21125 7.65581 6.88852 6.16768 5.21393
## 1997 9.53950 9.18448 8.32565 7.76247 6.98449 6.25361 5.28657
## 1998 9.68973 9.33762 8.48013 7.91690 7.13757 6.40394 5.43075
## 1999 9.60388 9.25213 8.39740 7.83629 7.06031 6.33032 5.36278
## 2000 9.67558 9.32120 8.46010 7.89479 7.11302 6.37758 5.40281
## 2001 9.58596 9.23486 8.38173 7.82167 7.04714 6.31851 5.35277
## Conc_1.11
## 1996 2.66718
## 1997 2.70434
## 1998 2.81326
## 1999 2.76662
## 2000 2.78728
## 2001 2.76146str.pt <- which(names(ffpk.setpts.out)=="Conc_1.1")
end.pt <- which(names(ffpk.setpts.out)=="Conc_1.11")
plot(t1, Observed, pch = "")
for(i in 501:1000){
lines(t1, ffpk.setpts.out[i,str.pt:end.pt], col="grey")
}
points(t1, Observed)
2.3 Multi-chains simulation
2.3.1 Input simulation
set.seed(2)
ffpk.mcmc.2 <- mcsim(model = "models/ffpk.model", input = "inputs/ffpk.mcmc.in")## * Creating executable program, pleas wait...## * Created executable program 'models/mcsim.ffpk.model'.## Executing...## * Create 'MCMC.check.out' from the last iteration.set.seed(3)
ffpk.mcmc.3 <- mcsim(model = "models/ffpk.model", input = "inputs/ffpk.mcmc.in")## * Creating executable program, pleas wait...## * Created executable program 'models/mcsim.ffpk.model'.## Executing...## * Create 'MCMC.check.out' from the last iteration.set.seed(4)
ffpk.mcmc.4 <- mcsim(model = "models/ffpk.model", input = "inputs/ffpk.mcmc.in")## * Creating executable program, pleas wait...## * Created executable program 'models/mcsim.ffpk.model'.## Executing...## * Create 'MCMC.check.out' from the last iteration.2.3.2 Posterior checks
Use mcmc_array to create the 3-D array (iterations * chains * [parameters+4]) to store all outputs.
sims <- mcmc_array(data = list(ffpk.mcmc.1, ffpk.mcmc.2, ffpk.mcmc.3, ffpk.mcmc.4))
dim(sims)## [1] 2001 4 8The bayesplot package can use to visualize output.
names(ffpk.mcmc.1)## [1] "iter" "Fgutabs.1." "kgutabs.1." "kelim.1." "Vdist.1."
## [6] "LnPrior" "LnData" "LnPosterior"parms_name <- c("Fgutabs.1.", "kgutabs.1.", "kelim.1.", "Vdist.1.")
mcmc_trace(sims, pars = parms_name, facet_args = list(ncol = 1, strip.position = "left"))
mcmc_dens_overlay(x = sims[1001:2000,,], pars = parms_name)
mcmc_pairs(sims[1001:2000,,], pars = parms_name, off_diag_fun = "hex")
monitor(sims[,,parms_name], digit=4) ## Inference for the input samples (4 chains: each with iter = 2001; warmup = 1000):
##
## Q5 Q50 Q95 Mean SD Rhat Bulk_ESS Tail_ESS
## Fgutabs.1. 0.819 0.920 0.994 0.915 0.055 1.047 82 237
## kgutabs.1. 3.445 4.225 5.055 4.248 0.484 1.007 472 682
## kelim.1. 0.053 0.054 0.057 0.054 0.001 1.019 458 568
## Vdist.1. 0.318 0.362 0.397 0.360 0.024 1.049 68 157
##
## For each parameter, Bulk_ESS and Tail_ESS are crude measures of
## effective sample size for bulk and tail quantities respectively (an ESS > 100
## per chain is considered good), and Rhat is the potential scale reduction
## factor on rank normalized split chains (at convergence, Rhat <= 1.05).