Parameters
Variance - Covariance (Least-square estimate)
Variance - Covariance (prior distribution)
Variance - Covariance (posterior distribution)
ui.R
Download
#-------------------------------------------------------------------------
# This application is governed by the CeCILL-B license.
# You can use, modify and/ or redistribute this code under the terms
# of the CeCILL license: http://www.cecill.info/index.en.html
#
# Marc Lavielle, Inria Saclay
# June 28th, 2015
#-------------------------------------------------------------------------
library(shinydashboard)
sidebar <- dashboardSidebar(
hr(),
sidebarMenu(id="tabs",
menuItem("Plots", icon = icon("line-chart"),
menuSubItem("data and predictions", tabName = "plot1", icon = icon("angle-right")), selected=TRUE,
menuSubItem("parameters distributions", tabName = "plot2", icon = icon("angle-right")),
menuSubItem("MCMC", tabName = "plot3", icon = icon("angle-right"))
),
menuItem("Parameters", tabName = "tabparam", icon = icon("table")),
menuItem("Codes", icon = icon("file-text-o"),
menuSubItem("ui.R", tabName = "ui", icon = icon("angle-right")),
menuSubItem("server.R", tabName = "server", icon = icon("angle-right"))
),
menuItem("ReadMe", tabName = "readme", icon = icon("mortar-board")),
menuItem("About", tabName = "about", icon = icon("question"))
),
hr(),
conditionalPanel("input.tabs=='plot1'",
sidebarMenu(
menuItem("Display", icon = icon("chevron-circle-right"),
fluidRow(
column(1),
column(10, checkboxInput("disp_data", "data", TRUE))
),
fluidRow(
column(1),
column(6, checkboxInput("disp_cond", "LS estimate", TRUE)),
column(4, checkboxInput("disp_pint1", "C.I. ", FALSE))
),
fluidRow(
column(1),
column(6, checkboxInput("disp_prior", "Prior ", TRUE)),
column(4, checkboxInput("disp_pint0", "P.I ", FALSE))
),
fluidRow(
column(1),
column(6, checkboxInput("disp_pmedian", "MAP", TRUE)),
column(4, checkboxInput("disp_pint2", "P.I ", FALSE))
),
fluidRow(
column(1),
# column(5, checkboxInput("disp_pint", "Intervals ", FALSE)),
column(5, numericInput("predint", "", min=1, max=99, value=90))
)
)
)
),
conditionalPanel("input.tabs=='plot2'",
sidebarMenu(
menuItem("Display", icon = icon("chevron-circle-right"),
fluidRow(
column(1),
column(10, checkboxInput("dens_prior", "prior densities", TRUE))
),
fluidRow(
column(1),
column(10, checkboxInput("dens_post", "posterior densities", TRUE))
)
)
)
),
conditionalPanel("input.tabs=='plot1' | input.tabs=='plot2' | input.tabs=='plot3'",
sidebarMenu(
menuItem("Distributions", icon = icon("chevron-circle-right"),
fluidRow(
column(1),
column(10, radioButtons("distribution","",c("prior" = "1","observations" = "2"), selected="2"))
)
),
menuItem("MCMC", icon = icon("chevron-circle-right"),
fluidRow(
column(1),
column(8, numericInput("niter", "Total iterations", 1000))
),
fluidRow(
column(1),
column(8, numericInput("nburn", "Burning iterations", 200))
),
fluidRow(
column(1),
column(10, actionButton("run", "Run"))
)
)
)
),
hr()
)
body <- dashboardBody(
tabItems(
tabItem(tabName = "plot1",
fluidRow(
box(width = 12, status = "primary",
plotOutput("plot1", height="400px",
click = "plot1_click",
brush = brushOpts(id = "plot1_brush")
)
)
)
),
tabItem(tabName = "plot2",
fluidRow(
box(width = 12, status = "primary",
plotOutput("plot2", height="400px")
)
)
),
tabItem(tabName = "plot3",
fluidRow(
box(width = 12, status = "primary",
plotOutput("plot3", height="400px")
)
)
),
tabItem(tabName = "tabparam",
box( width = NULL, status = "primary", solidHeader = TRUE, title="Parameters",
tableOutput("tablep")
),
box( width = NULL, status = "primary", solidHeader = TRUE, title="Variance - Covariance (Least-square estimate)",
tableOutput("tableC1")
),
box( width = NULL, status = "primary", solidHeader = TRUE, title="Variance - Covariance (prior distribution)",
tableOutput("tableC0")
),
box( width = NULL, status = "primary", solidHeader = TRUE, title="Variance - Covariance (posterior distribution)",
tableOutput("tableC2"),digits=5
)
),
tabItem(tabName = "ui",
box( width = NULL, status = "primary", solidHeader = TRUE, title="ui.R",
downloadButton('downloadData2', 'Download'),
br(),br(),
pre(includeText("ui.R"))
)
),
tabItem(tabName = "server",
box( width = NULL, status = "primary", solidHeader = TRUE, title="server.R",
downloadButton('downloadData3', 'Download'),
br(),br(),
pre(includeText("server.R"))
)
),
tabItem(tabName = "readme",
withMathJax(),
includeMarkdown("readMe.Rmd")
),
tabItem(tabName = "about",
includeMarkdown("about.Rmd")
)
),
conditionalPanel(condition = "input.distribution=='1' & (input.tabs=='plot1'|input.tabs=='plot2'|input.tabs=='plot3') ",
box(width = 4, status = "primary", solidHeader = FALSE,
sliderInput("ka", "ka*:", value = 1, min = 0, max = 2, step=0.1)
),
box(width = 4, status = "primary", solidHeader = FALSE,
sliderInput("V", "V*:", value = 1, min = 0, max = 2, step=0.1)
),
box(width = 4, status = "primary", solidHeader = FALSE,
sliderInput("Cl", "Cl*:", value = 0.05, min = 0, max = 0.1, step=0.005)
),
box(width = 4, status = "primary", solidHeader = FALSE,
sliderInput("oka", "omega_ka:", value = 0.2, min = 0, max = 1, step=0.05)
),
box(width = 4, status = "primary", solidHeader = FALSE,
sliderInput("oV", "omega_V:", value = 0.2, min = 0, max = 1, step=0.05)
),
box(width = 4, status = "primary", solidHeader = FALSE,
sliderInput("oCl", "omega_Cl:", value = 0.2, min = 0, max = 1, step=0.05)
)
),
conditionalPanel(condition = "input.distribution=='2' & (input.tabs=='plot1'|input.tabs=='plot2'|input.tabs=='plot3') ",
box(width = 4, status = "primary", solidHeader = FALSE,
sliderInput("sigma", "sigma:", value = 5, min = 0.5, max = 20, step=0.5)
),
box(width = 4, status = "primary", solidHeader = FALSE,
selectInput("N", label = h5("Number of data points"),choices = c(5,10,25,50,100),selected="25")
),
conditionalPanel(condition = "input.tabs=='plot1'",
box(width = 4, status = "primary", solidHeader = FALSE,
checkboxInput("disp_exclude", "display unused data", FALSE),
radioButtons("add","",c("add data" = "1","exclude data" = "2")),
actionButton("exclude_reset", "Reset")
# actionButton("run", "Run")
))
# box(width = 4, status = "primary", solidHeader = FALSE,
# actionButton("exclude_reset", "Reset")
# )
)
)
dashboardPage(
dashboardHeader(title = "Bayesian fitting"),
sidebar,
body
)
server.R
Download
#-------------------------------------------------------------------------
# This application is governed by the CeCILL-B license.
# You can use, modify and/ or redistribute this code under the terms
# of the CeCILL license: http://www.cecill.info/index.en.html
#
# Marc Lavielle, Inria Saclay
# June 28th, 2015
#-------------------------------------------------------------------------
library(ggplot2)
library(gridExtra)
# source("mlxmcmc1.R")
# source("model.R")
ggplotmlx <- function(...) {
ggplot(...) + theme_bw() + theme(plot.background = element_rect(fill=rgb(1,1,1)))
}
shinyServer(function(input, output) {
d <- read.csv(file="pkdata.txt", header=TRUE, sep="\t")
d.time <- d$time
n <- length(d$time)
c.time <- seq(0,100,by=0.25)
pk.prior <- list( name = c('ka','V','Cl'),
reference = c(1, 0.1, 0.01),
sd = c(0.2, 0.2, 0.2)*0.1)
vals <- reactiveValues(
keeprows = rep(TRUE, nrow(d))
)
r0 <- reactive({
muphi <- log(c(input$ka,input$V,input$Cl))
chol.omega <- c(input$oka,input$oV,input$oCl)
n.simul <- 1000
PSI = matrix(0,nrow=n.simul,ncol=length(muphi))
for (k in (1:n.simul)){
phi <- muphi + chol.omega*rnorm(3)
PSI[k,] <- exp(phi)
}
PSI <- data.frame(PSI)
names(PSI) <- c('ka','V','Cl')
r=list(P=PSI, p=exp(muphi))
return(r)
})
r1 <- reactive({
keep <- d[vals$keeprows, , drop = FALSE]
l0 <- log(c(0.1,1,0.04,2))
pk.nlm1 <- nlm(fmin1, l0, keep$y, keep$time, hessian="true")
muphi <- pk.nlm1$estimate[1:3]
solve.omega <- solve(pk.nlm1$hessian)*input$sigma^2
chol.omega <- chol(solve.omega[1:3,1:3])
n.simul <- 1000
PSI = matrix(0,nrow=n.simul,ncol=length(muphi))
for (k in (1:n.simul)){
phi <- muphi + chol.omega%*%rnorm(3)
PSI[k,] <- t(exp(phi))
}
PSI <- data.frame(PSI)
names(PSI) <- c('ka','V','Cl')
r=list(P=PSI, p=exp(muphi), C=solve.omega[1:3,1:3])
return(r)
})
r2 <- reactive({
input$run
keeprows <- vals$keeprows
niter <- isolate(input$niter)
pk.control <- list(niter=niter, nkernel=c(2,1))
keep <- d[keeprows, , drop = FALSE]
exclude <- d[!keeprows, , drop = FALSE]
pk.prior$reference <- (c(input$ka,input$V,input$Cl))
pk.prior$sd <- (c(input$oka,input$oV,input$oCl))
sigma <- (input$sigma)
if (length(which(pk.prior$sd>0))>0){
res <- mlxmcmc1(model="pk_1cpt",
data=keep,
xi=c(sigma,0),
prior=pk.prior,
control=pk.control)
pmap <- mlxoptim(model="pk_1cpt",
data=keep,
xi=c(sigma,0),
prior=pk.prior)
}else{
res <- matrix(rep(pk.prior$reference,niter),nrow=niter,byrow=TRUE)
pmap <- pk.prior$reference
}
p.table <- data.frame(res)
names(p.table) <- pk.prior$name
r=list(P=p.table, p=pmap)
return(r)
})
output$plot1 <- renderPlot({
pr <- 1 -input$predint/100
probs <- c(pr/2, 1-pr/2)
r <- r0()
S <- do.call("pkv_1cpt",list(r$P,c.time))
s <- apply(S,2,function(x) quantile(x,probs=probs))
bd0 <- data.frame(time=c.time,s1=s[1,],s2=s[2,])
s <- do.call("pk_1cpt",list(r$p,c.time))
C0 <- data.frame(time=c.time,C=s)
r <- r1()
S <- do.call("pkv_1cpt",list(r$P,c.time))
s <- apply(S,2,function(x) quantile(x,probs=probs))
bd1 <- data.frame(time=c.time,s1=s[1,],s2=s[2,])
s <- do.call("pk_1cpt",list(r$p,c.time))
C1 <- data.frame(time=c.time,C=s)
r <- r2()
r$P <- isolate(r$P[(input$nburn+1):input$niter,])
S <- do.call("pkv_1cpt",list(r$P,c.time))
s <- apply(S,2,function(x) quantile(x,probs=probs))
bd2 <- data.frame(time=c.time,s1=s[1,],s2=s[2,])
s <- do.call("pk_1cpt",list(r$p,c.time))
C2 <- data.frame(time=c.time,C=s)
keep <- d[vals$keeprows, , drop = FALSE]
exclude <- d[!vals$keeprows, , drop = FALSE]
vc=c(rgb(1,0.,0.),rgb(1,0.8,0.3),rgb(0.3,0.3,0.7),rgb(0.3,0.7,0.3),rgb(0.7,0.3,0.3))
names(vc)=letters[1:5]
lc <- c("used data", "unused data", "prior median",
"LS estimate", "posterior median")
vdisp <- rep(FALSE,5)
pl <- ggplotmlx()
if (input$disp_prior){
if (input$disp_pint0)
pl <- pl + geom_ribbon(data=bd0,aes(x=time,ymin=s1,ymax=s2), fill=rgb(0.8,0.8,0.9))
pl <- pl + geom_line(data=C0, aes(x=time, y=C, colour="c"))
vdisp[3] <- TRUE}
if (input$disp_cond){
if (input$disp_pint1)
pl <- pl + geom_ribbon(data=bd1,aes(x=time,ymin=s1,ymax=s2), fill=rgb(0.8,0.9,0.8))
pl <- pl + geom_line(data=C1, aes(x=time, y=C, colour="d"))
vdisp[4] <- TRUE}
if (input$disp_pmedian){
if (input$disp_pint2)
pl <- pl + geom_ribbon(data=bd2,aes(x=time,ymin=s1,ymax=s2), fill=rgb(0.9,0.8,0.8))
pl <- pl + geom_line(data=C2, aes(x=time, y=C, colour="e"))
vdisp[5] <- TRUE}
if (input$disp_data){
pl <- pl + geom_point(data=keep, aes(time, y, colour="a"), size=3)
vdisp[1] <- TRUE}
if (input$disp_exclude){
pl <- pl + geom_point(data=exclude, aes(time, y, colour="b") )
vdisp[2] <- TRUE}
if (any(vdisp)){
pl <- pl + scale_colour_manual(values=vc[vdisp], labels=lc[vdisp])
pl <- pl + theme(legend.position=c(0.8,0.85), legend.justification=c(0,1), legend.title=element_blank())
}
return(pl)
})
observeEvent(input$N, {
res = rep(FALSE, n)
L <- round(n/as.numeric(input$N))
i <- seq(floor(L/2),n,by=L)
res[i] <- TRUE
vals$keeprows <- res
})
observeEvent(input$plot1_click, {
res <- nearPoints(d, input$plot1_click, yvar="y", allRows=TRUE)
if (input$add=="1")
vals$keeprows[res$selected_] <- TRUE
else
vals$keeprows[res$selected_] <- FALSE
# vals$keeprows <- xor(vals$keeprows, res$selected_)
})
observeEvent(input$plot1_brush, {
res <- brushedPoints(d, input$plot1_brush, yvar="y", allRows=TRUE)
if (input$add=="1")
vals$keeprows[res$selected_] <- TRUE
else
vals$keeprows[res$selected_] <- FALSE
# vals$keeprows <- xor(vals$keeprows, res$selected_)
})
observeEvent(input$exclude_reset, {
res = rep(FALSE, n)
L <- round(n/as.numeric(input$N))
i <- seq(floor(L/2),n,by=L)
res[i] <- TRUE
vals$keeprows <- res
})
output$plot2 <- renderPlot({
pk.p <- log(c(input$ka,input$V,input$Cl))
pk.sd <- c(input$oka,input$oV,input$oCl)
pk.names <- c('ka','V','Cl')
p2 <- r2()$P
p2 <- isolate(p2[(input$nburn+1):input$niter,])
p0 <- r0()$P
p <- rbind(p0,p2)
x1 <- sapply(p,min)
x2 <- sapply(p,max)
n <- 200
pl <- list()
for (k in (1:3)){
plk <- ggplotmlx()
plk <- plk + ylab("") + xlab(pk.names[k])
if (input$dens_post){
dk <- density(p2[[k]],adjust=2)
dk <- data.frame(x=dk$x,y=dk$y)
plk <- plk + geom_line(data=dk,aes(x,y,colour="posterior"))
# plk <- plk + geom_line(aes(x,y),colour=rgb(0.7,0.3,0.3))
}
if (input$dens_prior){
xk <- seq(x1[k],x2[k],length.out=n)
dlk <- dlnorm(xk,meanlog=pk.p[k],sdlog=pk.sd[k])
dp <- data.frame(x=xk,y=dlk)
plk <- plk + geom_line(data=dp, aes(x=x,y=y,colour="prior"))
# plk <- plk + geom_line(data=dp, aes(x=x,y=y), colour=rgb(0.3,0.3,0.7))
}
if (k==1)
plk <- plk + theme(legend.position=c(0.9,0.6),legend.title=element_blank())
else
plk <- plk + theme(legend.position="none")
# plk <- plk + scale_colour_discrete(name ="distribution", labels=c("Prior", "Posterior"))
pl[[k]] <- plk
}
grid.arrange(pl[[1]],pl[[2]],pl[[3]])
})
output$plot3 <- renderPlot({
p <- r2()$P
pk.prior$sd <- c(input$oka,input$oV,input$oCl)
ind.eta = which(pk.prior$sd>0)
nb.etas <- length(ind.eta)
niter <- dim(p)[1]
p$iteration <- (1:niter)
# p1 <- qplot((1:niter),ka,data=p,colour="blue")
p1 <- ggplotmlx(data=p) + geom_line(aes(x=iteration,y=ka))
p2 <- ggplotmlx(data=p) + geom_line(aes(x=iteration,y=V))
p3 <- ggplotmlx(data=p) + geom_line(aes(x=iteration,y=Cl))
grid.arrange(p1,p2,p3)
})
output$tablep <- renderTable({
r <- matrix(ncol=3,nrow=3)
r[,1]=r1()$p
r[,2]=r0()$p
r[,3]=r2()$p
r <- data.frame(r)
names(r) <- c('LS','prior','MAP')
rownames(r) <- c('ka','V','Cl')
return(r)
},include.colnames=TRUE)
output$tableC1 <- renderTable({
r <- r1()$C
colnames(r) <- c('ka','V','Cl')
rownames(r) <- c('ka','V','Cl')
return(r)
},include.colnames=TRUE,digits=4)
output$tableC2 <- renderTable({
t2 <- r2()$P
r <- log(t2[(input$nburn+1):length(t2[[1]]),])
vr <- var(r)
return(vr)
},include.colnames=TRUE,digits=4)
output$tableC0 <- renderTable({
r <- matrix(ncol=3,nrow=3)
r[1,1]=input$oka^2
r[2,2]=input$oV^2
r[3,3]=input$oCl^2
r <- data.frame(r)
names(r) <- c('ka','V','Cl')
rownames(r) <- c('ka','V','Cl')
return(r)
},include.colnames=TRUE,digits=4)
})
mlxmcmc1 <- function(model=NULL,data=NULL,xi=NULL,prior=NULL,control=NULL)
{
if (is.null(control))
control=list()
if (is.null(control$niter))
control$niter = 100
if (is.null(control$nkernel))
control$nkernel = c(2,2)
if (is.null(control$stepsize.rw))
control$stepsize.rw=0.4
if (is.null(control$proba.mcmc))
control$proba.mcmc=0.4
if (is.null(control$nb.max))
control$nb.max=3
t.obs <- data$time
y.obs <- data$y
p0=exp(prior$reference)
y.pred <- do.call(model,list(p0,t.obs))
ind.eta = which(prior$sd>0)
nb.etas <- length(ind.eta)
if (nb.etas>1){
Omega.eta <- diag(prior$sd[ind.eta]^2)
chol.omega <- try(chol(Omega.eta))
solve.omega<- try(solve(Omega.eta))
domega <- diag(Omega.eta)*0.3;
}else{
Omega.eta <- prior$sd[ind.eta]^2
chol.omega <- prior$sd[ind.eta]
solve.omega<- 1/Omega.eta
domega <- Omega.eta*0.3;
}
# cat("Running Metropolis-Hastings algorithm\n")
# print(date())
mean.phi <- log(prior$reference)
psi <- prior$reference
phi <- log(psi)
f <- do.call(model,list(psi,t.obs))
g <- error.model(f,xi)
U.y <- sum(0.5 * ((y.obs - f)/g)^2 + log(g))
eta <- phi[ind.eta] - mean.phi[ind.eta]
phic <- phi
PSI = matrix(0,nrow=control$niter+1,ncol=length(psi))
PSI[1,]=psi
for (kiter in 1:control$niter)
{
#print(kiter)
if (control$nkernel[1] > 0)
{
for (u in 1:control$nkernel[1])
{
etac <- as.vector(chol.omega%*%rnorm(nb.etas))
phic[ind.eta] <- mean.phi[ind.eta] + etac
psic <- exp(phic)
f <- do.call(model,list(psic,t.obs))
g <- error.model(f,xi)
Uc.y <- sum(0.5 * ((y.obs - f)/g)^2 + log(g))
deltu <- Uc.y - U.y
if(deltu < (-1) * log(runif(1)))
{
eta <- etac
phi <- phic
U.y <- Uc.y
}
}
}
if (control$nkernel[2] > 0)
{
nb.max <- min(nb.etas,control$nb.max)
nbc2<-nt2 <- replicate(nb.etas,0)
U.eta <- 0.5 * eta %*% solve.omega %*% eta
for (u in 1:control$nkernel[2])
{
for (nrs2 in 1:nb.max)
{
for (j in 1:nb.etas)
{
jr <- sample(c(1:nb.etas), nrs2)
jr <- jr -jr[1] + j
vk2 <- jr%%nb.etas + 1
etac <- eta
etac[vk2] <- eta[vk2] + rnorm(nrs2)*domega[vk2]
phic[ind.eta] <- mean.phi[ind.eta] + etac
psic <- exp(phic)
f <- do.call(model,list(psic,t.obs))
g <- error.model(f,xi)
Uc.y <- sum(0.5 * ((y.obs - f)/g)^2 + log(g))
Uc.eta <- 0.5 * etac %*% solve.omega %*% etac
deltu <- Uc.y - U.y + Uc.eta - U.eta
if(deltu < (-1) * log(runif(1)))
{
eta <- etac
U.y <- Uc.y
U.eta <- Uc.eta
phi[ind.eta]=phic[ind.eta]
nbc2[vk2] <- nbc2[vk2]+1
}
nt2[vk2] <- nt2[vk2] + 1
}
}
}
domega <- domega*(1 + control$stepsize.rw*(nbc2/nt2 - control$proba.mcmc))
}
PSI[kiter+1,] <- exp(phi)
}
return(PSI)
}
#--------------------------------------------
mlxoptim <- function(model=NULL,data=NULL,xi=NULL,prior=NULL,initial=NULL)
{
errpred <- function(phi,model,y,t,muphi,ind.eta,xi,solve.omega){
psi <- exp(muphi)
psi[ind.eta] <- exp(phi)
eta <- phi-muphi[ind.eta]
f <- do.call(model,list(psi,t.obs))
g <- error.model(f,xi)
U.y <- sum(0.5 * ((y.obs - f)/g)^2 + log(g))
U.eta <- 0.5 * eta %*% solve.omega %*% eta
e <- U.y + U.eta
return(e)}
t.obs <- data[,1]
y.obs <- data[,2]
ind.eta = which(prior$sd>0)
nb.etas <- length(ind.eta)
if (nb.etas>1){
Omega.eta <- diag(prior$sd[ind.eta]^2)
solve.omega<- try(solve(Omega.eta))
}else{
Omega.eta <- prior$sd[ind.eta]^2
solve.omega<- 1/Omega.eta
}
l0 <-log(prior$reference[ind.eta])
muphi=log(prior$reference)
if (nb.etas>1){
r <- optim(l0,errpred,model=model,y=y.obs,t=t.obs,muphi=muphi,ind.eta=ind.eta,xi=xi,solve.omega=solve.omega,hessian=TRUE)
}else{
r <- optim(l0,errpred,model=model,y=y.obs,t=t.obs,muphi=muphi,ind.eta=ind.eta,xi=xi,solve.omega=solve.omega,
method="Brent",lower=-5, upper=5)
}
pest <- exp(r$par)
return(pest)
}
pk_1cpt <- function(x,t) {
ka <- x[1]
V<- x[2]
k <- x[3]/V
C <- 100*ka/(V*(ka-k))*(exp(-k*t)-exp(-ka*t))
return(C)
}
pkv_1cpt <- function(x,time) {
ka <- x[,1]
V<- x[,2]
k <- x[,3]/V
C <- 100*ka/(V*(k-ka))*(exp(-ka%*%t(time))-exp(-k%*%t(time)))
return(C)
}
error.model <- function(f,xi){
g <- xi[1] + xi[2]*f
return(g)
}
fmin1 <- function(x,y,t){
z=exp(x)
f=pk_1cpt(z[1:3],t)
g=z[4]
e=sum( ((y-f)/g)^2 + log(g^2))
return(e)
}
\[ \def\iid{{\sim}_{\rm i.i.d.}} \newcommand{\argmax}[1]{ {\rm arg}\max_{#1} } \newcommand{\argmin}[1]{ \operatorname*{arg\,min}_{#1} } \]
Objective
Fit a model to longitudinal data from a single individual using a Bayesian approach.
Let \(y=(y_1,y_2,\ldots,y_n)\) be the observed data at times \(t_1, t_2,\ldots,t_n\). We consider the following model for \(y_j\):
\[y_j = f(t_j,\phi) + \sigma\varepsilon_j\]
where
- \(f\) is the structural model,
- \(\phi\) is a vector of structural parameters. We assume that \(\phi\) is a random vector with prior distribution \(\pi\): \(\ \phi \sim \pi\).
- \((\varepsilon_j)\) is a sequence of independent and normally distributed residual errors with mean 0 and variance 1: \(\ \varepsilon_j \, \iid \, {\cal N}(0,1)\).
We aim to estimate the posterior distribution of \(\phi\) for a given vector of observations \(y\), a given residual error parameter \(\sigma\) and a given prior distribution \(\pi\). We also aim to estimate the conditional distribution of \(f(t,\phi)\):
The pharmacokinetics (PK) model
We will use in this example a one compartment PK model for oral administration with first order absorption and linear elimination. Assuming that a single dose \(D\) is given at time 0, the predicted concentration is thus given by:
\[f(t,\phi) = \frac{D \, ka}{V\,ka-Cl}\left(e^{-(Cl/V)\, t} - e^{-ka\, t} \right)\] where \(\phi=(ka,V,Cl)\) is a vector of PK parameters
We will use log-normal distributions as prior distributions:
\[ \begin{aligned} \log(ka) & \sim {\cal N}(\log(ka^{\star}), \omega_{ka}^2) \\ \log(V) & \sim {\cal N}(\log(V^{\star}), \omega_{V}^2) \\ \log(Cl) & \sim {\cal N}(\log(Cl^{\star}), \omega_{Cl}^2) \end{aligned} \]
Thus,the prior distribution of \(\phi\) depends on a vector of parameters \[ \zeta = (ka^{\star},V^{\star},Cl^{\star},\omega_{ka},\omega_{V},\omega_{Cl},\sigma)\] The joint model of \((y,\phi)\) is therefore a parametric model that can be described by the set of parameters \((\zeta,\sigma)\):
\[p(y,\phi ; \zeta,\sigma) = p(y|\phi;\sigma)p(\phi ; \zeta)\]
The posterior distribution of \(\phi\) combines the conditional distribution of the observations \(y\) and the prior distribution of \(\phi\):
\[p(\phi |y ; \zeta,\sigma) = \frac{p(y|\phi;\sigma)p(\phi ; \zeta)}{p(y;\zeta,\sigma)}\]
For any \(t\), \(f(t,\phi)\) is a random variable. We can then define the prior and posterior distributions of \(f(t,\phi)\).
Estimation of the prior and posterior distributions
Since \(f\) is non linear, the prior distribution of \(f(t,\phi)\) cannot be computed in a closed form. It can be estimated by Monte-Carlo simulations: \(N\) vectors of PK parameters \(\phi_1, \phi_2, \ldots, \phi_N\) are drawn with the prior distribution of \(\phi=(ka,V,Cl)\). Then, the prior distribution \(p(f(t,\phi) ; \zeta,\sigma)\) is estimated using \(f(t,\phi_1),\ldots, f(t,\phi_N)\).
We can estimate, for instance, the prior median of \(f(t,\phi)\) as well as any prediction interval.
Similarly, the posterior distribution \(p(f(t,\phi) | y ; \zeta,\sigma)\) can be estimated by Markov Chain Monte Carlo (MCMC).
Since the residual errors are normally distributed with a constant error model, the value of \(\phi\) that maximizes the conditional distribution of \(y\) is the least-square estimate:
\[ \begin{aligned} \phi_{\rm LS} & = \argmax{\phi} p(y|\phi;\sigma) \\ & = \argmin{\phi} \sum_{j=1}^n (y_j - f(t_j,\phi))^2 \end{aligned} \]
Menus
Plots
- data and predictions: display the data \(y=(y_1,\ldots y_n)\) and various predictions of \(f(t,\phi)\),
- parameters distributions display the prior and posterior distributions of \(ka\), \(V\) and \(Cl\),
- MCMC: display the sequences \((ka_k)\), \((V_k)\) and \((Cl_k)\).
Display
Select to display:
- the data \(y=(y_1,\ldots y_n)\),
- The Least-Square estimate \(f(t,\phi_{\rm LS})\) and a confidence interval
- The estimated prior median of \(f(t,\phi)\) and a prediction interval (estimated by Monte Carlo simulations) ,
- The estimated posterior median of \(f(t,\phi)\) and a prediction interval (estimated by MCMC).
Distributions
- prior: modify interactively the prior distribution of \(\phi\), i.e. the value of \(\zeta = (ka^{\star},V^{\star},Cl^{\star},\omega_{ka},\omega_{V},\omega_{Cl},\sigma)\).
- observations modify interactively the residual error parameter \(\sigma\) and the number of data points. You can also interactively add or remove data points by selecting isolated points (clicking) and regions (brushing, i.e. clicking and dragging a selection box).
MCMC
Modify the total number of MCMC iterations and the number of burning iterations (iterations which are not used for the estimation).
Tasks
Display the data and the Least-Square estimate of \(f\): this model provides the best fit for the observed data
Display the prior median of \(f\): this model only uses prior information and ignores the data.
Display the posterior median: this model combines the information given by the data with the prior.
Modify the value of \(\sigma\) and look at the impact on the posterior mean.
Modify the number of data points and look at the impact on the posterior mean.
Set \(\sigma=3\) and \(n=100\). Select exclude data to remove the first data points, either by clicking on brushing a region and look at the impact on the shape of the posterior mean during the absorption phase.
Set \(\sigma=5\) and \(n=25\). Select Distributions->prior and modify the values of the prior standard deviations \(\omega_{ka}\), \(\omega_V\) and \(\omega_Cl\).
Display the prediction and confidence intervals of \(f\) and repeat steps 1–7.
This application is governed by the CeCILL-B license.
You can use, modify and/or redistribute these codes under the terms of the
CeCILL license.
Marc Lavielle
Inria Saclay, Popix team
June 29th, 2015