mhw-flux.Rmd

---
title: "MHWs vs. heat flux"
author: "Robert Schlegel"
date: "2020-02-25"
output: workflowr::wflow_html
editor_options:
  chunk_output_type: console
csl: FMars.csl
bibliography: MHWflux.bib
---

```{r global_options, include = FALSE}
knitr::opts_chunk$set(fig.width = 8, fig.align = 'center',
                      echo = TRUE, warning = FALSE, message = FALSE, 
                      eval = TRUE, tidy = FALSE)
```

## Introduction

This vignette will walk through the thinking and the process for how to link physical variables to their potential effect on driving the onset and/or decline of MHWs. The primary source that inspired this work was @Chen2016. In this paper the authors were able to illustrate which parts of the heat budget were most likely driving the anomalous heat content in the surface of the ocean. What this analysis seeks to do is to build on this methodology by applying the fundamental concept to ALL of the MHWs detected in the NW Atlantic. Fundamentally we are running thousands of RMSE and correlation calculations between SST anomalies and the co-occurring anomalies for a range of physical variables. The lower the RMSE and the stronger the correlation (both positive and negative) the more of an indication this is to us that these phenomena are related. We will also be comparing the magnitude of change for T_Qnet and SSTa during the onset and decline of MHWs to see how much of these events can be attributed to heat flux.

```{r startup}
# All of the libraries and objects used in the project
# Note that this also loads the data we will be using in this vignette
source("code/functions.R")
```

## Relating SSTa to physical variables

We know when MHWs occurred, and our physical data are prepped, so what we need to do is run magnitude comparisons, RMSE, and correlations between SSTa from the start to peak and peak to end of each event for the full suite of variables. This will show us for each event which values related the best for the onset AND decline of the events. We will run RMSE and correlations on the full time series, but these were decided not to be included in the final manuscript as they don't add much other than to show how much better the calculations for onset or decline alone are.

The magnitude of change in SSTa and T_Qnet for the onset and decline of MHWs will be calculated by finding how much SSTa and T_Qnet increase from the onset to the peak of a MHW, and how much SSTa and T_Qnet decrease on the decline of the same event. The thinking here is that if SSTa exceeds the T_Qnet increase for the onset, this will indicate that some advective process is aiding the development of the SSTa, but if T_Qnet exceeds SSTa, this indicates that the heat flux term is responsible for the creation of the event. Likewise for the decline of events, if the T_Qnet value decreases more than SSTa this would mean that negative heat flux into the ocean was terminating the event, whereas a greater decrease in SSTa would imply that advection was primarily terminating the event.

Running correlations is useful to see how well certain variables change in relation to SSTa, but these do not show how similar those proportions of change are. In order to do this we need to calculate the RMSE between SSTa and the T_Qx variables. This may only be done with the T_Qx variables because they are in the same units as SST. One additional tweak we will need here is to add the cumulative T_Qx term for the following day to the SSTa from the first day of the MHW. This is done so that we may see more closely how the developing T_Qx term matches to the development of the SSTa. Effectively this tells us how much the T_Qx term may be driving SSTa.

```{r MHW-var-cor, eval=FALSE}
# Event index used for calculations
OISST_MHW_event_index <- OISST_MHW_event %>% 
  ungroup() %>% 
  mutate(row_index = 1:n())

# Run all the stats
ALL_cor <- plyr::ddply(OISST_MHW_event_index, c("row_index"), stats_all, .parallel = T) %>% 
  dplyr::select(-row_index) %>% 
  arrange(region, event_no) %>% 
  mutate(Parameter2 = factor(Parameter2))

# Save
saveRDS(ALL_cor, "data/ALL_cor.Rda")
saveRDS(ALL_cor, "shiny/ALL_cor.Rda")
```

What we have now is a long dataframe containing the magnitude of change, RMSE, and correlations of different variables with SSTa. It must be pointed out that the correlations for the non-Qx terms are for the same day, there is no time lag introduced, which may be important. Below we are going to visualise the range of correlations for each variable to see how much each distribution is skewed. This skewness could probably be quantified in a meaningful way... but let's look at the data first.

```{r shiny-histo}
# source("shiny/app.R")
# Or it is live here:
# https://robert-schlegel.shinyapps.io/MHWflux/
```

There are some really clear patterns coming through in the data. In particular SSS seems to be strongly related to the onset of MHWs. There are a lot of nuances in these data and so I think this is actually an example of where a shiny app is useful to interrogate the data.

In the shiny app it also comes out that the longer events tend not to correlate strongly with a single variable. This is to be expected and supports the argument that very persistent MHWs are supported by a confluence of variables. How to parse that out is an interesting challenge.

## 12 hour offset

While going through the data it was discovered that the OISST temperature values represent an average from midnight to midnight. This therefore puts them offset to the Qx variables, which were also being calculated over the same time period. But to compare the Qx terms to SST there should be a 12 hour offset, with the Qx variables preceding the SST, thereby representing the integral of time between daily SST values. So 12 hours were added to the heat flux variables before re-running the pipeline and recalculating all of the magnitude, RMSE, and correlation stats. In the code chunk below we compare the RMSE stats pre-12 hour shift and afterwards.

```{r shift-compare}
# Load old and new data
ALL_cor_old <- readRDS("data/ALL_cor_old.Rda") %>% 
  mutate(region = toupper(region)) %>% 
  mutate(run = "old") %>% 
  dplyr::select(-row_index)
ALL_cor <- readRDS("data/ALL_cor.Rda") %>% 
  dplyr::select(region:n_Obs, rmse) %>% 
  mutate(run = "new")
ALL_cor_stack <- rbind(ALL_cor, ALL_cor_old) %>% 
  filter(rmse > 0) %>% 
  mutate(Parameter2 = as.character(Parameter2),
         Parameter2 = case_when(Parameter2 == "qnet_budget" ~ "qnet",
                                Parameter2 == "lhf_budget" ~ "lhf",
                                Parameter2 == "shf_budget" ~ "shf",
                                Parameter2 == "lwr_budget" ~ "lwr",
                                Parameter2 == "swr_budget" ~ "swr",
                                Parameter2 == "qnet_mld_cum" ~ "qnet",
                                Parameter2 == "lhf_mld_cum" ~ "lhf",
                                Parameter2 == "shf_mld_cum" ~ "shf",
                                Parameter2 == "lwr_mld_cum" ~ "lwr",
                                Parameter2 == "swr_mld_cum" ~ "swr",
                                TRUE ~ Parameter2))

# Compare RMSE quantiles
ALL_cor_stack %>% 
  filter(rmse > 0) %>% 
  group_by(run, region, ts, Parameter2) %>% 
  summarise(q10 = quantile(rmse, 0.10), 
            q50 = quantile(rmse, 0.50),
            q90 = quantile(rmse, 0.90)) %>% 
  pivot_wider(names_from = run, values_from = q10:q90)

# Plot difference
ALL_cor_stack %>% 
  filter(rmse > 0) %>% 
  ggplot(aes(x = ts, y = rmse)) +
  geom_boxplot(aes(fill = run)) +
  facet_grid(Parameter2~region) #+
  # scale_y_continuous(limits = c(0, 2))
```

The new an improved method that correctly integrates Qx into SSTa creates much lower RMSE values for almost every variable and region. A smashing success.

## Small events and the shoaling MLD

An observation was made by one of the reviewers of the manuscript that there appear to be many small MHWs that occur close together. Events that barely peak out above the 90th percentile threshold before dropping back down. More importantly, some of these events appear to be close together in time when doing this, implying that some sort of preconditioning is at play, which was hypothesised to be a continuously shoaled MLD. While we acknowledge that the methods in this manuscript do not attempt to tackle the issue of preconditioning, it is certainly worth investigating the role of the MLD with these MHWs that appear perhaps as a skipping stone, hoping shallowly and rapidly along the curve of the 90th percentile threshold. In the following code we will find these shallow rapid events and relate them to the potential ongoing presence of an anomolously shallow (shoaled) MLD.

```{r MLD-shoal-1}
# Find the gaps in days between events
dist_days <- function(df){
  res_df <- data.frame()
  for(i in 1:nrow(df)-1){ # A for loop? I know right... is this even tidy code anymore!?
    df_sub <- df[c(i,i+1),]
    df_res <- data.frame(event_1 = df_sub$event_no[1],
                         event_2 = df_sub$event_no[2],
                         days = df_sub$date_start[2] - df_sub$date_end[1])
    res_df <- rbind(res_df, df_res)
  }
  res_df <- na.omit(res_df)
  return(res_df)
}
event_day_dist <- plyr::ddply(OISST_MHW_event, c("region"), dist_days)

# Plot the distribution of differences in days between events
event_day_dist %>% 
  group_by(region) %>% 
  mutate(count = n()) %>% 
  filter(days <= 366) %>% # Let's ignore events more than a year apart
  mutate(count_filter = n()) %>%
  ggplot() +
  geom_histogram(aes(x = as.integer(days), fill = region), 
                 bins = 12, show.legend = F) +
  geom_label(aes(x = 200, y = 10, label = paste0(count_filter,"/",count))) +
  facet_wrap(~region)
```

Of the 291 MHWs detected in the data for this study, 254 of them occur within 366 days of another MHW in the same region. As we may see from the histogram above, most of these events are occurring within two months of another event within the same region. Let's zoom in on this behaviour and visualise the distribution of gaps between events occurring not more than two months apart.

```{r MLD-shoal-2}
event_day_dist %>% 
  group_by(region) %>% 
  mutate(count = n()) %>% 
  filter(days <= 62) %>% # Let's ignore events more than a year apart
  mutate(count_filter = n()) %>%
  ggplot() +
  geom_histogram(aes(x = as.integer(days), fill = region), 
                 bins = 10, show.legend = F) +
  geom_label(aes(x = 30, y = 5, label = paste0(count_filter,"/",count))) +
  facet_wrap(~region)
```

The above histogram shows us that for the events that are occurring within two months of another month, there is a tendency for the events to occur within one or two weeks of another event. For our main goal here of relating temporally adjacent MHWs to a shoaled MLD we will use a two week filter. This limits the sample to 83 MHWs.

```{r MLD-shoal-3}
# Subset events
event_day_dist_sub <- event_day_dist %>% 
  filter(days <= 14) %>% 
  mutate(plyr_idx = 1:n()) # For parallel processing
  # dplyr::select(-days) %>% 
  # pivot_longer(event_1:event_2, names_to = "event_no")

# Function for calculating the MLD anomaly between neighbouring MHWs
MLD_neighbour <- function(df){
  
  # Filter out the two events being compared
  event_1 <- filter(OISST_MHW_event, 
                    region == df$region[1], event_no == df$event_1[1])
  event_2 <- filter(OISST_MHW_event, 
                    region == df$region[1], event_no == df$event_2[1])
  
  # Find the dates separating the events
  # We are making a conscious choice to also include the last and first days of the neighbouring events
  # This will account for the ocean state at the very end and beginning of the events
  event_gap <- seq(event_1$date_end, event_2$date_start, by = "day")
  
  # Extract the MLD data and calculated the cumulative anomaly
  MLD_cum_anom <- filter(ALL_ts_anom, var == "mld",
                         region == df$region[1], t %in% event_gap) %>% 
    summarise(count = n(),
              anom_cum = sum(anom),
              anom_mean = mean(anom))
  
  # Combine and exit
  res <- cbind(df, MLD_cum_anom)
  res$event_1_dur <- event_1$duration[1]
  res$event_2_dur <- event_2$duration[1]
  return(res)
}
MLD_res <- plyr::ddply(event_day_dist_sub, c("plyr_idx"), MLD_neighbour, .parallel = T)
```

What do the mean and cumulative MLD anomalies look like in the gaps between MHWs occurring within two weeks of each other??

```{r MLD-shoal-4}
# Cumulative values
hist_cum <- MLD_res %>% 
  ggplot() +
  geom_histogram(aes(x = anom_cum, fill = region), 
                 bins = 10, show.legend = F) +
  facet_wrap(~region)

# Mean values
hist_mean <- MLD_res %>% 
  ggplot() +
  geom_histogram(aes(x = anom_mean, fill = region), 
                 bins = 10, show.legend = F) +
  facet_wrap(~region)

# Combine intoa single plot
ggpubr::ggarrange(hist_cum, hist_mean, ncol = 1)
```

Above we have plotted histograms of the cumulative and mean MLD anomaly on days between MHWs that occur within two weeks of each other. We have shown the mean values because the cumulative values of events with larger gaps will likely tend to be larger than shorter gaps, even if the MLD may be shoaled more in those gaps. The reader is free to draw their conclusions from either result, but they both show a similar picture, which is that the MLD anomalies during the gaps between nearby MHWs skew towards positive (deeper) depth anomalies. The exception to this being the Mid-Atlantic Bight, and the Newfoundland Shelf to a lesser extent.

Perhaps there is a relationship between the length of the gap and the MLD anomaly?

```{r MLD-shoal-5}
MLD_res %>% 
  ggplot(aes(x = as.integer(days), y = anom_cum)) +
  geom_point(aes(colour = region)) +
  geom_smooth(method = "lm", se = F, aes(colour = region))
  # geom_histogram(aes(x = anom_cum, fill = region), 
                 # bins = 10, show.legend = F) +
  # facet_wrap(~region)
```

But this doesn't quite answer the issue that was posed. Remember that specifically the concern was about very short event occurring close in time to each other. So as a final step we will filter the results from above to only show when at least one of the two events was shorter than 10 days (double the absolute minimum length of five days).

```{r MLD-shoal-6}

```


## Results

### Magnitude of change summary

In this sub-section we will look at how the different magnitudes of change for onset and decline of the T_Qnet terms compare to the same changes in SSTa.

```{r mag-summary}
# Load magnitude data
ALL_mag <- readRDS("data/ALL_cor.Rda") %>% 
  dplyr::select(region:ts, Parameter2, n_Obs, mag) %>% 
  na.omit()

# Calculate the proportions of change in magnitudes
ALL_mag_prop <- ALL_mag %>% 
  filter(Parameter2 != "sst") %>% 
  left_join(filter(ALL_mag, Parameter2 == "sst"), 
            by = c("region", "season", "event_no", "ts", "n_Obs")) %>% 
  filter(ts != "full") %>% 
  dplyr::select(-Parameter2.y) %>% 
  dplyr::rename(var = Parameter2.x, mag_Qx = mag.x, mag_SSTa = mag.y) %>% 
  mutate(prop = mag_Qx/mag_SSTa) %>% 
  mutate(var = as.character(var),
         var = case_when(var == "lwr_budget" ~ "Qlw",
                         var == "swr_budget" ~ "Qsw",
                         var == "lhf_budget" ~ "Qlh",
                         var == "shf_budget" ~ "Qsh",
                         var == "qnet_budget" ~ "Qnet",
                         TRUE ~ var),
         var = factor(var, levels = c("Qnet", "Qlh", "Qsh", "Qlw", "Qsw")),
         region = toupper(region))

# Find an event with good onset and decline match
ALL_mag_prop %>% filter(var == "Qnet", prop >= 1, prop <= 2)

# Create schematic figure showing how magnitudes are calculated
p1 <- fig_mag_func(2, "GSL") +
  ggtitle("Good onset, bad decline")
p2 <- fig_mag_func(44, "GM") +
  ggtitle("Bad onset, good decline")
p3 <- fig_mag_func(32, "GSL") +
  ggtitle("Good onset, good decline")
mag_schematic <- ggpubr::ggarrange(p1, p2, p3, ncol = 3, nrow = 1, common.legend = TRUE)
# mag_schematic

# Scatterplot of SSTa vs SST_Qx (facets per variable)
mag_scatter <- ALL_mag_prop %>% 
  ggplot(aes(x = mag_SSTa, y = mag_Qx)) +
  geom_point(aes(colour = n_Obs, shape = ts)) +
  geom_smooth(aes(linetype = ts), method = "lm", colour = "black") +
  scale_colour_gradientn(colours = rainbow(10)) +
  facet_wrap(~var, nrow = 1) +
  labs(x = "ΔSSTa", y = "ΔSST_Qx", shape = "Time series", linetype = "Time series") +
  theme(legend.position = "top")
# mag_scatter

# Boxplots of proportions of change for SST_Qx for onset/decline
mag_box <- ALL_mag_prop %>% 
  ggplot(aes(x = ts, y = prop)) +
  geom_hline(aes(yintercept = 0), colour = "red") +
    geom_boxplot(aes(fill = season), outlier.shape = NA) +
  # geom_boxplot(aes(fill = region), outlier.shape = NA) +
  scale_y_continuous(limits = c(-2, 2)) +
  facet_wrap(~var, nrow = 1) +
  labs(x = NULL, y = "ΔSST_Qx / ΔSSTa")
# mag_box

# Combine for one massive display
mag_full <- ggpubr::ggarrange(mag_schematic, mag_scatter, mag_box, ncol = 1, nrow = 3, labels = c("A)", "B)", "C)"))
ggsave("output/magnitude_schematic.png", height = 8, width = 10)
ggsave("output/magnitude_schematic.pdf", height = 8, width = 10)
mag_full
```

The first thing we want to look at with these magnitude scores is how many MHWs appear to be driven or decayed by Qnet. We determine this if the proportion of change in Qnet is more than half that of the change in SSTa. But also, because the Qx variables can be greater than the sum (Qnet), we will see when individual variables are contributing more than half of the change in SSTa, too.

```{r mag-prop}
# Create TRUE/FALSE table when prop is over 0.5 for any variable
ALL_mag_TF <- ALL_mag_prop %>% 
  mutate(prop_TF = ifelse(prop > 0.5, TRUE, FALSE)) %>% 
  dplyr::select(region:n_Obs, prop_TF) %>% 
  pivot_wider(names_from = var, values_from = prop_TF)

# Total proportion of MHWs forced by Qnet
sum(ALL_mag_TF$Qnet)/nrow(ALL_mag_TF)

# Total number of events driven by Qnet
nrow(filter(ALL_mag_TF, ts == "onset", Qnet == TRUE))

# Total number of events decayed by Qnet
nrow(filter(ALL_mag_TF, ts == "decline", Qnet == TRUE))

# Proportion of MHWs with Qnet inset and decline
ALL_mag_TF %>% 
  dplyr::select(region:ts, Qnet) %>% 
  pivot_wider(names_from = ts, values_from = Qnet) %>% 
  filter(onset == TRUE, decline == TRUE)

# Counts by time series phase
ALL_mag_count_ts <- ALL_mag_prop %>% 
  group_by(ts, var) %>% 
  mutate(ts_count = n(),
         group = "Total") %>% 
  filter(prop > 0.5) %>% 
  mutate(filter_count = n(),
         filter_prop = round(filter_count/ts_count, 2)) %>% 
  dplyr::select(group, ts, var, filter_prop) %>% 
  distinct() %>% 
  mutate(var = factor(var, levels = c("Qnet", "Qlh", "Qsh", "Qlw", "Qsw"))) %>% 
  pivot_wider(names_from = var, values_from = filter_prop) %>% 
  dplyr::select(group, ts, Qnet, Qlh, Qsh, Qlw, Qsw)
# knitr::kable(ALL_mag_count_ts)

# Counts by time region phase
ALL_mag_count_region <- ALL_mag_prop %>% 
  group_by(region, ts, var) %>% 
  mutate(ts_count = n()) %>% 
  filter(prop > 0.5) %>% 
  mutate(filter_count = n(),
         filter_prop = round(filter_count/ts_count, 2)) %>% 
  dplyr::select(region, ts, var, filter_prop) %>% 
  distinct() %>% 
  pivot_wider(names_from = var, values_from = filter_prop) %>% 
  dplyr::select(region, ts, Qnet, Qlh, Qsh, Qlw, Qsw) %>% 
  arrange(region, ts) %>% 
  dplyr::rename(group = region)
# knitr::kable(ALL_mag_count_region)

# Counts by time series phase
ALL_mag_count_season <- ALL_mag_prop %>% 
  group_by(season, ts, var) %>% 
  mutate(ts_count = n()) %>% 
  filter(prop > 0.5) %>% 
  mutate(filter_count = n(),
         filter_prop = round(filter_count/ts_count, 2)) %>% 
  dplyr::select(season, ts, var, filter_prop) %>% 
  distinct() %>% 
  pivot_wider(names_from = var, values_from = filter_prop) %>% 
  dplyr::select(season, ts, Qnet, Qlh, Qsh, Qlw, Qsw) %>% 
  arrange(season, ts) %>% 
  dplyr::rename(group = season)
# knitr::kable(ALL_mag_count_season)
ALL_mag_count_ts_region_season <- rbind(ALL_mag_count_ts, ALL_mag_count_region, ALL_mag_count_season) %>% 
  dplyr::select(group:Qnet) %>% 
  pivot_wider(names_from = ts, values_from = Qnet) %>% 
  dplyr::select(group, onset, decline) %>% 
  mutate(onset = paste0(onset*100,"%"),
         decline = paste0(decline*100,"%"))
knitr::kable(ALL_mag_count_ts_region_season)
```

We also want a table showing the count of best magnitude match per variable per region/season.

```{r mag-table}
# Get the top results other than Qnet
ALL_mag_top <- ALL_mag_prop %>% 
  filter(var != "Qnet") %>% # Remove Qnet
  group_by(region, season, event_no, ts) %>% 
  filter(prop == max(prop))

# The top count by region
ALL_mag_region <- ALL_mag_top %>% 
  group_by(region, ts, var) %>% 
  summarise(count = n(), .groups = "drop") %>% 
  pivot_wider(names_from = var, values_from = count) %>% 
  mutate(Qlw = replace_na(Qlw, 0)) %>% 
  left_join(count_region, by = c("region", "ts")) %>% 
  dplyr::select(region, ts, count, Qlh:Qsw) %>% 
  mutate(Qlh = paste0(Qlh, " (",round((Qlh/count)*100),"%)"),
         Qsh = paste0(Qsh, " (",round((Qsh/count)*100),"%)"),
         Qlw = paste0(Qlw, " (",round((Qlw/count)*100),"%)"),
         Qsw = paste0(Qsw, " (",round((Qsw/count)*100),"%)")) %>% 
  dplyr::rename(group = region)
# knitr::kable(ALL_mag_region)

# The top count by season
ALL_mag_season <- ALL_mag_top %>% 
  group_by(season, ts, var) %>% 
  summarise(count = n(), .groups = "drop") %>% 
  pivot_wider(names_from = var, values_from = count) %>% 
  mutate(Qsh = replace_na(Qsh, 0)) %>%
  left_join(count_season, by = c("season", "ts")) %>% 
  dplyr::select(season, ts, count, Qlh:Qsw) %>% 
  mutate(Qlh = paste0(Qlh, " (",round((Qlh/count)*100),"%)"),
         Qsh = paste0(Qsh, " (",round((Qsh/count)*100),"%)"),
         Qlw = paste0(Qlw, " (",round((Qlw/count)*100),"%)"),
         Qsw = paste0(Qsw, " (",round((Qsw/count)*100),"%)")) %>% 
  dplyr::rename(group = season)
# knitr::kable(ALL_mag_season)

# Print table
ALL_mag_region_season <- rbind(ALL_mag_region, ALL_mag_season)
knitr::kable(ALL_mag_region_season)
```

These tables help to reiterate how much more important Qlh is for the onset and decline of MHWs over the other variables, and that Qsw often comes in second. For the seasons though we see a much more complex relationship that I won't go into here.

### RMSE summary

Here we will go about summarising the RMSE results. We will also create some boxplots to help with the process.

```{r RMSE-boxplot}
# The base RMSE results
ALL_RMSE <- ALL_cor %>% 
  filter(rmse > 0 ) %>% 
  dplyr::rename(var = Parameter2) %>% 
  dplyr::select(region:ts, var, n_Obs, rmse) %>% 
  mutate(var = as.character(var),
         var = case_when(var == "lwr_budget" ~ "Qlw",
                         var == "swr_budget" ~ "Qsw",
                         var == "lhf_budget" ~ "Qlh",
                         var == "shf_budget" ~ "Qsh",
                         var == "qnet_budget" ~ "Qnet",
                         TRUE ~ var),
         var = factor(var, levels = c("Qnet", "Qlh", "Qsh", "Qlw", "Qsw")),
         region = toupper(region)) %>% 
  # filter out events where Qnet is not a primary drivers
  left_join(ALL_mag_TF, by = c("region", "season", "event_no", "ts", "n_Obs")) %>% 
  filter(Qnet == TRUE)

# Summaries by region
ALL_RMSE_region <- ALL_RMSE %>% 
  dplyr::select(-season, -event_no, -n_Obs) %>% 
  group_by(region, ts, var) %>% 
  summarise(rmse_min = min(rmse),
           rmse_mean = mean(rmse),
           rmse_max = max(rmse), .groups = "drop")

# Boxplot of region summaries
ALL_RMSE %>% 
  dplyr::select(-season, -event_no, -n_Obs) %>% 
  ggplot(aes(x = var, y = rmse)) +
  geom_boxplot(aes(fill = ts)) +
  facet_wrap(~region) #+
  # theme(axis.text.x = element_text(angle = 45))

# Summaries by season
ALL_RMSE_season <- ALL_RMSE %>% 
  dplyr::select(-region, -event_no, -n_Obs) %>% 
  group_by(season, ts, var) %>% 
  summarise(rmse_min = min(rmse),
           rmse_mean = mean(rmse),
           rmse_max = max(rmse), .groups = "drop")

# Boxplot of season results
ALL_RMSE %>% 
  dplyr::select(-region, -event_no, -n_Obs) %>% 
  ggplot(aes(x = var, y = rmse)) +
  geom_boxplot(aes(fill = ts)) +
  facet_wrap(~season)

# Boxplot by variable
ALL_RMSE %>% 
  filter(var != "Qnet") %>% 
  ggplot(aes(x = region, y = rmse)) +
  geom_boxplot(aes(fill = ts)) +
  facet_wrap(~var)
ALL_RMSE %>% 
    filter(var != "Qnet") %>% 
  ggplot(aes(x = season, y = rmse)) +
  geom_boxplot(aes(fill = ts)) +
  facet_wrap(~var)
```

These boxplots help to tell a high level story about the importance of the Qx variables with SSTa, but I want a slightly more fine-grain level of results. To do this I am going to count which of the Qx variables has the lowest RMSE per MHW. These will then be grouped by region and season to provide an understanding for how this may vary.

```{r RMSE-table}
# Get the top results other than Qnet
ALL_RMSE_top <- ALL_RMSE %>% 
  filter(var != "Qnet") %>% # Remove Qnet
  group_by(region, season, event_no, ts) %>% 
  filter(rmse == min(rmse))

# The top count overall
ALL_RMSE_ts <- ALL_RMSE_top %>%
  group_by(ts) %>%
  mutate(total_count = n(),
         group = "Total") %>% 
  group_by(group, ts, total_count, var) %>% 
  summarise(count = n(), .groups = "drop") %>% 
  pivot_wider(names_from = var, values_from = count) %>% 
  replace_na(list(Qlh = 0, Qsh = 0, Qlw = 0, Qsw = 0)) %>% 
  dplyr::select(group, ts, total_count, Qlh:Qsw) %>% 
  filter(ts != "full") %>% 
  mutate(Qlh = paste0(Qlh, " (",round((Qlh/total_count)*100),"%)"),
         Qsh = paste0(Qsh, " (",round((Qsh/total_count)*100),"%)"),
         Qlw = paste0(Qlw, " (",round((Qlw/total_count)*100),"%)"),
         Qsw = paste0(Qsw, " (",round((Qsw/total_count)*100),"%)"))

# The top count by region
ALL_RMSE_region <- ALL_RMSE_top %>%
  group_by(region, ts) %>%
  mutate(total_count = n()) %>% 
  group_by(region, ts, total_count, var) %>% 
  summarise(count = n(), .groups = "drop") %>% 
  pivot_wider(names_from = var, values_from = count) %>% 
  replace_na(list(Qlh = 0, Qsh = 0, Qlw = 0, Qsw = 0)) %>% 
  dplyr::select(region, ts, total_count, Qlh:Qlw) %>% 
  filter(ts != "full") %>% 
  mutate(Qlh = paste0(Qlh, " (",round((Qlh/total_count)*100),"%)"),
         Qsh = paste0(Qsh, " (",round((Qsh/total_count)*100),"%)"),
         Qlw = paste0(Qlw, " (",round((Qlw/total_count)*100),"%)"),
         Qsw = paste0(Qsw, " (",round((Qsw/total_count)*100),"%)")) %>% 
  dplyr::rename(group = region)
# knitr::kable(ALL_RMSE_region)

# The top count by season
ALL_RMSE_season <- ALL_RMSE_top %>% 
  group_by(season, ts) %>%
  mutate(total_count = n()) %>% 
  group_by(season, ts, total_count, var) %>% 
  summarise(count = n(), .groups = "drop") %>% 
  pivot_wider(names_from = var, values_from = count) %>% 
  replace_na(list(Qlh = 0, Qsh = 0, Qlw = 0, Qsw = 0)) %>% 
  dplyr::select(season, ts, total_count, Qlh:Qlw) %>% 
  filter(ts != "full") %>% 
  mutate(Qlh = paste0(Qlh, " (",round((Qlh/total_count)*100),"%)"),
         Qsh = paste0(Qsh, " (",round((Qsh/total_count)*100),"%)"),
         Qlw = paste0(Qlw, " (",round((Qlw/total_count)*100),"%)"),
         Qsw = paste0(Qsw, " (",round((Qsw/total_count)*100),"%)")) %>% 
  dplyr::rename(group = season)
# knitr::kable(ALL_RMSE_season)

# Print table
ALL_RMSE_ts_region_season <- rbind(ALL_RMSE_ts, ALL_RMSE_region, ALL_RMSE_season) %>% 
  dplyr::rename(count = total_count)
knitr::kable(ALL_RMSE_ts_region_season)
```

I'm pretty happy with how these results are displayed. In the last chunk in this sub-section we will grab some specific summary stats that are used in the publication.

```{r RMSE-stats}
# Short event best fits
ALL_RMSE_short <- ALL_RMSE %>% 
  filter(n_Obs <= 14, var != "Qnet") %>% 
  group_by(region, season, event_no, ts) %>% 
  filter(rmse == min(rmse))

# Short event index
ALL_RMSE_long <- ALL_RMSE %>% 
  filter(n_Obs > 14, var != "Qnet") %>% 
  group_by(region, season, event_no, ts) %>% 
  filter(rmse == min(rmse))

# The top Qx for events of 14 days or less by region
ALL_RMSE_short_region <- ALL_RMSE_short %>%
  group_by(region, ts) %>% 
  mutate(count = n()) %>% 
  group_by(region, ts, count, var) %>% 
  summarise(count_Q = n(), .groups = "drop") %>% 
  pivot_wider(names_from = var, values_from = count_Q) %>% 
  replace_na(list(Qlh = 0, Qsh = 0, Qlw = 0, Qsw = 0)) %>% 
  dplyr::select(region, ts, count, Qlh:Qlw) %>% 
  mutate(Qlh = paste0(Qlh, " (",round((Qlh/count)*100),"%)"),
         Qsh = paste0(Qsh, " (",round((Qsh/count)*100),"%)"),
         Qlw = paste0(Qlw, " (",round((Qlw/count)*100),"%)"),
         Qsw = paste0(Qsw, " (",round((Qsw/count)*100),"%)"))
knitr::kable(ALL_RMSE_short_region)

# The top Qx for events of more than 14 days by region
ALL_RMSE_long_region <- ALL_RMSE_long %>%
  group_by(region, ts) %>% 
  mutate(count = n()) %>% 
  group_by(region, ts, count, var) %>% 
  summarise(count_Q = n(), .groups = "drop") %>% 
  pivot_wider(names_from = var, values_from = count_Q) %>% 
  replace_na(list(Qlh = 0, Qsh = 0, Qlw = 0, Qsw = 0)) %>% 
  dplyr::select(region, ts, count, Qlh:Qsh) %>% 
  mutate(Qlh = paste0(Qlh, " (",round((Qlh/count)*100),"%)"),
         Qsh = paste0(Qsh, " (",round((Qsh/count)*100),"%)"),
         Qlw = paste0(Qlw, " (",round((Qlw/count)*100),"%)"),
         Qsw = paste0(Qsw, " (",round((Qsw/count)*100),"%)"))
knitr::kable(ALL_RMSE_long_region)

# The top Qx for events of 14 days or less by season
ALL_RMSE_short_season <- ALL_RMSE_short %>%
  group_by(season, ts) %>% 
  mutate(count = n()) %>% 
  group_by(season, ts, count, var) %>% 
  summarise(count_Q = n(), .groups = "drop") %>% 
  pivot_wider(names_from = var, values_from = count_Q) %>% 
  replace_na(list(Qlh = 0, Qsh = 0, Qlw = 0, Qsw = 0)) %>% 
  dplyr::select(season, ts, count, Qlh:Qlw) %>% 
  mutate(Qlh = paste0(Qlh, " (",round((Qlh/count)*100),"%)"),
         Qsh = paste0(Qsh, " (",round((Qsh/count)*100),"%)"),
         Qlw = paste0(Qlw, " (",round((Qlw/count)*100),"%)"),
         Qsw = paste0(Qsw, " (",round((Qsw/count)*100),"%)"))
knitr::kable(ALL_RMSE_short_season)

# The top Qx for events of more than 14 days by season
ALL_RMSE_long_season <- ALL_RMSE_long %>%
  group_by(season, ts) %>% 
  mutate(count = n()) %>% 
  group_by(season, ts, count, var) %>% 
  summarise(count_Q = n(), .groups = "drop") %>% 
  pivot_wider(names_from = var, values_from = count_Q) %>% 
  replace_na(list(Qlh = 0, Qsh = 0, Qlw = 0, Qsw = 0)) %>% 
  dplyr::select(season, ts, count, Qlh:Qlw) %>% 
  mutate(Qlh = paste0(Qlh, " (",round((Qlh/count)*100),"%)"),
         Qsh = paste0(Qsh, " (",round((Qsh/count)*100),"%)"),
         Qlw = paste0(Qlw, " (",round((Qlw/count)*100),"%)"),
         Qsw = paste0(Qsw, " (",round((Qsw/count)*100),"%)"))
knitr::kable(ALL_RMSE_long_season)
```

We also want to see how RMSE scores change the shorter the phase of the MHW is.

```{r RMSE-short-change}
# Faceted scatter plot showing how RMSE changes based on ts length
ggplot(ALL_RMSE, aes(x = n_Obs, y = rmse)) +
  geom_point(aes(colour = ts)) +
  geom_smooth(aes(colour = ts), method = "lm") +
  facet_wrap(~region)

# Faceted scatter plot showing how RMSE changes based on ts length
ggplot(ALL_RMSE, aes(x = n_Obs, y = rmse)) +
  geom_point(aes(colour = ts)) +
  geom_smooth(aes(colour = ts), method = "lm") +
  facet_wrap(~season)
```

### Correlations: summary, regions, seasons

With the correlations calculated for the onset, decline, and full extent of each MHW, we want to know if any signals emerge from the regions and/or seasons of occurrence of these events. Is the relationship between SSS and MHW onset stronger in the winter? Stronger in certain region? Having manually looked through the shiny app it does look like there are some patterns. Below is the code used to determine the stats referred to in the paper.

```{r}
# Overall correlation results
cor_summary_all <- ALL_cor %>% 
  filter(Parameter1 == "sst",
         Parameter2 != "sst") %>% 
  group_by(ts, Parameter2) %>% 
  summarise(q10_r = quantile(r, 0.1),
            mean_r = mean(r, na.rm = T),
            q90_r = quantile(r, 0.9),
            mean_rmse = mean(rmse, na.rm = T))

# Overall positive only results
cor_positive_all <- ALL_cor %>% 
  filter(Parameter1 == "sst",
         Parameter2 != "sst",
         r > 0) %>% 
  group_by(ts, Parameter2) %>% 
  summarise(q10_r = quantile(r, 0.1),
            mean_r = mean(r, na.rm = T),
            q90_r = quantile(r, 0.9),
            mean_rmse = mean(rmse, na.rm = T))

# Overall negative only results
cor_negative_all <- ALL_cor %>% 
  filter(Parameter1 == "sst",
         Parameter2 != "sst",
         r < 0) %>% 
  group_by(ts, Parameter2) %>% 
  summarise(q10_r = quantile(r, 0.1),
            mean_r = mean(r, na.rm = T),
            q90_r = quantile(r, 0.9),
            mean_rmse = mean(rmse, na.rm = T))

# region correlation results
cor_summary_region <- ALL_cor %>% 
  filter(Parameter1 == "sst",
         Parameter2 != "sst") %>% 
  group_by(ts, region, Parameter2) %>% 
  summarise(q10_r = quantile(r, 0.1),
            mean_r = mean(r, na.rm = T),
            q90_r = quantile(r, 0.9),
            mean_rmse = mean(rmse, na.rm = T))

# Look at specific results
cor_summary_region %>% 
  filter(Parameter2 == "mld_cum")

# region correlation results
cor_summary_season <- ALL_cor %>% 
  filter(Parameter1 == "sst",
         Parameter2 != "sst") %>% 
  group_by(ts, season, Parameter2) %>% 
  summarise(q10_r = quantile(r, 0.1),
            mean_r = mean(r, na.rm = T),
            q90_r = quantile(r, 0.9),
            mean_rmse = mean(rmse, na.rm = T))

# Look at specific results
cor_summary_season %>% 
  filter(Parameter2 == "mld")

# Test for normality of r distributions

```

### Relationships

With patterns pulled out by region and season, we want to see if there are any relationships between MHWs that show strong correlations at onset/decline with a particular Qx variable and strong correlations at onset/decline with an air/sea variable. We will look for this within regions and seasons as well. For example, do MHWs that correlate well with an increase in SSS also correlate well with a decrease in long-wave radiation during the decline of the event? I'm not sure how best to go about this in a clean manner.

Another thing to consider would be if fast onset slow decline (and vice versa) events have different characteristics to slower evolving events. The same question could be posed to long vs short events and those with high intensities vs low. In order to begin this investigation we must join the MHW results to the correlation results. We will visualise these patterns with heatmaps.

```{r metrics_cor}
# Strong positive correlation at onset and strong positive and decline and vice versa
ALL_cor_wide <- readRDS("data/ALL_cor.Rda") %>% 
  dplyr::rename(var = Parameter2) %>% 
  ungroup() %>% 
  filter(Parameter1 == "sst") %>% 
  dplyr::select(region:ts, var, r, n_Obs) %>% 
  mutate(var = as.character(var)) %>% 
  mutate(var = case_when(var == "sst" ~ "SST",
                         var == "bottomT" ~ "SBT",
                         var == "sss" ~ "SSS",
                         var == "mld_cum" ~ "MLD_c",
                         var == "mld_1_cum" ~ "MLD_1_c",
                         var == "t2m" ~ "Air_temp",
                         var == "tcc_cum" ~ "Cloud_cover_c",
                         var == "p_e_cum" ~ "Precip_Evap_c",
                         var == "mslp_cum" ~ "MSLP_c",
                         var == "lwr_budget" ~ "Qlw",
                         var == "swr_budget" ~ "Qsw",
                         var == "lhf_budget" ~ "Qlh",
                         var == "shf_budget" ~ "Qsh",
                         var == "qnet_budget" ~ "Qnet",
                         TRUE ~ var)) %>% 
  filter(var %in% c("SST", "SSS", "SBT", "MLD_c", "MLD_1_c", "Air_temp", "Cloud_cover_c",
                           "Precip_Evap_c", "MSLP_c", "Qlw", "Qsw", "Qlh", "Qsh", "Qnet")) %>% 
  pivot_wider(values_from = r, names_from = var)

# Combine MHW metrics and correlation results
events_cor_prep <- OISST_MHW_event %>% 
  dplyr::select(region, season, event_no, duration, intensity_mean, intensity_max, 
                intensity_cumulative, rate_onset, rate_decline) %>% 
  left_join(ALL_cor_wide, by = c("region", "season", "event_no"))

# Heatmap showing average correlations by MHW duration
events_cor_prep %>% 
  mutate(duration = plyr::round_any(duration, 10)) %>% 
  group_by(ts, duration) %>% 
  mutate(count = n()) %>% 
  summarise_if(is.numeric, mean) %>% 
  pivot_longer(cols = Air_temp:Qsw) %>% 
  filter(name != "sst",
         ts != "full",
         duration <= 50) %>% # Only 1 event is longer than this 
  ggplot(aes(x = duration, y = name)) +
  geom_tile(aes(fill = value)) +
  geom_text(aes(label = count)) +
  facet_wrap(~ts) +
  scale_fill_gradient2(low = "blue", high = "red") +
  coord_cartesian(expand = F) +
  labs(y = NULL, x = "Duration (10 day steps)", fill = "r (mean)") +
  theme(legend.position = "bottom")

# Heatmap showing average correlations by MHW max intensity
events_cor_prep %>% 
  mutate(intensity_max = plyr::round_any(intensity_max, 0.5)) %>% 
  group_by(ts, intensity_max) %>%
  mutate(count = n()) %>% 
  summarise_if(is.numeric, mean) %>% 
  pivot_longer(cols = Air_temp:Qsw) %>% 
  filter(name != "sst",
         ts != "full") %>% 
  ggplot(aes(x = intensity_max, y = name)) +
  geom_tile(aes(fill = value)) +
  geom_text(aes(label = count)) +
  facet_wrap(~ts) +
  scale_fill_gradient2(low = "blue", high = "red") +
  coord_cartesian(expand = F) +
  labs(y = NULL, x = "Max Intensity (°C; 0.5° steps)", fill = "r (mean)") +
  theme(legend.position = "bottom")

# Heatmap showing average correlations by MHW rate onset
events_cor_prep %>% 
  mutate(rate_onset = round(rate_onset, 1)) %>% 
  group_by(ts, rate_onset) %>% 
  mutate(count = n()) %>% 
  summarise_if(is.numeric, mean) %>% 
  pivot_longer(cols = Air_temp:Qsw) %>% 
  filter(name != "sst",
         ts != "full",
         rate_onset <= 0.75) %>% # Only two is faster than this 
  ggplot(aes(x = rate_onset, y = name)) +
  geom_tile(aes(fill = value)) +
  geom_text(aes(label = count)) +
  facet_wrap(~ts) +
  scale_fill_gradient2(low = "blue", high = "red") +
  coord_cartesian(expand = F) +
  labs(y = NULL, x = "Rate of onset (°C; 0.1° steps)", fill = "r (mean)") +
  theme(legend.position = "bottom")

# Heatmap showing average correlations by MHW rate decline
events_cor_prep %>% 
  mutate(rate_decline = round(rate_decline, 1)) %>% 
  filter(rate_decline <= 0.6) %>% # nip off a couple of outliers
  group_by(ts, rate_decline) %>% 
  mutate(count = n()) %>% 
  summarise_if(is.numeric, mean) %>% 
  pivot_longer(cols = Air_temp:Qsw) %>% 
  filter(name != "sst",
         ts != "full") %>% 
  ggplot(aes(x = rate_decline, y = name)) +
  geom_tile(aes(fill = value)) +
  geom_text(aes(label = count)) +
  facet_wrap(~ts) +
  scale_fill_gradient2(low = "blue", high = "red") +
  coord_cartesian(expand = F) +
  labs(y = NULL, x = "Rate of decline (°C; 0.1° steps)", fill = "r (mean)") +
  theme(legend.position = "bottom")
```

### Linearity of events

Another observation I made while going through the results in the shiny app was how the onset and decline of events tend to have a better RMSE when the SSTa trend is more linear. My hypothesis is that this is because there is less advective pressure on the dispersion of the heat anomaly, therefore more of the heatflux is responsible for the temperature anomaly. This is a pretty straight forward statement, but I think it is useful in that one may be able to take linearity of SSTa as a proxy for the proportion of heat flux that is responsible for it. So in the chunk below I use simple linear models to create an R2 value for each onset and decline portion of an event. Then I find how close to 1.0 that R2 value is (meaning more linear SSTa) and I see how RMSE changes as R2 improves. I hypothesised that more linear SSTa onset will provide better (lower) RMSE values. The same is likely true for decline but I'm not certain.

```{r linearity}
# R2 for RMSE vs R2
ALL_cor_R2 <- readRDS("data/ALL_cor.Rda") %>% 
  ungroup() %>% 
  filter(Parameter1 == "sst", rmse > 0) %>% 
  nest_by(region, season, ts) %>% 
  summarise(broom::glance(lm(rmse ~ sst_R2, data = data)))

# Scatterplots
readRDS("data/ALL_cor.Rda") %>% 
  ungroup() %>% 
  filter(Parameter1 == "sst", rmse > 0) %>% 
  ggplot(aes(x = sst_R2, y = rmse)) +
  geom_point() +
  geom_smooth(method = "lm") +
  facet_grid(region ~ Parameter2)

# Boxplots
readRDS("data/ALL_cor.Rda") %>% 
  ungroup() %>% 
  filter(Parameter1 == "sst", rmse > 0) %>% 
  mutate(sst_R2_cut = cut(sst_R2, breaks = c(0, 0.25, 0.5, 0.75,1.0))) %>% 
  na.omit() %>% 
  ggplot(aes(x = sst_R2_cut, y = rmse)) +
  geom_boxplot() +
  facet_grid(region ~ Parameter2)
```

Nope. If there is any relationship it is incredibly noisy. I'll not be pursuing this further for this project, but I'm not entirely dissuaded from pursuing this avenue of research in the future with a deeper learning approach.

### Choice events

There are a lot of results to wade through and though it is clear there are important signals in the results, but it is proving difficult to distill them. One thought is that we don't need to look at all of the events, just the longest/most intense events with strong r values. This is first done by cutting out all Category I events. We then find strong correlations with long events. There should just be a few.

Once this has been done we group events by their strongest Qx relationship, then find their strongest relationship with the next level of variables (e.g. MLD, MSLP, and so on). Ideally one may find the top four flavours.

```{r choice-events}
# Filter out smol events
events_cor_cat <- events_cor_prep %>% 
  left_join(OISST_MHW_cats[,c("region", "event_no", "category")], by = c("region", "event_no")) %>% 
  filter(category != "I Moderate", duration >= 21)

# Events with high Qlw correlations at onset
events_cor_cat %>% 
  filter(ts == "onset", Qlw >= 0.7)

# Melt the data frame and find the q term with the highest correlation
# Those are then used to separate events into groups
events_q_onset <- events_cor_cat %>% 
  filter(ts == "onset") %>% 
  dplyr::select(region:event_no, ts, Qnet:Qsh) %>% 
  pivot_longer(cols = Qnet:Qsh, names_to = "var", values_to = "val") %>% 
  group_by(region, season, event_no) %>% 
  filter(abs(val) == max(abs(val)))

# With this guide one may then parse out the q groups
events_lhf_p_onset <- events_q_onset %>% 
  filter(val >= 0, var == "Qlh") %>% 
  left_join(events_cor_cat) %>% 
  ungroup() %>% 
  # unite(region, season, event_no, sep = "_")
  mutate(event_no = as.factor(event_no)) %>% 
  dplyr::select(region:ts, Air_temp:Qsw) %>% 
  pivot_longer(Air_temp:Qsw) %>% 
  ggplot(aes(x = event_no, y = name)) +
  geom_tile(aes(fill = value)) +
  # facet_wrap(~name, scales = "free") +
  scale_fill_gradient2(low = "blue", high = "red") +
  coord_cartesian(expand = F) +
  labs(y = NULL, fill = "r (mean)") +
  theme(legend.position = "bottom")
events_lhf_p_onset
```

### Table

In the following table a more concise summary of the results is presented.

```{r, echo=FALSE, message=FALSE}
# NB: This table was created manually by going through the Shiny app one variable at a time.
res_table <- read_csv("data/res_table.csv")
knitr::kable(res_table, caption = "Most of the variables that have been correlated against the temperature anomalies during the onset, decline, and full duration of MHWs. The cumulative heat flux terms were corrected for by the daily MLD (Q/(rho x Cp x hmld)) before the correlations were calculated. Correlations were also run on the cumulative flux terms without correcting for MLD, but there was little difference so the results are not itemised here. This table shows the full names of the variables, as well as the abbreviations used in the code. The 'onset' column describes (in shorthand) what the tendency of correlations for the MHWs is during the onset of events. This is repeated for the 'full' and 'decline' columns respectively. The 'season' column briefly states the most clear/noteworthy pattern(s) when looking at how the correlations are divided up by season. The same is done in the 'region' column. The last column, 'story', gives a TRUE/FALSE if I think the variable has a story to tell. Something worth pursuing further. Particularly to see if the variables relate strongly to other variables, not just temperature. This then could provide a framework for determining 'types' of MHWs (e.g. strong SSS change with strong latent heat flux).")
```

## References