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=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=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


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

  2. Compute the least-squares estimate of \((k_a, V, Cl)\) (by checking Estimation in the menu Outputs),

  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?


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