Observations and predictions
Predicted concentration
Parameters
Statistical results
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
# April 29th, 2015
#-------------------------------------------------------------------------
library(shinydashboard)
sidebar <- dashboardSidebar(
hr(),
sidebarMenu(id="tabs",
menuItem("Plots", icon = icon("line-chart"),
menuSubItem("observations vs time", tabName = "plot1", icon = icon("angle-right")), selected=TRUE,
menuSubItem("observations vs predictions", tabName = "plot2", icon = icon("angle-right")),
menuSubItem("residuals vs time", tabName = "plot3", icon = icon("angle-right")),
menuSubItem("residuals vs predictions", tabName = "plot4", icon = icon("angle-right"))
),
menuItem("Tables", icon = icon("table"), tabName="tabtable",
menuSubItem("parameters", tabName = "tabparam", icon = icon("angle-right")),
menuSubItem("observations and predictions", tabName = "tabobs", icon = icon("angle-right")),
menuSubItem("predicted concentration", tabName = "tabconc", icon = icon("angle-right")),
menuSubItem("statistics", tabName = "tabstat", icon = icon("angle-right"))
),
menuItem("Codes", icon = icon("file-text-o"),
menuSubItem("Mlxtran - first order", tabName = "pkmodel1", icon = icon("angle-right")),
menuSubItem("Mlxtran - zero order", tabName = "pkmodel0", icon = icon("angle-right")),
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' | input.tabs=='plot2' | input.tabs=='plot3' | input.tabs=='plot4'",
sidebarMenu(
menuItem("Absorption model", icon = icon("chevron-circle-right"),
fluidRow(
column(1),
column(10, radioButtons("absorption","",c("first order" = "1","zero order" = "0")))
)
),
menuItem("Outputs", icon = icon("chevron-circle-right"),
fluidRow(
column(1),
column(10,
# checkboxInput("checkObs", "Observations", value=TRUE),
checkboxInput("checkIni", "Initialization", value=TRUE),
checkboxInput("checkMLE", "Estimation", value=FALSE)
)
)
),
menuItem("Scale", icon = icon("chevron-circle-right"),
fluidRow(
column(1),
column(10, radioButtons("log","",c("linear" = "1","log" = "2")))
)
)
),
hr()
)
)
body <- dashboardBody(
tabItems(
tabItem(tabName = "plot1",
fluidRow(
box(width = 12, status = "primary",
plotOutput("plot1", height="400px")
)
)
),
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 = "plot4",
fluidRow(
box(width = 12, status = "primary",
plotOutput("plot4", height="400px")
)
)
),
tabItem(tabName = "tabobs",
box( width = NULL, status = "primary", solidHeader = TRUE, title="Observations and predictions",
tableOutput("tabley")
)
),
tabItem(tabName = "tabconc",
box( width = NULL, status = "primary", solidHeader = TRUE, title="Predicted concentration",
tableOutput("tablef")
)
),
tabItem(tabName = "tabparam",
box( width = NULL, status = "primary", solidHeader = TRUE, title="Parameters",
tableOutput("tablep")
)
),
tabItem(tabName = "tabstat",
box( width = NULL, status = "primary", solidHeader = TRUE, title="Statistical results",
tableOutput("tabler")
)
),
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 = "pkmodel0",
box( width = NULL, status = "primary", solidHeader = TRUE, title="pkmodel0.txt",
pre(includeText("pkmodel0.txt"))
)
),
tabItem(tabName = "pkmodel1",
box( width = NULL, status = "primary", solidHeader = TRUE, title="pkmodel1.txt",
pre(includeText("pkmodel1.txt"))
)
),
tabItem(tabName = "readme",
withMathJax(),
includeMarkdown("readMe.Rmd")
),
tabItem(tabName = "about",
includeMarkdown("../../about/about.Rmd")
)
),
conditionalPanel(condition = "(input.checkIni=='1') & (input.tabs=='plot1'|input.tabs=='plot2'|input.tabs=='plot3'|input.tabs=='plot4') ",
box(width = 4, status = "primary", solidHeader = FALSE,
conditionalPanel(condition = "input.absorption == '1'",
sliderInput("ka", "ka:", value = 0.8, min = 0, max = 1, step=0.05),
conditionalPanel(condition = "input.checkMLE == '1'",
tags$style(type='text/css', '#estka {color: red; font-size: 14px;}'),
verbatimTextOutput("estka")
)
),
conditionalPanel(condition = "input.absorption == '0'",
sliderInput("tk0", "duration Tk0:", value = 3, min = 0, max = 10, step=0.5),
conditionalPanel(condition = "input.checkMLE == '1'",
tags$style(type='text/css', '#estTk0 {color: red; font-size: 14px;}'),
verbatimTextOutput("estTk0")
)
)
),
box(width = 4, status = "primary", solidHeader = FALSE,
sliderInput("V", "V:", value = 15, min = 1, max = 20, step=1),
conditionalPanel(condition = "input.checkMLE == '1'",
tags$style(type='text/css', '#estV {color: red; font-size: 14px;}'),
verbatimTextOutput("estV")
)
),
box(width = 4, status = "primary", solidHeader = FALSE,
sliderInput("Cl", "Cl:", value = 5, min = 0.5, max = 10, step=0.5),
conditionalPanel(condition = "input.checkMLE == '1'",
tags$style(type='text/css', '#estCl {color: red; font-size: 14px;}'),
verbatimTextOutput("estCl")
)
)
)
)
dashboardPage(
dashboardHeader(title = "Fit a PK model"),
sidebar,
body
)
server.R
Download
library(mlxR)
shinyServer(function(input, output) {
d <- read.csv(file="individualFitting_data.txt", header=TRUE, sep="\t", quote="\"")
names(d)[2] <- "y"
Cc <- list(name='cc', time=seq(0,30,length.out=100))
ypred <- list(name='ypred', time=d$time)
adm <- list(time=0, amount=50)
d.param <- 4
n <- length(d$time)
r <- reactive({
if (input$absorption==1){
p.name <- c('ka','V','Cl')
pkmodel <- "pkmodel1.txt"
}else{
p.name <- c('Tk0','V','Cl')
pkmodel <- "pkmodel0.txt"
}
p <- list(name=p.name, value=c(1,2,1))
s <- list(load.design=TRUE, data.in=TRUE)
dataIn <- simulx(model=pkmodel,output=list(Cc,ypred),treatment=adm,parameter=p,settings=s)
if (input$absorption==1){
p.value <- c(input$ka,input$V,input$Cl)
}else{
p.value <- c(input$tk0,input$V,input$Cl)
}
dataIn$individual_parameters$value[1,]=p.value
s <- list(load.design=FALSE)
res = simulx(data=dataIn)
f <- data.frame(time=res[[1]][,1])
parameter <- data.frame(initial=p.value)
row.names(parameter) <- p.name
f$f_ini=res[[1]][,2]
d$ypred_ini=res[[2]][,2]
r=list()
r$sse_ini=sum((d$y-d$ypred_ini)^2)
r$textini <- paste("RSSE =",formatC(r$sse_ini,digits=3))
result <- rsse(r$sse_ini,n,d.param)
row.names(result)="initial"
if (input$checkMLE==TRUE){
p <- list(name=p.name, value=p.value)
rest <- mlxoptim(model=pkmodel,
output='cc',
treatment=adm,
data=d,
initial=p)
r$sse_est=rest$sse
r$textest <- paste("RSSE =",formatC(r$sse_est,digits=3))
s <- list(load.design=TRUE)
res <- simulx(model=pkmodel,output=list(Cc,ypred),treatment=adm,parameter=rest$parameter,settings=s)
d$ypred_est <- res[[2]][,2]
f$f_est <- res[[1]][,2]
parameter[,2] <-rest$parameter$value
names(parameter)[2]="estimated"
res.est <- rsse(rest$sse,n,d.param)
row.names(res.est) <- "estimate"
result[2,] <- res.est
}
r$parameter <- parameter
r$y <- d
r$f <- f
r$result <- t(result)
return(r)
})
szp <- 3.5
szl <- 1
dxini <- 0.9
dxest <- 0.9
dyini <- 1.07
dyest <- 1
output$plot1 <- renderPlot({
r=r()
pl <- ggplotmlx()
pl <- pl + geom_point(data=r$y, aes(x=time, y=y), color="blue", size=szp)
if (input$checkIni==TRUE){
pl <- pl + geom_line(data=r$f, aes(x=time, y=f_ini), size=szl)
}
if (input$checkMLE==TRUE){
pl <- pl + geom_line(data=r$f, aes(x=time, y=f_est), color="red", size=szl)
}
if (input$log==2){
pl <- pl + scale_y_log10()
}
my.ggp.yrange <- ggplot_build(pl)$panel$ranges[[1]]$y.range
my.ggp.xrange <- ggplot_build(pl)$panel$ranges[[1]]$x.range
xini <- my.ggp.xrange[1] + (my.ggp.xrange[2]-my.ggp.xrange[1])*dxini
xest <- my.ggp.xrange[1] + (my.ggp.xrange[2]-my.ggp.xrange[1])*dxest
yini <- my.ggp.yrange[1] + (my.ggp.yrange[2]-my.ggp.yrange[1])*dyini
yest <- my.ggp.yrange[1] + (my.ggp.yrange[2]-my.ggp.yrange[1])*dyest
if (input$checkIni==TRUE){
pl <- pl + annotate("text", label=r$textini, x=xini, y=yini, size=6, colour="black")
}
if (input$checkMLE==TRUE){
pl <- pl + annotate("text", label=r$textest, x=xest, y=yest, size=6, colour="red")
}
print(pl)
})
output$plot2 <- renderPlot({
r=r()
pl <- ggplotmlx()
if (input$checkIni==TRUE){
pl <- pl + geom_point(data=r$y, aes(x=ypred_ini, y=y), size=szp)
}
if (input$checkMLE==TRUE){
pl <- pl + geom_point(data=r$y, aes(x=ypred_est, y=y), color="red", size=szp)
}
pl <- pl + geom_abline(intercept = 0, slope = 1, color="blue", size=szl) + xlab("y_pred")
my.ggp.yrange <- ggplot_build(pl)$panel$ranges[[1]]$y.range
my.ggp.xrange <- ggplot_build(pl)$panel$ranges[[1]]$x.range
xini <- my.ggp.xrange[1] + (my.ggp.xrange[2]-my.ggp.xrange[1])*dxini
xest <- my.ggp.xrange[1] + (my.ggp.xrange[2]-my.ggp.xrange[1])*dxest
yini <- my.ggp.yrange[1] + (my.ggp.yrange[2]-my.ggp.yrange[1])*dyini
yest <- my.ggp.yrange[1] + (my.ggp.yrange[2]-my.ggp.yrange[1])*dyest
if (input$checkIni==TRUE){
pl <- pl + annotate("text", label=r$textini, x=xini, y=yini, size=6, colour="black")
}
if (input$checkMLE==TRUE){
pl <- pl + annotate("text", label=r$textest, x=xest, y=yest, size=6, colour="red")
}
print(pl)
})
output$plot3 <- renderPlot({
r=r()
pl <- ggplotmlx()
if (input$checkIni==TRUE){
r$y$e <- r$y$y - r$y$ypred_ini
pl <- pl + geom_point(data=r$y, aes(x=time, y=e), size=szp)
}
if (input$checkMLE==TRUE){
r$y$e <- r$y$y - r$y$ypred_est
pl <- pl + geom_point(data=r$y, aes(x=time, y=e), color="red", size=szp)
}
pl <- pl + geom_hline(aes(yintercept=0), color="blue", size=szl)
my.ggp.yrange <- ggplot_build(pl)$panel$ranges[[1]]$y.range
my.ggp.xrange <- ggplot_build(pl)$panel$ranges[[1]]$x.range
xini <- my.ggp.xrange[1] + (my.ggp.xrange[2]-my.ggp.xrange[1])*dxini
xest <- my.ggp.xrange[1] + (my.ggp.xrange[2]-my.ggp.xrange[1])*dxest
yini <- my.ggp.yrange[1] + (my.ggp.yrange[2]-my.ggp.yrange[1])*dyini
yest <- my.ggp.yrange[1] + (my.ggp.yrange[2]-my.ggp.yrange[1])*dyest
if (input$checkIni==TRUE){
pl <- pl + annotate("text", label=r$textini, x=xini, y=yini, size=6, colour="black")
}
if (input$checkMLE==TRUE){
pl <- pl + annotate("text", label=r$textest, x=xest, y=yest, size=6, colour="red")
}
print(pl)
})
output$plot4 <- renderPlot({
r=r()
pl <- ggplotmlx()
if (input$checkIni==TRUE){
r$y$e <- r$y$y - r$y$ypred_ini
pl <- pl + geom_point(data=r$y, aes(x=ypred_ini, y=e), size=szp) }
if (input$checkMLE==TRUE){
r$y$e <- r$y$y - r$y$ypred_est
pl <- pl + geom_point(data=r$y, aes(x=ypred_est, y=e), color="red", size=szp) }
pl <- pl + geom_hline(aes(yintercept=0), color="blue", size=szl)
my.ggp.yrange <- ggplot_build(pl)$panel$ranges[[1]]$y.range
my.ggp.xrange <- ggplot_build(pl)$panel$ranges[[1]]$x.range
xini <- my.ggp.xrange[1] + (my.ggp.xrange[2]-my.ggp.xrange[1])*dxini
xest <- my.ggp.xrange[1] + (my.ggp.xrange[2]-my.ggp.xrange[1])*dxest
yini <- my.ggp.yrange[1] + (my.ggp.yrange[2]-my.ggp.yrange[1])*dyini
yest <- my.ggp.yrange[1] + (my.ggp.yrange[2]-my.ggp.yrange[1])*dyest
if (input$checkIni==TRUE){
pl <- pl + annotate("text", label=r$textini, x=xini, y=yini, size=6, colour="black")
}
if (input$checkMLE==TRUE){
pl <- pl + annotate("text", label=r$textest, x=xest, y=yest, size=6, colour="red")
}
print(pl)
})
output$tabley <- renderTable({
r()$y
})
output$tablef <- renderTable({
r()$f
})
output$tablep <- renderTable({
r()$parameter
})
output$tabler <- renderTable({
r()$result
})
output$estka <- renderText({
r <- round(r()$parameter[1,],3)
pest=NULL
if (length(r)==2)
pest <- paste0("ka_est = ",r[2])
return(pest)
})
output$estTk0 <- renderText({
r <- round(r()$parameter[1,],3)
pest=NULL
if (length(r)==2)
pest <- paste0("Tk0_est = ",r[2])
return(pest)
})
output$estV <- renderText({
r <- round(r()$parameter[2,],3)
pest=NULL
if (length(r)==2)
pest <- paste0("V_est = ",r[2])
return(pest)
})
output$estCl <- renderText({
r <- round(r()$parameter[3,],3)
pest=NULL
if (length(r)==2)
pest <- paste0("Cl_est = ",r[2])
return(pest)
})
})
rsse <- function(sse,n,d)
{
deviance <- n*(1+ log(sse/n) +log(2*pi))
BIC <- deviance + log(n)*d
AIC <- deviance + 2*d
r <- data.frame(n=n,d=d,sse=sse,deviance=deviance,AIC,BIC)
return(r)
}
list <- structure(NA,class="result")
"[<-.result" <- function(x,...,value) {
args <- as.list(match.call())
args <- args[-c(1:2,length(args))]
length(value) <- length(args)
for(i in seq(along=args)) {
a <- args[[i]]
if(!missing(a)) eval.parent(substitute(a <- v,list(a=a,v=value[[i]])))
}
return(x)
}
mlxoptim <- function(model=NULL,output=NULL,treatment=NULL,data=NULL,initial=NULL)
{
errpred <- function(pp,y,dataIn){
dataIn$individual_parameters$value[1,]=exp(pp)
res = simulx(data=dataIn)
y.pred=res[[1]][,2]
e=sum((y.pred-y)^2)
return(e)}
t.obs <- data[,1]
y.obs <- data[,2]
dd <- simulx(model=model,
parameter=initial,
output=list(name=output,time=t.obs),
treatment=treatment,
settings=list(data.in=TRUE,load.design=TRUE))
l0=log(initial$value)
r <- optim(l0,errpred,y=y.obs,dataIn=dd)
n <- length(t.obs)
pest=initial
pest$value=exp(r$par)
final = list(parameter=pest,sse=r$value)
return(final)
}
pkmodel0.txt
[LONGITUDINAL]
input = {Tk0, V, Cl}
EQUATION:
cc = pkmodel(Tk0, V, Cl)
ypred = cc
pkmodel1.txt
[LONGITUDINAL]
input = {ka, V, Cl}
EQUATION:
cc = pkmodel(ka, V, Cl)
ypred = cc
Objective
Fit a PK model to PK data from a single individual.
Tasks
Using a first order absorption process, try to find some values of \((k_a, V, Cl)\) that fit well the observed PK data (blue dots),
Compute the least-squares estimate of \((k_a, V, Cl)\) (by checking
Estimationin the menuOutputs),Look at the diagnostic plots and compare the numerical results in the Table
Results,Repeat the exercice with a zero-order absorption process,
What is your conclusion?
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.
This Shiny application requires the mlxR package.
Marc Lavielle
Inria Saclay, Popix team
April 29th, 2015