individual fitting - PK model
















shinyUI(fluidPage(
  titlePanel(
    list(HTML('<p style="color:#4C0B5F; font-size:24px" fontsize=14>individual fitting - PK model</p>' )),
    windowTitle="individual fitting"),
  tabsetPanel(
    tabPanel("Plot",
             fluidRow(
               column(2,
                      br(),
                      br(),
                      br(),
                      br(),
                      selectInput("absorption","",c("first order absorption" = "1","zero order absorption" = "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, 
                             tabsetPanel( 
                               tabPanel("Fit", plotOutput("plot1")),
                               tabPanel("Obs vs Pred", plotOutput("plot2")), 
                               tabPanel("Residuals", 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")  

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") 
    
    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)) +  
        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")  +
        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)) +  
        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)
}


  1. 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),

  2. Compute the least-squares estimate of \((k_a, V, k)\) (by checking Estim.),

  3. Look at the diagnostic plots and compare the numerical results in the Table Results,

  4. Repeat the exercice with a zero-order absorption process,

  5. What is your conclusion?