individual fitting - comparing structural models I


















shinyUI(fluidPage(
  titlePanel(
    list(HTML('<p style="color:#4C0B5F; font-size:24px" fontsize=14>individual fitting - comparing structural models I</p>' )),
    windowTitle="individual fitting"),
  tabsetPanel(
    tabPanel("Plot",
             fluidRow(
               column(2,
                      br(),
                      br(),
                      radioButtons("model", "PK model", c("iv bolus" = "1","oral 1 cpt" = "2","oral 2 cpt" = "3")),                    
                      conditionalPanel(condition = "input.estim == 'FALSE'",
                                       conditionalPanel(condition = "input.model == '2'",
                                                        sliderInput("tk01", "Tk0:", value = 1, min = 0, max = 10, step=0.5)
                                       ),
                                       conditionalPanel(condition = "input.model == '3'",
                                                        sliderInput("tk02", "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),
                                       conditionalPanel(condition = "input.model == '3'",
                                                        sliderInput("k12", "k12:", value = 0.5, min = 0.1, max = 1, step=0.1),
                                                        sliderInput("k21", "k21:", value = 0.5, min = 0.1, max = 1, step=0.1)
                                       )
                      ),
                      conditionalPanel(condition = "input.estim == 'TRUE'",
                                       br(), tableOutput("tabler")
                      ),
                      br(),
                      br()
               ),
               column(10,
                      column(10, 
                             plotOutput("plot1")
                      ),
                      column(2, 
                             br(),br(),br(),
                             radioButtons("estim", "", c("initial" = FALSE,"estimation" = TRUE)),                    
                             br(),br(),br(),br(),br(),br(),
                             radioButtons("log", "", c("linear scale" = FALSE,"log scale" = TRUE))                    
                      )
               ))),
    tabPanel("Residuals", 
             column(3, br(),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({
    if (input$model==1){
      p.name <- c('V','k')
      p.value <- c(input$V,input$k)
      pkmodel <- "../pkiv_model.txt"
    }else if (input$model==2){
      p.name <- c('Tk0','V','k')
      p.value <- c(input$tk01,input$V,input$k)
      pkmodel <- "../pk0_model.txt"
    }else{
      p.name <- c('Tk0','V','k','k12','k21')
      p.value <- c(input$tk02,input$V,input$k,input$k12,input$k21)
      pkmodel <- "../pk02_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)
      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()
    pl <- ggplotmlx() 
    
    if (input$estim==FALSE){
      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{
      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$tabler <- renderTable({
    r()$result},
    digits = 4
  )
  
  output$tabley <- renderTable({ 
    r()$y
  })
  
})


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?