# Global sensitivity analysis plotting

#### Caroline Ring

#### 2017-03-02

This vignette includes the code to create the global sensitivity analysis plots included in the 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)
    }