# Ring et al. 2017 Global sensitivity analysis plotting

#### Caroline Ring

#### 2018-01-23

This vignette includes the code to create the global sensitivity analysis plots included in the HTTK-Pop paper.

    library(data.table)
    library(httk)
    library(ggplot2)
    library(scales)

First, read in the global sensitivity analysis data. Do this for all four combinations of the `poormetab` and `fup.censor` conditions.

    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'))

Next, add measured values of `Funbound.plasma` and `CLint` to the data table. Get these measured values from `httk`.

    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'))]

Next, set up a few things to prepare for plotting. First, change the factor levels for `poormetab` and `fup.censor` to be more informative than just TRUE and 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')

And set up a labeller function to more informatively label the plots by parameter.

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

Next, since we’re going to plot measured CLint on a log scale, we need to handle the case where measured CLint was 0. It’s a bit kludgy, but we’ll just plot CLint = 0 as CLint = 1e-8, so that those points will actually show up on the log-scale plot.

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

And define the colormap to use for the sensitivity values. I chose a multi-hue sequential palette from ColorBrewer, because the variations in hue make it easier to see where sensitivity index values fall into certain bins, and the yellow-green-blue palette shows up well against a gray plot background.

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

Finally, let’s actually make the plots.

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

We’d also like to plot the sensitivity analysis results for independent Monte Carlo – both showing the same data as for correlated Monte Carlo (sensitivity to the physiological parameters as a group), and showing sensitivity to each of the physiological parameters independently.

So first, we read in the independent MC data.

    #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

Now do the plotting.

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

And plot each physiological parameter independently.

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