individual fitting - The residual error model
















shinyUI(fluidPage(
  titlePanel(
    list(HTML('<p style="color:#4C0B5F; font-size:24px" fontsize=14>individual fitting - The residual error model</p>' )),
    windowTitle="individual fitting"),
  tabsetPanel(
    tabPanel("Plot",
             fluidRow(
               column(2,
                      br(),
                      conditionalPanel(condition = "input.estim == 'FALSE'",
                                       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()
                      ),
                      conditionalPanel(condition = "input.estim == 'TRUE'",
                                       br(), tableOutput("tabler")
                      ),
                      br(),
                      br()
               ),
               column(10,
                      column(10, 
                             plotOutput("plot1"),
                             column(1),
                             column(8,
                             radioButtons("log", "", c("linear scale" = FALSE,"log scale" = TRUE), inline=TRUE)
                             )
                      ),
                      column(2, 
                             br(),br(),
                             radioButtons("estim", "", c("initial" = FALSE,"estimation" = TRUE)),                    
                             br(),br(),br(),
                             radioButtons("model", "Error model", c("constant" = "1","proportional" = "2","combined" = "3","exponential" = "4")),                    
                             br()
                      )
               ))),
    tabPanel("Residuals", 
             column(3, 
                    br(),
                    radioButtons("plerr", "plot", c("residual vs time" = "1","residual vs prediction" = "2","observation vs prediction" = "3")),                    
                    br(),
                    radioButtons("table", "Table", c("hide" = FALSE,"see" = TRUE), inline=TRUE),
                    conditionalPanel(condition = "input.table == 'TRUE'",
                                     tableOutput("tabley"))
             ),
             column(9, plotOutput("plot3"))
    ),
    #     tabPanel("Results", br(), br(), 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)    
  
  r <- reactive({
    p.name <- c('Tk0','V','k')
    p.value <- c(input$tk0,input$V,input$k)
    pkmodel <- "../pk0_model.txt"
    p   <- list(name=p.name, value=p.value) 
    
    s <- list(load.design=TRUE, data.in=TRUE)
    dataIn <- simulx(model=pkmodel,output=list(Cc,ypred),treatment=adm,parameter=p,settings=s)
    s <- list(load.design=FALSE)
    res = simulx(data=dataIn)
    d$ypred <- res[[2]][,2]
    d$e <- d$y - d$ypred
    
    r=list()
    r$sse_ini=sum(d$e^2)
    if (input$estim==TRUE){
      p <- list(name=p.name, value=p.value) 
      rest <- mlxoptim(model=pkmodel,
                       output='cc',
                       treatment=adm,
                       data=d,
                       initial=p,
                       err.model=input$model)
      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 <- res[[2]][,2]
      d$e <- d$y - d$ypred
      r$result <- data.frame(estimate=rest$parameter$value)
      row.names(r$result) <- p.name
    }
    r$y <- d
    r$f <- data.frame(time=res[[1]][,1],f=res[[1]][,2])
    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), size=0.75)   
    }else{
      pl <- pl + geom_line(data=r$f, aes(x=time, y=f), color="red", size=0.75)  
    }
    
    if (input$log==TRUE){
      pl <- pl + scale_y_log10()
    }
    print(pl)
  })
  
  output$plot3 <- renderPlot({
    r=r()
    
    if (input$estim==FALSE){
      r.sse <- r$sse_ini
      r.color <- "black"
    }else{
      r.sse <- r$sse_est
      r.color <- "red"
    }
    pl <- ggplotmlx(data=r$y) 
    if (input$plerr==1){
      pl <- pl+ geom_point(aes(x=time, y=e), color=r.color) + geom_hline(aes(yintercept=0), color="blue")
    }else if (input$plerr==2){
      pl <- pl+ geom_point(aes(x=ypred, y=e), color=r.color) + geom_hline(aes(yintercept=0), color="blue")
    }else if (input$plerr==3){
      pl <- pl+ geom_point(aes(x=ypred, y=y), color=r.color) + geom_abline(intercept=0,slope=1, color="blue")
    }
    pl <- pl + ggtitle(paste("sum of squared errors of prediction =",formatC(r.sse,digits=3)))  +
      theme(plot.title = element_text(colour=r.color ))     
      
      print(pl)
  })
  
  output$tabler <- renderTable({
    r()$result},
    digits = 4
  )
  
  output$tabley <- renderTable({ 
    r()$y
  })
  
})


mlxoptim <- function(model=NULL,output=NULL,treatment=NULL,data=NULL,initial=NULL,err.model=1)
{
  errpred1 <- 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)}
  
  errpred2 <- function(pp,y,dataIn){
    dataIn$individual_parameters$value=exp(pp)
    res = simulx(data=dataIn)
    y.pred=res[[1]][,2]
    e=sum((1-y/y.pred)^2)
    return(e)}
  
  errpred3 <- function(pp,y,dataIn){
    dataIn$individual_parameters$value=exp(pp)
    res = simulx(data=dataIn)
    y.pred=res[[1]][,2]
    e=sum((1-y/y.pred)^2)
    return(e)}
  
  errpred4 <- function(pp,y,dataIn){
    dataIn$individual_parameters$value=exp(pp)
    res = simulx(data=dataIn)
    y.pred=res[[1]][,2]
    e=sum((log(y)-log(y.pred))^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)
  if (err.model==1){
    r <- optim(l0,errpred1,y=y.obs,dataIn=dd)
  }else if (err.model==2){
    r <- optim(l0,errpred2,y=y.obs,dataIn=dd)
  }else if (err.model==3){
    r <- optim(l0,errpred3,y=y.obs,dataIn=dd)
  }else{
    r <- optim(l0,errpred4,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?