Multiple Contaminants Highway Traffic Regression

## Script info ---------------------------
##
## Script name:"Highway Database Regression"
##
##
## Author: Christian Nilsen, Geosyntec Consultants
## Email: cnilsen@geosyntec.com
##
## Date Created: 2019-08-02
##
## Copyright (c) Geosyntec Consultants, 2019
##
## License---------------------------
##
##This program is free software: you can redistribute it and/or modify
##it under the terms of the GNU General Public License as published by
##the Free Software Foundation, version 3.0 
##
##This program is distributed in the hope that it will be useful,
##but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
##GNU General Public License for more details.
##For a copy of the GNU General Public License 
##  see <https://www.gnu.org/licenses/>.
##

library(tidyverse)

## -- Attaching packages ------------------------------------------------------------------------------------------------------------------ tidyverse 1.2.1 --

## v ggplot2 3.2.0     v purrr   0.3.2
## v tibble  2.1.3     v dplyr   0.8.3
## v tidyr   0.8.3     v stringr 1.4.0
## v readr   1.3.1     v forcats 0.4.0

## -- Conflicts --------------------------------------------------------------------------------------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

library(readr)
HwyDB <- read_csv("data/HwyDB_ADTvsEMC.csv")

## Parsed with column specification:
## cols(
##   ADT = col_double(),
##   `EMC Value` = col_double(),
##   `'tEMCQual'1` = col_character(),
##   Parameter = col_character(),
##   NonDetectFlag = col_logical()
## )

##Set local parameters - Non-PAH Contaminants
units <- 'ug/L'
T_Copper = 'Copper, water, unfiltered, recoverable, micrograms per liter'
D_Copper = 'Copper, water, filtered, micrograms per liter'
T_Zinc = 'Zinc'
Nitrite_Nitrate = 'Nitrite plus nitrate, water, unfiltered, milligrams per liter as nitrogen'
TKN = 'Ammonia plus organic nitrogen, water, unfiltered, milligrams per liter as nitrogen (TKN)'
T_Phosphorus = 'Phosphorus, water, unfiltered, milligrams per liter'
Suspended_Solids = 'Solids, suspended, water, milligrams per liter'

#Set PAH Contaminat parameters
Nine_H_Flourine = '9H-Fluorene, water, unfiltered, recoverable, micrograms per liter'
Acenaphthene = 'Acenaphthene, water, unfiltered, recoverable, micrograms per liter'
Acenaphthylene = 'Acenaphthylene, water, unfiltered, recoverable, micrograms per liter'
Anthracene = 'Anthracene, water, unfiltered, recoverable, micrograms per liter'
Anthraquinone = 'Anthraquinone, water, unfiltered, recoverable, micrograms per liter'
Benzo_a_anthracene = 'Benzo[a]anthracene, water, unfiltered, recoverable, micrograms per liter'
Benzo_a_pyrene = 'Benzo[a]pyrene, water, unfiltered, recoverable, micrograms per liter'
Benzo_b_fluoranthene = 'Benzo[b]fluoranthene, water, unfiltered, recoverable, micrograms per liter'
Benzo_k_fluoranthene = 'Benzo[k]fluoranthene, water, unfiltered, recoverable, micrograms per liter'
Chrysene = 'Chrysene, water, unfiltered, recoverable, micrograms per liter'
Dibenzo_a_h_anthracene = 'Dibenzo[a,h]anthracene, water, unfiltered, recoverable, micrograms per liter'
Flouranthene = 'Fluoranthene, water, unfiltered, recoverable, micrograms per liter'
Indeno_1_2_3_cd_pyrene = 'Indeno[1,2,3-cd]pyrene, water, unfiltered, recoverable, micrograms per liter'
Isoquinoline = 'Isoquinoline, water, unfiltered, recoverable, micrograms per liter'
Naphthalene = 'Naphthalene, water, unfiltered, recoverable, micrograms per liter'
Perylene = 'Perylene, water, unfiltered, recoverable, micrograms per liter'
Phenanthrene = 'Phenanthrene, water, unfiltered, recoverable, micrograms per liter'
Pyrene = 'Pyrene, water, unfiltered, recoverable, micrograms per liter'

#library(dplyr)
#Extract relevant data from overall data file

data_T_Copper <- HwyDB %>% 
  dplyr::filter(Parameter == T_Copper)
data_D_Copper <- HwyDB %>% 
  dplyr::filter(Parameter == D_Copper)
data_T_Zinc <- HwyDB %>% 
  dplyr::filter(Parameter == T_Zinc)
data_Nitrite_Nitrate <- HwyDB %>% 
  dplyr::filter(Parameter == Nitrite_Nitrate)
data_TKN <- HwyDB %>% 
  dplyr::filter(Parameter == TKN)
data_T_Phosphorus <- HwyDB %>% 
  dplyr::filter(Parameter == T_Phosphorus)
data_Suspended_Solids <- HwyDB %>% 
  dplyr::filter(Parameter == Suspended_Solids)

data_Nine_H_Flourine <- HwyDB %>% 
  dplyr::filter(Parameter == Nine_H_Flourine)
data_Acenaphthene <- HwyDB %>% 
  dplyr::filter(Parameter == Acenaphthene)
data_Acenaphthylene <- HwyDB %>% 
  dplyr::filter(Parameter == Acenaphthylene)
data_Anthracene <- HwyDB %>% 
  dplyr::filter(Parameter == Anthracene)
data_Anthraquinone <- HwyDB %>% 
  dplyr::filter(Parameter == Anthraquinone)
data_Benzo_a_anthracene <- HwyDB %>% 
  dplyr::filter(Parameter == Benzo_a_anthracene)
data_Benzo_a_pyrene <- HwyDB %>% 
  dplyr::filter(Parameter == Benzo_a_pyrene)
data_Benzo_b_fluoranthene <- HwyDB %>% 
  dplyr::filter(Parameter == Benzo_b_fluoranthene)
data_Benzo_k_fluoranthene <- HwyDB %>% 
  dplyr::filter(Parameter == Benzo_k_fluoranthene)
data_Chrysene <- HwyDB %>% 
  dplyr::filter(Parameter == Chrysene)
data_Dibenzo_a_h_anthracene <- HwyDB %>% 
  dplyr::filter(Parameter == Dibenzo_a_h_anthracene)
data_Flouranthene <- HwyDB %>% 
  dplyr::filter(Parameter == Flouranthene)
data_Indeno_1_2_3_cd_pyrene <- HwyDB %>% 
  dplyr::filter(Parameter == Indeno_1_2_3_cd_pyrene)
data_Isoquinoline <- HwyDB %>% 
  dplyr::filter(Parameter == Isoquinoline)
data_Naphthalene <- HwyDB %>% 
  dplyr::filter(Parameter == Naphthalene)
data_Perylene <- HwyDB %>% 
  dplyr::filter(Parameter == Perylene)
data_Phenanthrene <- HwyDB %>% 
  dplyr::filter(Parameter == Phenanthrene)
data_Pyrene <- HwyDB %>% 
  dplyr::filter(Parameter == Pyrene)

#Select lower threshold for ADT data
lowerBound <- 5000

#Filter data below ADT Threshold
data_T_Copper <- subset(data_T_Copper, ADT > lowerBound)
data_D_Copper <- subset(data_D_Copper, ADT > lowerBound)
data_T_Zinc <- subset(data_T_Zinc, ADT > lowerBound)
data_Nitrite_Nitrate <- subset(data_Nitrite_Nitrate, ADT > lowerBound)
data_TKN <- subset(data_TKN, ADT > lowerBound)
data_T_Phosphorus <- subset(data_T_Phosphorus, ADT > lowerBound)
data_Suspended_Solids <- subset(data_Suspended_Solids, ADT > lowerBound)
data_Nine_H_Flourine <- subset(data_Nine_H_Flourine, ADT > lowerBound)
data_Acenaphthene <- subset(data_Acenaphthene, ADT > lowerBound)
data_Acenaphthylene <- subset(data_Nine_H_Flourine, ADT > lowerBound)
data_Anthracene <- subset(data_Acenaphthylene, ADT > lowerBound)
data_Anthraquinone <- subset(data_Anthraquinone, ADT > lowerBound)
data_Benzo_a_anthracene <- subset(data_Benzo_a_anthracene, ADT > lowerBound)
data_Benzo_a_pyrene <- subset(data_Benzo_a_pyrene, ADT > lowerBound)
data_Benzo_b_fluoranthene <- subset(data_Benzo_b_fluoranthene, ADT > lowerBound)
data_Benzo_k_fluoranthene <- subset(data_Benzo_k_fluoranthene, ADT > lowerBound)
data_Chrysene <- subset(data_Nine_H_Flourine, ADT > lowerBound)
data_Dibenzo_a_h_anthracene <- subset(data_Dibenzo_a_h_anthracene, ADT > lowerBound)
data_Flouranthene <- subset(data_Flouranthene, ADT > lowerBound)
data_Indeno_1_2_3_cd_pyrene <- subset(data_Indeno_1_2_3_cd_pyrene, ADT > lowerBound)
data_Isoquinoline <- subset(data_Isoquinoline, ADT > lowerBound)
data_Naphthalene <- subset(data_Naphthalene, ADT > lowerBound)
data_Perylene <- subset(data_Perylene, ADT > lowerBound)
data_Phenanthrene <- subset(data_Phenanthrene, ADT > lowerBound)
data_Pyrene <- subset(data_Pyrene, ADT > lowerBound)

# Total Copper Plot
plot.title <- paste0("Total Copper Runoff ", units)
plot <- ggplot(data_T_Copper, aes(x=`ADT`,y=`EMC Value`)) + 
  geom_point()+
  ggtitle(plot.title) +
  theme(plot.title = element_text(hjust = 0.5))
plot

# Dissolved Copper Plot
plot.title <- paste0("Dissolved Copper Runoff ", units)
plot <- ggplot(data_D_Copper, aes(x=`ADT`,y=`EMC Value`)) + 
  geom_point()+
  ggtitle(plot.title) +
  theme(plot.title = element_text(hjust = 0.5))
plot

# Total Zinc Plot
plot.title <- paste0("Total Zinc Runoff ", units)
plot <- ggplot(data_T_Zinc, aes(x=`ADT`,y=`EMC Value`)) + 
  geom_point()+
  ggtitle(plot.title) +
  theme(plot.title = element_text(hjust = 0.5))
plot

# Nitrite and Nitrate Plot
plot.title <- paste0("Nitrite and Nitrate Runoff ", units)
plot <- ggplot(data_Nitrite_Nitrate, aes(x=`ADT`,y=`EMC Value`)) + 
  geom_point()+
  ggtitle(plot.title) +
  theme(plot.title = element_text(hjust = 0.5))
plot

# TKN Plot
plot.title <- paste0("Total Kjeldahl Nitrogen Runoff ", units)
plot <- ggplot(data_TKN, aes(x=`ADT`,y=`EMC Value`)) + 
  geom_point()+
  ggtitle(plot.title) +
  theme(plot.title = element_text(hjust = 0.5))
plot

# Total Phosphorus Plot
plot.title <- paste0("Total Phosphorus Runoff", units)
plot <- ggplot(data_T_Phosphorus, aes(x=`ADT`,y=`EMC Value`)) + 
  geom_point()+
  ggtitle(plot.title) +
  theme(plot.title = element_text(hjust = 0.5))
plot

#Suspended Solids Plot
plot.title <- paste0("Suspended Solids Runoff ", units)
plot <- ggplot(data_Suspended_Solids, aes(x=`ADT`,y=`EMC Value`)) + 
  geom_point()+
  ggtitle(plot.title) +
  theme(plot.title = element_text(hjust = 0.5))
plot

# Acemaphthene Plot
plot.title <- paste0("Acenaphthene Runoff ", units)
plot <- ggplot(data_Acenaphthene, aes(x=`ADT`,y=`EMC Value`)) + 
  geom_point()+
  ggtitle(plot.title) +
  theme(plot.title = element_text(hjust = 0.5))
plot

# Acenaphthylene Plot
plot.title <- paste0("Acenaphthylene Runoff ", units)
plot <- ggplot(data_Acenaphthylene, aes(x=`ADT`,y=`EMC Value`)) + 
  geom_point()+
  ggtitle(plot.title) +
  theme(plot.title = element_text(hjust = 0.5))
plot

# Anthracene Plot
plot.title <- paste0("Anthracene Runoff ", units)
plot <- ggplot(data_Anthracene, aes(x=`ADT`,y=`EMC Value`)) + 
  geom_point()+
  ggtitle(plot.title) +
  theme(plot.title = element_text(hjust = 0.5))
plot

# Anthraquinone Plot
plot.title <- paste0("Anthraquinone Runoff ", units)
plot <- ggplot(data_Anthraquinone, aes(x=`ADT`,y=`EMC Value`)) + 
  geom_point()+
  ggtitle(plot.title) +
  theme(plot.title = element_text(hjust = 0.5))
plot

# Benzo[a]anthracene Plot
plot.title <- paste0("Benzo[a]anthracene Runoff ", units)
plot <- ggplot(data_Benzo_a_anthracene, aes(x=`ADT`,y=`EMC Value`)) + 
  geom_point()+
  ggtitle(plot.title) +
  theme(plot.title = element_text(hjust = 0.5))
plot

# Benzo[a]pyrene Plot
plot.title <- paste0("Benzo[a]pyrene Runoff ", units)
plot <- ggplot(data_Benzo_a_pyrene, aes(x=`ADT`,y=`EMC Value`)) + 
  geom_point()+
  ggtitle(plot.title) +
  theme(plot.title = element_text(hjust = 0.5))
plot

# Benzo[b]fluoranthene Plot
plot.title <- paste0("Benzo[b]fluoranthene Runoff ", units)
plot <- ggplot(data_Benzo_b_fluoranthene, aes(x=`ADT`,y=`EMC Value`)) + 
  geom_point()+
  ggtitle(plot.title) +
  theme(plot.title = element_text(hjust = 0.5))
plot

# Benzo[k]fluoranthene Plot
plot.title <- paste0("Benzo[k]fluoranthene Runoff ", units)
plot <- ggplot(data_Benzo_k_fluoranthene, aes(x=`ADT`,y=`EMC Value`)) + 
  geom_point()+
  ggtitle(plot.title) +
  theme(plot.title = element_text(hjust = 0.5))
plot

# Chrysene Plot
plot.title <- paste0("Chrysene Runoff ", units)
plot <- ggplot(data_Chrysene, aes(x=`ADT`,y=`EMC Value`)) + 
  geom_point()+
  ggtitle(plot.title) +
  theme(plot.title = element_text(hjust = 0.5))
plot

# Dibenzo[a,h]anthracene Plot
plot.title <- paste0("Dibenzo[a,h]anthracene Runoff", units)
plot <- ggplot(data_Dibenzo_a_h_anthracene, aes(x=`ADT`,y=`EMC Value`)) + 
  geom_point()+
  ggtitle(plot.title) +
  theme(plot.title = element_text(hjust = 0.5))
plot

# Flouranthene Plot
plot.title <- paste0("Flouranthene Runoff", units)
plot <- ggplot(data_Flouranthene, aes(x=`ADT`,y=`EMC Value`)) + 
  geom_point()+
  ggtitle(plot.title) +
  theme(plot.title = element_text(hjust = 0.5))
plot

# Indeno[1,2,3-cd]pyrene Plot
plot.title <- paste0("Indeno[1,2,3-cd]pyrene Runoff ", units)
plot <- ggplot(data_Indeno_1_2_3_cd_pyrene, aes(x=`ADT`,y=`EMC Value`)) + 
  geom_point()+
  ggtitle(plot.title) +
  theme(plot.title = element_text(hjust = 0.5))
plot

# Isoquinoline Plot
plot.title <- paste0("Isoquinoline Runoff ", units)
plot <- ggplot(data_Isoquinoline, aes(x=`ADT`,y=`EMC Value`)) + 
  geom_point()+
  ggtitle(plot.title) +
  theme(plot.title = element_text(hjust = 0.5))
plot

# Naphthalene Plot
plot.title <- paste0("Naphthalene Runoff ", units)
plot <- ggplot(data_Naphthalene, aes(x=`ADT`,y=`EMC Value`)) + 
  geom_point()+
  ggtitle(plot.title) +
  theme(plot.title = element_text(hjust = 0.5))
plot

# 9H-Fluorene
plot.title <- paste0("9H-Fluorene Runoff ", units)
plot <- ggplot(data_Nine_H_Flourine, aes(x=`ADT`,y=`EMC Value`)) + 
  geom_point()+
  ggtitle(plot.title) +
  theme(plot.title = element_text(hjust = 0.5))
plot

# Perylene Plot
plot.title <- paste0("Perylene Runoff ", units)
plot <- ggplot(data_Perylene, aes(x=`ADT`,y=`EMC Value`)) + 
  geom_point()+
  ggtitle(plot.title) +
  theme(plot.title = element_text(hjust = 0.5))
plot

# Phenantrene Plot
plot.title <- paste0("Phenanthrene Runoff ", units)
plot <- ggplot(data_Phenanthrene, aes(x=`ADT`,y=`EMC Value`)) + 
  geom_point()+
  ggtitle(plot.title) +
  theme(plot.title = element_text(hjust = 0.5))
plot

# Pyrene Plot
plot.title <- paste0("Pyrene Runoff ", units)
plot <- ggplot(data_Pyrene, aes(x=`ADT`,y=`EMC Value`)) + 
  geom_point()+
  ggtitle(plot.title) +
  theme(plot.title = element_text(hjust = 0.5))
plot

Note: There is no relevant data for PAH’s: Anthraquinone, Isoquinoline, and Perylene. Therefore, regressions will not be run for these PAH’s.

library (NADA)

## Loading required package: survival

## 
## Attaching package: 'NADA'

## The following object is masked from 'package:stats':
## 
##     cor

with(data_T_Copper, cenxyplot(ADT, FALSE,`EMC Value`, NonDetectFlag, xlab = 'ADT', ylab = 'EMC', main = 'Total Copper Regression'))
reg = with(data_T_Copper, cenken(cen(`EMC Value`, NonDetectFlag) ~ ADT))
lines(reg)
a <- reg$slope
b <- reg$intercept
p <- reg$p
equation <- sprintf("%f * x", a)
equation <- sprintf("%s + %f", equation, b)
legend('topright', legend = c(sprintf("Equation = %s", equation), sprintf("P-value = %f", p)))

with(data_D_Copper, cenxyplot(ADT, FALSE,`EMC Value`, NonDetectFlag, xlab = 'ADT', ylab = 'EMC', main = 'Dissolved Copper Regression'))
reg = with(data_D_Copper, cenken(cen(`EMC Value`, NonDetectFlag) ~ ADT))
lines(reg)
a <- reg$slope
b <- reg$intercept
p <- reg$p
equation <- sprintf("%f * x", a)
equation <- sprintf("%s + %f", equation, b)
legend('topright', legend = c(sprintf("Equation = %s", equation), sprintf("P-value = %f", p)))

# with(data_T_Zinc, cenxyplot(ADT, FALSE,`EMC Value`, NonDetectFlag, xlab = 'ADT', ylab = 'EMC', main = 'Total Zinc Regression'))
# reg = with(data_T_Zinc, cenken(cen(`EMC Value`, NonDetectFlag) ~ ADT))
# lines(reg)

with(data_Nitrite_Nitrate, cenxyplot(ADT, FALSE,`EMC Value`, NonDetectFlag, xlab = 'ADT', ylab = 'EMC', main = 'Nitrate-Nitrite Regression'))
reg = with(data_Nitrite_Nitrate, cenken(cen(`EMC Value`, NonDetectFlag) ~ ADT))
lines(reg)
a <- reg$slope
b <- reg$intercept
p <- reg$p
equation <- sprintf("%f * x", a)
equation <- sprintf("%s + %f", equation, b)
legend('topright', legend = c(sprintf("Equation = %s", equation), sprintf("P-value = %f", p)))

with(data_TKN, cenxyplot(ADT, FALSE,`EMC Value`, NonDetectFlag, xlab = 'ADT', ylab = 'EMC', main = 'TKN Regression'))
reg = with(data_TKN, cenken(cen(`EMC Value`, NonDetectFlag) ~ ADT))
lines(reg)
a <- reg$slope
b <- reg$intercept
p <- reg$p
equation <- sprintf("%f * x", a)
equation <- sprintf("%s + %f", equation, b)
legend('topright', legend = c(sprintf("Equation = %s", equation), sprintf("P-value = %f", p)))

with(data_T_Phosphorus, cenxyplot(ADT, FALSE,`EMC Value`, NonDetectFlag, xlab = 'ADT', ylab = 'EMC', main = 'Total Phosphorus Regression'))
reg = with(data_T_Phosphorus, cenken(cen(`EMC Value`, NonDetectFlag) ~ ADT))
lines(reg)
a <- reg$slope
b <- reg$intercept
p <- reg$p
equation <- sprintf("%f * x", a)
equation <- sprintf("%s + %f", equation, b)
legend('topright', legend = c(sprintf("Equation = %s", equation), sprintf("P-value = %f", p)))

with(data_Suspended_Solids, cenxyplot(ADT, FALSE,`EMC Value`, NonDetectFlag, xlab = 'ADT', ylab = 'EMC', main = 'Suspended Solids Regression'))
reg = with(data_Suspended_Solids, cenken(cen(`EMC Value`, NonDetectFlag) ~ ADT))
lines(reg)
a <- reg$slope
b <- reg$intercept
p <- reg$p
equation <- sprintf("%f * x", a)
equation <- sprintf("%s + %f", equation, b)
legend('topright', legend = c(sprintf("Equation = %s", equation), sprintf("P-value = %f", p)))

with(data_Anthracene, cenxyplot(ADT, FALSE,`EMC Value`, NonDetectFlag, xlab = 'ADT', ylab = 'EMC', main = 'Anthracene Regression'))
reg = with(data_Anthracene, cenken(cen(`EMC Value`, NonDetectFlag) ~ ADT))
lines(reg)
a <- reg$slope
b <- reg$intercept
p <- reg$p
equation <- sprintf("%f * x", a)
equation <- sprintf("%s + %f", equation, b)
legend('topright', legend = c(sprintf("Equation = %s", equation), sprintf("P-value = %f", p)))

with(data_Acenaphthene, cenxyplot(ADT, FALSE,`EMC Value`, NonDetectFlag, xlab = 'ADT', ylab = 'EMC', main = 'Acenaphthene Regression'))
reg = with(data_Acenaphthene, cenken(cen(`EMC Value`, NonDetectFlag) ~ ADT))
lines(reg)
a <- reg$slope
b <- reg$intercept
p <- reg$p
equation <- sprintf("%f * x", a)
equation <- sprintf("%s + %f", equation, b)
legend('topright', legend = c(sprintf("Equation = %s", equation), sprintf("P-value = %f", p)))

with(data_Acenaphthylene, cenxyplot(ADT, FALSE,`EMC Value`, NonDetectFlag, xlab = 'ADT', ylab = 'EMC', main = 'Acenaphthylene Regression'))
reg = with(data_Acenaphthylene, cenken(cen(`EMC Value`, NonDetectFlag) ~ ADT))
lines(reg)
a <- reg$slope
b <- reg$intercept
p <- reg$p
equation <- sprintf("%f * x", a)
equation <- sprintf("%s + %f", equation, b)
legend('topright', legend = c(sprintf("Equation = %s", equation), sprintf("P-value = %f", p)))

# with(Benzo_a_anthracene, cenxyplot(ADT, FALSE,`EMC Value`, NonDetectFlag, xlab = 'ADT', ylab = 'EMC', main = 'Benzo[a]anthracene Regression'))
# reg = with(Benzo_a_anthracene, cenken(cen(`EMC Value`, NonDetectFlag) ~ ADT))
# lines(reg)

with(data_Benzo_a_pyrene, cenxyplot(ADT, FALSE,`EMC Value`, NonDetectFlag, xlab = 'ADT', ylab = 'EMC', main = 'Benzo[a]pyrene Regression'))
reg = with(data_Benzo_a_pyrene, cenken(cen(`EMC Value`, NonDetectFlag) ~ ADT))
lines(reg)
a <- reg$slope
b <- reg$intercept
p <- reg$p
equation <- sprintf("%f * x", a)
equation <- sprintf("%s + %f", equation, b)
legend('topright', legend = c(sprintf("Equation = %s", equation), sprintf("P-value = %f", p)))

with(data_Benzo_b_fluoranthene, cenxyplot(ADT, FALSE,`EMC Value`, NonDetectFlag, xlab = 'ADT', ylab = 'EMC', main = 'Benzo[b]fluoranthene Regression'))
reg = with(data_Benzo_b_fluoranthene, cenken(cen(`EMC Value`, NonDetectFlag) ~ ADT))
lines(reg)
a <- reg$slope
b <- reg$intercept
p <- reg$p
equation <- sprintf("%f * x", a)
equation <- sprintf("%s + %f", equation, b)
legend('topright', legend = c(sprintf("Equation = %s", equation), sprintf("P-value = %f", p)))

with(data_Benzo_k_fluoranthene, cenxyplot(ADT, FALSE,`EMC Value`, NonDetectFlag, xlab = 'ADT', ylab = 'EMC', main = 'Benzo[k]fluoranthene Regression'))
reg = with(data_Benzo_k_fluoranthene, cenken(cen(`EMC Value`, NonDetectFlag) ~ ADT))
lines(reg)
a <- reg$slope
b <- reg$intercept
p <- reg$p
equation <- sprintf("%f * x", a)
equation <- sprintf("%s + %f", equation, b)
legend('topright', legend = c(sprintf("Equation = %s", equation), sprintf("P-value = %f", p)))

with(data_Chrysene, cenxyplot(ADT, FALSE,`EMC Value`, NonDetectFlag, xlab = 'ADT', ylab = 'EMC', main = 'Chrysene Regression'))
reg = with(data_Chrysene, cenken(cen(`EMC Value`, NonDetectFlag) ~ ADT))
lines(reg)
a <- reg$slope
b <- reg$intercept
p <- reg$p
equation <- sprintf("%f * x", a)
equation <- sprintf("%s + %f", equation, b)
legend('topright', legend = c(sprintf("Equation = %s", equation), sprintf("P-value = %f", p)))

with(data_Dibenzo_a_h_anthracene, cenxyplot(ADT, FALSE,`EMC Value`, NonDetectFlag, xlab = 'ADT', ylab = 'EMC', main = 'Dibenzo[a,h]anthracene Regression'))
reg = with(data_Dibenzo_a_h_anthracene, cenken(cen(`EMC Value`, NonDetectFlag) ~ ADT))
lines(reg)
a <- reg$slope
b <- reg$intercept
p <- reg$p
equation <- sprintf("%f * x", a)
equation <- sprintf("%s + %f", equation, b)
legend('topright', legend = c(sprintf("Equation = %s", equation), sprintf("P-value = %f", p)))

with(data_Flouranthene, cenxyplot(ADT, FALSE,`EMC Value`, NonDetectFlag, xlab = 'ADT', ylab = 'EMC', main = 'Fluoranthene Regression'))
reg = with(data_Flouranthene, cenken(cen(`EMC Value`, NonDetectFlag) ~ ADT))
lines(reg)
a <- reg$slope
b <- reg$intercept
p <- reg$p
equation <- sprintf("%f * x", a)
equation <- sprintf("%s + %f", equation, b)
legend('topright', legend = c(sprintf("Equation = %s", equation), sprintf("P-value = %f", p)))

with(data_Indeno_1_2_3_cd_pyrene, cenxyplot(ADT, FALSE,`EMC Value`, NonDetectFlag, xlab = 'ADT', ylab = 'EMC', main = 'Indeno[1,2,3-cd]pyrene Regression'))
reg = with(data_Indeno_1_2_3_cd_pyrene, cenken(cen(`EMC Value`, NonDetectFlag) ~ ADT))
lines(reg)
a <- reg$slope
b <- reg$intercept
p <- reg$p
equation <- sprintf("%f * x", a)
equation <- sprintf("%s + %f", equation, b)
legend('topright', legend = c(sprintf("Equation = %s", equation), sprintf("P-value = %f", p)))

with(data_Naphthalene, cenxyplot(ADT, FALSE,`EMC Value`, NonDetectFlag, xlab = 'ADT', ylab = 'EMC', main = 'Napthalene Regression'))
reg = with(data_Naphthalene, cenken(cen(`EMC Value`, NonDetectFlag) ~ ADT))
lines(reg)
a <- reg$slope
b <- reg$intercept
p <- reg$p
equation <- sprintf("%f * x", a)
equation <- sprintf("%s + %f", equation, b)
legend('topright', legend = c(sprintf("Equation = %s", equation), sprintf("P-value = %f", p)))

with(data_Nine_H_Flourine, cenxyplot(ADT, FALSE,`EMC Value`, NonDetectFlag, xlab = 'ADT', ylab = 'EMC', main = '9H-Fluorine Regression'))
reg = with(data_Nine_H_Flourine, cenken(cen(`EMC Value`, NonDetectFlag) ~ ADT))
lines(reg)
a <- reg$slope
b <- reg$intercept
p <- reg$p
equation <- sprintf("%f * x", a)
equation <- sprintf("%s + %f", equation, b)
legend('topright', legend = c(sprintf("Equation = %s", equation), sprintf("P-value = %f", p)))

with(data_Phenanthrene, cenxyplot(ADT, FALSE,`EMC Value`, NonDetectFlag, xlab = 'ADT', ylab = 'EMC', main = 'Phenanthrene Regression'))
reg = with(data_Phenanthrene, cenken(cen(`EMC Value`, NonDetectFlag) ~ ADT))
lines(reg)
a <- reg$slope
b <- reg$intercept
p <- reg$p
equation <- sprintf("%f * x", a)
equation <- sprintf("%s + %f", equation, b)
legend('topright', legend = c(sprintf("Equation = %s", equation), sprintf("P-value = %f", p)))

with(data_Pyrene, cenxyplot(ADT, FALSE,`EMC Value`, NonDetectFlag, xlab = 'ADT', ylab = 'EMC', main = 'Pyrene Regression'))
reg = with(data_Pyrene, cenken(cen(`EMC Value`, NonDetectFlag) ~ ADT))
lines(reg)
a <- reg$slope
b <- reg$intercept
p <- reg$p
equation <- sprintf("%f * x", a)
equation <- sprintf("%s + %f", equation, b)
legend('topright', legend = c(sprintf("Equation = %s", equation), sprintf("P-value = %f", p)))

Note that the regressions were only significant for the following contaminants: Total Copper, Dissolved Copper, Nitrite-Nitrate, TKN, and Suspended Solids. As a result, quantile regressions will be run for those contaminants. For the other contaminants, the median EMC value will be extracted and assumed for ADT ranges.

# mycenfit = cenfit(data_T_Zinc$`EMC Value`, data_T_Zinc$NonDetectFlag)
# Value_T_Zinc = median(mycenfit)
# print(Value_T_Zinc)

Value_Suspended_Solids = median(cenfit(data_Suspended_Solids$`EMC Value`, data_Suspended_Solids$NonDetectFlag))
print('Value_Suspended Solids')

## [1] "Value_Suspended Solids"

print(Value_Suspended_Solids)

## [1] 72

Value_Acenaphthene = median(cenfit(data_Acenaphthene$`EMC Value`, data_Acenaphthene$NonDetectFlag))
print('Value_Acenaphthene')

## [1] "Value_Acenaphthene"

print(Value_Acenaphthene)

## [1] NA

Value_Acenaphthylene = median(cenfit(data_Acenaphthylene$`EMC Value`, data_Acenaphthylene$NonDetectFlag))
print('Value_Acenaphthylene')

## [1] "Value_Acenaphthylene"

print(Value_Acenaphthylene)

## [1] NA

Value_Anthracene = median(cenfit(data_Anthracene$`EMC Value`, data_Anthracene$NonDetectFlag))
print('Value_Anthracene')

## [1] "Value_Anthracene"

print(Value_Anthracene)

## [1] NA

Value_Benzo_a_anthracene = median(cenfit(data_Benzo_a_anthracene$`EMC Value`, data_Benzo_a_anthracene$NonDetectFlag))
print('Value_Benzo_a_anthracene')

## [1] "Value_Benzo_a_anthracene"

print(Value_Benzo_a_anthracene)

## [1] 0.02

Value_Benzo_a_pyrene = median(cenfit(data_Benzo_a_pyrene$`EMC Value`, data_Benzo_a_pyrene$NonDetectFlag))
print('Value_Benzo_a_pyrene')

## [1] "Value_Benzo_a_pyrene"

print(Value_Benzo_a_pyrene)

## [1] 0.03

Value_Benzo_b_fluoranthene = median(cenfit(data_Benzo_b_fluoranthene$`EMC Value`, data_Benzo_b_fluoranthene$NonDetectFlag))
print('Value_Benzo_b_fluoranthene')

## [1] "Value_Benzo_b_fluoranthene"

print(Value_Benzo_b_fluoranthene)

## [1] 0.04

Value_Benzo_k_fluoranthene = median(cenfit(data_Benzo_k_fluoranthene$`EMC Value`, data_Benzo_k_fluoranthene$NonDetectFlag))
print('Value_Benzo_k_fluoranthene')

## [1] "Value_Benzo_k_fluoranthene"

print(Value_Benzo_k_fluoranthene)

## [1] 0.02

Value_Chrysene = median(cenfit(data_Chrysene$`EMC Value`, data_Chrysene$NonDetectFlag))
print('Value_Chrysene')

## [1] "Value_Chrysene"

print(Value_Chrysene)

## [1] NA

Value_Dibenzo_a_h_anthracene = median(cenfit(data_Dibenzo_a_h_anthracene$`EMC Value`, data_Dibenzo_a_h_anthracene$NonDetectFlag))
print('Value_Dibenzo_a_h_anthracene')

## [1] "Value_Dibenzo_a_h_anthracene"

print(Value_Dibenzo_a_h_anthracene)

## [1] 0.0091

Value_Fluoranthene = median(cenfit(data_Flouranthene$`EMC Value`, data_Flouranthene$NonDetectFlag))
print('Value_Fluoranthene')

## [1] "Value_Fluoranthene"

print(Value_Fluoranthene)

## [1] 0.07

Value_Indeno_1_2_3_cd_pyrene = median(cenfit(data_Indeno_1_2_3_cd_pyrene$`EMC Value`, data_Indeno_1_2_3_cd_pyrene$NonDetectFlag))
print('Value_Indeno_1_2_3_cd_pyrene')

## [1] "Value_Indeno_1_2_3_cd_pyrene"

print(Value_Indeno_1_2_3_cd_pyrene)

## [1] 0.03

Value_Napthalene = median(cenfit(data_Naphthalene$`EMC Value`, data_Naphthalene$NonDetectFlag))
print('Value_Napthalene')

## [1] "Value_Napthalene"

print(Value_Napthalene)

## [1] 0.02

Value_Nine_H_Fluorine = median(cenfit(data_Nine_H_Flourine$`EMC Value`, data_Nine_H_Flourine$NonDetectFlag))
print('Value_Nine_H_Fluorine')

## [1] "Value_Nine_H_Fluorine"

print(Value_Nine_H_Fluorine)

## [1] NA

Value_Phenanthrene = median(cenfit(data_Phenanthrene$`EMC Value`, data_Phenanthrene$NonDetectFlag))
print('Value_Phenanthrene')

## [1] "Value_Phenanthrene"

print(Value_Phenanthrene)

## [1] 0.05

Value_Pyrene = median(cenfit(data_Pyrene$`EMC Value`, data_Pyrene$NonDetectFlag))
print('Value_Pyrene')

## [1] "Value_Pyrene"

print(Value_Pyrene)

## [1] 0.12

# 
library(quantreg)

## Loading required package: SparseM

## 
## Attaching package: 'SparseM'

## The following object is masked from 'package:base':
## 
##     backsolve

## 
## Attaching package: 'quantreg'

## The following object is masked from 'package:survival':
## 
##     untangle.specials

library(MASS)

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:dplyr':
## 
##     select

library(leaps)

# Quantile regression for Total Copper
obs = (data_T_Copper$`EMC Value`)
censored = data_T_Copper$NonDetectFlag
KM = Cen(obs, censored, type = 'left')#try without groups first
data_T_Copper$KM = KM@Surv

cen.quan.reg <- crq(KM ~ ADT, data = data_T_Copper, method = "Portnoy")
reg.summary <- summary(cen.quan.reg, tau = seq(25, 75, by=25)/100)
print(reg.summary)

## 
## tau: [1] 0.25
## 
## Coefficients:
##             Value   Lower Bd Upper Bd Std Error T Value Pr(>|t|)
## (Intercept) 4.97321 2.47274  9.65064  1.83113   2.71592 0.00661 
## ADT         0.00015 0.00010  0.00021  0.00003   5.45883 0.00000 
## 
## tau: [1] 0.5
## 
## Coefficients:
##             Value    Lower Bd Upper Bd Std Error T Value  Pr(>|t|)
## (Intercept) 12.66238 10.23962 15.40575  1.31792   9.60788  0.00000
## ADT          0.00021  0.00017  0.00025  0.00002   9.96358  0.00000
## 
## tau: [1] 0.75
## 
## Coefficients:
##             Value    Lower Bd Upper Bd Std Error T Value  Pr(>|t|)
## (Intercept) 24.67275 21.71428 29.99142  2.11156  11.68463  0.00000
## ADT          0.00031  0.00028  0.00036  0.00002  13.72152  0.00000

## [[1]]
## [[1]]$tau
## [1] 0.25
## 
## [[1]]$coefficients
##                    Value     Lower Bd     Upper Bd    Std Error  T Value
## (Intercept) 4.9732050410 2.472743e+00 9.6506355145 1.831129e+00 2.715923
## ADT         0.0001525388 9.979672e-05 0.0002093332 2.794349e-05 5.458833
##                 Pr(>|t|)
## (Intercept) 6.609123e-03
## ADT         4.792741e-08
## 
## [[1]]$NAs
## [1] 0
## 
## [[1]]$cov
##               [,1]          [,2]
## [1,]  2.150292e+00 -2.936753e-05
## [2,] -2.936753e-05  5.684182e-10
## 
## [[1]]$Brep
## [1] 100
## 
## [[1]]$bmethod
## [1] "jack"
## 
## 
## [[2]]
## [[2]]$tau
## [1] 0.5
## 
## [[2]]$coefficients
##                    Value     Lower Bd     Upper Bd    Std Error  T Value
## (Intercept) 1.266238e+01 1.023962e+01 15.405752988 1.317916e+00 9.607885
## ADT         2.102544e-04 1.660486e-04  0.000248768 2.110229e-05 9.963580
##             Pr(>|t|)
## (Intercept)        0
## ADT                0
## 
## [[2]]$NAs
## [1] 0
## 
## [[2]]$cov
##               [,1]          [,2]
## [1,]  2.378592e+00 -2.341408e-05
## [2,] -2.341408e-05  3.581638e-10
## 
## [[2]]$Brep
## [1] 100
## 
## [[2]]$bmethod
## [1] "jack"
## 
## 
## [[3]]
## [[3]]$tau
## [1] 0.75
## 
## [[3]]$coefficients
##                    Value     Lower Bd     Upper Bd    Std Error  T Value
## (Intercept) 2.467275e+01 2.171428e+01 2.999142e+01 2.111555e+00 11.68463
## ADT         3.126856e-04 2.753287e-04 3.646559e-04 2.278797e-05 13.72152
##             Pr(>|t|)
## (Intercept)        0
## ADT                0
## 
## [[3]]$NAs
## [1] 0
## 
## [[3]]$cov
##               [,1]          [,2]
## [1,]  1.084265e+01 -5.425886e-05
## [2,] -5.425886e-05  7.361972e-10
## 
## [[3]]$Brep
## [1] 100
## 
## [[3]]$bmethod
## [1] "jack"

seq <- c(1:3)
B0 <- c()

for (i in seq(1,3) ){
  B0[[i]] = reg.summary[i][[1]]$coefficients[1]
}
  
B0 <- reg.summary[2][[1]]$coefficients[1]
B1 <- reg.summary[2][[1]]$coefficients[2]
plot(data_T_Copper$ADT, (data_T_Copper$`EMC Value`), log = 'y', main = "Total Copper Regression", xlab = "ADT", ylab = "EMC")+abline(B0,B1, untf = TRUE)

## integer(0)

# Quantile regression for Dissolved Copper
obs = (data_D_Copper$`EMC Value`)
censored = data_D_Copper$NonDetectFlag
KM = Cen(obs, censored, type = 'left')#try without groups first
data_D_Copper$KM = KM@Surv

cen.quan.reg <- crq(KM ~ ADT, data = data_D_Copper, method = "Portnoy")
reg.summary <- summary(cen.quan.reg, tau = seq(25, 75, by=25)/100)
print(reg.summary)

## 
## tau: [1] 0.25
## 
## Coefficients:
##             Value   Lower Bd Upper Bd Std Error T Value Pr(>|t|)
## (Intercept) 3.36134 2.49973  4.71926  0.56621   5.93652 0.00000 
## ADT         0.00004 0.00003  0.00006  0.00001   4.04201 0.00005 
## 
## tau: [1] 0.5
## 
## Coefficients:
##             Value   Lower Bd Upper Bd Std Error T Value Pr(>|t|)
## (Intercept) 5.08611 4.07999  6.10605  0.51686   9.84040 0.00000 
## ADT         0.00008 0.00006  0.00010  0.00001   7.34124 0.00000 
## 
## tau: [1] 0.75
## 
## Coefficients:
##             Value    Lower Bd Upper Bd Std Error T Value  Pr(>|t|)
## (Intercept)  7.62817  5.98976  8.43476  0.62374  12.22980  0.00000
## ADT          0.00015  0.00010  0.00018  0.00002   7.13739  0.00000

## [[1]]
## [[1]]$tau
## [1] 0.25
## 
## [[1]]$coefficients
##                    Value     Lower Bd     Upper Bd    Std Error  T Value
## (Intercept) 3.361345e+00 2.499735e+00 4.719256e+00 5.662149e-01 5.936517
## ADT         3.681546e-05 2.581179e-05 6.151532e-05 9.108210e-06 4.042008
##                 Pr(>|t|)
## (Intercept) 2.911402e-09
## ADT         5.299546e-05
## 
## [[1]]$NAs
## [1] 0
## 
## [[1]]$cov
##               [,1]          [,2]
## [1,]  2.213650e-01 -3.144709e-06
## [2,] -3.144709e-06  7.923354e-11
## 
## [[1]]$Brep
## [1] 100
## 
## [[1]]$bmethod
## [1] "jack"
## 
## 
## [[2]]
## [[2]]$tau
## [1] 0.5
## 
## [[2]]$coefficients
##                    Value     Lower Bd     Upper Bd    Std Error  T Value
## (Intercept) 5.086112e+00 4.0799933767 6.1060496319 5.168606e-01 9.840395
## ADT         8.080929e-05 0.0000596888 0.0001028377 1.100758e-05 7.341240
##                 Pr(>|t|)
## (Intercept) 0.000000e+00
## ADT         2.116085e-13
## 
## [[2]]$NAs
## [1] 0
## 
## [[2]]$cov
##               [,1]          [,2]
## [1,]  2.328463e-01 -4.676212e-06
## [2,] -4.676212e-06  1.421323e-10
## 
## [[2]]$Brep
## [1] 100
## 
## [[2]]$bmethod
## [1] "jack"
## 
## 
## [[3]]
## [[3]]$tau
## [1] 0.75
## 
## [[3]]$coefficients
##                    Value     Lower Bd     Upper Bd    Std Error   T Value
## (Intercept) 7.6281701838 5.989757e+00 8.4347573732 6.237362e-01 12.229803
## ADT         0.0001489645 9.799951e-05 0.0001798124 2.087101e-05  7.137391
##                 Pr(>|t|)
## (Intercept) 0.000000e+00
## ADT         9.512391e-13
## 
## [[3]]$NAs
## [1] 0
## 
## [[3]]$cov
##               [,1]          [,2]
## [1,]  3.874989e-01 -7.013586e-06
## [2,] -7.013586e-06  3.133370e-10
## 
## [[3]]$Brep
## [1] 100
## 
## [[3]]$bmethod
## [1] "jack"

seq <- c(1:3)
B0 <- c()

for (i in seq(1,3) ){
  B0[[i]] = reg.summary[i][[1]]$coefficients[1]
}
  
B0 <- reg.summary[2][[1]]$coefficients[1]
B1 <- reg.summary[2][[1]]$coefficients[2]
plot(data_D_Copper$ADT, (data_D_Copper$`EMC Value`), log = 'y', main = "Dissolved Copper Regression", xlab = "ADT", ylab = "EMC")+abline(B0,B1, untf = TRUE)

## integer(0)

# # Quantile regression for Total Zinc
# obs = (data_T_Zinc$`EMC Value`)
# censored = data_T_Zinc$NonDetectFlag
# KM = Cen(obs, censored, type = 'left')#try without groups first
# data_T_Zinc$KM = KM@Surv
# 
# cen.quan.reg <- crq(KM ~ ADT, data = data_T_Zinc, method = "Portnoy")
# reg.summary <- summary(cen.quan.reg, tau = seq(25, 75, by=25)/100)
# 
# seq <- c(1:3)
# B0 <- c()
# 
# for (i in seq(1,3) ){
#   B0[[i]] = reg.summary[i][[1]]$coefficients[1]
# }
#   
# B0 <- reg.summary[2][[1]]$coefficients[1]
# B1 <- reg.summary[2][[1]]$coefficients[2]
# plot(data_T_Zinc$ADT, (data_T_Zinc$`EMC Value`), log = 'y', main = "Total Zinc Regression", xlab = "ADT", ylab = "EMC")+abline(B0,B1, untf = TRUE)

# Quantile regression for Nitrite-Nitrate
obs = (data_Nitrite_Nitrate$`EMC Value`)
censored = data_Nitrite_Nitrate$NonDetectFlag
KM = Cen(obs, censored, type = 'left')#try without groups first
data_Nitrite_Nitrate$KM = KM@Surv

cen.quan.reg <- crq(KM ~ ADT, data = data_Nitrite_Nitrate, method = "Portnoy")
reg.summary <- summary(cen.quan.reg, tau = seq(25, 75, by=25)/100)
print(reg.summary)

## 
## tau: [1] 0.25
## 
## Coefficients:
##             Value    Lower Bd Upper Bd Std Error T Value  Pr(>|t|)
## (Intercept)  0.34469  0.26026  0.45807  0.05046   6.83092  0.00000
## ADT          0.00000  0.00000  0.00000  0.00000  -0.72600  0.46784
## 
## tau: [1] 0.5
## 
## Coefficients:
##             Value    Lower Bd Upper Bd Std Error T Value  Pr(>|t|)
## (Intercept)  0.52331  0.39425  0.62300  0.05836   8.96755  0.00000
## ADT          0.00000  0.00000  0.00000  0.00000  -0.81998  0.41223
## 
## tau: [1] 0.75
## 
## Coefficients:
##             Value    Lower Bd Upper Bd Std Error T Value  Pr(>|t|)
## (Intercept)  0.74218  0.48742  1.07561  0.15005   4.94617  0.00000
## ADT          0.00000 -0.00001  0.00001  0.00001   0.65637  0.51158

## [[1]]
## [[1]]$tau
## [1] 0.25
## 
## [[1]]$coefficients
##                     Value      Lower Bd     Upper Bd    Std Error
## (Intercept)  3.446932e-01  2.602649e-01 4.580674e-01 5.046075e-02
## ADT         -1.569182e-06 -3.569827e-06 4.902744e-06 2.161410e-06
##                T Value     Pr(>|t|)
## (Intercept)  6.8309180 8.437251e-12
## ADT         -0.7259991 4.678393e-01
## 
## [[1]]$NAs
## [1] 0
## 
## [[1]]$cov
##               [,1]          [,2]
## [1,]  2.306061e-03 -6.720172e-08
## [2,] -6.720172e-08  3.069095e-12
## 
## [[1]]$Brep
## [1] 100
## 
## [[1]]$bmethod
## [1] "jack"
## 
## 
## [[2]]
## [[2]]$tau
## [1] 0.5
## 
## [[2]]$coefficients
##                     Value      Lower Bd     Upper Bd    Std Error
## (Intercept)  5.233131e-01  3.942476e-01 6.230000e-01 5.835628e-02
## ADT         -9.797392e-07 -2.744253e-06 1.939426e-06 1.194838e-06
##                T Value  Pr(>|t|)
## (Intercept)  8.9675543 0.0000000
## ADT         -0.8199766 0.4122294
## 
## [[2]]$NAs
## [1] 0
## 
## [[2]]$cov
##               [,1]          [,2]
## [1,]  5.382136e-03 -1.110617e-07
## [2,] -1.110617e-07  4.140320e-12
## 
## [[2]]$Brep
## [1] 100
## 
## [[2]]$bmethod
## [1] "jack"
## 
## 
## [[3]]
## [[3]]$tau
## [1] 0.75
## 
## [[3]]$coefficients
##                    Value      Lower Bd     Upper Bd    Std Error  T Value
## (Intercept) 7.421832e-01  4.874198e-01 1.075614e+00 1.500522e-01 4.946165
## ADT         4.055834e-06 -1.183005e-05 1.239175e-05 6.179142e-06 0.656375
##                 Pr(>|t|)
## (Intercept) 7.568975e-07
## ADT         5.115829e-01
## 
## [[3]]$NAs
## [1] 0
## 
## [[3]]$cov
##               [,1]          [,2]
## [1,]  2.573993e-02 -1.019547e-06
## [2,] -1.019547e-06  4.840667e-11
## 
## [[3]]$Brep
## [1] 100
## 
## [[3]]$bmethod
## [1] "jack"

seq <- c(1:3)
B0 <- c()

for (i in seq(1,3) ){
  B0[[i]] = reg.summary[i][[1]]$coefficients[1]
}
  
B0 <- reg.summary[2][[1]]$coefficients[1]
B1 <- reg.summary[2][[1]]$coefficients[2]
plot(data_Nitrite_Nitrate$ADT, (data_Nitrite_Nitrate$`EMC Value`), log = 'y', main = "Nitrite-Nitrate Regression", xlab = "ADT", ylab = "EMC")+abline(B0,B1, untf = TRUE)

## integer(0)

# Quantile regression for TKN
obs = (data_TKN$`EMC Value`)
censored = data_TKN$NonDetectFlag
KM = Cen(obs, censored, type = 'left')#try without groups first
data_TKN$KM = KM@Surv

cen.quan.reg <- crq(KM ~ ADT, data = data_TKN, method = "Portnoy")
reg.summary <- summary(cen.quan.reg, tau = seq(25, 75, by=25)/100)
print(reg.summary)

## 
## tau: [1] 0.25
## 
## Coefficients:
##             Value    Lower Bd Upper Bd Std Error T Value  Pr(>|t|)
## (Intercept)  0.67792  0.62165  0.77371  0.03879  17.47651  0.00000
## ADT          0.00000  0.00000  0.00001  0.00000   7.53086  0.00000
## 
## tau: [1] 0.5
## 
## Coefficients:
##             Value    Lower Bd Upper Bd Std Error T Value  Pr(>|t|)
## (Intercept)  1.25265  0.96270  1.35152  0.09919  12.62865  0.00000
## ADT          0.00001  0.00001  0.00001  0.00000   7.21030  0.00000
## 
## tau: [1] 0.75
## 
## Coefficients:
##             Value    Lower Bd Upper Bd Std Error T Value  Pr(>|t|)
## (Intercept)  2.51642  2.16839  3.01926  0.21706  11.59303  0.00000
## ADT          0.00001  0.00000  0.00002  0.00000   2.14569  0.03190

## [[1]]
## [[1]]$tau
## [1] 0.25
## 
## [[1]]$coefficients
##                    Value     Lower Bd     Upper Bd    Std Error   T Value
## (Intercept) 6.779248e-01 6.216549e-01 7.737114e-01 3.879063e-02 17.476511
## ADT         4.910239e-06 3.009770e-06 5.565625e-06 6.520158e-07  7.530858
##                 Pr(>|t|)
## (Intercept) 0.000000e+00
## ADT         5.040413e-14
## 
## [[1]]$NAs
## [1] 0
## 
## [[1]]$cov
##               [,1]          [,2]
## [1,]  1.785407e-03 -2.754425e-08
## [2,] -2.754425e-08  8.026084e-13
## 
## [[1]]$Brep
## [1] 100
## 
## [[1]]$bmethod
## [1] "jack"
## 
## 
## [[2]]
## [[2]]$tau
## [1] 0.5
## 
## [[2]]$coefficients
##                    Value     Lower Bd     Upper Bd    Std Error   T Value
## (Intercept) 1.252652e+00 9.626988e-01 1.351521e+00 9.919127e-02 12.628647
## ADT         7.305679e-06 5.467782e-06 9.439563e-06 1.013228e-06  7.210301
##                 Pr(>|t|)
## (Intercept) 0.000000e+00
## ADT         5.582201e-13
## 
## [[2]]$NAs
## [1] 0
## 
## [[2]]$cov
##               [,1]          [,2]
## [1,]  7.568315e-03 -6.285958e-08
## [2,] -6.285958e-08  1.096119e-12
## 
## [[2]]$Brep
## [1] 100
## 
## [[2]]$bmethod
## [1] "jack"
## 
## 
## [[3]]
## [[3]]$tau
## [1] 0.75
## 
## [[3]]$coefficients
##                    Value     Lower Bd     Upper Bd    Std Error   T Value
## (Intercept) 2.516424e+00 2.168390e+00 3.019263e+00 2.170634e-01 11.593033
## ADT         7.581754e-06 2.774913e-06 1.662589e-05 3.533479e-06  2.145691
##               Pr(>|t|)
## (Intercept) 0.00000000
## ADT         0.03189763
## 
## [[3]]$NAs
## [1] 0
## 
## [[3]]$cov
##               [,1]          [,2]
## [1,]  6.379429e-02 -5.422178e-07
## [2,] -5.422178e-07  9.384135e-12
## 
## [[3]]$Brep
## [1] 100
## 
## [[3]]$bmethod
## [1] "jack"

seq <- c(1:3)
B0 <- c()

for (i in seq(1,3) ){
  B0[[i]] = reg.summary[i][[1]]$coefficients[1]
}
  
B0 <- reg.summary[2][[1]]$coefficients[1]
B1 <- reg.summary[2][[1]]$coefficients[2]
plot(data_TKN$ADT, (data_TKN$`EMC Value`), log = 'y', main = "Total Kjeldahl Nitrogen Regression", xlab = "ADT", ylab = "EMC")+abline(B0,B1, untf = TRUE)

## integer(0)

# Quantile regression for Suspended Solids
obs = (data_Suspended_Solids$`EMC Value`)
censored = data_Suspended_Solids$NonDetectFlag
KM = Cen(obs, censored, type = 'left')#try without groups first
data_Suspended_Solids$KM = KM@Surv

cen.quan.reg <- crq(KM ~ ADT, data = data_Suspended_Solids, method = "Portnoy")
reg.summary <- summary(cen.quan.reg, tau = seq(25, 75, by=25)/100)
print(reg.summary)

## 
## tau: [1] 0.25
## 
## Coefficients:
##             Value    Lower Bd Upper Bd Std Error T Value  Pr(>|t|)
## (Intercept) 22.51298 14.37121 27.33081  3.30608   6.80957  0.00000
## ADT          0.00015  0.00011  0.00021  0.00002   6.02042  0.00000
## 
## tau: [1] 0.5
## 
## Coefficients:
##             Value    Lower Bd Upper Bd Std Error T Value  Pr(>|t|)
## (Intercept) 56.37397 46.50795 65.28205  4.78940  11.77058  0.00000
## ADT          0.00019  0.00010  0.00035  0.00006   2.97588  0.00292
## 
## tau: [1] 0.75
## 
## Coefficients:
##             Value     Lower Bd  Upper Bd  Std Error T Value   Pr(>|t|) 
## (Intercept) 165.50082 135.32035 214.82727  20.28275   8.15968   0.00000
## ADT          -0.00008  -0.00032   0.00006   0.00010  -0.82681   0.40835

## [[1]]
## [[1]]$tau
## [1] 0.25
## 
## [[1]]$coefficients
##                    Value     Lower Bd     Upper Bd    Std Error  T Value
## (Intercept) 2.251298e+01 1.437121e+01 2.733081e+01 3.306079e+00 6.809571
## ADT         1.482503e-04 1.144631e-04 2.109897e-04 2.462458e-05 6.020420
##                 Pr(>|t|)
## (Intercept) 9.789058e-12
## ADT         1.739650e-09
## 
## [[1]]$NAs
## [1] 0
## 
## [[1]]$cov
##               [,1]          [,2]
## [1,] 15.9068741875 -1.162102e-04
## [2,] -0.0001162102  1.119974e-09
## 
## [[1]]$Brep
## [1] 100
## 
## [[1]]$bmethod
## [1] "jack"
## 
## 
## [[2]]
## [[2]]$tau
## [1] 0.5
## 
## [[2]]$coefficients
##                    Value     Lower Bd     Upper Bd    Std Error  T Value
## (Intercept) 56.373969579 4.650795e+01 6.528205e+01 4.789398e+00 11.77058
## ADT          0.000193263 9.577606e-05 3.503485e-04 6.494314e-05  2.97588
##                Pr(>|t|)
## (Intercept) 0.000000000
## ADT         0.002921491
## 
## [[2]]$NAs
## [1] 0
## 
## [[2]]$cov
##               [,1]          [,2]
## [1,] 34.4681759572 -3.180831e-04
## [2,] -0.0003180831  4.302996e-09
## 
## [[2]]$Brep
## [1] 100
## 
## [[2]]$bmethod
## [1] "jack"
## 
## 
## [[3]]
## [[3]]$tau
## [1] 0.75
## 
## [[3]]$coefficients
##                     Value      Lower Bd     Upper Bd    Std Error
## (Intercept)  1.655008e+02  1.353204e+02 2.148273e+02 2.028275e+01
## ADT         -8.111251e-05 -3.245261e-04 6.003219e-05 9.810341e-05
##                T Value     Pr(>|t|)
## (Intercept)  8.1596842 4.440892e-16
## ADT         -0.8268062 4.083469e-01
## 
## [[3]]$NAs
## [1] 0
## 
## [[3]]$cov
##               [,1]          [,2]
## [1,] 355.994048682 -1.811046e-03
## [2,]  -0.001811046  1.171885e-08
## 
## [[3]]$Brep
## [1] 100
## 
## [[3]]$bmethod
## [1] "jack"

seq <- c(1:3)
B0 <- c()

for (i in seq(1,3) ){
  B0[[i]] = reg.summary[i][[1]]$coefficients[1]
}
  
B0 <- reg.summary[2][[1]]$coefficients[1]
B1 <- reg.summary[2][[1]]$coefficients[2]
plot(data_Suspended_Solids$ADT, (data_Suspended_Solids$`EMC Value`), log = 'y', main = "Suspended Solids Regression", xlab = "ADT", ylab = "EMC")+abline(B0,B1, untf = TRUE)

## integer(0)

# # # Quantile regression for Total Phosphorus
# # obs = (data_T_Phosphorus$`EMC Value`)
# # censored = data_T_Phosphorus$NonDetectFlag
# # KM = Cen(obs, censored, type = 'left')#try without groups first
# # data_T_Phosphorus$KM = KM@Surv
# # 
# # cen.quan.reg <- crq(KM ~ ADT, data = data_T_Phosphorus, method = "Portnoy")
# # reg.summary <- summary(cen.quan.reg, tau = seq(25, 75, by=25)/100)
# # print(reg.summary)
# # seq <- c(1:3)
# # B0 <- c()
# # 
# # for (i in seq(1,3) ){
# #   B0[[i]] = reg.summary[i][[1]]$coefficients[1]
# # }
# #   
# # B0 <- reg.summary[2][[1]]$coefficients[1]
# # B1 <- reg.summary[2][[1]]$coefficients[2]
# # plot(data_T_Phosphorus$ADT, (data_T_Phosphorus$`EMC Value`), log = 'y', main = "Total Phosphorus Regression", xlab = "ADT", ylab = "EMC")+abline(B0,B1, untf = TRUE)
# 
# # # Quantile regression for Acenaphthene
# # obs = (data_Acenaphthene$`EMC Value`)
# # censored = data_Acenaphthene$NonDetectFlag
# # KM = Cen(obs, censored, type = 'left')#try without groups first
# # data_Acenaphthene$KM = KM@Surv
# # 
# # 
# # cen.quan.reg <- crq(KM ~ ADT, data = data_Acenaphthene, method = "Portnoy")
# # reg.summary <- summary(cen.quan.reg, tau = seq(25, 75, by=25)/100)
# # 
# # #B0 = vector(,3)
# # seq <- c(1:3)
# # B0 <- c()
# # 
# # #xs <- sort((data_T_Copper$`EMC Value`))
# # for (i in seq(1,3) ){
# #   B0[[i]] = reg.summary[i][[1]]$coefficients[1]
# # }
# #   
# # B0 <- reg.summary[2][[1]]$coefficients[1]
# # B1 <- reg.summary[2][[1]]$coefficients[2]
# # plot(data_Acenaphthene$ADT, (data_Acenaphthene$`EMC Value`), log = 'y')+abline(B0,B1, untf = TRUE)
# 
# # 
# # # Quantile regression for Acenaphthylene
# # 
# # obs = (data_Acenaphthylene$`EMC Value`)
# # censored = data_Acenaphthylene$NonDetectFlag
# # KM = Cen(obs, censored, type = 'left')#try without groups first
# # data_Acenaphthylene$KM = KM@Surv
# # 
# # 
# # cen.quan.reg <- crq(KM ~ ADT, data = data_Acenaphthylene, method = "Portnoy")
# # reg.summary <- summary(cen.quan.reg, tau = seq(25, 75, by=25)/100)
# # 
# # #B0 = vector(,3)
# # seq <- c(1:3)
# # B0 <- c()
# # 
# # #xs <- sort((data_T_Copper$`EMC Value`))
# # for (i in seq(1,3) ){
# #   B0[[i]] = reg.summary[i][[1]]$coefficients[1]
# # }
# #   
# # B0 <- reg.summary[2][[1]]$coefficients[1]
# # B1 <- reg.summary[2][[1]]$coefficients[2]
# # plot(data_Acenaphthylene$ADT, (data_Acenaphthylene$`EMC Value`), log = 'y')+abline(B0,B1, untf = TRUE)
# 
# # # Quantile regression for Anthracene
# # 
# # obs = (data_Anthracene$`EMC Value`)
# # censored = data_Anthracene$NonDetectFlag
# # KM = Cen(obs, censored, type = 'left')#try without groups first
# # data_Anthracene$KM = KM@Surv
# # 
# # 
# # cen.quan.reg <- crq(KM ~ ADT, data = data_Anthracene, method = "Portnoy")
# # reg.summary <- summary(cen.quan.reg, tau = seq(25, 75, by=25)/100)
# # 
# # #B0 = vector(,3)
# # seq <- c(1:3)
# # B0 <- c()
# # 
# # #xs <- sort((data_T_Copper$`EMC Value`))
# # for (i in seq(1,3) ){
# #   B0[[i]] = reg.summary[i][[1]]$coefficients[1]
# # }
# #   
# # B0 <- reg.summary[2][[1]]$coefficients[1]
# # B1 <- reg.summary[2][[1]]$coefficients[2]
# # plot(data_Anthracene$ADT, (data_Anthracene$`EMC Value`), log = 'y')+abline(B0,B1, untf = TRUE)
# 
# # # Quantile regression for Benzo_a_anthracene
# # 
# # obs = (data_Benzo_a_anthracene$`EMC Value`)
# # censored = data_Benzo_a_anthracene$NonDetectFlag
# # KM = Cen(obs, censored, type = 'left')#try without groups first
# # data_Benzo_a_anthracene$KM = KM@Surv
# # 
# # 
# # cen.quan.reg <- crq(KM ~ ADT, data = data_Benzo_a_anthracene, method = "Portnoy")
# # reg.summary <- summary(cen.quan.reg, tau = seq(25, 75, by=25)/100)
# # 
# # #B0 = vector(,3)
# # seq <- c(1:3)
# # B0 <- c()
# # 
# # #xs <- sort((data_T_Copper$`EMC Value`))
# # for (i in seq(1,3) ){
# #   B0[[i]] = reg.summary[i][[1]]$coefficients[1]
# # }
# #   
# # B0 <- reg.summary[2][[1]]$coefficients[1]
# # B1 <- reg.summary[2][[1]]$coefficients[2]
# # plot(data_Benzo_a_anthracene$ADT, (data_Benzo_a_anthracene$`EMC Value`), log = 'y')+abline(B0,B1, untf = TRUE) 
# 
# # Quantile regression for Benzo[a]pyrene
# 
# obs = (data_Benzo_a_pyrene$`EMC Value`)
# censored = data_Benzo_a_pyrene$NonDetectFlag
# KM = Cen(obs, censored, type = 'left')#try without groups first
# data_Benzo_a_pyrene$KM = KM@Surv
# 
# 
# cen.quan.reg <- crq(KM ~ ADT, data = data_Benzo_a_pyrene, method = "Portnoy")
# reg.summary <- summary(cen.quan.reg, tau = seq(25, 75, by=25)/100)
# print(reg.summary)
# 
# #B0 = vector(,3)
# seq <- c(1:3)
# B0 <- c()
# 
# #xs <- sort((data_T_Copper$`EMC Value`))
# for (i in seq(1,3) ){
#   B0[[i]] = reg.summary[i][[1]]$coefficients[1]
# }
# 
# B0 <- reg.summary[2][[1]]$coefficients[1]
# B1 <- reg.summary[2][[1]]$coefficients[2]
# plot(data_Benzo_a_pyrene$ADT, (data_Benzo_a_pyrene$`EMC Value`), log = 'y', main = "Benzo[a]pyrene Regression", xlab = "ADT", ylab = "EMC")+abline(B0,B1, untf = TRUE)
# 
# # Quantile regression for Benzo[b]flouranthene
# 
# obs = (data_Benzo_b_fluoranthene$`EMC Value`)
# censored = data_Benzo_b_fluoranthene$NonDetectFlag
# KM = Cen(obs, censored, type = 'left')#try without groups first
# data_Benzo_b_fluoranthene$KM = KM@Surv
# 
# 
# cen.quan.reg <- crq(KM ~ ADT, data = data_Benzo_b_fluoranthene, method = "Portnoy")
# reg.summary <- summary(cen.quan.reg, tau = seq(25, 75, by=25)/100)
# print(reg.summary)
# 
# seq <- c(1:3)
# B0 <- c()
# 
# for (i in seq(1,3) ){
#   B0[[i]] = reg.summary[i][[1]]$coefficients[1]
# }
# 
# B0 <- reg.summary[2][[1]]$coefficients[1]
# B1 <- reg.summary[2][[1]]$coefficients[2]
# plot(data_Benzo_b_fluoranthene$ADT, (data_Benzo_b_fluoranthene$`EMC Value`), log = 'y', main = "Benzo[b]flouranthene Regression", xlab = "ADT", ylab = "EMC")+abline(B0,B1, untf = TRUE)
# 
# # # Quantile regression for Benzo[k]fluoranthene
# # 
# # obs = (data_Benzo_k_fluoranthene$`EMC Value`)
# # censored = data_Benzo_k_fluoranthene$NonDetectFlag
# # KM = Cen(obs, censored, type = 'left')#try without groups first
# # data_Benzo_k_fluoranthene$KM = KM@Surv
# # 
# # 
# # cen.quan.reg <- crq(KM ~ ADT, data = data_Benzo_k_fluoranthene, method = "Portnoy")
# # reg.summary <- summary(cen.quan.reg, tau = seq(25, 75, by=25)/100)
# # 
# # seq <- c(1:3)
# # B0 <- c()
# # 
# # for (i in seq(1,3) ){
# #   B0[[i]] = reg.summary[i][[1]]$coefficients[1]
# # }
# # 
# # B0 <- reg.summary[2][[1]]$coefficients[1]
# # B1 <- reg.summary[2][[1]]$coefficients[2]
# # plot(data_Benzo_k_fluoranthene$ADT, (data_Benzo_k_fluoranthene$`EMC Value`), log = 'y', main = "Benzo[k]flouranthene Regression", xlab = "ADT", ylab = "EMC")+abline(B0,B1, untf = TRUE)
# 
# # # Quantile regression for Chrysene
# # 
# # obs = (data_Chrysene$`EMC Value`)
# # censored = data_Chrysene$NonDetectFlag
# # KM = Cen(obs, censored, type = 'left')#try without groups first
# # data_Chrysene$KM = KM@Surv
# # 
# # 
# # cen.quan.reg <- crq(KM ~ ADT, data = data_Chrysene, method = "Portnoy")
# # reg.summary <- summary(cen.quan.reg, tau = seq(25, 75, by=25)/100)
# # 
# # #B0 = vector(,3)
# # seq <- c(1:3)
# # B0 <- c()
# # 
# # #xs <- sort((data_T_Copper$`EMC Value`))
# # for (i in seq(1,3) ){
# #   B0[[i]] = reg.summary[i][[1]]$coefficients[1]
# # }
# # 
# # B0 <- reg.summary[2][[1]]$coefficients[1]
# # B1 <- reg.summary[2][[1]]$coefficients[2]
# # plot(data_Chrysene$ADT, (data_Chrysene$`EMC Value`), log = 'y')+abline(B0,B1, untf = TRUE)
# 
# # # Quantile regression for Dibenzo[a,h]anthracene
# # 
# # obs = (data_Dibenzo_a_h_anthracene$`EMC Value`)
# # censored = data_Dibenzo_a_h_anthracene$NonDetectFlag
# # KM = Cen(obs, censored, type = 'left')#try without groups first
# # data_Dibenzo_a_h_anthracene$KM = KM@Surv
# # 
# # 
# # cen.quan.reg <- crq(KM ~ ADT, data = data_Dibenzo_a_h_anthracene, method = "Portnoy")
# # reg.summary <- summary(cen.quan.reg, tau = seq(25, 75, by=25)/100)
# # 
# # #B0 = vector(,3)
# # seq <- c(1:3)
# # B0 <- c()
# # 
# # #xs <- sort((data_T_Copper$`EMC Value`))
# # for (i in seq(1,3) ){
# #   B0[[i]] = reg.summary[i][[1]]$coefficients[1]
# # }
# # 
# # B0 <- reg.summary[2][[1]]$coefficients[1]
# # B1 <- reg.summary[2][[1]]$coefficients[2]
# # plot(data_Dibenzo_a_h_anthracene$ADT, (data_Dibenzo_a_h_anthracene$`EMC Value`), log = 'y')+abline(B0,B1, untf = TRUE)
# 
# # Quantile regression for Fluoranthene
# 
# obs = (data_Flouranthene$`EMC Value`)
# censored = data_Flouranthene$NonDetectFlag
# KM = Cen(obs, censored, type = 'left')#try without groups first
# data_Flouranthene$KM = KM@Surv
# 
# 
# cen.quan.reg <- crq(KM ~ ADT, data = data_Flouranthene, method = "Portnoy")
# reg.summary <- summary(cen.quan.reg, tau = seq(25, 75, by=25)/100)
# print(reg.summary)
# 
# #B0 = vector(,3)
# seq <- c(1:3)
# B0 <- c()
# 
# #xs <- sort((data_T_Copper$`EMC Value`))
# for (i in seq(1,3) ){
#   B0[[i]] = reg.summary[i][[1]]$coefficients[1]
# }
#   
# B0 <- reg.summary[2][[1]]$coefficients[1]
# B1 <- reg.summary[2][[1]]$coefficients[2]
# plot(data_Flouranthene$ADT, (data_Flouranthene$`EMC Value`), log = 'y', main = "Fluoranthene", xlab = "ADT", ylab = "EMC")+abline(B0,B1, untf = TRUE)
# 
# # Quantile regression for Indeno[1,2,3-cd]pyrene
# 
# obs = (data_Indeno_1_2_3_cd_pyrene$`EMC Value`)
# censored = data_Indeno_1_2_3_cd_pyrene$NonDetectFlag
# KM = Cen(obs, censored, type = 'left')#try without groups first
# data_Indeno_1_2_3_cd_pyrene$KM = KM@Surv
# 
# 
# cen.quan.reg <- crq(KM ~ ADT, data = data_Indeno_1_2_3_cd_pyrene, method = "Portnoy")
# reg.summary <- summary(cen.quan.reg, tau = seq(25, 75, by=25)/100)
# print(reg.summary)
# 
# seq <- c(1:3)
# B0 <- c()
# 
# for (i in seq(1,3) ){
#   B0[[i]] = reg.summary[i][[1]]$coefficients[1]
# }
# 
# B0 <- reg.summary[2][[1]]$coefficients[1]
# B1 <- reg.summary[2][[1]]$coefficients[2]
# plot(data_Indeno_1_2_3_cd_pyrene$ADT, (data_Indeno_1_2_3_cd_pyrene$`EMC Value`), log = 'y', main = "Indeno[1,2,3-cd]pyrene", xlab = "ADT", ylab = "EMC")+abline(B0,B1, untf = TRUE)
# 
# # # Quantile regression for Naphthalene
# # 
# # obs = (data_Naphthalene$`EMC Value`)
# # censored = data_Naphthalene$NonDetectFlag
# # KM = Cen(obs, censored, type = 'left')#try without groups first
# # data_Naphthalene$KM = KM@Surv
# # 
# # 
# # cen.quan.reg <- crq(KM ~ ADT, data = data_Naphthalene, method = "Portnoy")
# # reg.summary <- summary(cen.quan.reg, tau = seq(25, 75, by=25)/100)
# # 
# # #B0 = vector(,3)
# # seq <- c(1:3)
# # B0 <- c()
# # 
# # #xs <- sort((data_T_Copper$`EMC Value`))
# # for (i in seq(1,3) ){
# #   B0[[i]] = reg.summary[i][[1]]$coefficients[1]
# # }
# #   
# # B0 <- reg.summary[2][[1]]$coefficients[1]
# # B1 <- reg.summary[2][[1]]$coefficients[2]
# # plot(data_Naphthalene$ADT, (data_Naphthalene$`EMC Value`), log = 'y')+abline(B0,B1, untf = TRUE)
# 
# # # Quantile regression for 9H-Flourine
# # 
# # obs = (data_Nine_H_Flourine$`EMC Value`)
# # censored = data_Nine_H_Flourine$NonDetectFlag
# # KM = Cen(obs, censored, type = 'left')#try without groups first
# # data_Nine_H_Flourine$KM = KM@Surv
# # 
# # 
# # cen.quan.reg <- crq(KM ~ ADT, data = data_Nine_H_Flourine, method = "Portnoy")
# # reg.summary <- summary(cen.quan.reg, tau = seq(25, 75, by=25)/100)
# # 
# # #B0 = vector(,3)
# # seq <- c(1:3)
# # B0 <- c()
# # 
# # #xs <- sort((data_T_Copper$`EMC Value`))
# # for (i in seq(1,3) ){
# #   B0[[i]] = reg.summary[i][[1]]$coefficients[1]
# # }
# #   
# # B0 <- reg.summary[2][[1]]$coefficients[1]
# # B1 <- reg.summary[2][[1]]$coefficients[2]
# # plot(data_Nine_H_Flourine$ADT, (data_Nine_H_Flourine$`EMC Value`), log = 'y')+abline(B0,B1, untf = TRUE)
# 
# # Quantile regression for Phenanthrene
# 
# obs = (data_Phenanthrene$`EMC Value`)
# censored = data_Phenanthrene$NonDetectFlag
# KM = Cen(obs, censored, type = 'left')#try without groups first
# data_Phenanthrene$KM = KM@Surv
# 
# 
# cen.quan.reg <- crq(KM ~ ADT, data = data_Phenanthrene, method = "Portnoy")
# reg.summary <- summary(cen.quan.reg, tau = seq(25, 75, by=25)/100)
# print(reg.summary)
# 
# seq <- c(1:3)
# B0 <- c()
# 
# for (i in seq(1,3) ){
#   B0[[i]] = reg.summary[i][[1]]$coefficients[1]
# }
#   
# B0 <- reg.summary[2][[1]]$coefficients[1]
# B1 <- reg.summary[2][[1]]$coefficients[2]
# plot(data_Phenanthrene$ADT, (data_Phenanthrene$`EMC Value`), log = 'y', main = "Phenanthrene", xlab = "ADT", ylab = "EMC")+abline(B0,B1, untf = TRUE)
# 
# # Quantile regression for Pyrene
# 
# obs = (data_Pyrene$`EMC Value`)
# censored = data_Pyrene$NonDetectFlag
# KM = Cen(obs, censored, type = 'left')#try without groups first
# data_Pyrene$KM = KM@Surv
# 
# 
# cen.quan.reg <- crq(KM ~ ADT, data = data_Pyrene, method = "Portnoy")
# reg.summary <- summary(cen.quan.reg, tau = seq(25, 75, by=25)/100)
# print(reg.summary)
# 
# seq <- c(1:3)
# B0 <- c()
# 
# for (i in seq(1,3) ){
#   B0[[i]] = reg.summary[i][[1]]$coefficients[1]
# }
#   
# B0 <- reg.summary[2][[1]]$coefficients[1]
# B1 <- reg.summary[2][[1]]$coefficients[2]
# plot(data_Pyrene$ADT, (data_Pyrene$`EMC Value`), log = 'y', main = "Pyrene", xlab = "ADT", ylab = "EMC")+abline(B0,B1, untf = TRUE)