# Author: Shirley Pu
# Code used for the paper "Examining Environmental Matrix Effects on Quantitative Non-Targeted Analysis Estimates of Per- and Polyfluoroalkyl Substances"
# Contains code used for data analysis

# Load libraries
library(tidyverse)
library(investr)
library(ggridges)
library(ggpattern)
library(xgboost)
library(caret)
library(lme4)
library(lmerTest)
library(scales)

# Source custom functions
source("Pu_et_al_Env_Matrices_functions.R")

# Read in data
load("Matrices_Data_cleaned_final.RData")

# Setup of variables
uniq_media <- unique(Matrices_Data_cleaned %>% pull('Media'))
uniq_chem_media <- unique(Matrices_Data_cleaned %>% select(Chem,Media))

#### QUANTITATION ####
# Using just reconstitution method 2 data for Solvent, DW, and Sludge
Matrices_Data_cleaned_RM2 <- Matrices_Data_cleaned %>% 
  filter(Media %in% c("DW","Sludge","Solvent")) %>%
  mutate(Media = factor(Media,levels=c("Solvent","DW","Sludge")),
         Chem = gsub("FTS42","4:2 FTS",Chem) %>%
           gsub("FTS62","6:2 FTS",.) %>%
           gsub("FTS82","8:2 FTS",.))

uniq_media_RM2 <- unique(Matrices_Data_cleaned_RM2 %>% pull('Media') %>% as.character())

## Figure 1 ##
## Response factor distributions (Solvent, DW, Sludge RM2 only)
# Using joint bandwidth of 0.05 to alleviate multimodality in Sludge distributions

ggplot(Matrices_Data_cleaned_RM2 %>%
         rename(Matrix=Media) %>%
         group_by(Chem,Matrix) %>%
         dplyr::left_join(count(.)) %>%
         ungroup(),aes(x=RF_uM,y=reorder(Chem,-RT),color=Matrix,fill=Matrix))+
  scale_color_manual(values=c("tan","limegreen","purple"))+
  scale_fill_manual(values=c("tan","limegreen","purple"))+
  xlab(expression(paste("Response Factor (", mu, "M"^"-1",")")))+
  ylab("Chemical")+
  geom_density_ridges(alpha=.5,bandwidth=0.05)+
  scale_x_log10(breaks = scales::trans_breaks("log10", function(x) 10^x),
                labels = scales::trans_format("log10", scales::math_format(10^.x)))+
  #geom_text(aes(label=paste0("n=",n),group=Chem),hjust=0.5)+
  theme_bw(base_size=20)

### INVERSE PREDICTIONS ###
# Nested loop (loop over Media, loop over Chem in that Media)
invpred_results <- data.frame()
for(i in uniq_media_RM2){
  uniq_media_rows <- Matrices_Data_cleaned_RM2 %>% filter(Media==i)
  uniq_media_chem <- unique(uniq_media_rows %>% pull('Chem'))
  for(j in uniq_media_chem){
    invpred_results <- rbind(invpred_results,get_invest_summarize(j,i,Matrices_Data_cleaned_RM2,leaveout=TRUE,logspace=FALSE,more_cols=TRUE))
  }
}

### qNTA ###
## EXPERT-SELECTED SURROGATES BOOTSTRAP ##
# Expert surrogate list 
# Formatted as four columns; 1st column- Compound, 2nd column- Surrogate_1, 3rd column- Surrogate_2, 4th column- Surrogate_3
# Surrogates are selected based on homologous series of the chemical, and nearest eluting compound with the same ionizing group
expert_surrogates<-read.csv("Pu_et_al_Env_Matrices_Surrogate_Selection.csv")

# Loop over media
start_time <- Sys.time()
expert_boot_out <- data.frame()
set.seed(1)
for(i in uniq_media_RM2){
  curr_expert_boot <- Matrices_Data_cleaned_RM2[Matrices_Data_cleaned_RM2$Media==i,]
  curr_expert_cmps <- unique(curr_expert_boot$Chem)
  curr_expert_cmps_length <- length(curr_expert_cmps)
  curr_expert_surrogates <- expert_surrogates[expert_surrogates$Compound %in% curr_expert_cmps,] # Limit compounds to those in the current media
  # NOTE: it uses as many of the three surrogates as are present in the Media; if any surrogate compounds are missing, then the bootstrap will be using less than 3 surrogates
  # Media-specific surrogates may need to be chosen
  for(j in 1:nrow(curr_expert_surrogates)){
    expert_boot_out <- rbind(expert_boot_out,cbind(Chem=curr_expert_surrogates$Compound[j],Media=i,full_boot(curr_expert_boot[curr_expert_boot$Chem %in% curr_expert_surrogates[j,2:4],],runs=10000,med_out=FALSE,resample_size=curr_expert_cmps_length)))
  }  
}
elapsed <- Sys.time()-start_time
print(elapsed)
# took ~18.7 minutes to complete

# Get exceedances for expert-selected surrogates bootstrap
expert_boot_out_exceedances<-dplyr::left_join(invpred_results %>% select(Conc,Chem,Media,Rep,Area,Mass_Labeled_Chem,Mass_Labeled_Molar_Mass),expert_boot_out) %>%
  mutate(RF_0.975 = med_97.5th_percentile,
         RF_0.5 = med_50th_percentile,
         RF_0.025 = med_2.5th_percentile,
         estimate=Area/med_50th_percentile,
         lower=Area/med_97.5th_percentile,
         upper=Area/med_2.5th_percentile,
         estimate_ng_g=estimate*Mass_Labeled_Molar_Mass,
         lower_ng_g=lower*Mass_Labeled_Molar_Mass,
         upper_ng_g=upper*Mass_Labeled_Molar_Mass,
         AQ=10^log10(estimate_ng_g/Conc),
         AAQ=10^abs(log10(estimate_ng_g/Conc)),
         CLFR=upper_ng_g/lower_ng_g,
         LCQ=Conc/lower_ng_g,
         UCQ=upper_ng_g/Conc,
         lower_exceeded=Conc<lower_ng_g,
         upper_exceeded=Conc>upper_ng_g,
         either_exceeded=Conc<lower_ng_g|Conc>upper_ng_g) %>%
  select(-contains("perc"))

## GLOBAL SURROGATES BOOTSTRAP ## 
start_time <- Sys.time()
global_boot_out <- data.frame()
set.seed(1)
# do LOOCV bootstrap
for(i in uniq_media_RM2){
  curr_global_boot <- Matrices_Data_cleaned_RM2[Matrices_Data_cleaned_RM2$Media==i,]
  curr_global_cmps <- unique(curr_global_boot$Chem)
  global_boot_out <- rbind(global_boot_out,cbind(Media=i,Chem=curr_global_cmps,loocv_new(curr_global_boot,runs=10000,med_out=FALSE)))
}
elapsed <- Sys.time()-start_time
print(elapsed)
# Takes ~19.4  minutes to complete

# Get exceedances for global surrogates bootstrap
global_boot_out_exceedances <- dplyr::left_join(invpred_results %>% select(Conc,Chem,Media,Rep,Area,Mass_Labeled_Chem,Mass_Labeled_Molar_Mass),global_boot_out) %>%
  mutate(RF_0.975 = med_97.5th_percentile,
         RF_0.5 = med_50th_percentile,
         RF_0.025 = med_2.5th_percentile,
         estimate=Area/med_50th_percentile,
         lower=Area/med_97.5th_percentile,
         upper=Area/med_2.5th_percentile,
         estimate_ng_g=estimate*Mass_Labeled_Molar_Mass,
         lower_ng_g=lower*Mass_Labeled_Molar_Mass,
         upper_ng_g=upper*Mass_Labeled_Molar_Mass,
         AQ=10^log10(estimate_ng_g/Conc),
         AAQ=10^abs(log10(estimate_ng_g/Conc)),
         CLFR=upper_ng_g/lower_ng_g,
         LCQ=Conc/lower_ng_g,
         UCQ=upper_ng_g/Conc,
         lower_exceeded=Conc<lower_ng_g,
         upper_exceeded=Conc>upper_ng_g,
         either_exceeded=Conc<lower_ng_g|Conc>upper_ng_g) %>%
  select(-contains("perc"))

## IONIZATION EFFICIENCY REGRESSION - LINEAR MIXED MODELS ##
pred_IEs <- read.csv("Pu_et_al_Env_Matrices_predicted_IEs_RT.csv")

# For PFBA, its logIE_predicted_LOO is the same as its logIE_predicted because it was not in the training data
pred_IEs$Corresponding_Chemical[pred_IEs$ID=="PFBA"] <- pred_IEs$ID[pred_IEs$ID=="PFBA"]
pred_IEs$logIE_predicted_LOO[pred_IEs$ID=="PFBA"] <- pred_IEs$logIE_predicted[pred_IEs$ID=="PFBA"]

Matrices_Data_IE <- Matrices_Data_cleaned_RM2 %>%
  dplyr::left_join(pred_IEs %>% select(ID,logIE_predicted_LOO),by=c("Chem"="ID"))

# Calculated leverage for each chemical for per-media regressions (see "Pu_et_al_Env_Matrices_SI.R")
# Highest leverage is for PFBA- 0.0417 in Solvent, 0.0666 in DW, 0.554 in Sludge

# Data with PFBA removed (since it has high leverage)
Matrices_Data_IE_no_PFBA <- Matrices_Data_cleaned_RM2 %>%
  dplyr::left_join(pred_IEs %>% 
                     mutate(ID = gsub("FTS42","4:2 FTS",ID) %>%
                              gsub("FTS62","6:2 FTS",.) %>%
                              gsub("FTS82","8:2 FTS",.)) %>%
                     select(ID,logIE_predicted_LOO),by=c("Chem"="ID")) %>%
  filter(Chem != "PFBA")

# Load leave-one-out models
for(i in 1:19){
  model_name <- paste0("model",i)
  assign(model_name,readRDS(paste0("./loo_models/model_",i,".RData")))
}

# Load molecular descriptor dataframe
mol_descs <- read.csv("sepman_water_matrices_descs_calc.csv",check.names=FALSE)
mol_descs <- mol_descs %>%
  mutate(ID = gsub("FTS42","4:2 FTS",ID) %>%
           gsub("FTS62","6:2 FTS",.) %>%
           gsub("FTS82","8:2 FTS",.))

# Get predicted IEs from each leave-one-out model
loo_model_results <- data.frame()
for(i in 1:nrow(mol_descs)){
  curr_loo_results <- predict(get(paste0("model",i)),newdata=mol_descs)
  loo_model_results <- rbind(loo_model_results,data.frame(model = paste0("model",i),
                                                          loo_Chem = mol_descs$ID[i],
                                                          Chem = mol_descs$ID,
                                                          logIE_predicted = curr_loo_results))
}

## MIXED MODEL FITTING
# Separate models for each chemical (leave-one-out) and each matrix
RF_vs_IE_per_chem_matrix <- vector(mode="list",length=3L)
for(i in 1:length(uniq_media_RM2)){
  curr_media_chem <- unique(Matrices_Data_IE_no_PFBA[Matrices_Data_IE_no_PFBA$Media == uniq_media_RM2[i],c("Chem")])
  for(j in 1:nrow(curr_media_chem)){
    curr_RF_vs_IE_data <- Matrices_Data_IE_no_PFBA %>%
      filter(Media == uniq_media_RM2[i] & Chem != curr_media_chem[j,1][[1]]) %>%
      dplyr::left_join(loo_model_results %>%
                         filter(model == paste0("model",i+1))) #skip past model 1 because that's for PFBA
    RF_vs_IE_per_chem_matrix[[i]][[j]] <- lme4::lmer(LogRF_uM ~ logIE_predicted + (1|Chem),curr_RF_vs_IE_data)    
  }
}

## PARAMETRIC BOOTSTRAP
# Separated models per chemical and per matrix
start_time <- Sys.time()
RF_vs_IE_per_chem_matrix_boot <- vector(mode="list",length=length(uniq_media_RM2))
for(i in 1:length(RF_vs_IE_per_chem_matrix)){
  for(j in 1:18){
    curr_boot_data <- RF_vs_IE_per_chem_matrix[[i]][[j]]@frame
    # Function outputs predictions for all input data based on the fitted model
    # re.form=NULL indicates to use all random effects (see predict.merMod)
    bootMer_fun <- function(currLmer){
      predict(currLmer,newdata=curr_boot_data,re.form=NULL)
    }
    # Mixed model bootstrap
    # NOTE: this is pretty slow (because it is a n=10,000 iteration bootstrap)
    # 
    RF_vs_IE_per_chem_matrix_boot[[i]][[j]] <- bootMer(x= RF_vs_IE_per_chem_matrix[[i]][[j]],
                                                  FUN = bootMer_fun,
                                                  nsim=10000,
                                                  seed=1,
                                                  type='parametric',
                                                  use.u=FALSE)    
  }
}
elapsed <- Sys.time()-start_time
print(elapsed)
# Took ~4.1 hours

# Get 95% PI bounds
# same as sumBoot() from Groff et al 2022
# Iterate over the columns of bootMer_out$t, where each column is an observation from the data and each row is all observations from a bootstrap iteration 
# Get the 2.5%, 50%, and 97.5% of the given observation point from all the bootstrap simulations
boot_PI_fun <- function(bootMer_out){
  return(
    data.frame(fit=apply(bootMer_out$t,2,function(x) as.numeric(quantile(x,probs=.5,na.rm=TRUE))),
               lwr=apply(bootMer_out$t,2,function(x) as.numeric(quantile(x,probs=.025,na.rm=TRUE))),
               upr=apply(bootMer_out$t,2,function(x) as.numeric(quantile(x,probs=.975,na.rm=TRUE)))
    )
  )
}

# This is for RF space, we want concentration space
# Separate models for each chemical in each matrix
RF_vs_IE_per_chem_matrix_boot_PI <- vector(mode="list",length=length(uniq_media_RM2))
for(i in 1:length(RF_vs_IE_per_chem_matrix_boot)){
  for(j in 1:length(RF_vs_IE_per_chem_matrix_boot[[i]])){
    curr_IEs <- RF_vs_IE_per_chem_matrix[[i]][[j]]@frame$logIE_predicted
    curr_PI <- boot_PI_fun(RF_vs_IE_per_chem_matrix_boot[[i]][[j]])
    # if all values were the 2.5% value for the observations
    lwr_lm <- lm(formula=curr_PI$lwr ~ curr_IEs)
    lwr_smooth <- (lwr_lm$coefficients[[2]]*curr_IEs)+(lwr_lm$coefficients[[1]])
    # if all values were the 50% value for the observations
    fit_lm <- lm(formula=curr_PI$fit ~ curr_IEs)
    fit_smooth <-( fit_lm$coefficients[[2]]*curr_IEs)+(fit_lm$coefficients[[1]])
    # if all values were the 97.5% value for the observations
    upr_lm <- lm(formula=curr_PI$upr ~ curr_IEs)
    upr_smooth <- (upr_lm$coefficients[[2]]*curr_IEs)+(upr_lm$coefficients[[1]])
    RF_vs_IE_per_chem_matrix_boot_PI[[i]][[j]] <- list(PI=curr_PI,
                                             lwr_lm=lwr_lm,
                                             fit_lm=fit_lm,
                                             upr_lm=upr_lm)
  }
}

# Join the different per-chemical, per-matrix model output to one dataframe
RF_vs_IE_per_chem_matrix_boot_PI_df <- data.frame()
for(i in 1:length(RF_vs_IE_per_chem_matrix_boot)){
  for(j in 1:length(RF_vs_IE_per_chem_matrix_boot[[i]])){
    uniq_media_chem <- unique(Matrices_Data_IE_no_PFBA %>%
                                filter(Media == uniq_media_RM2[i]) %>%
                                pull('Chem'))
    curr_chem <- uniq_media_chem[j]
    curr_IEs <- data.frame(curr_IEs = pred_IEs %>%
                             mutate(ID = gsub("FTS42","4:2 FTS",ID) %>%
                                      gsub("FTS62","6:2 FTS",.) %>%
                                      gsub("FTS82","8:2 FTS",.)) %>%
                             filter(ID == curr_chem) %>%
                             pull('logIE_predicted_LOO'))
    RF_vs_IE_per_chem_matrix_boot_PI_df <- rbind(RF_vs_IE_per_chem_matrix_boot_PI_df,
                                                 cbind(Matrices_Data_IE_no_PFBA %>%
                                                         filter(Media == uniq_media_RM2[i] & Chem == curr_chem),
                                                       data.frame(lwr = predict(RF_vs_IE_per_chem_matrix_boot_PI[[i]][[j]]$lwr_lm,curr_IEs),
                                                                  fit = predict(RF_vs_IE_per_chem_matrix_boot_PI[[i]][[j]]$fit_lm,curr_IEs),
                                                                  upr = predict(RF_vs_IE_per_chem_matrix_boot_PI[[i]][[j]]$upr_lm,curr_IEs)
                                                       ))
                                                 )
  }
}

# This is for RF space, we want concentration space
# Separate models for each chemical in each matrix
RF_vs_IE_per_chem_matrix_boot_PI <- vector(mode="list",length=length(uniq_media_RM2))
for(i in 1:length(RF_vs_IE_per_chem_matrix_boot)){
  curr_chem <- Matrices_Data_IE_no_PFBA %>% filter(Media == uniq_media_RM2[i]) %>% pull('Chem') %>% unique()
  for(j in 1:length(RF_vs_IE_per_chem_matrix_boot[[i]])){
    curr_IEs <- RF_vs_IE_per_chem_matrix_boot[[i]][[j]][["data"]]$logIE_predicted
    curr_PI <- boot_PI_fun(RF_vs_IE_per_chem_matrix_boot[[i]][[j]])
    # if all values were the 2.5% value for the observations
    lwr_lm <- lm(formula=curr_PI$lwr ~ curr_IEs)
    lwr_smooth <- (lwr_lm$coefficients[[2]]*curr_IEs)+(lwr_lm$coefficients[[1]])
    # if all values were the 50% value for the observations
    fit_lm <- lm(formula=curr_PI$fit ~ curr_IEs)
    fit_smooth <-( fit_lm$coefficients[[2]]*curr_IEs)+(fit_lm$coefficients[[1]])
    # if all values were the 97.5% value for the observations
    upr_lm <- lm(formula=curr_PI$upr ~ curr_IEs)
    upr_smooth <- (upr_lm$coefficients[[2]]*curr_IEs)+(upr_lm$coefficients[[1]])
    RF_vs_IE_per_chem_matrix_boot_PI[[i]][[j]] <- list(PI=curr_PI,
                                                       lwr_lm=lwr_lm,
                                                       fit_lm=fit_lm,
                                                       upr_lm=upr_lm)
    # Plot the LOO mixed model IE regressions
    RF_vs_IE_plot_data <- cbind(RF_vs_IE_per_chem_matrix_boot[[i]][[j]]$data,
                                lwr_smooth,
                                fit_smooth,
                                upr_smooth)
    loo_chem_data <- Matrices_Data_IE_no_PFBA %>% 
      filter(Media == uniq_media_RM2[i] & Chem == curr_chem[j])
    # Figure S3 uses PFOS plots
    ggplot(data=RF_vs_IE_per_chem_matrix_boot[[i]][[j]]$data,aes(x=logIE_predicted,y=LogRF_uM))+
      geom_point()+
      geom_point(data=loo_chem_data,aes(x=logIE_predicted_LOO,y=LogRF_uM),color="red")+
      geom_line(aes(x=logIE_predicted,y=lwr_smooth))+
      geom_line(aes(x=logIE_predicted,y=fit_smooth))+
      geom_line(aes(x=logIE_predicted,y=upr_smooth))+
      xlab("Log(IE)")+
      ylab("Log(RF)")+
      theme_bw(base_size=18)+
      theme(axis.text=element_text(size=18))+
      ggtitle(paste0("Log(RF) vs Log(IE) regression ",uniq_media_RM2[i]," ",curr_chem[j]))
    #ggsave(paste0("RF_vs_IE_per_chem_matrix_",uniq_media_RM2[i],"_",gsub(":2 ","2",curr_chem[j]),".png"),width=8,height=8,units="in")
  }
}

## Join prediction intervals to the original data
# By chemical, by matrix
Matrices_Data_RF_vs_IE_per_chem_matrix <- RF_vs_IE_per_chem_matrix_boot_PI_df %>%
  mutate(estimate=Area/(10^fit),
         lower=Area/(10^upr),
         upper=Area/(10^lwr),
         estimate_ng_g=estimate*Mass_Labeled_Molar_Mass,
         lower_ng_g=lower*Mass_Labeled_Molar_Mass,
         upper_ng_g=upper*Mass_Labeled_Molar_Mass,
         AQ=10^log10(estimate_ng_g/Conc),
         AAQ=10^abs(log10(estimate_ng_g/Conc)),
         CLFR=upper_ng_g/lower_ng_g,
         LCQ=Conc/lower_ng_g,
         UCQ=upper_ng_g/Conc,
         lower_exceeded=Conc<lower_ng_g,
         upper_exceeded=Conc>upper_ng_g,
         either_exceeded=Conc<lower_ng_g|Conc>upper_ng_g)
# Group per-chemical to get exceedance percentages (for LRP, URP, ORP)
Matrices_Data_RF_vs_IE_bychem_matrix <- Matrices_Data_RF_vs_IE_per_chem_matrix %>%
  group_by(Chem,Media) %>%
  dplyr::left_join(count(.)) %>%
  mutate(perc_above=sum(upper_exceeded)/n,
         perc_below=sum(lower_exceeded)/n)
Matrices_Data_RF_vs_IE_bychem_matrix_uniq <- Matrices_Data_RF_vs_IE_bychem_matrix %>%
  select(Chem,Media,logIE_predicted_LOO,perc_above,perc_below) %>%
  unique()

# Remove some columns before joining with other data
Matrices_Data_RF_vs_IE_all_approaches <- Matrices_Data_RF_vs_IE_per_chem_matrix %>%
  select(-c(lwr,fit,upr))

### Combine results into one dataframe (for Supplementary Material)
# Wide has different columns for results from each approach, long has one column called 'approach'
all_approaches_results_wide <- invpred_results %>%
  rename_with(~paste0(.x,"_","cal_curve"),.cols=c(estimate,lower,upper,estimate_ng_g,lower_ng_g,upper_ng_g,AQ,AAQ,CLFR,LCQ,UCQ,lower_exceeded,upper_exceeded,either_exceeded))
approach_list <- list(expert_boot=expert_boot_out_exceedances,global_boot=global_boot_out_exceedances,IE_regression=Matrices_Data_RF_vs_IE_all_approaches)
for(i in 1:length(approach_list)){
  curr_approach <- approach_list[[i]]
  all_approaches_results_wide <- dplyr::left_join(all_approaches_results_wide,curr_approach %>% rename_with(~paste0(.x,"_",names(approach_list)[i]),.cols=c(RF_0.975,RF_0.5,RF_0.025,estimate,lower,upper,estimate_ng_g,lower_ng_g,upper_ng_g,AQ,AAQ,CLFR,LCQ,UCQ,lower_exceeded,upper_exceeded,either_exceeded)))
}

all_approaches_results_long <- rbind(invpred_results %>% select(-RT) %>% select(-contains(c("Nominal","Free_Acid","_uM","Percent","Log","Rsquared","RF"))) %>% mutate(approach="cal_curve"),
                                     global_boot_out_exceedances %>% mutate(approach="global_boot"),
                                     expert_boot_out_exceedances %>% mutate(approach="expert_boot"),
                                     Matrices_Data_RF_vs_IE_all_approaches %>% select(-RT) %>% select(-contains(c("Nominal","Free_Acid","Log","_uM","logIE"))) %>% mutate(approach="IE_regression"))

# Change Rep column to integer
all_approaches_results_wide <- all_approaches_results_wide %>% 
  mutate(Rep = as.integer(Rep))
all_approaches_results_long <- all_approaches_results_long %>% 
  mutate(Rep = as.integer(Rep))

# Change Conc to numeric and sort by Conc, Chem, Media, Rep
all_approaches_results_wide <- all_approaches_results_wide %>% 
  mutate(Conc = as.numeric(Conc)) %>%
  arrange(across(c(Conc,Chem,Media,Rep)))

# Replace NA with empty string
# Rename
all_approaches_results_SI <- all_approaches_results_wide %>%
  replace(is.na(.),"") %>%
  dplyr::left_join(pred_IEs %>%
                     pivot_longer(Solvent_RT:Sludge_RT,
                                  names_to="Media",
                                  values_to="RT_input") %>%
                     mutate(Media =gsub("_RT","",Media)) %>%
                     rename(Chem = ID) %>%
                     select(Chem,Media,RT_input)) %>%
  mutate(Media=factor(Media,levels=c("Solvent","DW","Sludge"))) %>%
  arrange(Media,RT_input,Conc,Rep) %>%
  select(-RT_input) %>%
  mutate(Mass_Labeled_Chem = Mass_Labeled_Chem %>%
           gsub("FTS42","4:2 FTS",.) %>%
           gsub("FTS62","6:2 FTS",.) %>%
           gsub("FTS82","8:2 FTS",.))

###
## Make plots of mixed-effect models for predicted IE

### STATISTICAL COMPARISON OF APPROACHES ### 
## CALIBRATION CURVE METRICS (SLOPE AND RSQUARED)
# Load calibration curve metrics
load("Pu_et_al_Env_Matrices_cal_curve_params.RData")

# Keep only RM2 DW, Sludge, Solvent
cal_curve_params_RM2 <- cal_curve_params %>%
  filter(Media %in% c("DW","Sludge","Solvent"))

# Median of slopes for each media
cal_curve_params_RM2 %>% 
  group_by(Media) %>%
  summarize_at("slope",median,na.rm=TRUE)

# Range of slopes for each media
cal_curve_params_RM2 %>%
  group_by(Media) %>%
  mutate(min_slope = min(slope,na.rm=TRUE),
         max_slope = max(slope,na.rm=TRUE)) %>% 
  select(Media,min_slope,max_slope) %>%
  unique()

slope_binom_tests <- data.frame()

# Sign test (median of 1)
for(i in uniq_media){
  curr_slopes <- cal_curve_params %>% filter(Media==i) %>% pull('slope')
  curr_slopes_sub <- curr_slopes - 1
  slope_binom_tests <- rbind(slope_binom_tests,data.frame(Media=i,p_val=binom.test(sum(curr_slopes_sub<0),n=length(curr_slopes_sub),p=0.5)$p.value))
}

# Is AQ significantly different from 1?
# AQ median
AQ_wilcx_tests <- data.frame()
approaches_AQ_median<-all_approaches_results_wide %>%
  select(Chem,Conc,Media,Rep,starts_with("AQ")) %>%
  group_by(Chem,Media) %>%
  summarize(across(starts_with("AQ_"), ~median(.x,na.rm=TRUE))) %>%
  ungroup()
for(i in uniq_media_RM2){
  approaches_AQ_median_media <- approaches_AQ_median %>%
    filter(Media == i)
  for(j in c("_cal_curve","_expert_boot","_global_boot","_IE_regression")){
    AQ_wilcx_tests <- rbind(AQ_wilcx_tests,data.frame(Media=i,
                                                      appraoch=j,
                                                      p_val=wilcox.test(approaches_AQ_median_media %>% 
                                                                                     pull(paste0('AQ',j)),mu=1)$p.value))
  }
}

## OVERALL STATISTICS (AAQ, CLFR, UCQ)
# By approach and media
all_approaches_results_median <- all_approaches_results_long %>%
  mutate(approach = factor(approach, levels=c("cal_curve","expert_boot","global_boot","IE_regression"))) %>%
  group_by(Media,approach) %>%
  summarize_at(c("AQ","AAQ","CLFR","UCQ"),median)

all_approaches_results_ORP <- all_approaches_results_long %>%
  mutate(approach = factor(approach, levels=c("cal_curve","expert_boot","global_boot","IE_regression"))) %>%
  group_by(Media,approach) %>%
  mutate(n=n()) %>%
  group_by(Media,approach,n) %>%
  summarize_at(c("either_exceeded"),sum) %>%
  mutate(ORP=1-(either_exceeded/n))

all_approaches_results_ranges <- all_approaches_results_long %>%
  mutate(approach = factor(approach, levels=c("cal_curve","expert_boot","global_boot","IE_regression"))) %>%
  group_by(Media,approach) %>%
  mutate(AQ_min = min(AQ,na.rm=TRUE),
         AQ_max = max(AQ,na.rm=TRUE),
         AAQ_min = min(AAQ,na.rm=TRUE),
         AAQ_max = max(AAQ,na.rm=TRUE),
         CLFR_min = min(CLFR,na.rm=TRUE),
         CLFR_max = max(CLFR,na.rm=TRUE)) %>% 
  select(Media,approach,AQ_min:CLFR_max) %>%
  unique()

# For all media
all_approaches_results_all_median <- all_approaches_results_long %>% 
  mutate(approach = factor(approach, levels=c("cal_curve","expert_boot","global_boot","IE_regression"))) %>%
  group_by(approach) %>%
  summarize_at(c("AQ","AAQ","CLFR","UCQ"),median)

all_approaches_results_all_ORP <- all_approaches_results_long %>% 
  mutate(approach = factor(approach, levels=c("cal_curve","expert_boot","global_boot","IE_regression"))) %>%
  group_by(approach) %>%
  mutate(ORP = 1-(sum(either_exceeded)/n())) %>%
  select(approach,ORP) %>% 
  unique()

 all_approaches_results_all_ranges <- all_approaches_results_long %>%
    mutate(approach = factor(approach, levels=c("cal_curve","expert_boot","global_boot","IE_regression"))) %>%
    group_by(approach) %>%
    mutate(AQ_min = min(AQ,na.rm=TRUE),
           AQ_max = max(AQ,na.rm=TRUE),
           AAQ_min = min(AAQ,na.rm=TRUE),
           AAQ_max = max(AAQ,na.rm=TRUE),
           CLFR_min = min(CLFR,na.rm=TRUE),
           CLFR_max = max(CLFR,na.rm=TRUE)) %>% 
    select(approach,AQ_min:CLFR_max) %>%
    unique()

# LIST OF ALL APPROACH COMBINATIONS (USED FOR COMPARISON)
all_approach_list <- c("cal_curve","expert_boot","global_boot","IE_regression")
approach_combns <- combn(x=all_approach_list,m=2) %>% t()

# Accuracy Quotient (AQ)
# AQ (test for median = 0 for each approach in each media, not approach comparison)
AQ_median_test <- data.frame()

for(i in uniq_media_RM2){
  for(j in all_approach_list){
    curr_AQs <- all_approaches_results_long %>% filter(Media==i & approach==j) %>% pull('AQ')
    curr_AQs_sub <- curr_AQs-median(curr_AQs)
    AQ_median_test <- rbind(AQ_median_test,data.frame(Media=i,approach=j,median_AQ=median(curr_AQs),p_val=binom.test(sum(curr_AQs_sub<0),n=length(curr_AQs_sub),p=0.5)$p.value))
  }
}

# Absolute Accuracy Quotient (AAQ)
# Pairwise comparison of approaches with Wilcoxon 
approaches_AAQ_median<-all_approaches_results_wide %>%
  select(Chem,Conc,Media,Rep,starts_with("AAQ")) %>%
  group_by(Chem,Media) %>%
  summarize(across(contains("AAQ_"), ~median(.x,na.rm=TRUE))) %>%
  ungroup()

approaches_AAQ_wilcox<-vector(mode="list",length=nrow(approach_combns))
for(i in 1:length(uniq_media_RM2)){
  approaches_AAQ_median_media <- approaches_AAQ_median %>%
    filter(Media == uniq_media_RM2[i])
  for(j in 1:nrow(approach_combns)){
    pair1<-approach_combns[j,1]
    pair2<-approach_combns[j,2]
    AAQ_1<-paste0("AAQ_",pair1)
    AAQ_2<-paste0("AAQ_",pair2)
    AAQ_diff<-(approaches_AAQ_median_media %>% pull(AAQ_2) %>% log10())-(approaches_AAQ_median_media %>% pull(AAQ_1) %>% log10())
    AAQ_diff_name<-paste0("log_AAQ_diff_",pair2,"_",pair1)
    approaches_AAQ_median_media<-approaches_AAQ_median_media %>% 
      mutate(!!AAQ_diff_name := AAQ_diff)
    print(length(AAQ_diff[!is.na(AAQ_diff)]))
    approach_AAQ_wilcox<-wilcox.test(approaches_AAQ_median_media %>% pull(AAQ_diff_name),conf.int=TRUE)
    # the estimate is the pseudomedian
    # we will try using the mean of the difference in log10(AAQ) as the estimate
    AAQ_wilcox_mean<-mean(approaches_AAQ_median_media %>% pull(AAQ_diff_name),na.rm=TRUE)
    AAQ_wilcox_GM_ratio<-10^(AAQ_wilcox_mean)
    curr_list_num <- ((i-1)*nrow(approach_combns))+j
    approaches_AAQ_wilcox[[curr_list_num]]<-c(approach_AAQ_wilcox,mean=AAQ_wilcox_mean,GM_ratio=AAQ_wilcox_GM_ratio)
    names(approaches_AAQ_wilcox)[[curr_list_num]]<-paste0(pair2,"_VS_",pair1)
  }
}

approach_AAQ_wilcox_df<-sapply(approaches_AAQ_wilcox,function(x) unlist(x)) %>% t()

# Wilcoxon using results from all media
approaches_AAQ_median<-all_approaches_results_wide %>%
  select(Chem,Conc,Media,Rep,starts_with("AAQ")) %>%
  group_by(Chem) %>%
  summarize(across(contains("AAQ_"), ~median(.x,na.rm=TRUE))) %>%
  ungroup()

approaches_AAQ_wilcox<-vector(mode="list",length=nrow(approach_combns))
for(i in 1:nrow(approach_combns)){
  pair1<-approach_combns[i,1]
  pair2<-approach_combns[i,2]
  AAQ_1<-paste0("AAQ_",pair1)
  AAQ_2<-paste0("AAQ_",pair2)
  AAQ_diff<-(approaches_AAQ_median %>% pull(AAQ_2) %>% log10())-(approaches_AAQ_median %>% pull(AAQ_1) %>% log10())
  AAQ_diff_name<-paste0("log_AAQ_diff_",pair2,"_",pair1)
  approaches_AAQ_median<-approaches_AAQ_median %>% 
    mutate(!!AAQ_diff_name := AAQ_diff)
  approach_AAQ_wilcox<-wilcox.test(approaches_AAQ_median %>% pull(AAQ_diff_name),conf.int=TRUE)
  # the estimate is the pseudomedian
  # we will try using the mean of the difference in log10(AAQ) as the estimate
  AAQ_wilcox_mean<-mean(approaches_AAQ_median %>% pull(AAQ_diff_name),na.rm=TRUE)
  AAQ_wilcox_GM_ratio<-10^(AAQ_wilcox_mean)
  approaches_AAQ_wilcox[[i]]<-c(approach_AAQ_wilcox,mean=AAQ_wilcox_mean,GM_ratio=AAQ_wilcox_GM_ratio)
  names(approaches_AAQ_wilcox)[[i]]<-paste0(pair2,"_VS_",pair1)
}

approach_AAQ_wilcox_df<-sapply(approaches_AAQ_wilcox,function(x) unlist(x)) %>% t()

## CONFIDENCE LIMIT FOLD RANGE (CLFR)
# Pairwise comparison
# Wilcoxon signed-rank test (per media)
approaches_CLFR_median<-all_approaches_results_wide %>%
  select(Chem,Conc,Media,Rep,starts_with("CLFR")) %>%
  group_by(Chem,Media) %>%
  summarize(across(contains("CLFR_"), ~median(.x,na.rm=TRUE))) %>%
  ungroup()

approaches_CLFR_wilcox<-vector(mode="list",length=nrow(approach_combns))
for(i in 1:length(uniq_media_RM2)){
  approaches_CLFR_median_media <- approaches_CLFR_median %>%
    filter(Media==uniq_media_RM2[i])
  for(j in 1:nrow(approach_combns)){
    pair1<-approach_combns[j,1]
    pair2<-approach_combns[j,2]
    CLFR_1<-paste0("CLFR_",pair1)
    CLFR_2<-paste0("CLFR_",pair2)
    CLFR_diff<-(approaches_CLFR_median_media %>% pull(CLFR_2) %>% log10())-(approaches_CLFR_median_media %>% pull(CLFR_1) %>% log10())
    CLFR_diff_name<-paste0("log_CLFR_diff_",pair2,"_",pair1)
    approaches_CLFR_median_media<-approaches_CLFR_median_media %>% 
      mutate(!!CLFR_diff_name := CLFR_diff)
    print(nrow(approaches_CLFR_median_media %>% filter(!is.na(!!!CLFR_diff_name))))
    approach_CLFR_wilcox<-wilcox.test(approaches_CLFR_median_media %>% pull(CLFR_diff_name),conf.int=TRUE)
    # the estimate is the pseudomedian
    # we will try using the mean of the difference in log10(CLFR) as the estimate
    CLFR_wilcox_mean<-mean(approaches_CLFR_median_media %>% pull(CLFR_diff_name),na.rm=TRUE)
    CLFR_wilcox_GM_ratio<-10^(CLFR_wilcox_mean)
    curr_list_num <- ((i-1)*nrow(approach_combns))+j
    approaches_CLFR_wilcox[[curr_list_num]]<-c(approach_CLFR_wilcox,mean=CLFR_wilcox_mean,GM_ratio=CLFR_wilcox_GM_ratio)
    names(approaches_CLFR_wilcox)[[curr_list_num]]<-paste0(pair2,"_VS_",pair1)
  }
}

approach_CLFR_wilcox_df<-sapply(approaches_CLFR_wilcox,function(x) unlist(x)) %>% t()

# Wilcoxon signed-rank tests (non-parametric) for all media
approaches_CLFR_median<-all_approaches_results_wide %>%
  select(Chem,Conc,Media,Rep,starts_with("CLFR")) %>%
  group_by(Chem) %>%
  summarize(across(contains("CLFR_"), ~median(.x,na.rm=TRUE)))

approaches_CLFR_wilcox<-vector(mode="list",length=nrow(approach_combns))
for(i in 1:nrow(approach_combns)){
  pair1<-approach_combns[i,1]
  pair2<-approach_combns[i,2]
  CLFR_1<-paste0("CLFR_",pair1)
  CLFR_2<-paste0("CLFR_",pair2)
  CLFR_diff<-(approaches_CLFR_median %>% pull(CLFR_2) %>% log10())-(approaches_CLFR_median %>% pull(CLFR_1) %>% log10())
  CLFR_diff_name<-paste0("log_CLFR_diff_",pair2,"_",pair1)
  approaches_CLFR<-approaches_CLFR_median %>% 
    mutate(!!CLFR_diff_name := CLFR_diff)
  approach_CLFR_wilcox<-wilcox.test(approaches_CLFR %>% pull(CLFR_diff_name),conf.int=TRUE)
  # the estimate is the pseudomedian
  # we will try using the mean of the difference in log10(CLFR) as the estimate
  CLFR_wilcox_mean<-mean(approaches_CLFR %>% pull(CLFR_diff_name),na.rm=TRUE)
  CLFR_wilcox_GM_ratio<-10^(CLFR_wilcox_mean)
  approaches_CLFR_wilcox[[i]]<-c(approach_CLFR_wilcox,mean=CLFR_wilcox_mean,GM_ratio=CLFR_wilcox_GM_ratio)
  names(approaches_CLFR_wilcox)[[i]]<-paste0(pair2,"_VS_",pair1)
}

approach_CLFR_wilcox_df<-sapply(approaches_CLFR_wilcox,function(x) unlist(x)) %>% t()

## RELIABILITY (EXCEEDANCE PERCENTAGE)
#ORP
all_approaches_results_long %>% 
  group_by(Media,approach) %>% 
  summarize_at(c("lower_exceeded","upper_exceeded"),sum) %>% 
  mutate(lower_exceeded_perc=lower_exceeded/nrow(all_approaches_results_wide),
         upper_exceeded_perc=upper_exceeded/nrow(all_approaches_results_wide)) %>%
  mutate(URP=1-upper_exceeded_perc,
         ORP=1-lower_exceeded_perc-upper_exceeded_perc)

## Statistical tests for performance between environmental matrices
media_combns <- combn(x=uniq_media_RM2,m=2) %>% t()
uniq_approaches <- c("cal_curve","expert_boot","global_boot","IE_regression")

all_approaches_permedia <- all_approaches_results_long %>%
  select(Conc:Rep,AQ:CLFR,approach) %>%
  pivot_wider(names_from=Media,
              values_from=AQ:CLFR)

# Wilcoxon signed-rank tests (non-parametric) for all media
approaches_AAQ_median<-all_approaches_permedia %>%
  select(Chem,Conc,Rep,approach,starts_with("AAQ")) %>%
  group_by(Chem,approach) %>%
  summarize(across(contains("AAQ_"), ~median(.x,na.rm=TRUE))) %>%
  ungroup()

approaches_AAQ_wilcox<-vector(mode="list",length=nrow(media_combns)*length(uniq_approaches))
for(i in 1:length(uniq_approaches)){
  for(j in 1:nrow(media_combns)){
    pair1<-media_combns[j,1]
    pair2<-media_combns[j,2]
    AAQ_1<-paste0("AAQ_",pair1)
    AAQ_2<-paste0("AAQ_",pair2)
    curr_AAQ_data <- approaches_AAQ_median %>%
      ungroup() %>%
      filter(approach == uniq_approaches[i])
    AAQ_diff<-(curr_AAQ_data %>% pull(AAQ_2) %>% log10())-(curr_AAQ_data %>% pull(AAQ_1) %>% log10())
    AAQ_diff_name<-paste0("log_AAQ_diff_",pair2,"_",pair1)
    curr_AAQ_data<-curr_AAQ_data %>% 
      mutate(!!AAQ_diff_name := AAQ_diff)
    print(nrow(curr_AAQ_data %>% filter(!is.na(!!AAQ_diff_name))))
    approach_AAQ_wilcox<-wilcox.test(curr_AAQ_data %>% pull(AAQ_diff_name),conf.int=TRUE)
    # the estimate is the pseudomedian
    # we will try using the mean of the difference in log10(CLFR) as the estimate
    curr_list_num <- ((i-1)*nrow(media_combns))+j
    AAQ_wilcox_mean<-mean(curr_AAQ_data %>% pull(AAQ_diff_name),na.rm=TRUE)
    AAQ_wilcox_GM_ratio<-10^(AAQ_wilcox_mean)
    approaches_AAQ_wilcox[[curr_list_num]]<-c(approach_AAQ_wilcox,mean=AAQ_wilcox_mean,GM_ratio=AAQ_wilcox_GM_ratio)
    names(approaches_AAQ_wilcox)[[curr_list_num]]<-paste0(pair2,"_VS_",pair1,"_",uniq_approaches[i])
  }  
}
approach_AAQ_wilcox_df<-sapply(approaches_AAQ_wilcox,function(x) unlist(x)) %>% t()

# CLFR
# Wilcoxon signed-rank tests (non-parametric)
approaches_CLFR_median<-all_approaches_permedia %>%
  select(Chem,Conc,Rep,approach,starts_with("CLFR")) %>%
  group_by(Chem,approach) %>%
  summarize(across(contains("CLFR_"), ~median(.x,na.rm=TRUE)))

approaches_CLFR_wilcox<-vector(mode="list",length=nrow(media_combns)*length(uniq_approaches))
for(i in 1:length(uniq_approaches)){
  for(j in 1:nrow(media_combns)){
    pair1<-media_combns[j,1]
    pair2<-media_combns[j,2]
    CLFR_1<-paste0("CLFR_",pair1)
    CLFR_2<-paste0("CLFR_",pair2)
    curr_CLFR_data <- approaches_CLFR_median %>%
      ungroup() %>%
      filter(approach == uniq_approaches[i])
    CLFR_diff<-(curr_CLFR_data %>% pull(CLFR_2) %>% log10())-(curr_CLFR_data %>% pull(CLFR_1) %>% log10())
    CLFR_diff_name<-paste0("log_CLFR_diff_",pair2,"_",pair1)
    curr_CLFR_data<-curr_CLFR_data %>% 
      mutate(!!CLFR_diff_name := CLFR_diff)
    print(nrow(curr_CLFR_data %>% filter(!is.na(!!CLFR_diff_name))))
    approach_CLFR_wilcox<-wilcox.test(curr_CLFR_data %>% pull(CLFR_diff_name),conf.int=TRUE)
    # the estimate is the pseudomedian
    # we will try using the mean of the difference in log10(CLFR) as the estimate
    curr_list_num <- ((i-1)*nrow(media_combns))+j
    CLFR_wilcox_mean<-mean(curr_CLFR_data %>% pull(CLFR_diff_name),na.rm=TRUE)
    CLFR_wilcox_GM_ratio<-10^(CLFR_wilcox_mean)
    approaches_CLFR_wilcox[[curr_list_num]]<-c(approach_CLFR_wilcox,mean=CLFR_wilcox_mean,GM_ratio=CLFR_wilcox_GM_ratio)
    names(approaches_CLFR_wilcox)[[curr_list_num]]<-paste0(pair2,"_VS_",pair1,"_",uniq_approaches[i])
  }  
}
approach_CLFR_wilcox_df<-sapply(approaches_CLFR_wilcox,function(x) unlist(x)) %>% t()

# Solvent-based predictions for DW and sludge
# Using solvent training data for quantification of DW and Sludge
Matrices_Data_solvent <- Matrices_Data_cleaned_RM2 %>% filter(Media == "Solvent")
uniq_solvent_chem <- Matrices_Data_solvent %>% pull('Chem') %>% unique()

solvent_cals <- list(length(uniq_solvent_chem))
for(i in 1:length(uniq_solvent_chem)){
  plot_lm<-lm(data=Matrices_Data_solvent %>% filter(Chem == uniq_solvent_chem[i]),LogArea~LogConc_uM)
  solvent_cals[[i]] <- plot_lm
}
names(solvent_cals) <- uniq_solvent_chem

## Using solvent training data for quantification of DW and Sludge
# Inverse predictions using solvent calibration curves
solvent_invpred_results <- data.frame()
for(i in c("DW","Sludge")){
  uniq_media_rows <- Matrices_Data_cleaned_RM2 %>% filter(Media==i)
  uniq_media_chem <- unique(uniq_media_rows %>% pull('Chem'))
  for(j in uniq_media_chem){
    curr_data <- uniq_media_rows %>% filter(Chem == j)
    if(j %in% uniq_solvent_chem){
      curr_LogArea <- curr_data %>% pull('LogArea')
      curr_invpred <- sapply(1:nrow(curr_data), function(x) calibrate(solvent_cals[[j]],y0=curr_LogArea[x]))%>% t()
      solvent_invpred_results <- rbind(solvent_invpred_results,cbind(curr_data,curr_invpred))
    }
  }
}

solvent_invpred_results <- solvent_invpred_results %>%
  mutate(across(estimate:upper,as.numeric)) %>%
  select(-interval) %>%
  mutate(estimate=(10^estimate),
         lower=(10^lower),
         upper=(10^upper),
         estimate_ng_g=estimate*Mass_Labeled_Molar_Mass,
         lower_ng_g=lower*Mass_Labeled_Molar_Mass,
         upper_ng_g=upper*Mass_Labeled_Molar_Mass,
         AQ=10^log10(estimate_ng_g/Conc),
         AAQ=10^abs(log10(estimate_ng_g/Conc)),
         CLFR=upper_ng_g/lower_ng_g,
         LCQ=Conc/lower_ng_g,
         UCQ=upper_ng_g/Conc,
         lower_exceeded=Conc<lower_ng_g,
         upper_exceeded=Conc>upper_ng_g,
         either_exceeded=Conc<lower_ng_g|Conc>upper_ng_g)

# Solvent expert bootstrap
expert_boot_out_solvent <- expert_boot_out %>% filter(Media == 'Solvent')

solvent_expert_boot_preds <- Matrices_Data_cleaned_RM2 %>%
  filter(Media %in% c("DW","Sludge")) %>%
  dplyr::left_join(expert_boot_out_solvent %>%
                     select(-Media) %>%
                     select(Chem,contains("med"))) %>%
  mutate(estimate=Area/med_50th_percentile,
         lower=Area/med_97.5th_percentile,
         upper=Area/med_2.5th_percentile,
         estimate_ng_g=estimate*Mass_Labeled_Molar_Mass,
         lower_ng_g=lower*Mass_Labeled_Molar_Mass,
         upper_ng_g=upper*Mass_Labeled_Molar_Mass,
         AQ=10^log10(estimate_ng_g/Conc),
         AAQ=10^abs(log10(estimate_ng_g/Conc)),
         CLFR=upper_ng_g/lower_ng_g,
         LCQ=Conc/lower_ng_g,
         UCQ=upper_ng_g/Conc,
         lower_exceeded=Conc<lower_ng_g,
         upper_exceeded=Conc>upper_ng_g,
         either_exceeded=Conc<lower_ng_g|Conc>upper_ng_g) %>%
  select(-contains("perc"))

# Solvent global bootstrap
global_boot_out_solvent <- global_boot_out %>% filter(Media == 'Solvent')

solvent_global_boot_preds <- Matrices_Data_cleaned_RM2 %>%
  filter(Media %in% c("DW","Sludge")) %>%
  dplyr::left_join(global_boot_out_solvent %>%
                     select(-Media) %>%
                     select(Chem,contains("med"))) %>%
  mutate(estimate=Area/med_50th_percentile,
         lower=Area/med_97.5th_percentile,
         upper=Area/med_2.5th_percentile,
         estimate_ng_g=estimate*Mass_Labeled_Molar_Mass,
         lower_ng_g=lower*Mass_Labeled_Molar_Mass,
         upper_ng_g=upper*Mass_Labeled_Molar_Mass,
         AQ=10^log10(estimate_ng_g/Conc),
         AAQ=10^abs(log10(estimate_ng_g/Conc)),
         CLFR=upper_ng_g/lower_ng_g,
         LCQ=Conc/lower_ng_g,
         UCQ=upper_ng_g/Conc,
         lower_exceeded=Conc<lower_ng_g,
         upper_exceeded=Conc>upper_ng_g,
         either_exceeded=Conc<lower_ng_g|Conc>upper_ng_g) %>%
  select(-contains("perc"))

# Join the different per-chemical, per-matrix model output to one dataframe
RF_vs_IE_per_chem_matrix_boot_PI_df <- data.frame()
for(i in 2:length(RF_vs_IE_per_chem_matrix_boot)){
  for(j in 1:length(RF_vs_IE_per_chem_matrix_boot[[i]])){
    uniq_media_chem <- unique(Matrices_Data_IE_no_PFBA %>%
                                filter(Media == uniq_media_RM2[i]) %>%
                                pull('Chem'))
    curr_chem <- uniq_media_chem[j]
    curr_IEs <- data.frame(curr_IEs = pred_IEs %>%
                             filter(ID == curr_chem) %>%
                             pull('logIE_predicted_LOO'))
    RF_vs_IE_per_chem_matrix_boot_PI_df <- rbind(RF_vs_IE_per_chem_matrix_boot_PI_df,
                                                 cbind(Matrices_Data_IE_no_PFBA %>%
                                                         filter(Media == uniq_media_RM2[i] & Chem == curr_chem),
                                                       data.frame(lwr = predict(RF_vs_IE_per_chem_matrix_boot_PI[[1]][[j]]$lwr_lm,curr_IEs),
                                                                  fit = predict(RF_vs_IE_per_chem_matrix_boot_PI[[1]][[j]]$fit_lm,curr_IEs),
                                                                  upr = predict(RF_vs_IE_per_chem_matrix_boot_PI[[1]][[j]]$upr_lm,curr_IEs)
                                                       ))
    )
  }
}

Matrices_Data_RF_vs_IE <- RF_vs_IE_per_chem_matrix_boot_PI_df %>%
  mutate(estimate=Area/(10^fit),
         lower=Area/(10^upr),
         upper=Area/(10^lwr),
         estimate_ng_g=estimate*Mass_Labeled_Molar_Mass,
         lower_ng_g=lower*Mass_Labeled_Molar_Mass,
         upper_ng_g=upper*Mass_Labeled_Molar_Mass,
         AQ=10^log10(estimate_ng_g/Conc),
         AAQ=10^abs(log10(estimate_ng_g/Conc)),
         CLFR=upper_ng_g/lower_ng_g,
         LCQ=Conc/lower_ng_g,
         UCQ=upper_ng_g/Conc,
         lower_exceeded=Conc<lower_ng_g,
         upper_exceeded=Conc>upper_ng_g,
         either_exceeded=Conc<lower_ng_g|Conc>upper_ng_g)

# Save solvent-based results
solvent_all_approaches_results_wide <- solvent_invpred_results %>%
  select(Conc:LogRF_uM)
approach_list <- list(cal_curve=solvent_invpred_results,global_boot=solvent_global_boot_preds,expert_boot=solvent_expert_boot_preds,IE_regression=Matrices_Data_RF_vs_IE)
for(i in 1:length(approach_list)){
  curr_approach <- approach_list[[i]]
  solvent_all_approaches_results_wide <- dplyr::left_join(solvent_all_approaches_results_wide,curr_approach %>% rename_with(~paste0(.x,"_",names(approach_list)[i]),.cols=c(estimate,lower,upper,estimate_ng_g,lower_ng_g,upper_ng_g,AQ,AAQ,CLFR,LCQ,UCQ,lower_exceeded,upper_exceeded,either_exceeded)))
}

solvent_all_approaches_results_long <- rbind(solvent_invpred_results %>% select(-RT) %>% select(-contains(c("Nominal","Free_Acid","_uM","Percent","Log","Rsquared","RF"))) %>% mutate(approach="cal_curve"),
                                     solvent_global_boot_preds %>% mutate(approach="global_boot") %>% select(-c("RT","LogArea","Nominal","Free_Acid","Conc_uM","RF_uM","LogConc_uM","LogRF_uM")),
                                     solvent_expert_boot_preds %>% mutate(approach="expert_boot") %>% select(-c("RT","LogArea","Nominal","Free_Acid","Conc_uM","RF_uM","LogConc_uM","LogRF_uM")),
                                     Matrices_Data_RF_vs_IE %>% select(-RT) %>% select(-contains(c("Nominal","Free_Acid","Log","_uM","logIE","lwr","fit","upr"))) %>% mutate(approach="IE_regression"))

# Change Conc to numeric and sort by Conc, Chem, Media, Rep
solvent_all_approaches_results_wide <- solvent_all_approaches_results_wide %>% 
  mutate(Conc = as.numeric(Conc)) %>%
  arrange(across(c(Conc,Chem,Media,Rep)))

# Replace NA with empty string
# Rename
solvent_all_approaches_results_SI <- solvent_all_approaches_results_wide %>%
  replace(is.na(.),"") %>%
  dplyr::left_join(pred_IEs %>%
                     pivot_longer(Solvent_RT:Sludge_RT,
                                  names_to="Media",
                                  values_to="RT_input") %>%
                     mutate(Media = gsub("_RT","",Media)) %>%
                     rename(Chem = ID) %>%
                     select(Chem,Media,RT_input)) %>%
  arrange(Media,RT_input,Conc,Rep) %>%
  select(-RT_input)

# Figure 2 - One randomly selected point for each chemical. Comparing matrix-based vs solvent-based results for DW
random_per_chem <- data.frame()
DW_data <- all_approaches_results_long %>% filter(Media == "DW")
uniq_chem_DW <- DW_data %>% pull('Chem') %>% unique()

set.seed(3)
for(i in 1:length(uniq_chem_DW)){
  per_chem_rows <- DW_data %>% filter(Chem == uniq_chem_DW[i]) %>% drop_na()
  random_per_chem <- rbind(random_per_chem,per_chem_rows[sample(nrow(per_chem_rows),size=1),])
}

# Figure 2
# Facet by Training_Data (qNTA approaches only)
random_per_chem_solvent_vs_matrix_data <- dplyr::left_join(random_per_chem %>% select(Conc,Chem,Rep),all_approaches_results_wide) %>%
  select(Conc,Chem,Rep,Media,estimate_expert_boot:either_exceeded_IE_regression) %>%
  select(-contains("RF_0")) %>%
  filter(Chem != "PFBA" & Media == "DW") %>%
  mutate(Training_Data = "Matrix-matched") %>%
  rbind(dplyr::left_join(random_per_chem %>% select(Conc,Chem,Rep),solvent_all_approaches_results_wide %>%
                           mutate(Rep=as.character(Rep))) %>%
          select(Conc,Chem,Rep,Media,estimate_global_boot:either_exceeded_IE_regression) %>%
          select(-c(lwr:upr)) %>%
          relocate(estimate_global_boot:either_exceeded_global_boot,.after=Media) %>%
          filter(Chem != "PFBA" & Media == "DW") %>%
          mutate(Training_Data = "Solvent-based"))

ggplot(data=random_per_chem_solvent_vs_matrix_data) + 
         geom_point(aes(x=reorder(Chem,Conc),y=estimate_ng_g_expert_boot,color="Expert surrogates bootstrap"),position=position_nudge(x=-0.3))+
         geom_linerange(aes(x=reorder(Chem,Conc),ymin=lower_ng_g_expert_boot,ymax=upper_ng_g_expert_boot),position=position_nudge(x=-0.3),color="blue")+
         geom_point(aes(x=reorder(Chem,Conc),y=estimate_ng_g_global_boot,color="Global surrogates bootstrap"))+
         geom_linerange(aes(x=reorder(Chem,Conc),ymin=lower_ng_g_global_boot,ymax=upper_ng_g_global_boot),color="black")+
         geom_point(aes(x=reorder(Chem,Conc),y=estimate_ng_g_IE_regression,color="IE Regression"),position=position_nudge(x=0.3))+
         geom_linerange(aes(x=reorder(Chem,Conc),ymin=lower_ng_g_IE_regression,ymax=upper_ng_g_IE_regression),position=position_nudge(x=0.3),color="purple")+
         geom_segment(aes(x=as.numeric(reorder(Chem,Conc))-0.4,xend=as.numeric(reorder(Chem,Conc))+0.4,y=Conc,yend=Conc),color="red")+
         scale_y_log10(breaks = scales::trans_breaks("log10", function(y) 10^y),
                       labels = scales::trans_format("log10", scales::math_format(10^.x)))+
         theme_bw(base_size = 24)+
         theme(axis.text.x=element_text(angle=90,vjust=0.5,hjust=1),text=element_text(size=24))+
         ylab("Concentration (ng/mL)")+
         xlab("Chemical")+
         scale_color_manual(limits=c("Expert surrogates bootstrap","Global surrogates bootstrap","IE Regression","True concentration"),
                            values=c("blue","black","purple","red"),
                            labels=c(expression(BRF[expert]),expression(BRF[global]),expression(IE[reg]),"True concentration"),
                            name="Legend")+
         guides(color=guide_legend(override.aes=list(linetype=c(0,0,0,0))))+
         theme(legend.title = element_blank(),
               legend.position = c(0.35, 0.2),
               legend.background = element_rect(colour = 'black', fill = 'white', linetype='solid'))+
         facet_grid(cols=vars(Training_Data))

# Statistical comparison of solvent-based vs. matrix-matched qNTA estimates
# Statistical comparison of solvent-based and matrix-matched predictions
solvent_vs_matrix_data <- all_approaches_results_wide %>%
  select(Conc,Chem,Rep,Media,estimate_expert_boot:either_exceeded_IE_regression) %>%
  select(-contains("RF_0")) %>%
  mutate(Training_Data = "Matrix-matched") %>%
  rbind(solvent_all_approaches_results_wide %>%
          select(Conc,Chem,Rep,Media,estimate_global_boot:either_exceeded_IE_regression) %>%
          select(-c(lwr:upr)) %>%
          relocate(estimate_global_boot:either_exceeded_global_boot,.after=Media) %>%
          mutate(Training_Data = "Solvent")) %>%
  filter(Media %in% c("DW","Sludge"))

solvent_vs_matrix_AAQ_median <- solvent_vs_matrix_data %>%
  select(Chem,Conc,Media,Rep,Training_Data,starts_with("AAQ")) %>%
  group_by(Chem,Media,Training_Data) %>%
  summarize(across(contains("AAQ_"), ~median(.x,na.rm=TRUE))) %>%
  ungroup()

qNTA_approaches <- c("expert_boot","global_boot","IE_regression")
solvent_vs_matrix_AAQ_wilcox_list<-vector(mode="list",length=3L)
n <- 1
for(i in c("DW","Sludge")){
  solvent_vs_matrix_AAQ_median_media <- solvent_vs_matrix_AAQ_median %>%
    filter(Media == i)
  for(j in 1:3){
    solvent_vs_matrix_AAQ_median_media_approach <- solvent_vs_matrix_AAQ_median_media %>%
      select(Chem,Media,Training_Data,paste0("AAQ_",qNTA_approaches[j])) %>% 
      pivot_wider(names_from=Training_Data,values_from=paste0("AAQ_",qNTA_approaches[j]))
    AAQ_diff<-(solvent_vs_matrix_AAQ_median_media_approach %>% pull('Solvent') %>% log10())-(solvent_vs_matrix_AAQ_median_media_approach %>% pull('Matrix-matched') %>% log10())
    AAQ_diff_name<-paste0("log_AAQ_diff_",qNTA_approaches[j],"_",i)
    print(length(AAQ_diff[!is.na(AAQ_diff)]))
    solvent_vs_matrix_AAQ_wilcox<-wilcox.test(AAQ_diff,conf.int=TRUE)
    # the estimate is the pseudomedian
    # we will try using the mean of the difference in log10(AAQ) as the estimate
    AAQ_wilcox_mean<-mean(AAQ_diff,na.rm=TRUE)
    AAQ_wilcox_GM_ratio<-10^(AAQ_wilcox_mean)
    solvent_vs_matrix_AAQ_wilcox_list[[n]]<-c(solvent_vs_matrix_AAQ_wilcox,mean=AAQ_wilcox_mean,GM_ratio=AAQ_wilcox_GM_ratio)
    names(solvent_vs_matrix_AAQ_wilcox_list)[[n]]<-paste0(AAQ_diff_name)
    n <- n+1
  }
}

solvent_vs_matrix_CLFR_median <- solvent_vs_matrix_data %>%
  select(Chem,Conc,Media,Rep,Training_Data,starts_with("CLFR")) %>%
  group_by(Chem,Media,Training_Data) %>%
  summarize(across(contains("CLFR_"), ~median(.x,na.rm=TRUE))) %>%
  ungroup()

solvent_vs_matrix_CLFR_wilcox_list<-vector(mode="list",length=3L)
n <- 1
for(i in c("DW","Sludge")){
  solvent_vs_matrix_CLFR_median_media <- solvent_vs_matrix_CLFR_median %>%
    filter(Media == i)
  for(j in 1:3){
    solvent_vs_matrix_CLFR_median_media_approach <- solvent_vs_matrix_CLFR_median_media %>%
      select(Chem,Media,Training_Data,paste0("CLFR_",qNTA_approaches[j])) %>% 
      pivot_wider(names_from=Training_Data,values_from=paste0("CLFR_",qNTA_approaches[j]))
    CLFR_diff<-(solvent_vs_matrix_CLFR_median_media_approach %>% pull('Solvent') %>% log10())-(solvent_vs_matrix_CLFR_median_media_approach %>% pull('Matrix-matched') %>% log10())
    CLFR_diff_name<-paste0("log_CLFR_diff_",qNTA_approaches[j],"_",i)
    print(length(CLFR_diff[!is.na(CLFR_diff)]))
    solvent_vs_matrix_CLFR_wilcox<-wilcox.test(CLFR_diff,conf.int=TRUE)
    # the estimate is the pseudomedian
    # we will try using the mean of the difference in log10(AAQ) as the estimate
    CLFR_wilcox_mean<-mean(CLFR_diff,na.rm=TRUE)
    CLFR_wilcox_GM_ratio<-10^(CLFR_wilcox_mean)
    solvent_vs_matrix_CLFR_wilcox_list[[n]]<-c(solvent_vs_matrix_CLFR_wilcox,mean=CLFR_wilcox_mean,GM_ratio=CLFR_wilcox_GM_ratio)
    names(solvent_vs_matrix_CLFR_wilcox_list)[[n]]<-paste0(CLFR_diff_name)
    n <- n+1
  }
}

# Save results to RData Object
save(invpred_results,
     global_boot_out_exceedances,
     expert_boot_out_exceedances,
     Matrices_Data_IE,
     Matrices_Data_RF_vs_IE_per_chem_matrix,
     all_approaches_results_wide,
     all_approaches_results_long,
     all_approaches_results_SI,
     solvent_all_approaches_results_wide,
     solvent_all_approaches_results_long,
     solvent_all_approaches_results_SI,
     file="Pu_et_al_qNTA_Env_Matrices_final_results.RData")