---
title: "Get Started with ssdtools"
author: "ssdtools Team"
date: "`r Sys.Date()`"
bibliography: references.bib
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Get Started with ssdtools}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 4,
  fig.height = 4
)
```

## Introduction

`ssdtools` is an R package to fit Species Sensitivity Distributions (SSDs) using Maximum Likelihood and model averaging.

SSDs are cumulative probability distributions that are used to estimate the percent of species that are affected by a given concentration of a chemical.
The concentration that affects 5% of the species is referred to as the 5% Hazard Concentration (HC).
For more information on SSDs the reader is referred to @posthuma_species_2001.

In order to use `ssdtools` you need to install R (see below) or use the Shiny [app](https://bcgov-env.shinyapps.io/ssdtools).
The shiny app includes a user guide. 
This vignette is a user manual for the R package.

## Philosophy

`ssdtools` provides the key functionality required to fit SSDs using Maximum Likelihood and model averaging in R.
It is intended to be used in conjunction with [tidyverse](https://www.tidyverse.org) packages such as `readr` to input data, `tidyr` and `dplyr` to group and manipulate data and `ggplot2` [@ggplot2] to plot data.
As such it endeavours to fulfill the tidyverse [manifesto](https://tidyverse.tidyverse.org/articles/manifesto.html).

## Installing

In order to install R [@r] the appropriate binary for the users operating system should be downloaded from [CRAN](https://cran.r-project.org) and then installed.

Once R is installed, the `ssdtools` package can be installed (together with the tidyverse) by executing the following code at the R console
```{r, eval = FALSE}
install.packages("ssdtools")
install.packages("tidyverse")
```

The `ssdtools` package (and key packages) can then be loaded into the current session using
```{r, message = FALSE}
library(ssdtools)
library(readr)
library(ggplot2)
```

## Getting Help

To get additional information on a particular function just type `?` followed by the name of the function at the R console.
For example `?ssd_gof` brings up the R documentation for the `ssdtools` goodness of fit function.

For more information on using R the reader is referred to [R for Data Science](https://r4ds.had.co.nz) [@wickham_r_2016].

If you discover a bug in `ssdtools` please file an issue with a [reprex](https://reprex.tidyverse.org/articles/reprex-dos-and-donts.html) (repeatable example) at <https://github.com/bcgov/ssdtools/issues>.

## Inputting Data

Once the `ssdtools` package has been loaded the next task is to input some data.
An easy way to do this is to save the concentration data for a *single* chemical as a column called `Conc` in a comma separated file (`.csv`). 
Each row should be the sensitivity concentration for a separate species.
If species and/or group information is available then this can be saved as `Species` and `Group` columns.
The `.csv` file can then be read into R using the following
```{r, eval = FALSE}
data <- read_csv(file = "path/to/file.csv")
```
For the purposes of this manual we use the CCME dataset for boron which is provided with the `ssdtools` package.
```{r}
boron_data <- ssdtools::boron_data
print(boron_data)
```

## Fitting Distributions

The function `ssd_fit_dists()` inputs a data frame and fits one or more distributions.
The user can specify a subset of the 

- [gamma](https://en.wikipedia.org/wiki/Gamma_distribution) (`gamma`), 
- [Gompertz](https://en.wikipedia.org/wiki/Gompertz_distribution) (`gompertz`),
- [log-Gumbel](https://en.wikipedia.org/wiki/Gumbel_distribution) (`lgumbel`),
- [log-logistic](https://en.wikipedia.org/wiki/Log-logistic_distribution) (`llog`),
- [log-normal](https://en.wikipedia.org/wiki/Log-normal_distribution) (`lnorm`) and
- [Weibull](https://en.wikipedia.org/wiki/Weibull_distribution) (`weibull`) 

distributions using the `dists` argument.
```{r}
boron_dists <- ssd_fit_dists(boron_data, dists = c("llogis", "lnorm", "gamma"))
```

The user can also specify one or more [custom distributions](https://bcgov.github.io/ssdtools/articles/faqs.html#how-do-i-fit-a-custom-distribution).

## Coefficients

The coefficients can be extracted using the `coef` function.
However, in and off themselves the coefficients are not that helpful.
```{r}
coef(boron_dists)
```

## Plots

It is generally much more informative to plot the fits using the `autoplot` generic function.
As `autoplot` returns a `ggplot` object it can be modified prior to plotting (printing) to make it look prettier.
```{r, fig.width = 5}
theme_set(theme_bw()) # set plot theme
gp <- autoplot(boron_dists)
gp <- gp + ggtitle("Species Sensitivity Distributions for Boron")
print(gp)
```

## Selecting One Distribution

Given multiple distributions the user is faced with choosing the best fitting distribution (or as discussed below averaging the results weighted by the fit).

```{r}
boron_gof <- ssd_gof(boron_dists)
boron_gof[order(boron_gof$delta), ]
```

The `ssd_gof()` function returns several goodness of fit measures that can be used to select the best distribution including three statistics 

- [Anderson-Darling](https://en.wikipedia.org/wiki/Anderson–Darling_test) (`ad`) statistic,
- [Kolmogorov-Smirnov](https://en.wikipedia.org/wiki/Kolmogorov–Smirnov_test) (`ks`) statistic and
- [Cramer-von Mises](https://en.wikipedia.org/wiki/Cramér–von_Mises_criterion) (`cvm`) statistic

and three information criteria

- Akaike's Information Criterion (`aic`),
- Akaike's Information Criterion corrected for sample size (`aicc`) and
- Bayesian Information Criterion (`bic`)

Following @burnham_model_2002 we recommend the `aicc` for model selection.
The best fitting model is that with the lowest `aicc` (indicated by the model with a `delta` value of 0.000 in the goodness of fit table).
In the current example the best fitting model is the gamma distribution but the lnorm distribution has some support.

For further information on the advantages of an information theoretic approach in the context of selecting SSDs the reader is referred to @schwarz_improving_2019

## Averaging Multiple Distributions

Often other distributions will fit the data almost as well as the best distribution as evidenced by `delta` values < 2 [@burnham_model_2002]. 
In this situation the recommended approach is to estimate the average fit based on the relative weights of the distributions  [@burnham_model_2002].
The `aicc` based weights are indicated by the `weight` column in the goodness of fit table.
In the current example, the gamma and log-normal distributions have `delta` values < 2.

## Estimating the Fit

The `predict` function can be used to generate estimates model-averaged (or if `average = FALSE` individual) estimates.
By default model averaging is based on `aicc`.

```{r, eval = FALSE}
set.seed(99)
boron_pred <- predict(boron_dists, ci = TRUE)
```

The resultant object is a data frame of the estimated concentration (`est`) with standard error (`se`) and lower (`lcl`) and upper (`ucl`) 95% confidence limits by percent of species affected (`percent`).
The confidence limits are estimated using parametric bootstrapping.
```{r}
boron_pred
```
The data frame of the estimates can then be plotted together with the original data using the `ssd_plot()` function to summarize an analysis.
Once again the returned object is a `ggplot` object which can be customized prior to plotting.
```{r, fig.height = 5, fig.width = 6}
gp <- ssd_plot(boron_data, boron_pred,
  color = "Group", label = "Species",
  xlab = "Concentration (mg/L)", ribbon = TRUE
)
gp <- gp + expand_limits(x = 5000) + # to ensure the species labels fit
  scale_color_manual(values = c(
    "Amphibian" = "Black", "Fish" = "Blue",
    "Invertebrate" = "Red", "Plant" = "Brown"
  )) +
  ggtitle("Species Sensitivity for Boron")
print(gp)
```

In the above plot the model-averaged 95% confidence interval is indicated by the shaded band and the model-averaged 5% Hazard Concentration ($HC_5$) by the dotted line.
Hazard concentrations are discussed below.

## Hazard Concentrations

The 5% hazard concentration ($HC_5$) is the concentration that affects 5% of the species tested.

```{r, eval = FALSE}
set.seed(99)
boron_hc5 <- ssd_hc(boron_dists, ci = TRUE)
```

```{r}
print(boron_hc5)
```

## Plotting

The `ssdtools` package provides three ggplot geoms to allow you construct your own plots.

The first is `geom_ssd()` which plots species sensitivity data 
```{r}
ggplot(boron_data) +
  geom_ssd(aes_string(x = "Conc"))
``` 

The second is `geom_xribbon()` which plots species sensitivity confidence intervals
```{r}
ggplot(boron_pred) +
  geom_xribbon(aes_string(xmin = "lcl", xmax = "ucl", y = "percent/100"))
```  

And the third is `geom_hcintersect()` which plots hazard concentrations
```{r}
ggplot() +
  geom_hcintersect(xintercept = c(1, 2, 3), yintercept = c(5, 10, 20) / 100)
```    

They can be combined together as follows
```{r}
gp <- ggplot(boron_pred, aes_string(x = "est")) +
  geom_xribbon(aes_string(xmin = "lcl", xmax = "ucl", y = "percent/100"), alpha = 0.2) +
  geom_line(aes_string(y = "percent/100")) +
  geom_ssd(data = boron_data, aes_string(x = "Conc")) +
  scale_y_continuous("Species Affected (%)", labels = scales::percent) +
  expand_limits(y = c(0, 1)) +
  xlab("Concentration (mg/L)")
print(gp + geom_hcintersect(xintercept = boron_hc5$est, yintercept = 5 / 100))
```

To log the x-axis add the following code.
```{r}
gp <- gp + coord_trans(x = "log10") +
  scale_x_continuous(
    breaks = scales::trans_breaks("log10", function(x) 10^x),
    labels = comma_signif
  )
print(gp + geom_hcintersect(xintercept = boron_hc5$est, yintercept = 5 / 100))
```

The most recent plot can be saved as a file using `ggsave()`, which also allows the user to set the resolution.
```r
ggsave("file_name.png", dpi = 600)
```

### Censored Data

Censored data is that for which only a lower and/or upper limit is known for a particular species. 
If the `right` argument in `ssd_fit_dists()` is different to the `left` argument then the data are considered to be censored. 
`fluazinam` is a censored data set from the `fitdistrplus` package.

```{r}
data(fluazinam, package = "fitdistrplus")
head(fluazinam)
```

There are less goodness-of-fit statistics available for
fits to censored data (currently just `aic` and `bic`).
The `delta` values are calculated using `aic`.

As the sample size `n` is undefined for censored data, `aicc` cannot be calculated.
However, if all the models have the same number of parameters, the `aic` `delta` values are identical to those for `aicc`.
For this reason, `ssdtools` only permits the analysis of censored data using two-parameter models.

```{r, eval = FALSE}
fluazinam_dists <- ssd_fit_dists(fluazinam, left = "left", right = "right")
ssd_gof(fluazinam_dists)
```
The model-averaged predictions (and hazard concentrations complete with 95% confidence limits) can be calculated using `aic`
```{r, eval = FALSE}
set.seed(99)
fluazinam_pred <- predict(fluazinam_dists, ci = TRUE)
```
and the results plotted complete with arrows indicating the censorship.
```{r}
ssd_plot(fluazinam, fluazinam_pred,
  left = "left", right = "right",
  xlab = "Concentration (mg/L)"
)
```

-----

<a rel="license" href="https://creativecommons.org/licenses/by/4.0/"><img alt="Creative Commons Licence"
style="border-width:0" src="https://i.creativecommons.org/l/by/4.0/80x15.png" /></a><br /><span
xmlns:dct="http://purl.org/dc/terms/" property="dct:title">ssdtools</span> by <span
xmlns:cc="http://creativecommons.org/ns#" property="cc:attributionName">the Province of British Columbia
</span> is licensed under a <a rel="license" href="https://creativecommons.org/licenses/by/4.0/">
Creative Commons Attribution 4.0 International License</a>.

## References

