README.Rmd

---
title: "pspforecast"
output: github_document
---


Shellfish toxicity forecast serving package


## Requirements

  + [R v4+](https://www.r-project.org/)
  
  + [rlang](https://CRAN.R-project.org/package=rlang)
  
  + [dplyr](https://CRAN.R-project.org/package=dplyr)
  
  + [readr](https://CRAN.R-project.org/package=readr)
  
  + [tidyr](https://CRAN.R-project.org/package=tidyr)
  
  + [httr](https://CRAN.R-project.org/package=httr)

## Installation

```
remotes::install_github("BigelowLab/pspforecast")
```

```{r message=FALSE, warning=FALSE, echo=FALSE}
suppressPackageStartupMessages({
  library(dplyr)
  library(ggplot2)
  library(plotly)
  
  library(pspdata)
  library(pspforecast)
  
  library(patchwork)
})
```

## Reading the forecast database 

### Variables:
 + version - the version/configuration of the model used to make the prediction
 
 + ensemble_n - number of ensemble members used to generate prediction
 
 + location - the sampling station the forecast is for
 
 + date - the date the forecast was made on
 
 + name - site name
 
 + lat - latitude
 
 + lon - longitude
 
 + class_bins - the bins used to classify shellfish total toxicity (i.e. 0: 0-10, 1: 10-30, 2: 30-80, 3: >80)
 
 + forecast_date - the date the forecast is valid for (i.e. one week ahead of when it was made)
 
 + predicted_class - the predicted classification at the location listed on the forecast_date (in this case 0-3)
 + p_0 - class 0 probability
 + p_1 - class 1 probability
 + p_2 - class 2 probability
 + p_3 - class 3 probability
 + p3_sd - class 3 probability standard deviation
 + p_3_min - class 3 minimum probability (from ensemble run)
 + p_3_max - class 3 maximum probability (from ensemble run)
 + predicted_class - the predicted classification


```{r echo=FALSE, message=FALSE, warning=FALSE}
psp <- read_psp_data()
```


```{r}
predictions <- read_forecast(year = "2024") |>
  distinct()

glimpse(predictions)
```

## 2024 Season Results

```{r echo=FALSE}
pred_w_results <- add_forecast_results(predictions, toxin_measurements = psp) 
```

```{r cm24, echo=FALSE, cahce=FALSE}
num_levels <- 4
levels <- seq(from=0, to=(num_levels-1))


cm <- as.data.frame(table(predicted = factor(pred_w_results$predicted_class, levels), actual = factor(pred_w_results$class, levels))) |> 
  mutate(frac = round(Freq/sum(Freq)*100)) |> 
  mutate(frac = sapply(.data$frac, function(x) if (x == "0") {x = "<1"} else {x}))
  #mutate(frac = sapply(.data$frac, function(x) if (x == "0") {x = "<1"} else {x}))

confusion_matrix_24 <- ggplot2::ggplot(data = cm,
                                    mapping = ggplot2::aes(x = .data$predicted, y = .data$actual)) +
  ggplot2::geom_tile(ggplot2::aes(fill = log(.data$Freq+1))) +
  ggplot2::geom_text(ggplot2::aes(label = sprintf("%1.0f", .data$Freq)), size=8) +
  #ggplot2::geom_text(ggplot2::aes(label = sprintf("%1.0f%%", .data$frac)), vjust = -1, hjust=1, size=5) +
  ggplot2::geom_text(ggplot2::aes(label = paste(.data$frac, "%", sep="")), vjust = 3, size=4)+
  ggplot2::scale_fill_gradient(low = "white", 
                               high = "blue") +
  ggplot2::labs(x = "Predicted Classifications", 
                y = "Actual Classifications") +
  ggplot2::theme_linedraw() +
  ggplot2::theme(axis.text=  ggplot2::element_text(size=14),
                 axis.title= ggplot2::element_text(size=14,face="bold"),
                 title =     ggplot2::element_text(size = 14, face = "bold"),
                 legend.position = "none") +
  ggplot2::geom_rect(aes(xmin=0.5, xmax=3.5, ymin=0.5, ymax=3.5), alpha=0) +
  ggplot2::geom_rect(aes(xmin=3.5, xmax=4.5, ymin=3.5, ymax=4.5), alpha=0)


confusion_matrix_24

```

```{r scatter24, echo=FALSE, cache=FALSE}

p3_v_tox_24 <- ggplot2::ggplot(data = pred_w_results, ggplot2::aes(x=.data$p_3, y=.data$toxicity)) +
  ggplot2::geom_point(alpha=0.7, size=3) +
  ggplot2::labs(x = "Predicted Probability of Closure-level Toxicity",
                y = "Measured Toxicity") +
  ggplot2::geom_hline(yintercept=80, linetype="dashed")

p3_v_tox_24
```


### Metrics 

 + tp - The model predicted class 3 and the following week's measurement was class 3
 + fp - The model predicted class 3 and the following week's measurement was not class 3
 + tn - The model predicted class 0,1,2 and the following week's measurement was in class 0,1,2
 + fn - The model predicted class 0,1,2 and the following week's measurement was class 3
 + accuracy - Measure of how many correct classifications were predicted
 + cl_accuracy - Considering predictions are those that correctly predicted toxicity above or below the closure limit or not
 + precision - TP/(TP+FP)
 + sensitivity - TP/(TP+FN)
 + specificity - TN/(TN+FP)
 + f_1


```{r echo=FALSE}
forecast_metrics(pred_w_results, measured_col = "class")
```


## 2023 Season Results

```{r}
predictions <- read_forecast(year = "2023")
```


```{r echo=FALSE}
pred_w_results <- add_forecast_results(predictions, toxin_measurements = psp)
```

### Confusion Matrix

```{r cm23, echo=FALSE, cahce=FALSE}
num_levels <- 4
levels <- seq(from=0, to=(num_levels-1))


cm <- as.data.frame(table(predicted = factor(pred_w_results$predicted_class, levels), actual = factor(pred_w_results$class, levels))) |> 
  mutate(frac = round(Freq/sum(Freq)*100)) |> 
  mutate(frac = sapply(.data$frac, function(x) if (x == "0") {x = "<1"} else {x}))
  #mutate(frac = sapply(.data$frac, function(x) if (x == "0") {x = "<1"} else {x}))

confusion_matrix_23 <- ggplot2::ggplot(data = cm,
                                    mapping = ggplot2::aes(x = .data$predicted, y = .data$actual)) +
  ggplot2::geom_tile(ggplot2::aes(fill = log(.data$Freq+1))) +
  ggplot2::geom_text(ggplot2::aes(label = sprintf("%1.0f", .data$Freq)), size=8) +
  #ggplot2::geom_text(ggplot2::aes(label = sprintf("%1.0f%%", .data$frac)), vjust = -1, hjust=1, size=5) +
  ggplot2::geom_text(ggplot2::aes(label = paste(.data$frac, "%", sep="")), vjust = 3, size=4)+
  ggplot2::scale_fill_gradient(low = "white", 
                               high = "blue") +
  ggplot2::labs(x = "Predicted Classifications", 
                y = "Actual Classifications") +
  ggplot2::theme_linedraw() +
  ggplot2::theme(axis.text=  ggplot2::element_text(size=14),
                 axis.title= ggplot2::element_text(size=14,face="bold"),
                 title =     ggplot2::element_text(size = 14, face = "bold"),
                 legend.position = "none") +
  ggplot2::geom_rect(aes(xmin=0.5, xmax=3.5, ymin=0.5, ymax=3.5), alpha=0) +
  ggplot2::geom_rect(aes(xmin=3.5, xmax=4.5, ymin=3.5, ymax=4.5), alpha=0)


confusion_matrix_23

```

### Probability of Closure-level Toxicity vs Measured Toxicity

```{r scatter23, echo=FALSE, cache=FALSE}

p3_v_tox_23 <- ggplot2::ggplot(data = pred_w_results, ggplot2::aes(x=.data$p_3, y=.data$toxicity)) +
  ggplot2::geom_point(alpha=0.7, size=3) +
  ggplot2::labs(x = "Predicted Probability of Closure-level Toxicity",
                y = "Measured Toxicity") +
  ggplot2::geom_hline(yintercept=80, linetype="dashed")

p3_v_tox_23
```

### Metrics 

```{r echo=FALSE}
forecast_metrics(pred_w_results, measured_col = "class")
```



## 2022 Season Results


```{r echo=FALSE}
predictions <- read_forecast(year = "2022")
```


```{r echo=FALSE}
pred_w_results <- add_forecast_results(predictions, toxin_measurements = psp)
```

### Confusion Matrix

```{r cm22, echo=FALSE, cahce=FALSE}
num_levels <- 4
levels <- seq(from=0, to=(num_levels-1))


cm <- as.data.frame(table(predicted = factor(pred_w_results$predicted_class, levels), actual = factor(pred_w_results$class, levels))) |> 
  mutate(frac = round(Freq/sum(Freq)*100)) |> 
  mutate(frac = sapply(.data$frac, function(x) if (x == "0") {x = "<1"} else {x}))
  #mutate(frac = sapply(.data$frac, function(x) if (x == "0") {x = "<1"} else {x}))

confusion_matrix_22 <- ggplot2::ggplot(data = cm,
                                    mapping = ggplot2::aes(x = .data$predicted, y = .data$actual)) +
  ggplot2::geom_tile(ggplot2::aes(fill = log(.data$Freq+1))) +
  ggplot2::geom_text(ggplot2::aes(label = sprintf("%1.0f", .data$Freq)), size=8) +
  #ggplot2::geom_text(ggplot2::aes(label = sprintf("%1.0f%%", .data$frac)), vjust = -1, hjust=1, size=5) +
  ggplot2::geom_text(ggplot2::aes(label = paste(.data$frac, "%", sep="")), vjust = 3, size=4)+
  ggplot2::scale_fill_gradient(low = "white", 
                               high = "blue") +
  ggplot2::labs(x = "Predicted Classifications", 
                y = "Actual Classifications") +
  ggplot2::theme_linedraw() +
  ggplot2::theme(axis.text=  ggplot2::element_text(size=14),
                 axis.title= ggplot2::element_text(size=14,face="bold"),
                 title =     ggplot2::element_text(size = 14, face = "bold"),
                 legend.position = "none")

confusion_matrix_22
```

### Probability of Closure-level Toxicity vs Measured Toxicity

```{r scatter22, echo=FALSE, cache=FALSE}

p3_v_tox_22 <- ggplot2::ggplot(data = pred_w_results, ggplot2::aes(x=.data$p_3, y=.data$toxicity)) +
  ggplot2::geom_point(alpha=0.7, size=3) +
  ggplot2::labs(x = "Predicted Probability of Closure-level Toxicity",
                y = "Measured Toxicity") +
  ggplot2::geom_hline(yintercept=80, linetype="dashed")

p3_v_tox_22
```

### Metrics 

```{r metrics22, echo=FALSE}
forecast_metrics(pred_w_results, measured_col = "class")
```


### Timing of initial closure-level predictions

```{r echo=FALSE, message=FALSE, warning=FALSE}
closures <- find_closures(psp, predictions)
```


```{r echo=FALSE, warning=FALSE, message = FALSE}
ggplot(closures, aes(x = predicted_week, y = actual_week)) +
  #geom_point() +
  geom_jitter() +
  geom_abline(intercept = 0, linetype="dashed") +
  labs(x = "First predicted closure-level toxicity (week of year)",
       y = "First measured closure-level toxicity (week of year)")
```



## 2021 Season Results

```{r echo=FALSE}
predictions <- read_forecast(year = "2021")
```


```{r echo=FALSE}
pred_w_results <- add_forecast_results(predictions, toxin_measurements = psp)
```

### Confusion Matrix

```{r cm21, echo=FALSE, cahce=FALSE}
num_levels <- 4
levels <- seq(from=0, to=(num_levels-1))


cm <- as.data.frame(table(predicted = factor(pred_w_results$predicted_class, levels), actual = factor(pred_w_results$class, levels))) |> 
  mutate(frac = round(Freq/sum(Freq)*100)) |> 
  mutate(frac = sapply(.data$frac, function(x) if (x == "0") {x = "<1"} else {x}))
  #mutate(frac = sapply(.data$frac, function(x) if (x == "0") {x = "<1"} else {x}))

confusion_matrix_21 <- ggplot2::ggplot(data = cm,
                                    mapping = ggplot2::aes(x = .data$predicted, y = .data$actual)) +
  ggplot2::geom_tile(ggplot2::aes(fill = log(.data$Freq+1))) +
  ggplot2::geom_text(ggplot2::aes(label = sprintf("%1.0f", .data$Freq)), size=8) +
  #ggplot2::geom_text(ggplot2::aes(label = sprintf("%1.0f%%", .data$frac)), vjust = -1, hjust=1, size=5) +
  ggplot2::geom_text(ggplot2::aes(label = paste(.data$frac, "%", sep="")), vjust = 3, size=4)+
  ggplot2::scale_fill_gradient(low = "white", 
                               high = "blue") +
  ggplot2::labs(x = "Predicted Classifications", 
                y = "Actual Classifications") +
  ggplot2::theme_linedraw() +
  ggplot2::theme(axis.text=  ggplot2::element_text(size=14),
                 axis.title= ggplot2::element_text(size=14,face="bold"),
                 title =     ggplot2::element_text(size = 14, face = "bold"),
                 legend.position = "none")

confusion_matrix_21
```

### Probability of Closure-level Toxicity vs Measured Toxicity

```{r scatter21, echo=FALSE, cache=FALSE}

p3_v_tox_21 <- ggplot2::ggplot(data = pred_w_results, ggplot2::aes(x=.data$prob_3, y=.data$toxicity)) +
  ggplot2::geom_point(alpha=0.7, size=3) +
  ggplot2::labs(x = "Predicted Probability of Closure-level Toxicity",
                y = "Measured Toxicity") +
  ggplot2::geom_hline(yintercept=80, linetype="dashed")

p3_v_tox_21
```

### Metrics 

```{r metrics21, echo=FALSE}
forecast_metrics(pred_w_results, measured_col = "class")
```


### Closure-level accuracy


### Timing of initial closure-level predictions

```{r echo=FALSE, message=FALSE, warning=FALSE}
closures <- find_closures(psp, predictions)
```


```{r echo=FALSE, warning=FALSE, message = FALSE}
ggplot(closures, aes(x = predicted_week, y = actual_week)) +
  geom_point() +
  geom_abline(intercept = 0, linetype="dashed") +
  labs(x = "First predicted closure-level toxicity (week of year)",
       y = "First measured closure-level toxicity (week of year)")
```


### Possible manuscript plot(s)

```{r echo=FALSE}
pred_w_results <- read_all_results()

num_levels <- 4
levels <- seq(from=0, to=(num_levels-1))


cm <- as.data.frame(table(predicted = factor(pred_w_results$predicted_class, levels), actual = factor(pred_w_results$class, levels), year=factor(pred_w_results$year, levels=2021:2024))) |> 
  dplyr::mutate(frac = round(Freq/sum(Freq)*100)) |> 
  dplyr::mutate(frac = sapply(.data$frac, function(x) if (x == "0") {x = "<1"} else {x}))

ggplot2::ggplot(data = cm, ggplot2::aes(x = .data$predicted, y = .data$actual)) +
  ggplot2::geom_tile(ggplot2::aes(fill = log(.data$Freq+1))) +
  ggplot2::geom_text(ggplot2::aes(label = sprintf("%1.0f", .data$Freq)), size=8) +
  ggplot2::facet_grid(cols=vars(.data$year)) +
  #ggplot2::geom_text(ggplot2::aes(label = paste(.data$frac, "%", sep="")), vjust = 3, size=4) +
  ggplot2::scale_fill_gradient(low = "white", 
                               high = "blue") +
  ggplot2::labs(x = "Predicted Classifications", 
                y = "Actual Classifications") +
  ggplot2::theme_linedraw() +
  ggplot2::theme(axis.text=  ggplot2::element_text(size=14),
                 axis.title= ggplot2::element_text(size=14,face="bold"),
                 title =     ggplot2::element_text(size = 14, face = "bold"),
                 legend.position = "none") +
  ggplot2::geom_rect(aes(xmin=0.5, xmax=3.5, ymin=0.5, ymax=3.5), alpha=0) +
  ggplot2::geom_rect(aes(xmin=3.5, xmax=4.5, ymin=3.5, ymax=4.5), alpha=0)
```


```{r echo=FALSE}
ggplot2::ggplot(data = pred_w_results, ggplot2::aes(x=.data$p_3, y=.data$toxicity)) +
  ggplot2::geom_point(alpha=0.7, size=3) +
  ggplot2::facet_grid(cols=vars(.data$year)) +
  ggplot2::labs(x = "Predicted Probability of Closure-level Toxicity",
                y = "Measured Toxicity") +
  ggplot2::geom_hline(yintercept=80, linetype="dashed") +
  ggplot2::theme_bw()

```

### Last Updated

```{r echo=FALSE}
Sys.Date()
```