---
title: "Literature Review"
output: html_document
date: "2026-04-06"
author: "Olivia Cronin-Golomb"
contact: "olivia.cronin-golomb@whoi.edu"
---

# Set Up Environment
## Load libraries
```{r}
rm(list=ls()) #clear environment

library(dplyr)
library(ggplot2)
library(RColorBrewer)
library(tidyverse)
library(ggnewscale)
library(cowplot)
library(scales)
library(ggpattern)
library(extrafont)
library(boot)
library(ggtext)
```

## Set input and output directories
```{r}
#Input directory
wd <- "C:/Users/ocron/OneDrive - Massachusetts Institute of Technology/Documents/EPA/Temporal_Mismatch_DataPackage1/Data/Processed/LitReview/"

#Output directory
Paper_figures_filepath <- "C:/Users/ocron/OneDrive - Massachusetts Institute of Technology/Documents/EPA/Temporal_Mismatch_DataPackage1/Outputs/Lit_review/"

Extra_figures_filepath <- "C:/Users/ocron/OneDrive - Massachusetts Institute of Technology/Documents/EPA/Temporal_Mismatch_DataPackage1/Outputs/Extra_Figures/"
```

# Read in lit review tables
## Temporal Mismatches
```{r}
MatchUp <- read.csv(paste0(wd, "Sentinel2_LitReview_Temporal_Mismatch.csv"))

#Remove rows with a 'NaN' index
MatchUp <- MatchUp[!is.na(MatchUp$Index),]

#Only keep rows with studies that used field data (not synthetic)
MatchUp <- MatchUp[MatchUp$Relevant. == "Y",]

#Rename any NaNs in caterogy columns as 'Unreported'
MatchUp <- MatchUp %>%
    replace_na(list(Category_jDay = 'Unreported'))
MatchUp <- MatchUp %>%
    replace_na(list(Category = 'Unreported'))
row.names(MatchUp) <- NULL
```

## MAE
```{r}
Algorithm <- read.csv(paste0(wd, "Sentinel2_LitReview_MAE.csv"))

#Filter out NaNs
Algorithm <- Algorithm[!is.na(Algorithm$MAE),]

# Keep 'Index', 'MAE' and 'MAE_type' columns
Algorithm <- Algorithm[,c(1,11,12)]
```

# Plot counts of temporal mismatch categories
## Get counts
```{r}
sortCat <- c("1 hour", "2 hours","3 hours", "4 hours", "5 hours", "6 hours", "Calendar Day","1 day", "2 days", "3 days", "4 days", "5 days", "6 days", "7 days", "2 weeks", "3 weeks", "4 weeks", "2 years", "3 years", "4 years", "5 years", "Unreported")

sortCatDF <- as.data.frame(sortCat)
sortCatDF$value <- seq(1, nrow(sortCatDF),1)
names(sortCatDF) <- c("Category_jDay", "value")

countPivot <- MatchUp %>%
  group_by(Category_jDay) %>%
  summarize(count_cat = n())
names(countPivot) <- c("Category_jDay", "Count")

pivot <- merge(sortCatDF, countPivot, by = "Category_jDay", all = T)
pivot <- pivot[order(pivot$value),]
pivot$Count[is.na(pivot$Count)] <- 0
row.names(pivot) <- NULL

pivot$Category_Day <- c(rep("1 day", 8), sortCat[9:22])

#Assign grey scale colors to eachcategory
colorList <- c(brewer.pal(9,"Greys"))
colorList <- c(colorList, rep(colorList[9], 12), "black")
pivot$Color <- colorList

#Break df into one day, greater than one day, and unreported for plotting purposes
day_one <- pivot[c(1:8),]
day_one$Index <- seq(nrow(day_one), 1, -1)
day_one <- day_one[order(day_one$Index),]
row.names(day_one) <- NULL
day_one$Category_jDay = factor(day_one$Category_jDay, levels = rev(c("1 hour", "2 hours","3 hours", "4 hours", "5 hours", "6 hours", "Calendar Day","1 day")))

day_rest <- pivot[c(10:nrow(pivot)-1),]
unreported <- pivot[c(nrow(pivot)),]
```

## Create barplot 
```{r}
# Example data for the new bar
new_data <- data.frame(
  Category_Day = "Unreported",  
  Count = 15                    
)

gg_count <- ggplot() + 
  # Internal divisions in "1 day" column with black
  geom_bar(data = day_one, aes(x = Category_Day, y = Count, fill = Category_jDay), 
           position = "stack", 
           stat = "identity", 
           width = 1, 
           color = "grey15",  
           linewidth = 0.5) + 
  # External border for "1 day" column with light grey
  geom_bar(data = day_one, aes(x = Category_Day, y = sum(Count)), 
           position = "stack", 
           stat = "identity", 
           width = 1, 
           color = "white",  
           fill = NA,
           linewidth = 0.75) + 
  # Bars for other days
  geom_bar(data = day_rest, aes(x = Category_jDay, y = Count), 
           color = "white",
           stat = "identity", 
           width = 1,
           fill = "black", 
           linewidth = 0.75) +
  # Add the "Unreported" bar with a pattern
  geom_col_pattern(data = unreported, 
                   aes(x = Category_jDay, y = Count, fill = Category_jDay), 
                   color = "black",
                   fill = "white",
                   pattern_fill = "black",
                   pattern_angle = 45,
                   pattern_density = 0.1,
                   pattern_spacing = 0.025,
                   pattern_key_scale_factor = 0.6,
                   linewidth = 0.1) +
  # Customize theme and other elements
  theme_classic() + 
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1, size = 16),
        axis.text.y = element_text(size = 16),
        axis.title.y = element_text(size = 18),
        axis.title.x = element_blank(),
        plot.title = element_text(hjust = 0.5, size = 20),
        panel.grid.minor = element_blank(),
        panel.grid.major.x = element_blank(),
        legend.text = element_text(size = 16),
        legend.title = element_text(size = 18),
        text = element_text(family="Times New Roman"),
        legend.position = c(0.87, 0.74),  
        legend.background = element_rect(fill='transparent'),
        panel.border = element_rect(color = "black", fill = NA, linewidth = 0.5)) +
  guides(fill = guide_legend(byrow = TRUE)) +
  # Scale for the legend, only including `day_one` categories
  scale_fill_manual(drop = FALSE, 
                    name = "",
                    breaks = c(day_one$Category_jDay),
                    labels = rev(c("1 hour", "2 hours", "3 hours", "4 hours", "5 hours", "6 hours", "Calendar day", "1 day")),
                    values = c(rev(pivot$Color[1:8]), pivot$Color[9:21])) +
  scale_y_continuous(breaks = seq(0, 40, 10), limits = c(0, 40)) +
  scale_x_discrete(limits = c(pivot$Category_jDay[8:21], "Unreported"), drop = FALSE) + 
  geom_segment(aes(x = 0.52, xend = 1.47, y = 0.19, yend = 0.19), 
               color = "grey15", linewidth = 0.5) +
  geom_segment(aes(x = 0.52, xend = 1.47, y = 37.9, yend = 37.9), 
               color = "grey15", linewidth = 0.5) +
  geom_segment(aes(x = 0.52, xend = 0.52, y = 0.12, yend = 37.98), 
             color = "grey15", linewidth = 0.5) +
  geom_segment(aes(x = 1.47, xend = 1.47, y = 0.12, yend = 37.98), 
             color = "grey15", linewidth = 0.5)

gg_count

ggsave(path = Paper_figures_filepath, filename = "2_LitReview_Temporal_Mismatch_Counts.png", device = "png", plot = gg_count, width = 7.98, height = 5, dpi = 600, units = "in")
```

# Plot MAE
## Only studies with additive MAE
```{r}
#Select only temporal mismatch studies with MAE report
alg_match <- merge(Algorithm, MatchUp, by = "Index")
alg_match <- alg_match[order(alg_match$Index, -abs(alg_match$MAE) ), ]

alg_match <- alg_match %>%
  mutate(Method_Group = case_when(
    is.na(Chl.a.method) ~ "Unspecified",
    grepl("Sonde|YSI", Chl.a.method, ignore.case = TRUE) ~ "Fluorescence",
    grepl("Spec|HPLC|ethanol", Chl.a.method, ignore.case = TRUE) ~ "Extraction",
    grepl("Radiometer|LiDAR|WV3", Chl.a.method, ignore.case = TRUE) ~ "Remote",
    TRUE ~ "Other"
  ))

#Only select additive MAE
vals <- c("|estimated - measured|/n", "|estimated - measured|/n *notExplicit", "|estimated - measured|/n *noSum")
alg_match <- alg_match[ alg_match$MAE_type %in% vals, ]
alg_match <- alg_match[alg_match$Category_nHours != 1534, ]

alg_match_clean <- na.omit(alg_match[,c("Index" ,'MAE', "Category_nHours","Method_Group")])

```

## Compute Spearmans Rank
### Correlation Coefficient
```{r}
cor_test <- cor.test(alg_match_clean$Category_nHours, alg_match_clean$MAE, method = "spearman")
print(cor_test$estimate)
print(cor_test$p.value)
```

### Confidence Intervals
```{r}
# Function to compute Spearman's rho for bootstrapping
spearman_rho <- function(data, indices) {
  d <- data[indices, ]
  return(cor(d$Category_nHours, d$MAE, method = "spearman"))
}

# Bootstrap (2000 replicates)
set.seed(123)
boot_result <- boot(data = alg_match_clean, statistic = spearman_rho, R = 2000)

# Compute 95% CI
boot.ci(boot_result, type = "perc")
```
### Spearman result summary
```{r}
rho_val <- round(cor_test$estimate,2)
ci_lower <- round(boot.ci(boot_result, type="perc")$percent[4],2)
ci_upper <- round(boot.ci(boot_result, type="perc")$percent[5],2)
p_val <- "< 0.001"

# Create text label
spearman_text <- paste0("ρ = ", rho_val, 
                        "\n95% CI: ", ci_lower, "\u2013", ci_upper, 
                        "\np ", p_val)
```

## MAE Scatterplot
### Regular axes
```{r}
gg_mae <- ggplot(alg_match_clean, aes(x = Category_nHours, y = MAE, color = Method_Group)) +
  geom_point(size = 3, alpha = 0.8) +
  theme_classic() +
  labs(
    x = "No. of days since sampling",
    y = "Mean absolute error (µg L⁻¹)",
    color = "In situ chl-a \nprocessing method"
  ) +
  scale_color_manual(values = c(
    "Fluorescence" = "steelblue",
    "Extraction" = "darkgreen",
    "Remote" = "darkorange",
    "Unspecified" = "darkorange",
    "Other" = "purple"
  )) +
  scale_x_continuous(breaks = seq(0, 30, 5), limits = c(0, 30)) +
  theme(
    axis.text.x = element_text(size = 16),
    axis.text.y = element_text(size = 16),
    axis.title.x = element_text(size = 18),
    axis.title.y = element_text(size = 18),
    plot.title = element_text(hjust = 0.5, size = 20),
    legend.text = element_text(size = 16),
    legend.title = element_text(size = 17),
    legend.spacing.y = unit(0.25, 'cm'),
    text = element_text(family = "Times New Roman"),
    legend.position = c(0.84, 0.8),  
    panel.border = element_rect(color = "black", fill = NA, linewidth = 0.5),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    legend.background = element_rect(fill = "transparent", color = "grey80", linewidth = 0.4),

  ) + 
  annotate("text", x = 0, y = 41, label = spearman_text, 
         hjust = 0, vjust = 1, size = 5.7, family = "Times New Roman")

gg_mae

ggsave(path = Extra_figures_filepath, filename = "3b_MAE_Scatter.png", device = "png", plot = gg_mae, width = 8, height = 5, dpi = 300, units = "in")
```

### Log-transformed y-axis
```{r}
gg_mae_log <- ggplot(alg_match_clean, aes(x = Category_nHours, y = MAE, color = Method_Group)) +
  geom_point(size = 3, alpha = 0.8) +
  theme_classic() +
  labs(
    x = "No. of days since sampling",
    y = expression(log[10]*"(mean absolute error) ("*mu*"g L⁻¹)"),
    color = "In situ chl-a \nprocessing method"
  ) +
  scale_color_manual(values = c(
    "Fluorescence" = "steelblue",
    "Extraction" = "darkgreen",
    "Remote" = "darkorange",
    "Unspecified" = "darkorange",
    "Other" = "purple"
  )) +
  scale_y_log10() +
  scale_x_continuous(breaks = seq(0, 30, 5), limits = c(0, 30)) +
  theme(
    axis.text.x = element_text(size = 16),
    axis.text.y = element_text(size = 16),
    axis.title.x = element_text(size = 18),
    axis.title.y = element_text(size = 18),
    plot.title = element_text(hjust = 0.5, size = 20),
    legend.text = element_text(size = 16),
    legend.title = element_text(size = 17),
    legend.spacing.y = unit(0.25, 'cm'),
    text = element_text(family = "Times New Roman"),
    legend.position = c(0.84, 0.2),
    panel.border = element_rect(color = "black", fill = NA, linewidth = 0.5),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    legend.background = element_rect(fill = "transparent", color = "white", linewidth = 0.4)) + 
  annotate("text", x = 12, y = 0.5, label = spearman_text, 
         hjust = 0, vjust = 1, size = 5.7, family = "Times New Roman")

gg_mae_log

ggsave(path = Paper_figures_filepath, filename = "3_MAE_ScatterPlot_logY.png", device = "png", plot = gg_mae_log, width = 8, height = 5, dpi = 600, units = "in")

#PDF for cleaning for publication in adobe illustartor later
#ggsave(
#  filename = file.path(Figure_filepath, "2b_MAE_ScatterPlot_logY.pdf"),
#  plot = gg_mae_log,
#  device = cairo_pdf,
#  width = 8, height = 5, units = "in")

```
