## ----configure_knitr, eval = TRUE--------------------------------------------- knitr::opts_chunk$set(collapse = TRUE, comment = '#>') options(rmarkdown.html_vignette.check_title = FALSE) ## ----clear_memory, eval = TRUE------------------------------------------------ rm(list=ls()) ## ----runchunks, eval = TRUE--------------------------------------------------- # Set whether or not the following chunks will be executed (run): execute.vignette <- FALSE ## ----load_packages, eval = execute.vignette----------------------------------- # # # # # # Setup # # # # # library(httk) # run Armitage.R in vitro distribution model # library(tidyverse) # data wrangling # library(data.table) # data table functionality # library(ggplot2) # plotting # library(ggrepel) # jitter labels in plots # library(ggpubr) # arrange plots for publication # library(viridis) # colorblind palette # library(scales) # required for label_log # # #save figures and tables to working directory: yes/no # save_output = FALSE # # #for saving outputs # data.date <- "16JUL2025" # path_out <- "C:/Users/mscherer/OneDrive - Environmental Protection Agency (EPA) # /Profile/Desktop/IVDManuscriptData/manuscript figures 16JUL25/" # ## ----rmsle_func, eval = execute.vignette-------------------------------------- # #function to calculate RMSLE # sle <- function (actual, predicted) # { # log.act <- log10(actual) # log.pred <- log10(predicted) # good <- !is.na(log.act) & is.finite(log.act) & # !is.na(log.pred) & is.finite(log.pred) # log.act <- log.act[good] # log.pred <- log.pred[good] # return((log.act - log.pred)^2) # } # msle <- function (actual, predicted) # { # mean(sle(actual, predicted)) # } # rmsle <-function (actual, predicted) # { # sqrt(msle(actual, predicted)) # } ## ----load_IVD_Manuscript_data, eval = execute.vignette------------------------ # # IVD_Manuscript_data <- httk::Scherer2025.IVD # # cat(paste("Summarized data from ", # length(unique(IVD_Manuscript_data$Citation)), # " studies covering ", # length(unique(IVD_Manuscript_data$ChemicalName)), # " unique chemicals ", # "with ", # length(IVD_Manuscript_data$CellConcentration_uM), # " observations total.", # sep="") # ) # # #14 studies, 43 chemicals, 153 observations # # #number each row # IVD_Manuscript_data[,X:=1:length(IVD_Manuscript_data$ChemicalName)] # ## ----eval = execute.vignette-------------------------------------------------- # # ### Create assay_component_endpoint_name columns for assay parameters ### # # #create assay_name column # IVD_Manuscript_data[, assay_name := # gsub(", ", "", IVD_Manuscript_data$Citation)] # # #create the assay_component_endpoint_name column # IVD_Manuscript_data[, assay_component_endpoint_name := # paste("IVD", gsub(", ", "", assay_name), sep="")] # # #tack on the "interesting" part for each assay to create the endpoint_name # IVD_Manuscript_data[assay_component_endpoint_name=="IVDBellwon2015b" | # assay_component_endpoint_name == "IVDBroeders2013" | # assay_component_endpoint_name == "IVDWei2010", # assay_component_endpoint_name := # paste(assay_component_endpoint_name, CellType, sep="_")] # #these two have multiple concentrations or time points # IVD_Manuscript_data[assay_component_endpoint_name=="IVDMundy2004" | # assay_component_endpoint_name=="IVDLundquist2014", # assay_component_endpoint_name := # paste(paste(assay_component_endpoint_name, # paste(ExposureDuration_hrs, "hr", sep = ""), # sep="_"), paste(FBSf, "FBSf", sep = ""), # sep = "_")] # # #create a copy with the assay parameters for the options table (Table S7) # options.IVD_Manuscript_data<-copy(IVD_Manuscript_data) # # ### Pull in httk assay parameter table ### # invitro.assay.params <- copy(invitro.assay.params) # # #copy the experimental data frame to pull out the parameter lines # parameters_only<-copy(IVD_Manuscript_data) # # ### Set up assay parameters to match with invitro.assay.params.RData ### # # #rename columns to align with the names in invitro.assay.params.RData # parameters_only[, timepoint_hr := ExposureDuration_hrs] %>% # .[, media_serum := SerumType] %>% # .[, cell_type := CellType] # # # get the list of columns we need # ivpnames<-names(invitro.assay.params) # # #add on the extra columns to match invitro.assay.params and parameters_only # if(!all(ivpnames%in%names(parameters_only))){ # parameters_only[,ivpnames[!(ivpnames %in% names(parameters_only))]] <- # as.double(NA)} # # #filter parameters_only to match up with ivpnames and only the unique rows # parameters_only<-parameters_only %>% select(as.character(ivpnames)) %>% unique() # # # reassign own parameters table to pass to armitage and kramer functions # user_assay_parameters<-parameters_only # # should have 29 + 1649 = 1678 rows # # #remove parameter columns since we will be pulling from "user_assay_parameters" # ivpnames[2]<-NA #remove endpoint_name bc we need to be able to reference it # IVD_Manuscript_data[,ivpnames[(ivpnames %in% names(IVD_Manuscript_data))]]<- # NULL # ## ----run_IVD_models, eval = execute.vignette---------------------------------- # # # run Armitage model, output concentrations in umol/L (e.g. uM) # armitageOutput_data.dt <- armitage_eval(tcdata=IVD_Manuscript_data, # surface.area.switch = TRUE, # user_assay_parameters = # user_assay_parameters, # restrict.ion.partitioning = TRUE) # # # # # run Kramer model, output concentrations in umol/L (e.g. uM) # kramerOutput_data.dt <- kramer_eval(tcdata=IVD_Manuscript_data, # surface.area.switch = TRUE, # user_assay_parameters = # user_assay_parameters, # restrict.ion.partitioning = TRUE) # ## ----manipulate_for_plotting, eval = execute.vignette------------------------- # # #delete nomconc column # Armitage.dt<-subset(armitageOutput_data.dt, select=-c(nomconc)) # # #pull ioc types from armitage since kramer doesn't generate them # IOC_Type_output<-Armitage.dt[,.(ChemicalName, IOC_Type)] %>% distinct() # # #delete nomconc column and add in the IOC type column for Kramer # Kramer.dt_subset<-subset(kramerOutput_data.dt, select=-c(nomconc)) # Kramer.dt<-Kramer.dt_subset[IOC_Type_output, on = "ChemicalName"] # # #delete cellconc column and add in the IOC type column for Nominal # Nominal.dt_subset<-subset(kramerOutput_data.dt, # select = -c(concentration_cells)) # Nominal.dt<-Nominal.dt_subset[IOC_Type_output, on = "ChemicalName"] # # #add tags # Armitage.dt[,label:= "Armitage"] # Kramer.dt[,label:= "Kramer"] # Nominal.dt[,label:= "Nominal"] # # #combine data for plotting (datasets are side by side) # combodata<-merge(Armitage.dt, kramerOutput_data.dt, by = "X", all.x = TRUE) # #kramer output still has nomconc values # # if(save_output){ # write.csv(combodata, file = paste0(path_out,"ManuscriptOutput_", # data.date,"_",Sys.Date(),".csv")) # } # ## ----ggplot_set_up, eval = execute.vignette----------------------------------- # # #set default theme elements for visually consistent manuscript plots # personal_theme <- function(){ # theme_bw() + # theme(axis.line = element_blank(), # panel.border = element_rect(colour = "black", linewidth=1), # axis.title = element_text(size = rel(1.25)), # axis.text = element_text(size = rel(1.15)), # plot.title = element_text(size = rel(1.5), # hjust = 0.5), # panel.grid = element_blank(), # legend.position = "none") # } # # #set plot colors to correspond with each chemica # plotcolors<-scales::viridis_pal()(length( # unique(armitageOutput_data.dt$ChemicalName))) # names(plotcolors)<-armitageOutput_data.dt$ChemicalName %>% unique() %>% sort() # ## ----joint_intracel_plot_fig2, eval = execute.vignette------------------------ # # #ARMITAGE: # # #Armitage RMSLE: # arm_rmsle_long<-rmsle((Armitage.dt$CellConcentration_uM),(Armitage.dt$ccells)) # arm_rmsle<-round(arm_rmsle_long, 2) # #1.12 # # #Armitage r^2: # ml_arm_overall = lm(CellConcentration_uM~ccells, data = Armitage.dt) # summary(ml_arm_overall)$r.squared # #0.33 # # #plot set up # yrng <- range(Armitage.dt$CellConcentration_uM) # xrng <- range(Armitage.dt$ccells) # n_chems<-n_distinct(Nominal.dt$ChemicalName) # n_data<-nrow(Nominal.dt) # # # now plot # Armitage_ccell<-ggplot(Armitage.dt) + # geom_point(aes(ccells, CellConcentration_uM, color = ChemicalName)) + # geom_abline(intercept=0,slope=1, linetype = "dashed") + # xlab(expression("Armitage Predicted "*italic("C"["cell"])~"(\u03BCM)")) + # ylab(expression("Experimental "*italic(" C"["cell"])~"(\u03BCM)")) + # scale_x_log10(lim = c(0.00001, 400000), breaks = c(10^-4, 10^0, 10^4), # labels = label_log(digits = 2)) + # scale_y_log10(lim = c(0.00001, 400000), breaks = c(10^-4, 10^0, 10^4), # labels = label_log(digits = 2)) + # annotate(geom = "text", x = xrng[2], y = yrng[1], label = parse(text = paste( # "N[data] ", "~'='~", n_data)), hjust = 1, vjust = 0, size = 4) + # annotate(geom = "text", x = xrng[2], y = yrng[1], label = parse(text = paste( # "N[chemicals] ", "~'='~", n_chems)), hjust = 1, vjust = 1.2, size = 4)+ # annotate(geom = "text", x = xrng[2], y = yrng[1], # label = paste("RMSLE =", arm_rmsle), # hjust = 1, vjust = 3, size = 4) + # personal_theme() + # scale_color_manual(values=plotcolors) # # #KRAMER: # # #Kramer RMSLE: # kram_rmsle_long<- # rmsle((Kramer.dt$CellConcentration_uM),(Kramer.dt$concentration_cells)) # kram_rmsle<-format(round(kram_rmsle_long, digits=2), nsmall = 2) # #1.30 # # #Kramer r^2: # ml_kram_overall = lm(CellConcentration_uM~concentration_cells, data = Kramer.dt) # summary(ml_kram_overall)$r.squared # #0.03 # # #plot # Kramer_ccell<-ggplot(Kramer.dt)+ # geom_point(aes(concentration_cells, CellConcentration_uM, # colour = ChemicalName)) + # geom_abline(intercept=0,slope=1, linetype = "dashed") + # scale_x_log10(lim = c(0.00001, 400000), breaks = c(10^-4, 10^0, 10^4), # labels = label_log(digits = 2)) + # scale_y_log10(lim = c(0.00001, 400000), breaks = c(10^-4, 10^0, 10^4), # labels = label_log(digits = 2)) + # xlab(expression("Kramer Predicted "*italic("C"["cell"])~"(\u03BCM)")) + # ylab(expression("Experimental "*italic(" C"["cell"])~"(\u03BCM)")) + # annotate(geom = "text", x = xrng[2], y = yrng[1], label = parse(text = paste( # "N[data] ", "~'='~", n_data)), hjust = 1, vjust = 0, size = 4) + # annotate(geom = "text", x = xrng[2], y = yrng[1], label = parse(text = paste( # "N[chemicals] ", "~'='~", n_chems)), hjust = 1, vjust = 1.2, size = 4)+ # annotate(geom = "text", x = xrng[2], y = yrng[1], # label = paste("RMSLE =", kram_rmsle), hjust = 1, vjust = 3, size = 4) + # personal_theme() + # scale_color_manual(values=plotcolors) # # #NOMINAL: # # #Nominal RMSLE: # nom_rmsle_long<-rmsle((Nominal.dt$CellConcentration_uM),(Nominal.dt$nomconc)) # nom_rmsle<-round(nom_rmsle_long, 2) # #1.45 # # Nominal_ccell<-ggplot(Nominal.dt) + # geom_point(aes(nomconc, CellConcentration_uM, colour = ChemicalName))+ # geom_abline(intercept=0,slope=1, linetype = "dashed") + # scale_x_log10(lim = c(0.00001, 400000), breaks = c(10^-4, 10^0, 10^4), # labels = label_log(digits = 2)) + # scale_y_log10(lim = c(0.00001, 400000), breaks = c(10^-4, 10^0, 10^4), # labels = label_log(digits = 2)) + # xlab(expression(italic("C"["nominal"])~"(\u03BCM)")) + # ylab(expression("Experimental "~italic("C"["cell"])~"(\u03BCM)")) + # labs(colour = "Chemical")+ # annotate(geom = "text", x = xrng[2], y = yrng[1], label = parse(text = paste( # "N[data] ", "~'='~", n_data)), hjust = 1, vjust = 0, size = 4) + # annotate(geom = "text", x = xrng[2], y = yrng[1], label = parse(text = paste( # "N[chemicals] ", "~'='~", n_chems)), hjust = 1, vjust = 1.2, size = 4) + # annotate(geom = "text", x = xrng[2], y = yrng[1], # label = paste("RMSLE =", nom_rmsle), # hjust = 1, vjust = 3, size = 4) + # personal_theme() + # theme(legend.position = "bottom", # legend.key.size = unit(2, 'mm'), # legend.spacing = unit(10, "pt"), # legend.text = element_text(size = 6))+ # guides(col = guide_legend(ncol = 7)) + # scale_color_manual(values=plotcolors) # # #Combine the three plots w RMSLE # Ccells_comboplot<-ggarrange(Nominal_ccell, Armitage_ccell, Kramer_ccell, # ncol=3 , legend = "bottom", common.legend = TRUE) # # #save # if(save_output){ # ggsave(file = paste0("Figure2_",data.date,"_",Sys.Date(),".png"), # path = path_out, # Ccells_comboplot, width = 10, height = 5, dpi = 300) # } # # print(Ccells_comboplot) # ## ----multiobs_intracel_plot_fig3, eval = execute.vignette--------------------- # # #Armitage.multobs.dt %>% count(ChemicalName, NominalDose_uM, Citation) # #some of these have super close nomconcs because they are experimentally # #measured - not really relevant for this comparison because they are replicates # #chemicals to remove: # #Hexachlorobenzene, Malathion, Pentachlorophenol, Propiconazole # #(all of the chemicals from Stadnicka, 2014) # # multiObsChemnameList<-Armitage.dt %>% # dplyr::filter(Citation != "Stadnicka, 2014") %>% # count(ChemicalName, NominalDose_uM) %>% # count(ChemicalName) %>% dplyr::filter(n>=2) %>% # select(ChemicalName) %>% distinct() # #10 chemicals with multiple nominal doses measured # # #filter measured data to just these chemicals # Armitage.multobs.dt_og <- Armitage.dt %>% # dplyr::filter(ChemicalName %in% multiObsChemnameList$ChemicalName) # # Kramer.multobs.dt_og <- Kramer.dt %>% # dplyr::filter(ChemicalName %in% multiObsChemnameList$ChemicalName) # # Nominal.multobs.dt_og <- Nominal.dt %>% # dplyr::filter(ChemicalName %in% multiObsChemnameList$ChemicalName) # #each have 89 observations # # #merge the lines with multiple observations for the same nomconc (ie tox21) # Armitage.multobs.dt<-Armitage.multobs.dt_og %>% # group_by(ChemicalName, NominalDose_uM) %>% # mutate(mean_measuredccell = mean(CellConcentration_uM), # mean_predictedccell = mean(ccells)) %>% # select(ChemicalName, NominalDose_uM, mean_measuredccell, # mean_predictedccell) %>% distinct() # # Kramer.multobs.dt<-Kramer.multobs.dt_og %>% # group_by(ChemicalName, NominalDose_uM) %>% # mutate(mean_measuredccell = mean(CellConcentration_uM), # mean_predictedccell = mean(concentration_cells)) %>% # select(ChemicalName, NominalDose_uM, mean_measuredccell, # mean_predictedccell) %>% distinct() # # Nominal.multobs.dt<-Nominal.multobs.dt_og %>% # group_by(ChemicalName, NominalDose_uM) %>% # mutate(mean_measuredccell = mean(CellConcentration_uM), # mean_predictedccell = mean(NominalDose_uM)) %>% # select(ChemicalName, NominalDose_uM, mean_measuredccell, # mean_predictedccell) %>% distinct() # #down to 50 observations because the duplicates have been averaged # # ### Plotting ### # # #plot set up # yrng_mult <- range(Armitage.multobs.dt$mean_measuredccell) # xrng_mult <- range(Armitage.multobs.dt$mean_predictedccell) # n_chems_mult<-n_distinct(Armitage.multobs.dt$ChemicalName) # n_data_mult<-nrow(Armitage.multobs.dt) # # #ARMITAGE # # #Armitage RMSLE: # arm_mult_rmsle_long<-rmsle((Armitage.multobs.dt$mean_measuredccell), # (Armitage.multobs.dt$mean_predictedccell)) # arm_mult_rmsle<-round(arm_mult_rmsle_long, 2) # #1.27 # # #r squared # ml_arm <- lm(mean_measuredccell~mean_predictedccell, data = Armitage.multobs.dt) # rsq_armitage_mult <- summary(ml_arm)$r.squared # #0.87 # # # now plot # Armitage_ccell_mult<-ggplot(Armitage.multobs.dt) + # geom_point(aes(mean_predictedccell, mean_measuredccell, # colour = ChemicalName)) + # geom_smooth(aes(mean_predictedccell, mean_measuredccell, # colour = ChemicalName), method=lm, se = FALSE) + # geom_abline(intercept=0,slope=1, linetype = "dashed") + # scale_x_log10(lim = c(0.0001, 400000), breaks = c(10^-4, 10^0, 10^4), # labels = label_log(digits = 2)) + # scale_y_log10(lim = c(0.0001, 400000), breaks = c(10^-4, 10^0, 10^4), # labels = label_log(digits = 2)) + # xlab(expression("Armitage Predicted "*italic("C"["cell"])~"(\u03BCM)")) + # ylab(expression("Experimental "*italic(" C"["cell"])~"(\u03BCM)")) + # labs(colour = "Chemical")+ # annotate(geom = "text", x = xrng_mult[2], y = yrng_mult[1], # label = parse(text = paste("N[data] ", "~'='~", n_data_mult)), # hjust = 1, vjust = 0, size = 4) + # annotate(geom = "text", x = xrng_mult[2], y = yrng_mult[1], # label = parse(text = paste("N[chemicals] ", "~'='~", n_chems_mult)), # hjust = 1, vjust = 1.2, size = 4)+ # annotate(geom = "text", x = xrng_mult[2], y = yrng_mult[1], # label = paste("RMSLE =", arm_mult_rmsle), hjust = 1, vjust = 3, size = 4) + # personal_theme() + # scale_color_manual(values=plotcolors) # # #KRAMER # # #Kramer RMSLE: # kram_mult_rmsle_long<-rmsle((Kramer.multobs.dt$mean_measuredccell), # (Kramer.multobs.dt$mean_predictedccell)) # kram_mult_rmsle<-round(kram_mult_rmsle_long, 2) # # 1.58 # # #r squared # ml_kram <- lm(mean_measuredccell~mean_predictedccell, data = Kramer.multobs.dt) # rsq_kramer_mult <- summary(ml_kram)$r.squared # #0.87 # # # now plot # Kramer_ccell_mult<-ggplot(Kramer.multobs.dt) + # geom_point(aes(mean_predictedccell, mean_measuredccell, # colour = ChemicalName)) + # geom_smooth(aes(mean_predictedccell, mean_measuredccell, # colour = ChemicalName), method=lm, se = FALSE) + # geom_abline(intercept=0,slope=1, linetype = "dashed") + # scale_x_log10(lim = c(0.0001, 400000), breaks = c(10^-4, 10^0, 10^4), # labels = label_log(digits = 2)) + # scale_y_log10(lim = c(0.0001, 400000), breaks = c(10^-4, 10^0, 10^4), # labels = label_log(digits = 2)) + # xlab(expression("Kramer Predicted "*italic("C"["cell"])~"(\u03BCM)")) + # ylab(expression("Experimental "*italic(" C"["cell"])~"(\u03BCM)")) + # labs(colour = "Chemical")+ # annotate(geom = "text", x = xrng_mult[2], y = yrng_mult[1], # label = parse(text = paste("N[data] ", "~'='~", n_data_mult)), # hjust = 1, vjust = 0, size = 4) + # annotate(geom = "text", x = xrng_mult[2], y = yrng_mult[1], # label = parse(text = paste("N[chemicals] ", "~'='~", n_chems_mult)), # hjust = 1, vjust = 1.2, size = 4) + # annotate(geom = "text", x = xrng_mult[2], y = yrng_mult[1], # label = paste("RMSLE =", kram_mult_rmsle), hjust = 1, vjust = 3, size = 4) + # personal_theme() + # scale_color_manual(values=plotcolors) + # theme(legend.key.size = unit(2, 'mm'), # legend.text = element_text(size = 6)) # # # #arrange # Ccells_comboplot_mult<-ggarrange(Armitage_ccell_mult, Kramer_ccell_mult, ncol=2, # legend = "bottom", common.legend = TRUE) # # if(save_output){ # ggsave(file = paste0("Figure3_",data.date,"_",Sys.Date(),".png"), # path = path_out, # Ccells_comboplot_mult, width = 8, height = 4, dpi = 300) # } # # print(Ccells_comboplot_mult) # ## ----chemicalPart_intracel_plot_fig4, eval = execute.vignette----------------- # # #Cwat_uM -------- Cwat_R # Cwater_comparison<-ggplot(combodata) + # geom_point(aes(concentration_medium, cwat_s, colour = ChemicalName.x)) + # geom_abline(intercept=0,slope=1, linetype = "dashed") + # ggtitle("Medium Concentration") + # xlab(expression("Kramer Predicted "*italic("C"["medium"])~"(\u03BCM)")) + # ylab(expression("Armitage Predicted "*italic(" C"["medium"])~"(\u03BCM)")) + # scale_x_log10(labels = label_log(digits = 2)) + # scale_y_log10(labels = label_log(digits = 2)) + # personal_theme() + # scale_color_manual(values=plotcolors) # # #Cair_uM -------- Cair_R # Cair_comparison<-ggplot(combodata) + # geom_point(aes(concentration_air, cair, colour = ChemicalName.x)) + # geom_abline(intercept=0,slope=1, linetype = "dashed") + # ggtitle("Air Concentration") + # xlab(expression("Kramer Predicted "*italic("C"["air"])~"(\u03BCM)")) + # ylab(expression("Armitage Predicted "*italic(" C"["air"])~"(\u03BCM)")) + # scale_x_log10(labels = label_log(digits = 2)) + # scale_y_log10(labels = label_log(digits = 2)) + # personal_theme() + # scale_color_manual(values=plotcolors) # # #C_cells_uM ----- Ccells # Ccells_comparison<-ggplot(combodata) + # geom_point(aes(concentration_cells, ccells, colour = ChemicalName.x)) + # geom_abline(intercept=0,slope=1, linetype = "dashed") + # ggtitle("Cell Concentration") + # xlab(expression("Kramer Predicted "*italic("C"["cell"])~"(\u03BCM)")) + # ylab(expression("Armitage Predicted "*italic(" C"["cell"])~"(\u03BCM)")) + # scale_x_log10(lim = c(10^-4, 10^5), labels = label_log(digits = 2), # breaks = c(10^-2, 10^0, 10^2, 10^4)) + # scale_y_log10(lim = c(10^-4, 10^5), labels = label_log(digits = 2), # breaks = c(10^-2, 10^0, 10^2, 10^4)) + # personal_theme() + # scale_color_manual(values=plotcolors) # # #Aplastic_uM_m2 - (Cplastic_R / Sarea_R) # Aplastic_comparison<-ggplot(combodata) + # geom_point(aes(concentration_plastic, cplastic, colour = ChemicalName.x)) + # geom_abline(intercept=0,slope=1, linetype = "dashed") + # ggtitle("Plastic Concentration") + # xlab(expression("Kramer Predicted "*italic("C"["plastic"])~"(\u03BCM)")) + # ylab(expression("Armitage Predicted "*italic(" C"["plastic"])~"(\u03BCM)")) + # labs(colour = "Chemical")+ # scale_x_log10(labels = label_log(digits = 2), # breaks = c(10^-6, 10^-4, 10^-2, 10^0)) + # scale_y_log10(labels = label_log(digits = 2)) + # personal_theme() + # scale_color_manual(values=plotcolors) + # theme(legend.key.size = unit(1, 'mm'), # legend.spacing = unit(2, "pt"), # legend.text = element_text(size = 6)) # # # compartment_plots<-ggarrange(Ccells_comparison, Aplastic_comparison, # Cwater_comparison, Cair_comparison, # common.legend = TRUE, legend = "none") # # if(save_output){ # ggsave(file = paste0("Figure4_",data.date,"_",Sys.Date(),".png"), # path = path_out, # compartment_plots, width = 8, height = 6, dpi= 300) # } # # print(compartment_plots) # ## ----curated_values_plot_fig5, eval = execute.vignette------------------------ # # #load curated data info # original_curated.data.dt <- httk::Dimitrijevic.IVD # # #first - run using the default/dashboard values # default.dt<-copy(original_curated.data.dt %>% # select(-c( "Arnot_pka", "Arnot_pkb", "Chemaxon_pKa", # "Chemaxon_pKb", "log.KOW.N", "log.KAW.N", # "Predicted_Ccell_µM", "FoA"))) # # # run the model, output concentrations in umol/L (e.g. uM) # armitageOutput_default.dt <- armitage_eval(tcdata=default.dt, # restrict.ion.partitioning = TRUE, # surface.area.switch = FALSE) # # #next- run using curated values (EAS-E Suite) # curated.data.dt<-copy(original_curated.data.dt) # # #overwrite using curated physchem properties # curated.data.dt[, pKa_Donor := as.character(Arnot_pka)] %>% #acidic # .[, pKa_Accept := as.character(Arnot_pkb)] %>% #basic # .[, gkow := log.KOW.N] %>% # .[, gkaw_n := log.KAW.N] # # # run the model, output concentrations in umol/L (e.g. uM) # armitageOutput_curated.data.dt <- armitage_eval(tcdata=curated.data.dt, # restrict.ion.partitioning = TRUE, # surface.area.switch = FALSE) # # #save output with these values (Supplemental Materials T9) # defaultvals<- armitageOutput_default.dt %>% # select(Name, casrn, gkow_n, gkaw_n, pKa_Donor, pKa_Accept) # curatedvals<-armitageOutput_curated.data.dt %>% # select(Name, casrn, gkow_n, gkaw_n, pKa_Donor, pKa_Accept) # #tie the two tables together for the supplement # s9_table<-rbind(curatedvals, defaultvals) # # #calculate rmsle # # #default values # # #individual # armitageOutput_default.dt[,rmsle_ccell_default:= # rmsle(Reported_Ccell_µM, ccells), by = Name] # # #total # RMSLE_default.dt<-rmsle(armitageOutput_default.dt$Reported_Ccell_µM, # armitageOutput_default.dt$ccells) # RMSLE_default.dt<-round(RMSLE_default.dt, 3) # #0.586 # # #curated values # # #individual # armitageOutput_curated.data.dt[,rmsle_ccell_curated:= # rmsle(Reported_Ccell_µM, ccells), by = Name] # # #total # RMSLE_curated.data.dt<-rmsle(armitageOutput_curated.data.dt$Reported_Ccell_µM, # armitageOutput_curated.data.dt$ccells) # RMSLE_curated.data.dt<-round(RMSLE_curated.data.dt, 3) # #0.570 # # #improvement when using curated data: # RMSLE_default.dt-RMSLE_curated.data.dt # #0.016 # # ## difference between the two on a per-chemical basis ## # curated_data.dt<-merge(armitageOutput_default.dt, # armitageOutput_curated.data.dt, by = "Name") # # #calculate rmsle difference for each chemical # curated_data.dt[,rmsle_difference:= # rmsle_ccell_curated-rmsle_ccell_default, by = Name] # # #which chemicals had the greatest magnitude of difference # #largest improvement when using curated data # head(curated_data.dt %>% # select(rmsle_difference, Name) %>% # arrange(rmsle_difference), 2) # # #largest improvement when using default data # tail(curated_data.dt %>% # select(rmsle_difference, Name) %>% # arrange(rmsle_difference), 2) # # #plot # curated_difference_httk_plot<-ggplot(curated_data.dt) + # geom_point(aes(ccells.y, Reported_Ccell_µM.y, colour = Name), size = 4) + # geom_label_repel(aes(ccells.y, Reported_Ccell_µM.y, # label = round(rmsle_difference, 2))) + # geom_abline(intercept=0,slope=1, linetype = "dashed") + # xlab(expression("Armitage Predicted "*italic("C"["cell"])~"(\u03BCM)")) + # ylab(expression("Experimental "*italic(" C"["cell"])~"(\u03BCM)")) + # labs(color = "Chemical")+ # scale_x_log10(lim = c(10, 10000), labels = label_log(digits = 2)) + # scale_y_log10(lim = c(10, 10000), labels = label_log(digits = 2)) + # personal_theme() + # scale_color_viridis(discrete=TRUE)+ # theme(legend.position = "bottom", # legend.key.size = unit(2, 'mm'), # legend.spacing = unit(10, "pt"), # legend.text = element_text(size = 10)) # # if(save_output){ # ggsave(file = paste0("Figure5_",data.date,"_",Sys.Date(),".png"), # path = path_out, # curated_difference_httk_plot, width = 8, height = 5, dpi = 300) # # write.csv(s9_table, file = paste0(path_out,"table_s9_",data.date,"_", # Sys.Date(),".csv")) # } # # ## ----fold_change_calculation, eval = execute.vignette------------------------- # # # The Armitage model resulted in a slight over prediction of Ccell whereas # # the Kramer model had the greatest spread but most evenly over and # # underpredicted Ccell (Figure 2). # range(Armitage.dt$ccells) # #1.555754e-02 4.845393e+04 # range(Kramer.dt$concentration_cells) # #5.838149e-04 2.663947e+04 # # #Kramer has a larger spread of predictions # # ##create datatable with both outputs stacked # stacked.data<-rbind(Armitage.dt, Kramer.dt, fill=TRUE) # # #overpredict or underpredict for each model compared to experimental # stacked.data_HL<-copy(stacked.data) # stacked.data_HL[ccells>CellConcentration_uM, # ARMITAGE_ccellcomp:="ARMITAGE higher"] %>% # .[ccells% # .[concentration_cells>CellConcentration_uM, # KRAMER_ccellcomp:="KRAMER higher"] %>% # .[concentration_cells% count(ARMITAGE_ccellcomp, KRAMER_ccellcomp) # #Armitage model more often over predicts # #Kramer is pretty even split of over and under predictions # # #pick out where armitage prediction or kramer prediction is larger # sidebyside.data_HL<-copy(combodata) # sidebyside.data_HL[ccells>concentration_cells, ARMITAGE:="ARMITAGE higher"] %>% # .[ccells% select(ChemicalName.x, ARMITAGE, KRAMER) %>% # distinct() %>% count(ARMITAGE, KRAMER) # #armitage: higher for 38 chemicals # #kramer: higher for 5 chemicals # #38+5 = 43 chemicals overall # ## ----experimental_nomconc, eval = execute.vignette---------------------------- # # #pull observations that did experimentally measure the nominal dose # exp_nom.dt<-copy(IVD_Manuscript_data[!(is.na(ExperimentalDose_uM)),]) # #6 citations did measure experimental dose # #42 observations # # #run armitage with nomconc as NominalDose_uM.x (not experimentally measured) # arm_nomconc<-copy(exp_nom.dt[, nomconc := NominalDose_uM]) # # # run Armitage model, output concentrations in umol/L (e.g. uM) # armitageOutput_arm_nomconc <- armitage_eval(tcdata=arm_nomconc, # surface.area.switch = TRUE, # user_assay_parameters = # user_assay_parameters, # restrict.ion.partitioning = TRUE) # # #run armitage with nomconc as ExperimentalDose_uM.x (experimentally measured) # arm_expconc<-copy(exp_nom.dt[, nomconc := ExperimentalDose_uM]) # # armitageOutput_arm_expconc <- armitage_eval(tcdata=arm_expconc, # surface.area.switch = TRUE, # user_assay_parameters = # user_assay_parameters, # restrict.ion.partitioning = TRUE) # # #Armitage NOMCONC RMSLE: # arm_nomconc_rmsle_long<-rmsle((armitageOutput_arm_nomconc$CellConcentration_uM), # (armitageOutput_arm_nomconc$ccells)) # arm_nomconc_rmsle<-round(arm_nomconc_rmsle_long, 2) # #1.62 # # #Armitage EXPCONC RMSLE: # arm_expconc_rmsle_long<-rmsle((armitageOutput_arm_expconc$CellConcentration_uM), # (armitageOutput_arm_expconc$ccells)) # arm_expconc_rmsle<-round(arm_expconc_rmsle_long, 2) # #1.53 # # #improvement between the experimental and nomconc rmsles # arm_nomconc_rmsle- arm_expconc_rmsle # #0.09 # ## ----physchem_data_tableS1, eval = execute.vignette--------------------------- # # #create datatable to hold the info # table1.dt<-copy(Nominal.dt) # # # pull out the columns we need # data.dt<-table1.dt[, c("Citation", "casrn", "ChemicalName", "IOC_Type", # "nomconc", "CellConcentration_uM"), with = FALSE] # # # # pull dtxsids, log P, henry's law, and pka: # # #dtxsids # cheminfo.dt<-as.data.table(get_chem_id(chem.cas=unique(data.dt$casrn))) # # cheminfo.dt[, casrn:=chem.cas] %>% #rename # .[,chem.cas:=NULL] %>% #delete # .[,chem.name:=NULL] #delete # # #add dtxsid info to the storage file # data.dt2<-data.dt[cheminfo.dt,on="casrn"] # # #collect log P, henry's law, and pka values # data.dt2[, c("logP","logHenry", "pKa_Donor", "pKa_Accept") := # as.data.frame(get_physchem_param(param = c("logP","logHenry", # "pKa_Donor", # "pKa_Accept"), # chem.cas = casrn))] # # # # calculate median CellConcentration_uM:cnom ratio # # (for each citation if multiple) and number of observations per chemical # # #get ccell : cnom ratio # data.dt2[,cellconc_cnom_ratio:=CellConcentration_uM/nomconc] # # #for each chemical/citation combination, calculate the median ccell:cnom ratio # data.dt2[, median_ratio := median(cellconc_cnom_ratio), # by=list(Citation, ChemicalName)] # # #count the number of observations by citation and chemname # data.dt2[, number_of_obs:=.N, by=list(Citation, ChemicalName)] # # #filter down to one observation per citation/casrn combo # final_output<-data.dt2 %>% # mutate(logP=round(logP, 2), # logHenry=round(logHenry, 2), # median_ratio=round(median_ratio, 3)) %>% # select(-c("casrn", "nomconc", "CellConcentration_uM", # "cellconc_cnom_ratio")) %>% #get rid of extra rows # distinct(Citation, ChemicalName, .keep_all = TRUE) # #1 obs per citation/chem pair # # #save the output # if(save_output){ # write.csv(final_output, file = paste0(path_out,"table_s1_",data.date,"_", # Sys.Date(),".csv")) # } # ## ----assay_descriptors_tableS2, eval = execute.vignette----------------------- # # #copy table to manipulate # parameter_table<-copy(Nominal.dt) # # table_s2<-parameter_table[, c("Citation", "ChemicalName", "nomconc", # "NominalDose_uM", "v_total", "v_working", # "sarea", "well_number", "cell_yield", "FBSf", # "FigureNumber", "CellType")] #subset columns # # #save the output # if(save_output){ # write.csv(table_s2, # file = paste0(path_out,"table_s2_raw_",data.date,"_", # Sys.Date(),".csv")) # } # ## ----multobs_RMSLE_tableS3, eval = execute.vignette--------------------------- # # #Armitage # # #copy for editing # multobs_rmsle_armitage<-copy(as.data.table(Armitage.multobs.dt)) # # #group by chemname, calc RMSLE, r^2, and number of obs # multobs_rmsle_armitage[, RMSLE_armitage_mult:= # rmsle(mean_measuredccell, mean_predictedccell), # by = "ChemicalName"] %>% #rmsle # .[, rsq_armitage_mult:= # summary(lm(mean_measuredccell~mean_predictedccell))$r.squared, # by = "ChemicalName"] %>% #r^2 # .[, numberofpoints:= .N, by = "ChemicalName"] # # #numberofpoints = number of datapoints considered (includes averaged points) # # #repeat for kramer # # #copy for editing # multobs_rmsle_kramer<-copy(as.data.table(Kramer.multobs.dt)) # # #group by chemname, calc RMSLE, r^2, and number of obs # multobs_rmsle_kramer[, RMSLE_kramer_mult:= # rmsle(mean_measuredccell, mean_predictedccell), # by = "ChemicalName"] %>% #rmsle # .[, rsq_kramer_mult:= # summary(lm(mean_measuredccell~mean_predictedccell))$r.squared, # by = "ChemicalName"] %>% #r^2 # .[, numberofpoints:= .N, by = "ChemicalName"] # # # #cut down both to the columns we need and only 1 obs per chemical # mult.dt<-merge(multobs_rmsle_armitage[, -c("NominalDose_uM", # "mean_measuredccell", # "mean_predictedccell")] %>% # distinct() , # multobs_rmsle_kramer[, -c("NominalDose_uM", # "mean_measuredccell", # "mean_predictedccell")] %>% # distinct()) # # #for higher/lower rmsle comparison # mult_HL<-copy(mult.dt[RMSLE_kramer_mult>RMSLE_armitage_mult, # higherorlower:="Kramer higher"] %>% # .[RMSLE_kramer_mult% count(higherorlower) # #armitage predictions have lower error for 6 chemicals # #kramer predictions are lower error for the other 4 # # #save the output # if(save_output){ # write.csv(mult.dt, # file = paste0(path_out,"table_s3_",data.date,"_",Sys.Date(),".csv")) # } # # # ## ----ionization_RMSLE_tableS4, eval = execute.vignette------------------------ # # #copy for editing # ionization_rmsle<-copy(stacked.data) # # #calculate RMSLSE and r^2 for each model and ionization type # ionization_rmsle_output<- copy(ionization_rmsle[label=="Armitage", # RMSLE_armitage_ionization:= # rmsle(CellConcentration_uM, ccells), # by = "IOC_Type"] %>% #RMSLE # .[label=="Armitage", rsq_armitage_ionization:= # summary(lm(CellConcentration_uM~ccells))$r.squared, # by = "IOC_Type"] %>% #r^2 # .[label=="Kramer", RMSLE_kramer_ionization:= # rmsle(CellConcentration_uM, concentration_cells), # by = "IOC_Type"] %>% #RMSLE # .[label=="Kramer", rsq_kramer_ionization:= # summary(lm(CellConcentration_uM~concentration_cells))$r.squared, # by = "IOC_Type"] %>% #r^2 # .[,numberofpoints:= .N, by = "IOC_Type"] %>% # .[, c("IOC_Type", "RMSLE_armitage_ionization", # "rsq_armitage_ionization", "RMSLE_kramer_ionization", # "rsq_kramer_ionization", "numberofpoints")] %>% #select columns we want # distinct()) # # print(ionization_rmsle_output) # #number of points is overall, not for each calculation # # #save the output # if(save_output){ # write.csv(ionization_rmsle_output, # file = paste0(path_out,"table_s4_",data.date,"_",Sys.Date(),".csv")) # } # ## ----logP_RMSLE_tableS5, eval = execute.vignette------------------------------ # # #cut up into chunks # gkowbins<-stacked.data %>% mutate(LogP_bins = # cut(gkow, breaks = c(-Inf, 1.30, 3.36, Inf), # labels = c("Hydrophilic", "Moderate Lipophilic", "Lipophilic"))) # # #calculate rmsle for each model and logP bin combo # gkowbins_rmsle_output<-suppressWarnings(copy( # gkowbins[label == "Armitage", Armitage_RMSLE := rmsle( # CellConcentration_uM, ccells), by = LogP_bins] %>% # .[label == "Kramer", Kramer_RMSLE := rmsle( # CellConcentration_uM, concentration_cells), by = LogP_bins] %>% # .[, numberofpoints := .N, by = list(LogP_bins, label)] %>% # .[, c("LogP_bins", "Armitage_RMSLE", "Kramer_RMSLE", # "numberofpoints", "label")] %>% # distinct())) # # print(gkowbins_rmsle_output) # #number of points is overall, not for each calculation # # #save the output # if(save_output){ # write.csv(gkowbins_rmsle_output, # file = paste0(path_out,"table_s5_",data.date,"_",Sys.Date(),".csv")) # } # ## ----kaw_RMSLE_tableS6, eval = execute.vignette------------------------------- # # #source for chunks is table 2 Volatility class in water: # #https://www.epa.gov/pesticide-science-and-assessing-pesticide-risks/ # #guidance-reporting-environmental-fate-and-transport#:~:text=is%20very%20low.-, # #Volatility%20from%20Water,the%20reciprocal%20of%20KAW.&text=GMW%20=%20gram%20 # #molecular%20weight%20in%20g/mole.&text=This%20value%20is%20the%20inverse,as%20 # #Cwater/Cair.&text=Note%20that%20all%20chemicals%20may,volatility%20potential%20 # #is%20very%20low. # # #cut up kaw into relevant bins # kawbins<-stacked.data %>% # mutate(Kaw_bins = cut(DR_kaw, breaks = c(Inf, 10^-2, 10^-3, 10^-5, -Inf), # labels = c("Rapidly lost from a water surface", # "Volatile from a water surface", # "Slightly volatile from a water surface", "Non-volatile"))) # # #calculate rmsle for each model and logP bin combo # kaw_rmsle_output<-suppressWarnings(kawbins[label == "Armitage", # Armitage_RMSLE := rmsle( # CellConcentration_uM, ccells), by = Kaw_bins] %>% # .[label == "Kramer", Kramer_RMSLE := rmsle( # CellConcentration_uM, concentration_cells), by = Kaw_bins] %>% # .[, numberofpoints := .N, by = list(Kaw_bins, label)] %>% # .[, c("Kaw_bins", "Armitage_RMSLE", "Kramer_RMSLE", # "numberofpoints", "label")] %>% # distinct()) # # # #save the output # if(save_output){ # write.csv(kaw_rmsle_output, # file = paste0(path_out,"table_s6_",data.date,"_",Sys.Date(),".csv")) # } # # ## ----options_tableS7, eval = execute.vignette--------------------------------- # # # #list of options that can be changed in the armitage model: # # # - this.option.kbsa2 Use alternative bovine-serum-albumin partitioning # # - this.option.swat2 Use alternative water solubility correction # # - this.option.kpl2 Use alternative plastic-water partitioning model # # - option.plastic == TRUE (default) gives nonzero surface area # # - option.bottom == TRUE (default) includes surface area of the bottom of # # the well in determining sarea # # - restrict.ion.partitioning # # § False (default): Treat entire chemical as neutral and allow it all # # to partition normally as such # # § True: Restrict so only the neutral fraction of the chemical is # # available to partition normally but the ionized portion is scaled # # # ##filter to just the ones we can run with the surface area calculation since # #option.plastic and option.bottom require that function: # options.IVD_Manuscript_data_2<- # options.IVD_Manuscript_data[!is.na(well_number), ] %>% # .[,sarea:=NULL] #delete sarea so the surface area function will run # #126 out of the original 153 observations # # length(options.IVD_Manuscript_data$Citation %>% unique()) #14 citations # # length(options.IVD_Manuscript_data_2$Citation %>% unique()) #12 citations # # 153-126 data points; 27 data points total removed # # # create table to append the results to # options.table <- NULL # # # run basic armitage, only setting the surface area switch # option.basic <- armitage_eval(tcdata=options.IVD_Manuscript_data_2, # surface.area.switch = TRUE, # user_assay_parameters = user_assay_parameters) # # option.basic_RMSLE<-rmsle((option.basic$CellConcentration_uM), # (option.basic$ccells)) # #1.148 # # #append to table # options.table<-rbind(options.table, c(round(option.basic_RMSLE ,3), # "All Default")) # # # run with this.option.kbsa2 == TRUE (default is false) # option.kbsa <- armitage_eval(tcdata=options.IVD_Manuscript_data_2, # this.option.kbsa2 = TRUE, # surface.area.switch = TRUE, # user_assay_parameters = user_assay_parameters) # # option.kbsa_RMSLE<-rmsle((option.kbsa$CellConcentration_uM), # (option.kbsa$ccells)) # #1.150 # # #append to table # options.table<-rbind(options.table, c(round(option.kbsa_RMSLE, 3), # "option.kbsa=TRUE")) # # # # run with this.option.swat2 == TRUE (default is false) # option.swat <- armitage_eval(tcdata=options.IVD_Manuscript_data_2, # this.option.swat2 = TRUE, # surface.area.switch = TRUE, # user_assay_parameters = user_assay_parameters) # # option.swat_RMSLE<-rmsle((option.swat$CellConcentration_uM), # (option.swat$ccells)) # #3.585 # # #append to table # options.table<-rbind(options.table, c(round(option.swat_RMSLE, 3), # "option.swat2=TRUE")) # # # run with this.option.kpl2 == TRUE (Fischer) (default is false/kramer) # option.kpl <- armitage_eval(tcdata=options.IVD_Manuscript_data_2, # this.option.kpl2 = TRUE, # surface.area.switch = TRUE, # user_assay_parameters = user_assay_parameters) # # option.kpl_RMSLE<-rmsle((option.kpl$CellConcentration_uM),(option.kpl$ccells)) # #1.101 # # #append to table # options.table<-rbind(options.table, c(round(option.kpl_RMSLE, 3), # "option.kpl2=TRUE")) # # # run with option.plastic == FALSE (default is true) # opt.plast_input<- copy(options.IVD_Manuscript_data_2) # opt.plast_input[,option.plastic := FALSE] # # option.plast <- armitage_eval(tcdata=opt.plast_input, # surface.area.switch = TRUE, # user_assay_parameters = user_assay_parameters) # # option.plast_RMSLE<-rmsle((option.plast$CellConcentration_uM), # (option.plast$ccells)) # #1.178 # # #append to table # options.table<-rbind(options.table, c(round(option.plast_RMSLE, 3), # "option.plastic=FALSE")) # # # run with option.bottom == FALSE (default is true) # option.bottom_input<- copy(options.IVD_Manuscript_data_2) # option.bottom_input[,option.bottom := FALSE] # # option.bottom <- armitage_eval(tcdata=option.bottom_input, # surface.area.switch = TRUE, # user_assay_parameters = user_assay_parameters) # # option.bottom_RMSLE<-rmsle((option.bottom$CellConcentration_uM), # (option.bottom$ccells)) # #3.956 # # #append to table # options.table<-rbind(options.table, c(round(option.bottom_RMSLE,3), # "option.bottom=FALSE")) # # # run with restrict.ion.partitioning == TRUE (default is FALSE) # option.restion <- armitage_eval(tcdata=options.IVD_Manuscript_data_2, # restrict.ion.partitioning = TRUE, # surface.area.switch = TRUE, # user_assay_parameters = user_assay_parameters) # # option.restion_RMSLE<-rmsle((option.restion$CellConcentration_uM), # (option.restion$ccells)) # #1.152 # # #append to table # options.table<-rbind(options.table, c(round(option.restion_RMSLE,3), # "restrict.ion.partitioning = TRUE")) # # #save the output # if(save_output){ # write.csv(options.table, # file = paste0(path_out,"table_s7_",data.date,"_",Sys.Date(),".csv")) # } # ## ----intracellular_histogram_plotS1, eval = execute.vignette------------------ # # #calculate RMSLE for each row of data # combodata[, RMSLE_armitage_by_chem:= rmsle(CellConcentration_uM.x, ccells), # by = "ChemicalName.x"] # combodata[, RMSLE_kramer_by_chem:= # rmsle(CellConcentration_uM.x, concentration_cells), # by = "ChemicalName.x"] # combodata[, RMSLE_nominal_by_chem:= rmsle(CellConcentration_uM.x, nomconc), # by = "ChemicalName.x"] # # #cut down to one observation per chemical # histo.data<-copy(combodata %>% # select(ChemicalName.x, RMSLE_armitage_by_chem, # RMSLE_kramer_by_chem, RMSLE_nominal_by_chem) %>% # unique()) # # #plot histograms # nominal_histogram_by_chem <-ggplot(histo.data, aes(x=RMSLE_nominal_by_chem)) + # geom_histogram(binwidth = 0.5, boundary = 0) + # labs(x = "Nominal RMSLE", # y = "Number of Observations") + # xlim(0, 5) + ylim(0, 17) + # theme_classic() + # theme(axis.text=element_text(size=12), # axis.title=element_text(size=14)) # # armitage_histogram_by_chem <-ggplot(histo.data, aes(x=RMSLE_armitage_by_chem)) + # geom_histogram(binwidth = 0.5, boundary = 0) + # labs(x = "Armitage RMSLE", # y = "Number of Observations") + # xlim(0, 5) + ylim(0, 17) + # theme_classic() + # theme(axis.text=element_text(size=12), # axis.title=element_text(size=14)) # # kramer_histogram_by_chem <- ggplot(histo.data, aes(x=RMSLE_kramer_by_chem)) + # geom_histogram(binwidth = 0.5, boundary = 0) + # labs(x = "Kramer RMSLE", # y = "Number of Observations") + # xlim(0, 5) + ylim(0, 17) + # theme_classic()+ # theme(axis.text=element_text(size=12), # axis.title=element_text(size=14)) # # histogram_comboplot_by_chem<-ggarrange(nominal_histogram_by_chem + # theme(axis.ticks.y = element_blank(), # plot.margin = margin(r = 1)), # armitage_histogram_by_chem + # theme(axis.text.y = element_blank(), # axis.ticks.y = element_blank(), # axis.title.y = element_blank(), # plot.margin = margin(r = 1, # l = 1)), # kramer_histogram_by_chem + # theme(axis.text.y = element_blank(), # axis.ticks.y = element_blank(), # axis.title.y = element_blank(), # plot.margin = margin(r = 1, # l = 1)), # ncol=3) # # if(save_output){ # ggsave(file = paste0("Figure_S1_by_chem_",data.date,"_",Sys.Date(),".png"), # path = path_out, # histogram_comboplot_by_chem, width = 10, height = 5, dpi = 300) # } # # print(histogram_comboplot_by_chem) #