---
title: "CalcThemAll_PRM"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{CalcThemAll_PRM}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

## The Pesticide Risk Metric (PRM)

This package provides code for estimating the combined toxicity of measured pesticides using their Species Sensitivity Distributions (SSDs) in a multi-substance, potentially affected fraction approach. The resulting value is an estimate of the percentage of species impacted by the measured pesticides, a value which is referred to in this package as the **Pesticide Risk Metric (PRM)**. This is based on the Queensland Department of Environment and Science Water Quality and Investigation's [method](https://www.publications.qld.gov.au/dataset/method-development-pesticide-risk-metric-baseline-condition-of-waterways-to-gbr/resource/c65858f9-d7ba-4aef-aa4f-e148f950220f) (referred to in this package as the WQI Method) used in the Reef 2050 Water Quality Improvement Plan. The main difference is the ability to add different pesticides to be measured provided they have one of the following SSD distributions:

* Burr_Type_III
* Gamma
* Inverse Weibull
* Log Gumbel
* Log Logistic
* Log Logistic Log Logistic
* Log Normal
* Log Normal Log Normal 

These PRM values can be compared to the three categories of ecosystem condition, as defined in the Australian and New Zealand Guidelines for Fresh & Marine Water Quality seen in the table below.

```{r out.width = '80%', echo = FALSE}
knitr::include_graphics("figures/guidelines_table.png")
```

# How To Calc Them All

## Adding Your Own Pesticides

One of the key functions of this package is **the ability to add your own pesticides** to the risk metric, provided that SSD information is available. This can be done by editing the `pesticide_info` table which is included in this package and is used as a look-up table for many of the later functions. You can add pesticides provided you know the SSD distribution type, scale parameters, shape/location parameters, and weight parameter if required. Depending on how/if LOR replacement is needed relative LOR replacement values for each new pesticide will be required, and if the impact of specific pesticide types is of interest a pesticide type parameter should be provided. The `add_your_own_pesticide()` function appends new pesticide information to the `pesticide_info` table as seen below.

```{r add_your_own_pesticide, echo=TRUE, warning=FALSE}
library(CalcThemAll.PRM)
library(DT)
pesticide_info <- CalcThemAll.PRM::pesticide_info #the original 22 pesticides
datatable(pesticide_info, options = list(pageLength = 10,lengthMenu = c(5, 10, 15, 20), scrollX = T))

pesticide_info <- add_your_own_pesticide(pesticides = "Poison", #adding one new pesticide
                                         relative_LORs = 0.023, pesticide_types = "Poison",
                                         distribution_types = "Log-Normal", scales = 0.09,
                                         shape_locations = 0.014)
datatable(pesticide_info, options = list(pageLength = 10,lengthMenu = c(5, 10, 15, 20), scrollX = T))

pesticide_info <- add_your_own_pesticide(pesticides = #adding multiple new pesticides
                                           c("Poison", "Acid", "Sludge"),
                                         relative_LORs = c(0.03, 0.01, 0.5), 
                                         pesticide_types = c("Ghost", "Bug", "Poison"),
                                         distribution_types = c("Log-Normal", "Log-Logistic
                                                                Log-Logistic", "Burr Type III"),
                                         scales = c(0.3, 0.002, 2),
                                         scale_2s = c(NA, 0.04, NA), 
                                         shape_locations = c(1, 0.07, 3),
                                         shape_location_2s = c(NA, 0.14, 2.3),
                                         weights = c(NA, 0.08, NA))
datatable(pesticide_info, options = list(pageLength = 10,lengthMenu = c(5, 10, 15, 20), scrollX = T))

```

It is **not necessary to remove unmeasured pesticides** from the `pesticide_info` table due to the independent action assumption in the methods, however **if they are not removed they will still require an empty column** in the concentration data set and **which pesticides are used in the metric should be clearly documented** for transparency. 

## Treating Limit Of Reporting Values (LORs)

The treatment of non-detect (left-censored) data is a common issue for the statistical analysis of lab concentration measurements for risk assessment or reporting purposes. Environmental monitoring programs such as those analysing pesticides, metals, or other contaminants, usually report non detects as being below the analytical limit of reporting (LOR) of the analysis. Most laboratories report the detection of a chemical at a concentration less than the LOR of as “< LOR”, for example “< 0.02 µg/L” if the limit of reporting is 0.02 µg/L. All that is known about the concentration of such data is that the ‘true’ value must lie somewhere between zero and the LOR. Proper treatment of such data is imperative for risk assessment, with common approaches being to report the data at the LOR (simply by removing the “<” sign) or as a fraction of the LOR (e.g. half or some other fraction). This package provides two ways to treat LORs using the `treat_LORs_all_data()` function. The first is the method created by the Department of Environment and Science Queensland Water Quality Investigations outlined in the methods document above and is referred to as the **"WQI"** method. This method looks for the first above LOR (detect) concentration recorded for each sampling year or other recurring time interval and replaces every <LOR value before this with a negligible value of 10^-11. Following the first detect, it is assumed there is a very small amount of the contaminant in the waterway, so the <LOR values are replaced with the absolute LOR (less than sign removed) value multiplied by the relative toxicity value defined in the `pesticide_info` table. This standardizes the LOR replacement method so that all contaminants are treated equally according to toxicity; and avoids the introduction of an artificial toxicity signature. The second method is the replacing all LOR values with zero and is referred to as the **"zero"** method. The `treat_LORs_all_data()` function requires a concentration data set with **a column that matches each pesticide name in the** `pesticide_info` **table and a "Site Name" and "Date" column**, and how to use it can be seen below with the example `Canto_pesticides` concentration data provided in the package.

```{r treat_LORs_all_data, echo=TRUE, warning=FALSE}

datatable(Canto_pesticides, options = list(pageLength = 10,lengthMenu = c(5, 10, 15, 20), scrollX = T)) #Canto pesticide concentrations before LOR treatment

Canto_pesticides_LOR_treated <- treat_LORs_all_data(raw_data = Canto_pesticides, #this is the pesticide concentration data set to be treated
pesticide_info = CalcThemAll.PRM::pesticide_info, #this specifies the pesticide info look-up table
treatment_method = "WQI") #this selects the LOR treatment method

datatable(Canto_pesticides_LOR_treated, options = list(pageLength = 10,lengthMenu = c(5, 10, 15, 20), scrollX = T)) #Canto pesticide concentrations after treatment, LORs should be replaced with either 0.0000001 or LOR replacement value

```


## Calculate Daily Average PRM

Once the `pesticide_info` table has all the required pesticides and tbe LOR values have been treated **we can begin the calculating daily average PRM values**, the main focus of this package! The `calculate_daily_average_PRM()` function calculates a daily average PRM value for each day that pesticide concentration data is provided in a "Total PRM" column. It also calculates daily average PRM values for each pesticide_type in the `pesticide_info` table in a column named after the type followed by " PRM", for example "PSII Herbicide PRM". This function is designed to run with a data frame exported from the `treat_LORs_all_data()` function. However, if you skipped the treat LORs step the input data frame needs all the same columns and an additional column for "Sampling Year". An example of the function is shown below.

```{r calculate_daily_average_PRM, echo=TRUE, warning=FALSE}
head(Canto_pesticides_LOR_treated) #Canto pesticide concentrations after LOR treatment

#calculate daily PRM
Canto_daily_PRM <- calculate_daily_average_PRM(LOR_treated_data = Canto_pesticides_LOR_treated)
head(Canto_daily_PRM)

```

### Plot Daily Average PRM

This package also contains a function to **plot daily PRM for a single site and sampling year**. `plot_daily_PRM()` creates an interactive plotly plot for a single PRM group/type, colour coded to match the PRM legend table shown at the beginning of this vignette.

```{r plot_daily_PRM, echo=TRUE, warning=FALSE, out.width = '100%'}
#filter daily PRM data for a single site and sampling year
Violet_Town_2017_2018_PRM <- Canto_daily_PRM %>%
 dplyr::filter(.data$`Sampling Year` ==  "2017-2018" &  .data$`Site Name` == "Violet Town")

plot_daily_PRM(daily_PRM_data = Violet_Town_2017_2018_PRM,
               title = F, #this toggles the title on and off
               wet_season_start = "2017-10-02", #start date of the wet season or high-risk window
                                                #this is optional and can be removed with = NULL
               wet_season_length = 182, #length of wet season or high-risk window
               PRM_group = "PSII Herbicide PRM") #PRM group to plot, for all PRM = "Total PRM"

```

## PRM Data Table

This package also contains a function to create a **Data Table colour coded by PRM risk level**. `PRM_DT()` creates an interactive data table with the ability to choose if colour is the cell fill or text.

```{r PRM_DT, echo=TRUE, warning=FALSE, out.width = '100%'}
PRM_DT(PRM_data = Canto_daily_PRM, fill_cols = "Total PRM",
 colour_cols = c("PSII Herbicide PRM", "Other Herbicide PRM", "Insecticide PRM"))
```

## Calculate Wet Season Average PRM

The final part of this package is for those interested in **estimating a PRM value for a reoccurring large window of time like a wet season** or a period of high pollution risk. It involves using the calculated daily average PRM values within a wet season window or high-risk window to generate a single mean value using **multiple imputation** to fill in non-sampling days. More detail about this can be found in the WQI methods document. For example if looking at wet seasons, if 4 sampling year's worth of concentration data is provided there will be a single wet season PRM for each year. This is done using the `calculate_wet_season_average_PRM()` function and at this time is only capable of calculating one pesticide_type at a time or "Total PRM" being all pesticide groups together. It is important to note that the multiple imputation method, like any statistical method, has its own underlying assumptions and limitations, the primary one being that the observed data are representative and that there are no unmeasured relationships that are causing the data to be missing. An example of this would be if a monitoring program only samples during wet weather events when runoff and pesticide risk is high. In this example, the multiple imputation method will replicate this data to fill in missing sampled days, and thereby the final average will be reflective of wet weather events and not the entire (low flow and high flow) risk profile. The multiple imputation approach can still be used, but consideration should be given to communication of the final average and what it means. 

```{r calculate_wet_season_PRM, echo=TRUE, warning=FALSE}
#This calculation can take a few minutes so one site & sampling year is used in this example
Celestial_City_2019_2020_daily_PRM <- Canto_daily_PRM %>% 
  dplyr::filter(`Site Name` == "Celestial City" & `Sampling Year` == "2019-2020")

CC2019_2020_wet_season_Total_PRM <- calculate_wet_season_average_PRM(daily_PRM_data = Celestial_City_2019_2020_daily_PRM, PRM_group = "Total PRM")
#this calculates the wet season average PRM for all pesticide groups in one total value
#to calculate for a specific group define it in "PRM_group ="

CC2019_2020_wet_season_Total_PRM

```

## Disclaimer

Information is from several sources and, as such, does not necessarily represent government or departmental policy. While every care is taken to ensure the accuracy of this information, the Department of Environment and Science makes no representations or warranties relating to  accuracy, reliability, completeness, currency or suitability for any particular purpose and disclaims all responsibility and all liability (including without limitation, liability in negligence) for all expenses, losses, damages (including indirect or consequential damage) and costs that might be incurred as a result of any use or of reliance on the information and calculated data in any way and for any reason.

## Citation

**R Package:**

*Bezzina A, Neelamraju C, Strauss J, Kaminski H, Roberts C, Glen J, Dias F. 2022. CalcThemAll.PRM: Pesticide Risk Metric Calculations. R package. Water Quality Monitoring & Investigations, Department of Environment and Science, Queensland Government. https://github.com/AlexWaterboyBezzina/CalcThemAll.PRM*

**Methods Behind Pesticide Risk Metric:**

*Warne MStJ, Neelamraju C, Strauss J, Smith RA, Turner RDR, Mann RM. 2020. Development of a method for estimating the toxicity of pesticide mixtures and a Pesticide Risk Baseline for the Reef 2050 Water Quality Improvement Plan. Brisbane: Department of Environment and Science, Queensland Government.*
