## ---- include=FALSE------------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = '#>')

## ----load_packages-------------------------------------------------------
library(data.table)
library(httk)
library(ggplot2) 
library(scales)

## ----read_sens_data, eval=FALSE------------------------------------------
#  model <- '3compartmentss'
#  css.method <- 'analytic'
#  
#  #All combinations of poormetab and fup.censor
#  pmfc <- expand.grid(poormetab=c(TRUE, FALSE),
#                      fup.censor=c(TRUE, FALSE),
#                      fuptofub=TRUE)
#  #Initialize list of data.tables
#  dt.list <- vector(mode="list", length=nrow(pmfc))
#  
#  for (i in 1:nrow(pmfc)){
#    dt.list[[i]] <- readRDS(paste0('data/',
#                     'sens_glenisaacs_nhanes_',
#                     'allchems_',
#                     'allgroups_',
#                     model,
#                     '_',
#                     css.method,
#                     '_fup_censor_',
#                     pmfc[i, 'fup.censor'],
#                     '_poormetab_',
#                     pmfc[i, 'poormetab'],
#                     ifelse(pmfc[i, 'fuptofub'], "_FuptoFub", ""),
#                     '.Rdata'))
#    dt.list[[i]][, poormetab:=pmfc[i, 'poormetab']]
#    dt.list[[i]][, fup.censor:=pmfc[i, 'fup.censor']]
#    dt.list[[i]][, fuptofub:=pmfc[i, 'fuptofub']]
#  }
#  
#  dat <- rbindlist(dt.list)
#  rm(dt.list)
#  setnames(dat, c('chemcas', 'rn'),
#           c('CAS', 'param'))

## ----get_measured_vals, eval=FALSE---------------------------------------
#  cheminfo.dt <- as.data.table(httk::get_cheminfo(model='3compartmentss',
#                                            info=c('CAS','Compound',
#                                                   'Funbound.plasma'),
#                                            exclude.fub.zero=FALSE))
#  setnames(cheminfo.dt,
#           'Human.Funbound.plasma',
#           'Funbound.plasma'
#  )
#  
#  cheminfo.dt[, Clint:=sapply(CAS,
#                              function(x) parameterize_steadystate(chem.cas=x,
#                                                                   species='Human')$Clint)]
#  
#  dat <- merge(dat,
#                    cheminfo.dt[, .(CAS,
#                                    Compound,
#                                    Funbound.plasma,
#                                    Clint)],
#                    by='CAS')
#  
#    setnames(dat, 'Clint', 'CLint')
#    dat[param=='Clint', param:='CLint']
#    dat[param=='Funbound.plasma', param:='Fub']
#    dat[, param:=factor(param,
#                                levels=c('CLint',
#                                         'Fub',
#                                         'phys.par'))]

## ----change_factor_levels, eval=FALSE------------------------------------
#  
#    #change poormetab factor levels to be informative
#    dat[, poormetab:=factor(poormetab,
#                                levels=c(FALSE, TRUE),
#                                labels=c('excluded',
#                                         'included'))]
#    #change fup.censor factor levels to be informative
#    dat[, fup.censor:=factor(fup.censor,
#                                 levels=c(FALSE, TRUE),
#                                 labels=c('<lod=lod/2',
#                                          '<lod censored'))]
#    #change fup.censor name is be more informative
#    setnames(dat, 'fup.censor', 'Fub')

## ----my_labeller_fun, eval=FALSE-----------------------------------------
#  my_labeller_fun <- function(DF){
#   if (variable=='param'){
#     return(paste('i =', value))
#   }else{
#     return(paste(variable,value))
#   }
#  }

## ----handle_clint_zero, eval=FALSE---------------------------------------
#  dat[CLint==0, CLint:=1e-8] #plot them as 1e-8 but label them as zero
#  #Hack the log-scaled x-axis so that CLint==0 points have some buffer space for display
#  CLint.breaks <- c(1e-8, 1e-6, 1e-4, 1e-2, 1e0, 1e2, 1e4)
#  #Label the CLint breaks
#  CLint.labels <- c(expression(0^{paste(' ')}), expression(10^{-6}),
#                    expression(10^{-4}),
#                    expression(10^{-2}), expression(10^0),
#                    expression(10^2), expression(10^4))

## ----define_sens_colormap, eval=FALSE------------------------------------
#  ylgnbu <- RColorBrewer::brewer.pal(n=9, name="YlGnBu")

## ----global_sens_plots, eval=FALSE---------------------------------------
#  for (ecg in c('Total',
#                    'Age.6.11',
#                    'Age.12.19',
#                    'Age.20.65',
#                    'Age.GT65',
#                    'BMIgt30',
#                    'BMIle30',
#                    'Males',
#                    'Females',
#                    'ReproAgeFemale')){
#    for (fpb in c(TRUE)){
#      dt.tmp <- dat[param %in% c('CLint',
#                                         'Fub',
#                                         'phys.par') &
#                           ExpoCast.group==ecg &
#                           fuptofub==fpb, ]
#      setorder(dt.tmp,Si)
#      p<-ggplot(data=dt.tmp) +
#          geom_point(aes(x=CLint,
#                         y=Funbound.plasma,
#                         color=Si,
#                         fill=Si,
#                         order=Si),
#                     shape=21)+
#          scale_color_gradientn(colours=ylgnbu,
#                                guide='colourbar',
#                                limits=c(0,1), #specify limits, breaks, and labels to get the colorbar tick marks where I want them
#                                breaks=c(0,0.25,0.5,0.75,1),
#                                labels=as.character(c(0,0.25,0.5,0.75,1)),
#                                oob=squish,
#                                name='First-order \nCss sensitivity \nto param i')+ #there may be some very small negative indexes; just treat them as zero
#          scale_fill_gradientn(colours=ylgnbu,
#                               guide=FALSE, #don't plot separate colorbar for fill
#                               na.value=NA, #don't fill points for NAs
#                               limits=c(0,1), #specify limits, breaks, and labels to get the colorbar tick marks where I want them
#                               breaks=c(0,0.25,0.5,0.75,1),
#                               labels=as.character(c(0,0.25,0.5,0.75,1)),
#                               oob=squish)+
#          scale_x_continuous(trans='log10', #use the hack to get CLint==0 points to not be clipped
#                             breaks=CLint.breaks,
#                             labels=CLint.labels) +
#          scale_y_continuous(limits=c(0,1),
#                             oob=squish) + #Funbound.plasma shouldn't be >1 -- if it is, just pretend it's 1
#  #         facet_grid(poormetab*Fub~param,
#  #                    labeller=my_labeller_fun)+
#    facet_grid(poormetab*Fub~param,
#               labeller=labeller(poormetab=function(x) paste("CLint: PM ", x),
#                                          Fub = label_both,
#                                          param=function(x) paste("i = ", x))) +
#          labs(x='Measured CLint (uL/min/million cells)',
#               y='Measured Fub',
#               title=ecg) +
#          theme(legend.title=element_text(size=12), #bump up the legend text size
#                legend.text=element_text(size=12),
#                strip.background=element_blank())
#    ggsave(filename=paste0('pdf_figures/',
#                    paste('Si',
#                          'httkpop',
#                          'group',
#                          ecg,
#                          'model',
#                          model,
#                          "FuptoFub",
#                          fpb,
#                          sep='_'),
#                    '.pdf'),
#           plot=p,
#           height=8.5,
#           width=11)
#  print(p)
#      }
#  }

## ---- eval=FALSE---------------------------------------------------------
#  #Initialize list of data.tables
#  dt.list <- vector(mode="list", length=nrow(pmfc))
#  
#  for (i in 1:nrow(pmfc)){
#    dt.list[[i]] <- readRDS(paste0('data/',
#                     'sens_glenisaacs_nhanes_',
#                     'allchems_',
#                     'indepMC_',
#                     model,
#                     '_',
#                     css.method,
#                     '_fup_censor_',
#                     pmfc[i, 'fup.censor'],
#                     '_poormetab_',
#                     pmfc[i, 'poormetab'],
#                      ifelse(pmfc[i, 'fuptofub'], "_FuptoFub", ""),
#                     '.Rdata'))
#    dt.list[[i]][, poormetab:=pmfc[i, 'poormetab']]
#    dt.list[[i]][, fup.censor:=pmfc[i, 'fup.censor']]
#    dt.list[[i]][, fuptofub:=pmfc[i, 'fuptofub']]
#  }
#  
#  dat_indep <- rbindlist(dt.list)
#  setnames(dat_indep, c('chemcas', 'rn'),
#           c('CAS', 'param'))
#   dat_indep <- merge(dat_indep,
#                    cheminfo.dt[, .(CAS,
#                                    Compound,
#                                    Funbound.plasma,
#                                    Clint)],
#                    by='CAS')
#  
#    setnames(dat_indep, 'Clint', 'CLint')
#    dat_indep[param=='Clint', param:='CLint']
#    dat_indep[param=='Funbound.plasma', param:='Fub']
#    dat_indep[, param:=factor(param,
#                            levels=c('CLint',
#                                     'Fub',
#                                     'phys.par',
#                                     unique(param)[!(unique(param) %in%
#                                                       c('CLint',
#                                                         'Fub',
#                                                         'phys.par'))]))]
#    #change poormetab factor levels to be informative
#    dat_indep[, poormetab:=factor(poormetab,
#                                levels=c(FALSE, TRUE),
#                                labels=c('excluded',
#                                         'included'))]
#    #change fup.censor factor levels to be informative
#    dat_indep[, fup.censor:=factor(fup.censor,
#                                 levels=c(FALSE, TRUE),
#                                 labels=c('<lod=lod/2',
#                                          '<lod censored'))]
#    #change fup.censor name is be more informative
#    setnames(dat_indep, 'fup.censor', 'Fub')
#  dat_indep[CLint==0, CLint:=1e-8] #plot them as 1e-8 but label them as zero

## ----indep_params, eval=FALSE--------------------------------------------
#  for (fpb in c(TRUE)){
#    dt.tmp <- dat_indep[param %in% c('CLint',
#                                         'Fub',
#                                         'phys.par') &
#                             fuptofub==fpb, ]
#    setorder(dt.tmp, Si)
#  p<-ggplot(data=dt.tmp) +
#          geom_point(aes(x=CLint,
#                         y=Funbound.plasma,
#                         color=Si,
#                         fill=Si,
#                         order=Si),
#                     shape=21)+
#          scale_color_gradientn(colours=ylgnbu,
#                                guide='colourbar',
#                                limits=c(0,1), #specify limits, breaks, and labels to get the colorbar tick marks where I want them
#                                breaks=c(0,0.25,0.5,0.75,1),
#                                labels=as.character(c(0,0.25,0.5,0.75,1)),
#                                oob=squish,
#                                name='First-order \nCss sensitivity \nto param i')+ #there may be some very small negative indexes; just treat them as zero
#          scale_fill_gradientn(colours=ylgnbu,
#                               guide=FALSE, #don't plot separate colorbar for fill
#                               na.value=NA, #don't fill points for NAs
#                               limits=c(0,1), #specify limits, breaks, and labels to get the colorbar tick marks where I want them
#                               breaks=c(0,0.25,0.5,0.75,1),
#                               labels=as.character(c(0,0.25,0.5,0.75,1)),
#                               oob=squish)+
#          scale_x_continuous(trans='log10', #use the hack to get CLint==0 points to not be clipped
#                             breaks=CLint.breaks,
#                             labels=CLint.labels) +
#          scale_y_continuous(limits=c(0,1),
#                             oob=squish) + #Funbound.plasma shouldn't be >1 -- if it is, just pretend it's 1
#  #         facet_grid(poormetab*Fub~param,
#  #                    labeller=my_labeller_fun)+
#    facet_grid(poormetab*Fub~param,
#               labeller=labeller(poormetab=function(x) paste("CLint: PM ", x),
#                                          Fub = label_both,
#                                          param=function(x) paste("i = ", x))) +
#          labs(x='Measured CLint (uL/min/million cells)',
#               y='Measured Fub',
#               title='Independent MC') +
#          theme(legend.title=element_text(size=12), #bump up the legend text size
#                legend.text=element_text(size=12),
#                strip.background=element_blank())
#    ggsave(filename=paste0('pdf_figures/',
#                    paste('Si',
#                          'indepMC',
#                          'model',
#                          model,
#                          "FuptoFub",
#                          fpb,
#                          sep='_'),
#                    '.pdf'),
#           plot=p,
#           height=8.5,
#           width=11)
#  print(p)
#  }

## ----indep_phys, eval=FALSE----------------------------------------------
#  for (fpb in c(TRUE)){
#    dt.tmp <- dat_indep[param %in% c('million.cells.per.gliver',
#                                          'Qgfrc',
#                                          'Qtotal.liverc',
#                                          'Vliverc') &
#                             fuptofub==fpb, ]
#    setorder(dt.tmp, Si)
#  p<-ggplot(data=dt.tmp) +
#          geom_point(aes(x=CLint,
#                         y=Funbound.plasma,
#                         color=Si,
#                         fill=Si,
#                         order=Si),
#                     shape=21)+
#          scale_color_gradientn(colours=ylgnbu,
#                                guide='colourbar',
#                                limits=c(0,1), #specify limits, breaks, and labels to get the colorbar tick marks where I want them
#                                breaks=c(0,0.25,0.5,0.75,1),
#                                labels=as.character(c(0,0.25,0.5,0.75,1)),
#                                oob=squish,
#                                name='First-order \nCss sensitivity \nto param i')+ #there may be some very small negative indexes; just treat them as zero
#          scale_fill_gradientn(colours=ylgnbu,
#                               guide=FALSE, #don't plot separate colorbar for fill
#                               na.value=NA, #don't fill points for NAs
#                               limits=c(0,1), #specify limits, breaks, and labels to get the colorbar tick marks where I want them
#                               breaks=c(0,0.25,0.5,0.75,1),
#                               labels=as.character(c(0,0.25,0.5,0.75,1)),
#                               oob=squish)+
#          scale_x_continuous(trans='log10', #use the hack to get CLint==0 points to not be clipped
#                             breaks=CLint.breaks,
#                             labels=CLint.labels) +
#          scale_y_continuous(limits=c(0,1),
#                             oob=squish) + #Funbound.plasma shouldn't be >1 -- if it is, just pretend it's 1
#  #         facet_grid(poormetab*Fub~param,
#  #                    labeller=my_labeller_fun)+
#    facet_grid(poormetab*Fub~param,
#               labeller=labeller(poormetab=function(x) paste("CLint: PM ", x),
#                                          Fub = label_both,
#                                          param=function(x) paste("i = ", x))) +
#          labs(x='Measured CLint (uL/min/million cells)',
#               y='Measured Fub',
#               title='Independent MC') +
#          theme(legend.title=element_text(size=12), #bump up the legend text size
#                legend.text=element_text(size=12),
#                strip.background=element_blank())
#    ggsave(filename=paste0('pdf_figures/',
#                    paste('Si',
#                          'indepMC',
#                          'physparams',
#                          'model',
#                          model,
#                          "FuptoFub",
#                          fpb,
#                          sep='_'),
#                    '.pdf'),
#           plot=p,
#           height=8.5,
#           width=11)
#  print(p)
#  }

