PK modelling: individual fitting
shinyUI(fluidPage(
titlePanel(
list(HTML('<p style="color:#4C0B5F; font-size:24px" fontsize=14>PK modelling: individual fitting</p>' )),
windowTitle="PK modelling"),
tabsetPanel(
tabPanel("Plot",
fluidRow(
column(2,
br(),
br(),
br(),
br(),
# selectInput("absorption","",c("first order absorption" = "1","zero order absorption" = "0")),
radioButtons("absorption", "absorption", c("first order" = '1',"zero order" = '0')),
conditionalPanel(condition = "input.absorption == '1'",
sliderInput("ka", "ka:", value = 0.8, min = 0, max = 2, step=0.05)
),
conditionalPanel(condition = "input.absorption == '0'",
sliderInput("tk0", "Tk0:", value = 1, min = 0, max = 10, step=0.5)
),
br(),
sliderInput("V", "V:", value = 15, min = 1, max = 20, step=1),
br(),
sliderInput("k", "k:", value = 0.5, min = 0.05, max = 1, step=0.05),
br(),
br()
),
column(10,
column(10,
br(),
tabsetPanel(
tabPanel("Fit", br(), plotOutput("plot1")),
tabPanel("Obs vs Pred", br(), plotOutput("plot2")),
tabPanel("Residuals", br(), plotOutput("plot3"))
)),
column(2,
br(),br(),br(),br(),br(),br(),
radioButtons("estim", "", c("initial" = FALSE,"estimation" = TRUE)),
br(),
radioButtons("log", "", c("linear scale" = FALSE,"log scale" = TRUE))
)
))),
tabPanel("Tables",
tabsetPanel(
tabPanel("Parameters", tableOutput("tablep")),
tabPanel("Observations", tableOutput("tabley")),
tabPanel("Concentration", tableOutput("tablef")),
tabPanel("Results", tableOutput("tabler"))
)),
tabPanel("ui.R", pre(includeText("ui.R"))),
tabPanel("server.R", pre(includeText("server.R"))),
tabPanel("ReadMe", withMathJax(), includeMarkdown("ReadMe.Rmd"))
)
))
source("../../initMlxR.R")
r.size <- 3.5
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','k')
pkmodel <- "pk1_model.txt"
}else{
p.name <- c('Tk0','V','k')
pkmodel <- "pk0_model.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$k)
}else{
p.value <- c(input$tk0,input$V,input$k)
}
# if (input$estim==FALSE){
dataIn$individual_parameters$value=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)
result <- rsse(r$sse_ini,n,d.param)
row.names(result)="initial"
# }
if (input$estim==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
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)
r
})
output$plot1 <- renderPlot({
r <- r()
pl <- ggplotmlx() + geom_point(data=r$y, aes(x=time, y=y), color="blue", size=r.size)
if (input$estim==FALSE){
pl <- pl + geom_line(data=r$f, aes(x=time, y=f_ini), size=0.75) +
ggtitle(paste("sum of squared errors of prediction =",formatC(r$sse_ini,digits=3)))
}else{
pl <- pl + geom_line(data=r$f, aes(x=time, y=f_est), color="red", size=0.75) +
ggtitle(paste("sum of squared errors of prediction =",formatC(r$sse_est,digits=3))) +
theme(plot.title = element_text(colour="red" ))
}
if (input$log==TRUE){
pl <- pl + scale_y_log10()
}
print(pl)
})
output$plot2 <- renderPlot({
r=r()
pl <- ggplotmlx()
if (input$estim==FALSE){
pl <- pl + geom_point(data=r$y, aes(x=ypred_ini, y=y), size=r.size) +
ggtitle(paste("sum of squared errors of prediction =",formatC(r$sse_ini,digits=3)))
}else{
pl <- pl + geom_point(data=r$y, aes(x=ypred_est, y=y), color="red", size=r.size) +
ggtitle(paste("sum of squared errors of prediction =",formatC(r$sse_est,digits=3))) +
theme(plot.title = element_text(colour="red" ))
}
pl <- pl + geom_abline(intercept = 0, slope = 1, color="blue") + xlab("y_pred")
print(pl)
})
output$plot3 <- renderPlot({
r=r()
pl <- ggplotmlx()
if (input$estim==FALSE){
r$y$e <- r$y$y - r$y$ypred_ini
pl <- pl + geom_point(data=r$y, aes(x=time, y=e), size=r.size) +
ggtitle(paste("sum of squared errors of prediction =",formatC(r$sse_ini,digits=3)))
}else{
r$y$e <- r$y$y - r$y$ypred_est
pl <- pl + geom_point(data=r$y, aes(x=time, y=e), color="red") +
ggtitle(paste("sum of squared errors of prediction =",formatC(r$sse_est,digits=3))) +
theme(plot.title = element_text(colour="red" ))
}
pl <- pl + geom_hline(aes(yintercept=0), color="blue")
print(pl)
})
output$tabley <- renderTable({
r()$y
})
output$tablef <- renderTable({
r()$f
})
output$tablep <- renderTable({
r()$parameter
})
output$tabler <- renderTable({
r()$result
})
})
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)
}
mlxoptim <- function(model=NULL,output=NULL,treatment=NULL,data=NULL,initial=NULL)
{
errpred <- function(pp,y,dataIn){
dataIn$individual_parameters$value=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)
}
Using a first order absorption process, try to find some values of \((k_a, V, k)\) that fit well the observed PK data (blue dots),
Compute the least-squares estimate of \((k_a, V, k)\) (by checking Estim.),
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?