jonahMethods.Rmd

---
title: 'Methods - Jonah Crab'
date-modified: '2023-07-28'
execute: 
  message: false
  warning: false
editor_options: 
  chunk_output_type: console
editor: 
  markdown: 
    wrap: 72
---

# Intro
Our approach here will be to utilize the multiple spatial scales and
different metrics of abundance (weight/biomass, in kg, and catch per
tow) to identify recurring patterns/trends in the dynamics of Jonah crab
in the Gulf of Maine. We start with univariate analysis to assess the
dimensionality, nonlinearity, and predictability of the system, and to
evaluate if these characteristics differ across the coastline or at
different depths. If the univariate analysis indicates that the system
has appropriately low-dimensional, nonlinear dynamics, we move on to
exploring the potential drivers of Jonah crab abundance, specifically
interspecies interactions with Atlantic rock crabs (primarily a
competitor with Jonah crabs) and Atlantic sea scallops (primarily a prey
species). The results for weight generally (but not always) are similar to those for the catch data, so here we focus only on catch.

### Load packages

```{r}
#Load packages
library(tidyverse)
library(lubridate) #date formatting
library(patchwork) #combining plots
library(viridis)
library(rEDM) #EDM
library(sf) #for spatial data
library(tseries)
```

### Import and clean data

```{r}
df_tows<-read.csv("data/Maine_inshore_trawl/MEtows.csv") #tow data
df_s_cat<- read.csv("data/Maine_inshore_trawl/MEscallopCatch.csv") #scallop catch
df_r_cat<- read.csv("data/Maine_inshore_trawl/MErockCatch.csv") #rock crab catch
df_j_cat<- read.csv("data/Maine_inshore_trawl/MEjonahCatch.csv") #jonah crab catch

surveyGrid <-st_read("~/Downloads/lab_notebook/Maine/MaineDMR_-_Inshore_Trawl_Survey_Grid") #CRS: WGS 84/EPSG 4326
```

Rename and add variables to survey grid

```{r}
surveyGrid <- surveyGrid %>% 
  mutate(Region = region_id,
         Stratum = depth_stra,
         GridID = grid_id, .keep="unused", .before=last_surve)

surveyGrid$area <- as.numeric(paste0(surveyGrid$Region, surveyGrid$Stratum))
```

Define function cleanCatch and apply to the three catch data files

```{r}
cleanCatch <- function(x) {
  full_join(x, df_tows) %>%
    arrange(Survey, Tow_Number) %>% 
    select(-c("Stratum", "Subsample_Weight_kg", "Subsample_Weight_kg_2", "Male_Wt_kg", "Female_Wt_kg","Date", "Surface_WaterTemp_DegC", "Surface_Salinity", "End_Latitude","End_Longitude", "Air_Temp", "Tow_Time")) %>%
    mutate(Number_Caught = replace_na(Number_Caught,0),
           Weight_kg = replace_na(Weight_kg,0),
           Expanded_Catch = replace_na(Expanded_Catch,0),
           Expanded_Weight_kg = replace_na(Expanded_Weight_kg,0)) %>% 
    mutate(Stratum = Depth_Stratum, Date = date(ymd_hms(Start_Date)), .keep="unused") %>%
    mutate(area = as.numeric(paste0(Region, Stratum)),.before= Survey)
}

s_cat_clean_seasons <- cleanCatch(df_s_cat) %>% 
mutate(Common_Name = "Scallop")

r_cat_clean_seasons <- cleanCatch(df_r_cat) %>% 
  mutate(Common_Name = "Rock")

j_cat_clean_seasons <- cleanCatch(df_j_cat) %>% 
  mutate(Common_Name = "Jonah")

```

Reorder columns so the data frames are easier to work with

```{r}
#Reorder columns
colOrder<-c("area", "Survey", "Tow_Number", "Region", "Stratum", "Expanded_Catch", 
            "Expanded_Weight_kg", "Date", "Common_Name", "Number_Caught", "Weight_kg",
            "Start_Latitude", "Start_Longitude","Season",
            "Year","Grid", "Start_Depth_fathoms", "End_Depth_fathoms",
            "Bottom_WaterTemp_DegC", "Bottom_Salinity")

j_cat_clean_seasons <- j_cat_clean_seasons %>% select(all_of(colOrder))
r_cat_clean_seasons <- r_cat_clean_seasons %>% select(all_of(colOrder))
s_cat_clean_seasons <- s_cat_clean_seasons %>% select(all_of(colOrder))
```

Define a function, summaryCatch, to summarize over tows for each area
(area = region/stratum combination) and apply to all three species.

```{r}
summaryCatch <- function(df) {
  df %>% group_by(area, Season, Year, Region, Stratum) %>%
    summarise(avgCatch = mean(Expanded_Catch),
              avgWt = mean(Expanded_Weight_kg)) 
}

#computes averages for each study area (area = region-stratum combination)
j_cat_sum_seasons <- summaryCatch(j_cat_clean_seasons)
r_cat_sum_seasons <- summaryCatch(r_cat_clean_seasons)
s_cat_sum_seasons <- summaryCatch(s_cat_clean_seasons)
```

Combine the three separate data frames

```{r}
catch_seasons <- s_cat_sum_seasons %>% left_join(j_cat_sum_seasons, by=c("area", "Season", "Region", "Stratum", "Year"), suffix = c("_s", "_j"))

catch_seasons <- catch_seasons %>% left_join(r_cat_sum_seasons, by=c("area", "Season", "Region", "Stratum", "Year")) %>% 
  mutate(avgCatch_r = avgCatch,avgWt_r = avgWt, .keep="unused")
```

Now make it "tidy" and convert area, Species, Season, Region, and
Stratum to factors.

```{r}
catchTidy_seasons <- pivot_longer(catch_seasons, 
                          cols = 6:ncol(catch_seasons)) %>% 
  mutate(Type = case_when(
    startsWith(name, "avgCatch_") ~"catch",
    startsWith(name,"avgWt_") ~"wt")) %>% 
  mutate(Species = case_when(
    endsWith(name, "s") ~"scallop",
    endsWith(name, "r") ~"rock",
    endsWith(name, "j") ~"jonah"))

catchTidy_seasons <- catchTidy_seasons %>% 
  mutate(area = as.factor(area), Species = as.factor(Species),Season = as.factor(Season),Region = as.factor(Region), Stratum = as.factor(Stratum)) %>% 
  select(-name)
```

The complete function will turn implicit missing values into explicit
ones. In this case, Stratum 4 was not surveyed until 2003, so this will
add in the appropriate NA values for 2000-2002. We do this for both the
tidy version and the original/long version. Turning Fall/Spring and Year
into specific dates makes for easier plotting, although they may be
slightly different than the original dates that were associated with
individual tows.

```{r}
catch_complete <- complete(data=catch_seasons %>% ungroup(), Region, Stratum, Season, Year) %>% 
  mutate(area = paste(Region, Stratum)) %>% 
  mutate(date=paste(Year, case_when(Season== "Fall" ~ "-11-01", Season =="Spring" ~"-05-01"), sep = ""), .before=Region) %>% 
  filter(date != "2000-05-01")

catchTidy_complete<- complete(data = catchTidy_seasons %>% ungroup() %>%  filter(Type=="catch"),  Region, Stratum, Season, Year) %>% 
        mutate(area = as.numeric(paste0(Region, Stratum))) %>% 
        mutate(date=paste(Year, case_when(Season== "Fall" ~ "-11-01", Season =="Spring" ~"-05-01"), sep = ""), .before=Region) %>% filter(date != "2000-05-01")

```

Parse the date column

```{r}
catch_complete <- catch_complete %>% mutate(date = lubridate::ymd(date))
catchTidy_complete <- catchTidy_complete %>% mutate(date = lubridate::ymd(date))
```

# General trends

First, we look at trends in overall abundance:

```{r}
#line graph of abundance over time by season, no spatial distinction
ggplot(data = catchTidy_seasons %>% filter(Type=="catch", Species=="jonah") %>% group_by(Year, Season) %>% summarise(value = mean(value)))+geom_line(aes(x=Year, y=value))+facet_wrap(~Season)+theme_classic()+labs(y="Abundance (catch/tow)")
```

For comparison, we also plot the catch time series of Jonah crabs
alongside rock crabs, both overall and separated by season:

```{r}
ggplot(data = catchTidy_complete %>% filter(Type == "catch", Species != "scallop") %>% group_by(date, Species) %>% 
         summarise(avg = mean(value, na.rm = TRUE)), aes(x=date, y=avg))+geom_line()+facet_wrap(~Species)+labs(y="catch")+theme_classic()
#a lot more variation in the jonah crabs


ggplot(data = catchTidy_complete %>% filter(Type == "catch", Species != "scallop") %>% group_by(Year, Season, Species) %>% 
          summarise(avg = mean(value, na.rm = TRUE)), aes(x=Year, y=avg))+geom_line()+facet_grid(Season~Species)+labs(y="catch")+theme_classic()
#a lot more variation in the jonah crabs
```

Fluctuations in catch of rock crabs appear smaller and more regular than
for Jonah crabs.

Now we'll incorporate the spatial aspect of the data. We start by
looking at abundance of Jonah and rock crabs averaged from 2000-2022.

```{r}
regionsGrid_orig <- surveyGrid %>% group_by(area) %>% summarise(num = n_distinct(GridID))
regionsGrid <- left_join(regionsGrid_orig, catchTidy_complete %>% filter(Type=="catch",) %>% group_by(area, date, Species) %>% summarise(avg = mean(value)))

ggplot(data=regionsGrid %>% filter(Species != "scallop") %>% group_by(Species, area) %>% summarise(avg = mean(avg)))+geom_sf(aes(fill=avg))+facet_wrap(~Species)+scale_fill_viridis_c()

```

We can also separate by seasons:

```{r}
regionsGrid_seasons <- left_join(regionsGrid_orig, catchTidy_complete %>% filter(Type=="catch",) %>% group_by(area, date, Species, Season) %>% summarise(avg = mean(value)))

#Jonah and rock crabs
ggplot(data=regionsGrid_seasons %>% filter(Species != "scallop") %>% group_by(area, Species, Season) %>% summarize(avg = mean(avg, na.rm=TRUE)))+geom_sf(aes(fill=avg))+facet_grid(Season~Species)+scale_fill_viridis_c(option = "F", name="avg catch")

# Jonah crab only
#ggplot(data=regionsGrid_seasons %>% filter(Species == "jonah") %>% group_by(area, Season) %>% summarize(avg = mean(avg)))+geom_sf(aes(fill=avg))+facet_wrap(~Season)+scale_fill_viridis_c()

```

We look at seasonal movement by year:

```{r}
jonahCatch <- left_join(regionsGrid_orig, catchTidy_complete %>% filter(Species=="jonah"))

jonahCatchFall <- jonahCatch %>% filter(Season == "Fall", Type == "catch") %>% rename(valueFall = value) %>% select(-c("Season", "date"))
jonahCatchSpring <- jonahCatch %>% filter(Season == "Spring", Type == "catch") %>% st_drop_geometry() %>% rename(valueSpring = value) %>% select(-c("Season", "date"))

jonahCatchDiff <- left_join(jonahCatchFall, jonahCatchSpring) %>% arrange(area, Year)
jonahCatchDiff <- jonahCatchDiff %>% mutate(diff = valueFall -valueSpring)

#difference fall vs spring by area
ggplot(data = jonahCatchDiff %>% group_by(area) %>% summarise(avg = mean(diff, na.rm=TRUE)))+geom_sf(aes(fill=avg))+
scale_fill_viridis_b(name="fall catch minus spring catch", option="G")

#difference fall vs spring by year and area
ggplot(data = jonahCatchDiff %>% filter(Year != 2000 & Year != 2020))+geom_sf(aes(fill=diff))+
scale_fill_viridis_b(name="fall catch minus spring catch")+facet_wrap(~Year)

seasonalDiff_by_area<- jonahCatchDiff %>% st_drop_geometry() %>% group_by(area) %>% summarise(avg = mean(diff, na.rm=TRUE), sd = sd(diff, na.rm=TRUE))
seasonalDiff_by_area

```

# Second differencing
For EDM, we want the time series to be as long as possible, so instead
of considering seasons separately, we utilize the entire time series but
take second-differences (i.e., x(t) = x(t)-x(t-2)) to remove the
seasonal effect. Unfortunately, because of COVID, we are also going to
get rid of the data after the missed Spring 2020 survey, since we can't
take second-differences with NA values and we want to avoid linear
interpolation that may obscure nonlinear system dynamics.

```{r}
lag2 <- function(x) {
  x_lagged <- (x - lag(x, 2))
  return(x_lagged)
}  # test with lag2(c(1, 3, 3, 5, 6, 9, 12))

catch_complete_diff <- catch_complete %>% arrange(date) %>% group_by(area) %>% 
  mutate(across(where(is.double) & !date, lag2)) %>% 
  arrange(area) %>% 
  filter(date != "2000-11-01" & date != "2001-05-01") %>%  filter(date < as.Date("2020-05-01"))

#Tidy it up
complete_tidy_diff <- pivot_longer(catch_complete_diff,cols = 7:ncol(catch_complete)) %>% 
  mutate(Type = case_when(
    startsWith(name, "avgCatch_") ~"catch",
    startsWith(name,"avgWt_") ~"wt",
    startsWith(name,"avgLogWt") ~"logWt",
    startsWith(name,"avgLogCatch") ~"logCatch")) %>% 
  mutate(Species = case_when(
    endsWith(name, "s") ~"scallop",
    endsWith(name, "r") ~"rock",
    endsWith(name, "j") ~"jonah")) %>%
  mutate(area = as.factor(area), Species = as.factor(Species),
         Region = as.factor(Region), Type = as.factor(Type),
         Stratum = as.factor(Stratum)) %>% 
  select(-name)
```

Let's look at the differenced data:

```{r}
#All areas on one graph, split by species
ggplot(data = complete_tidy_diff %>% 
         filter(Type == "catch", Species != "scallop"), aes(x=date, y=value, color=area))+geom_line()+facet_wrap(~Species) +labs(y="2nd-differenced catch", x="Year")

#Averaged across areas, split by species
ggplot(data = complete_tidy_diff %>% filter(Type == "catch", Species != "scallop") %>% group_by(date, Species) %>% 
    summarise(avg = mean(value, na.rm = TRUE)), aes(x=date, y=avg))+geom_line()+facet_wrap(~Species)+theme_classic()+labs(y="2nd-differenced catch", x="Year")

#Colored by species, split by area
ggplot(data = complete_tidy_diff %>% filter(Type == "catch", Species != "scallop"), aes(x=date, y=value, color=Species))+geom_line()+facet_grid(Region~Stratum)+labs(x="Depth stratum", y="Region")+theme(axis.title.x = element_text(margin = margin(t = 10, r = 0, b = 0, l = 0)))


```

As we noted earlier, there appears to be much more variation in Jonah
crab abundance compared to the rock crabs, although not in every area.

# E and θ (aggregate)

We start our EDM analyses by looking at the dynamics of the system in
aggregate, averaged across all areas. First we look at embedding
dimension using Simplex projection with leave-one-out cross-validation:

```{r}
EmbedDimension(dataFrame=complete_tidy_diff %>% filter(Species=="jonah", Type=="catch") %>% group_by(date) %>% 
                 summarise(avg = mean(value, na.rm = TRUE)) %>% 
                 ungroup() %>% select(date, avg),  columns ="avg", target="avg", lib = "1 37", pred="1 37")


EmbedDimension(dataFrame=complete_tidy_diff %>% filter(Species=="jonah", Type=="wt") %>% group_by(date) %>% 
                 summarise(avg = mean(value, na.rm = TRUE)) %>% 
                 ungroup() %>% select(date, avg),  columns ="avg", target="avg", lib = "1 37", pred="1 37")
```

The optimal embedding dimension is the one with the highest predictive
skill. For both catch and weight, E=2 has the highest rho value,
although neither measure of abundance shows particularly high rho
values.

We will also use S-mapping to evaluate nonlinearity, with E=2:

```{r}
PredictNonlinear(dataFrame=complete_tidy_diff %>% filter(Species=="jonah", Type=="catch") %>% group_by(date) %>% 
                   summarise(avg = mean(value, na.rm = TRUE)) %>% 
                   ungroup() %>% select(date, avg),  columns ="avg", target="avg", lib = "1 37", pred="1 37", E=2)

PredictNonlinear(dataFrame=complete_tidy_diff %>% filter(Species=="jonah", Type=="wt") %>% group_by(date) %>% 
                   summarise(avg = mean(value, na.rm = TRUE)) %>% 
                   ungroup() %>% select(date, avg),  columns ="avg", target="avg", lib = "1 37", pred="1 37", E=2)
```

We can see that both catch is better represented by a
nonlinear model (theta >0) than a linear one. In addition, the
predictive skill of the s-map projections is much higher than the Simplex
models.

Now we want to see how these characteristics change at different spatial
scales. We will first define some useful functions to more efficiently
repeat the previous tests:

```{r}
############ Find E and rho - vector input -------------------------------------------------

findE_v <- function(v, maxE = 7) {
  lib_vec <- paste(1, length(v))
  indices <- c(1:length(v))
  df <- data.frame(indices,v)
  colnames(df)<-c("index", "value")
  rho_E<- EmbedDimension(dataFrame = df, lib = lib_vec, pred = lib_vec, columns = "value",target = "value", maxE = maxE)
  E_out<-rho_E[which.max(rho_E$rho),"E"][1]
  return(E_out)
}

findErho_v <- function(v) {
  lib_vec <- paste(1, length(v))
  indices <- c(1:length(v))
  df <- data.frame(indices,v)
  colnames(df)<-c("index", "value")
  rho_E<- EmbedDimension(dataFrame = df, lib = lib_vec, pred = lib_vec, columns = "value",target = "value", maxE = 7)
  rho_out<-rho_E[which.max(rho_E$rho),"rho"][1]
  return(rho_out)
}

############ Find Theta and rho - vector input -------------------------------------------------

findTheta_v <- function(v, E) {
  lib_vec <- paste(1, length(v))
  indices <- c(1:length(v))
  df <- data.frame(indices,v)
  colnames(df)<-c("index", "value")
  rho_Theta<- PredictNonlinear(dataFrame = df, lib = lib_vec, pred = lib_vec, columns = "value",target = "value", E=E)
  Theta_out<-rho_Theta[which.max(rho_Theta$rho),"Theta"][1]
  return(Theta_out)
}

findThetaRho_v <- function(v, E) {
  lib_vec <- paste(1, length(v))
  indices <- c(1:length(v))
  df <- data.frame(indices,v)
  colnames(df)<-c("index", "value")
  rho_Theta<- PredictNonlinear(dataFrame = df, lib = lib_vec, pred = lib_vec, columns = "value",target = "value", E=E)
  Rho_out<-rho_Theta[which.max(rho_Theta$rho),"rho"][1]
  return(Rho_out)
}

############ findSpeciesE & findSpeciesErho -------------------------------------------------------

#returns a tibble with the optimal embedding dimension for time series from each region/stratum combination
#Find species E
findSpeciesE <- function(df, season=NULL, type) {
  if (is.null(season)) {
    df_out <- df %>% 
      filter(Type == type) %>% 
      group_by(Region, Stratum) %>% 
      select(Year, value) %>%
      summarise(E_opt = findE_v(value)) %>%
      pivot_wider(names_from = Stratum, values_from = E_opt) %>% 
      ungroup() %>% 
      select(-Region)
  }
  else {
  
  df_out <- df %>% 
    filter(Type == type, Season == season) %>% 
    group_by(Region, Stratum) %>% 
    select(Year, value) %>%
    summarise(E_opt = findE_v(value)) %>%
    pivot_wider(names_from = Stratum, values_from = E_opt) %>% 
    ungroup() %>% 
    select(-Region) }
  
  return(df_out)
}

#Find predictive skill for Simplex with optimal E
findSpeciesErho <- function(df, season=NULL, type) {
  if (is.null(season)) {
  df_out <- df %>% 
    filter(Type == type) %>% 
    group_by(Region, Stratum) %>% 
    select(Year, value) %>%
    summarise(E_opt_rho = findErho_v(value)) %>%
    pivot_wider(names_from = Stratum, values_from = E_opt_rho) %>% 
    ungroup() %>% 
    select(-Region) }
  
  else {
    df_out <- df %>% 
      filter(Type == type, Season == season) %>% 
      group_by(Region, Stratum) %>% 
      select(Year, value) %>%
      summarise(E_opt_rho = findErho_v(value)) %>%
      pivot_wider(names_from = Stratum, values_from = E_opt_rho) %>% 
      ungroup() %>% 
      select(-Region)
    
  }
  return(df_out)
}

############ findSpeciesTheta -------------------------------------------------------

findSpeciesTheta <- function(df, season=NULL, type) {
  df_E <- findSpeciesE(df=df, season=season, type=type)
  
  findE <- function (reg, strat) {
     E <- as.integer(df_E %>% slice(reg) %>% pull(strat))
    return(E) 
    }
  
  if (is.null(season)) {
    df_out <- df %>%
      filter(Type == type) %>%
      mutate(E_row = as.integer(Region), E_col = as.integer(Stratum)) %>% 
      rowwise() %>% 
      mutate(E = findE(reg=E_row, strat=E_col)) %>%
      group_by(Region, Stratum) %>%
      summarise(Theta_opt =  findTheta_v(value, E[1])) %>%
      pivot_wider(names_from = Stratum, values_from = Theta_opt) %>%
      ungroup() %>%
      select(-Region) }
  
  else {
      df_out <- df %>%
        filter(Type == type, Season==season) %>%
        mutate(E_row = as.integer(Region), E_col = as.integer(Stratum)) %>% 
        rowwise() %>% 
        mutate(E = findE(reg=E_row, strat=E_col)) %>%
        group_by(Region, Stratum) %>%
        summarise(Theta_opt =  findTheta_v(value, E[1])) %>%
        pivot_wider(names_from = Stratum, values_from = Theta_opt) %>%
        ungroup() %>%
        select(-Region) }
  
  return(df_out)
}
#Find predictive skill for s-map with optimal E and theta
findSpeciesTheta_rho <- function(df, season=NULL, type) {
  df_E <- findSpeciesE(df=df, season=season, type=type)
  
  findE <- function (reg, strat) {
    E <- as.integer(df_E %>% slice(reg) %>% pull(strat))
    return(E) 
  }
  
  if (is.null(season)) {
    df_out <- df %>%
      filter(Type == type) %>%
      mutate(E_row = as.integer(Region), E_col = as.integer(Stratum)) %>% 
      rowwise() %>% 
      mutate(E = findE(reg=E_row, strat=E_col)) %>%
      group_by(Region, Stratum) %>%
      summarise(Theta_opt =  findThetaRho_v(value, E[1])) %>%
      pivot_wider(names_from = Stratum, values_from = Theta_opt) %>%
      ungroup() %>%
      select(-Region) }
  
  else {
    df_out <- df %>%
      filter(Type == type, Season==season) %>%
      mutate(E_row = as.integer(Region), E_col = as.integer(Stratum)) %>% 
      rowwise() %>% 
      mutate(E = findE(reg=E_row, strat=E_col)) %>%
      group_by(Region, Stratum) %>%
      summarise(Theta_opt =  findThetaRho_v(value, E[1])) %>%
      pivot_wider(names_from = Stratum, values_from = Theta_opt) %>%
      ungroup() %>%
      select(-Region) }
  
  return(df_out)
}

```

# E and θ (by area)

Now we apply the functions to the catch data:

```{r, fig.show='hide'}
jonah_catchE <- findSpeciesE(complete_tidy_diff %>% filter(Species=="jonah") %>% na.omit(), type="catch")
jonah_catchE_rho<- findSpeciesErho(complete_tidy_diff %>% filter(Species=="jonah") %>% na.omit(), type="catch")%>% round(digits=3)

jonah_catch_theta<- findSpeciesTheta(complete_tidy_diff %>% filter(Species=="jonah") %>% na.omit(), type="catch")
jonah_catch_theta_rho<- findSpeciesTheta_rho(complete_tidy_diff %>% filter(Species=="jonah") %>% na.omit(), type="catch")
```

And repeat for the weight data:
```{r, fig.show='hide'}
jonah_wtE <- findSpeciesE(complete_tidy_diff %>% filter(Species=="jonah") %>% na.omit(), type="wt") 
jonah_wtE_rho<- findSpeciesErho(complete_tidy_diff %>% filter(Species=="jonah") %>% na.omit(), type="wt")%>% round(digits=3)

jonah_wt_theta<- findSpeciesTheta(complete_tidy_diff %>% filter(Species=="jonah") %>% na.omit(), type="wt")
jonah_wt_theta_rho<- findSpeciesTheta_rho(complete_tidy_diff %>% filter(Species=="jonah") %>% na.omit(), type="wt")

```

Let's take a look:

```{r}
jonah_catchE
jonah_catchE_rho

jonah_catch_theta
jonah_catch_theta_rho

jonah_wtE
jonah_wtE_rho

jonah_wt_theta
jonah_wt_theta_rho
```

There is clear spatial heterogeneity: some areas have an optimal
embedding dimension of only 1, while others reach 6; some areas display
essentially linear dynamics, with theta=0.01, but most areas appear to
be highly nonlinear.

The predictive skill also ranges widely. For example, Region 2, Stratum
4 appears especially low (catch s-map is 0.14, weight Simplex is 0.03,
and catch Simplex and weight s-map have negative rho), whereas Region 1,
Stratum 1 appears to be among the most predictable (all 4 models have
rho \>0.7). We note that the predictability for individual areas is
generally higher than the aggregate results we computed earlier.

The lowered-number regions and strata (further south and closer inshore,
respectively) appear to be generally more predictable than the outer
strata and more northern regions. This seems reminiscent of the
geographical area (in square miles) and number of tows allocated to each
region/stratum combination. Let's test this by first creating data
frames with the numbers from the ME-NH Inshore Trawl Survey protocols
and procedures manual.

```{r}
sq_miles <- data.frame(c(253.27, 279.63, 259.62, 205.3, 138.54), c(214.22, 191.23, 262.9, 206.12, 220.49), c(227.35, 211.66, 280.03, 310.49, 365.04), c(225.65, 263.49, 183.69, 170.72, 196.11))
colnames(sq_miles)<- colnames(jonah_catchE)

tows_per_area <- data.frame(c(6, 7, 6, 5, 4), c(6, 5, 7, 5, 6), c(6, 6, 7, 8, 9), c(5, 5, 4, 4, 4))
colnames(tows_per_area)<- colnames(jonah_catchE)

sq_miles_v <- c(as.matrix(sq_miles))
tows_v <- c(as.matrix(tows_per_area))

jonah_catchE_v <- c(as.matrix(jonah_catchE))
jonah_catchE_rho_v <- c(as.matrix(jonah_catchE_rho))
jonah_catch_theta_v <- c(as.matrix(jonah_catch_theta))
jonah_catch_theta_rho_v <- c(as.matrix(jonah_catch_theta_rho))
jonah_wtE_v <- c(as.matrix(jonah_wtE))
jonah_wtE_rho_v <- c(as.matrix(jonah_wtE_rho))
jonah_wt_theta_v <- c(as.matrix(jonah_wt_theta))
jonah_wt_theta_rho_v <- c(as.matrix(jonah_wt_theta_rho))

corDf <- data.frame(sq_miles_v, tows_v, jonah_catchE_v, jonah_catchE_rho_v, jonah_catch_theta_v, jonah_catch_theta_rho_v, jonah_wtE_v, jonah_wtE_rho_v, jonah_wt_theta_v, jonah_wt_theta_rho_v)
areaList <- c(11, 12, 13, 14, 21, 22, 23, 24, 31, 32, 33, 34, 41, 42, 43, 44, 51, 52, 53, 54)

library(corrplot)
corrplot(cor(corDf), method = "circle", order = 'hclust', type="lower", diag=FALSE)

testRes = cor.mtest(corDf, conf.level = 0.95, method="spearman")

corrplot(cor(corDf), method = "circle", order = 'hclust', type="lower", diag=FALSE, p.mat = testRes$p, sig.level = c(0.001, 0.01, 0.05), pch.cex = 0.9,
         insig = 'label_sig', pch.col = 'grey20')
```

Key takeaways:

1.  Weight E and catch E are positively correlated, as are weight theta
    and catch theta

2.  Areas that were highly predictable using Simplex are more likely to
    be highly predictable using s-mapping

3.  Within an abundance metric (catch or weight), the predictability of
    a system is not significantly correlated with its complexity
    (dimensionality) or linearity

4.  The size of an area is not significantly correlated (using
    Spearman's method) with any of the other metrics besides number of
    tows, which is to be expected because tows were specifically
    allocated to ensure even density in all areas.

5.  Number of tows per area appeared correlated with Simplex rho for
    catch data (positive correlation) and theta for weight data
    (negative correlation), but the strength of these correlations was
    relatively low, and it is likely that this is a case of statistical
    significance without practical/biological significance.

We will also look at these results when we average across regions or
strata.

# E and θ (by region/stratum)

First, we'll define a function that takes in the name of a group and
returns EDM stats on a data frame grouped appropriately:

```{r}
findSpeciesGroups_both<- function(df, species, type, g) {
  df_out <- df %>% na.omit() %>% 
    filter(Type == type, Species == species) %>% 
    group_by(!!sym(g), date) %>% 
    summarise(avg = mean(value)) %>% 
    group_by(!!sym(g))  %>%
    summarise(E_opt = findE_v(avg),
              rho_E = findErho_v(avg),
              Theta = findTheta_v(avg, E_opt),
              rho_theta = findThetaRho_v(avg, E_opt))
  return(df_out)
}
```

Now we apply to the catch and weight data, looking at regional
differences:

```{r, fig.show='hide'}
jonah_catch_regions_stats<- findSpeciesGroups_both(complete_tidy_diff, type="catch", g="Region", species="jonah")
jonah_wt_regions_stats <- findSpeciesGroups_both(complete_tidy_diff, type="wt", g="Region", species="jonah")
jonah_catch_regions_stats
jonah_wt_regions_stats
```

Dimensionality appears relatively consistent between regions, but
predictability and nonlinearity appear to vary widely. As when we looked
at individual areas earlier, the lower-numbered regions (further south)
seem more predictable than the more northern regions when modeled with
the appropriate nonlinear parameter (i.e., s-map with optimal theta). Is
there a negative correlation between region and s-map predictability?

We use Spearman's method rather than the default Pearson's because
Region is an ordinal variable and Pearson's requires at least interval
data.

```{r}
#Catch - s-map predictive skill
cor.test(jonah_catch_regions_stats$rho_theta, as.numeric(jonah_catch_regions_stats$Region), method = "spearman", alternative = "less")

# rho = -0.7, p=0.117

#Catch - nonlinearity
cor.test(jonah_catch_regions_stats$Theta, as.numeric(jonah_catch_regions_stats$Region), method = "spearman", alternative = "less")
# rho = -0.7, p=0.117

#Weight - s-map predictive skill
cor.test(jonah_wt_regions_stats$rho_theta, as.numeric(jonah_wt_regions_stats$Region), method = "spearman", alternative = "less")
# rho = -0.9, p=0.042

#Catch - show corr coefficients
corrplot(cor(data.frame(jonah_catch_regions_stats %>% mutate(Region = as.integer(Region))), method="spearman"), method="number")
# Catch - label significant correlations
corrplot(cor(data.frame(jonah_catch_regions_stats %>% mutate(Region = as.integer(Region)))), p.mat = 
           cor.mtest(jonah_catch_regions_stats %>% mutate(Region = as.integer(Region)), method="spearman", alternative="less")$p, insig = 'label_sig')

#Weight - show corr coefficients
corrplot(cor(data.frame(jonah_wt_regions_stats %>% mutate(Region = as.integer(Region))), method="spearman"), method="number")
# Weight - label significant correlations
corrplot(cor(data.frame(jonah_wt_regions_stats %>% mutate(Region = as.integer(Region)))), p.mat = 
           cor.mtest(jonah_wt_regions_stats %>% mutate(Region = as.integer(Region)), method="spearman", alternative="less")$p, insig = 'label_sig')

#both weight and catch, to increase sample size
corrplot(cor(data.frame(rbind(jonah_wt_regions_stats, jonah_catch_regions_stats) %>% mutate(Region = as.integer(Region))), method="spearman"), method="number")
corrplot(cor(data.frame(rbind(jonah_wt_regions_stats, jonah_catch_regions_stats) %>% mutate(Region = as.integer(Region)))), p.mat = 
           cor.mtest(rbind(jonah_wt_regions_stats, jonah_catch_regions_stats) %>% mutate(Region = as.integer(Region)), method="spearman", alternative="less")$p, insig = 'label_sig')
```

S-map predictability is indeed strongly negatively correlated with
region. Although the combined correlation plot also shows a significant
negative correlation between E and theta, this is only driven by the
catch data and is not present in the weight data. The only clear strong
association present in both sets of data is the negative relationship
between region and nonlinear predictive skill.

On to strata:

```{r, fig.show='hide'}
jonah_catch_strat_stats<- findSpeciesGroups_both(complete_tidy_diff, type="catch", g="Stratum", species="jonah")
jonah_wt_strat_stats <- findSpeciesGroups_both(complete_tidy_diff, type="wt", g="Stratum", species="jonah")
jonah_catch_strat_stats
jonah_wt_strat_stats
```

For each stratum, the values computed from the weight data were very
similar to those computed from the catch data. Simplex projection using
the disparate abundance metrics identified identical embedding dimension
values, rho values generally within \<0.1 of each other, and concurring
measures of nonlinearity: theta between 1-1.5 for strata 1, 3, and 4,
with theta=9 for stratum 2. This provides a reassuring indication that
the method is robust to at least small amounts of observation noise
and/or slightly different ways of "viewing" the attractor manifold of
the system.

After our analysis of the univariate dynamics of the Jonah crab fishery
at multiple spatial scales, we can feel confident that it has the
low-dimensional, nonlinear dynamics that make CCM an appropriate method
of determining causality.

# Convergent cross-mapping

As with the univariate analysis, we will start by averaging across all
spatial scales. Previously, we identified E=2 to be the optimal
embedding dimension for the overall Jonah crab population for both catch
and weight, so we will use that value for cross-mapping.

However, since we will also be attempting to predict the target species
(scallops or rock crabs) from Jonah crab abundance, we will need to
identify the appropriate E and theta for the target populations as well.
The CCM function automatically computes cross-map skill in both
directions, but we will be repeating the function calls twice: once with
E optimized for Jonah and once with E optimized for the target.

First, we'll define a function that adds a new column with the direction
for which we optimized the cross-mapping, and a vector containing the
different combinations we'll be trying.

```{r}
addDirection <- function(df) {
  df_out <- df %>% 
    mutate(xmap = case_when(prey=="scallop" & direction=="jonah -> prey" ~ "jonah -> scallop",
                            prey=="scallop" & direction=="prey -> jonah" ~ "scallop -> jonah",
                            prey=="rock" & direction == "jonah -> prey" ~ "jonah -> rock",
                            prey=="rock" & direction == "prey -> jonah" ~ "rock -> jonah")) %>% 
    select(-c(prey, jonah, direction))
  return(df_out)
}

combos <- c("jonah:scallop", "scallop:jonah", "jonah:rock","rock:jonah")
```

First, filter the second-differenced data frame to create the more
specific data frames needed for CCM analysis:

```{r}
# Catch only
catchCCMdf <- catch_complete_diff %>% ungroup() %>% na.omit() %>% 
  select(date, Region, Stratum, area, avgCatch_s, avgCatch_r, avgCatch_j) %>% 
  rename(rock = avgCatch_r , scallop = avgCatch_s , jonah= avgCatch_j)  %>% 
  mutate(area = as.integer(paste0(Region, Stratum)))

# Weight only
wtCCMdf <- catch_complete_diff %>% ungroup() %>% na.omit() %>% 
  select(date, Region, Stratum, area, avgWt_s, avgWt_r, avgWt_j) %>% 
  rename(rock = avgWt_r , scallop = avgWt_s , jonah= avgWt_j)  %>% 
  mutate(area = as.integer(paste0(Region, Stratum)))

catchCCMdf_agg <- catchCCMdf %>% group_by(date) %>% summarise(across(scallop:jonah, mean))
wtCCMdf_agg <- wtCCMdf %>% group_by(date) %>% summarise(across(scallop:jonah, mean))
```

We also create a data frame to hold the combinations that we want to
test:

```{r}
params_ccm_combos <- data.frame(jonah=c("jonah", "jonah"), prey=c("scallop", "rock"))
```

Now for the actual cross-mapping:
## Aggregate CCM - catch

```{r}
par(mfrow=c(2,2), mar=c(1,1,1,1))

RESULTS_ccm_aggregate <- pmap_dfr(params_ccm_combos,function(jonah,prey){
  lib_vec <- paste(1, nrow(catchCCMdf_agg))
  
  #find optimal E for predicting scallop/rock from jonah (i.e., jonah -> scallop/rock)
  rho_E_1<- EmbedDimension(dataFrame = catchCCMdf_agg, lib = lib_vec, pred = lib_vec, 
                           columns = jonah,target = prey, maxE = 7, showPlot = FALSE) 
  E_out_1<-rho_E_1[which.max(rho_E_1$rho),"E"][1] #store E
  
  #Run CCM - jonah to rock/scallop
  out_1 <- CCM(dataFrame= catchCCMdf_agg, columns=jonah, target=prey, E = E_out_1, Tp=1,
               libSizes = paste(E_out_1+2, nrow(catchCCMdf_agg) - E_out_1, "1",sep=" "), sample=100, verbose=FALSE, showPlot = TRUE) %>%
    mutate(jonah=jonah,
           prey=prey,
           direction= paste("jonah","->","prey"),
           E = E_out_1)
  
  #find optimal E for predicting jonah from scallop or rock (i.e., scallop/rock -> jonah)
  rho_E_2<- EmbedDimension(dataFrame = catchCCMdf_agg, lib = lib_vec, pred = lib_vec, 
                           columns = prey,target = jonah, maxE = 7, showPlot = FALSE)
  E_out_2<-rho_E_2[which.max(rho_E_1$rho),"E"][1] #store E
  
  #Run CCM - rock/scallop to jonah
  out_2 <- CCM(dataFrame= catchCCMdf_agg, columns=prey, target=jonah, E = E_out_2, Tp=1, 
               libSizes = paste(E_out_2+2, nrow(catchCCMdf_agg)-E_out_2, "1",sep=" "), 
               sample=100, verbose=FALSE, showPlot = TRUE) %>%
    mutate(jonah=jonah,
           prey=prey,
           direction= paste("prey","->","jonah"),
           E = E_out_2)
  
  bind_rows(out_1,out_2)
  
}) %>% addDirection()

```

Both directions gave the same system dimensionality. Jonah and scallop
may have bidirectional causality, with faster convergence for
scallop:jonah compared to jonah:scallop. Rock crab abundance could be
effectively cross-mapped from Jonah crabs but not vice versa, suggesting
unidirectional influence from rock crabs to Jonah crabs. We will repeat
with the weight data to see if the results are the same.

## Aggregate CCM - weight
```{r}
par(mfrow=c(2,2), mar=c(2,2,2,2))

RESULTS_ccm_wt_aggregate <- pmap_dfr(params_ccm_combos,function(jonah,prey){
  lib_vec <- paste(1, nrow(wtCCMdf_agg))
  
  #find optimal E for predicting scallop or rock from jonah (jonah -> scallop/rock)
  rho_E_1<- EmbedDimension(dataFrame = wtCCMdf_agg, lib = lib_vec, pred = lib_vec, 
                           columns = jonah,target = prey, maxE = 7, showPlot = FALSE) 
  E_out_1<-rho_E_1[which.max(rho_E_1$rho),"E"][1] #store E
  
  #Run CCM - jonah to rock/scallop
  out_1 <- CCM(dataFrame= wtCCMdf_agg, columns=jonah, target=prey, E = E_out_1, Tp=1,
               libSizes = paste(E_out_1+2, nrow(wtCCMdf_agg) - E_out_1, "1",sep=" "), sample=100, verbose=FALSE, showPlot = TRUE) %>%
    mutate(jonah=jonah,
           prey=prey,
           direction= paste("jonah","->","prey"),
           E = E_out_1)
  
  #find optimal E for predicting jonah from scallop or rock (scallop/rock -> jonah)
  rho_E_2<- EmbedDimension(dataFrame = wtCCMdf_agg, lib = lib_vec, pred = lib_vec, 
                           columns = prey,target = jonah, maxE = 7, showPlot = FALSE)
  E_out_2<-rho_E_2[which.max(rho_E_1$rho),"E"][1] #store E
  
  #Run CCM - rock/scallop to jonah
  out_2 <- CCM(dataFrame= wtCCMdf_agg, columns=prey, target=jonah, E = E_out_2, Tp=1, 
               libSizes = paste(E_out_2+2, nrow(wtCCMdf_agg)-E_out_2, "1",sep=" "), 
               sample=100, verbose=FALSE, showPlot = TRUE) %>%
    mutate(jonah=jonah,
           prey=prey,
           direction= paste("prey","->","jonah"),
           E = E_out_2)
  
  bind_rows(out_1,out_2)
  
}) %>% addDirection()
```

The Jonah and rock crab weight results are very similar to the catch
results: same E and convergence in the direction of Jonah crabs --\>
rock crabs but not vice versa, suggesting that rock crabs have a causal
influence on Jonah crabs. However, cross-mapping with scallops had very
low predictive skill in both directions, with neither direction
exhibiting clear convergence.

We can visually compare all four cross-mappings with the following:

```{r}
#Overall xmap skill - catch
ggplot(data = RESULTS_ccm_aggregate %>% ungroup() %>% group_by(LibSize) %>% summarise(across(all_of(combos), ~mean(.x, na.rm = TRUE))) %>% pivot_longer(cols = 2:5)) +geom_line(aes(x=LibSize, y=value, color=name))

#Overall xmap skill - weight
ggplot(data = RESULTS_ccm_wt_aggregate %>% ungroup() %>% group_by(LibSize) %>% summarise(across(all_of(combos), ~mean(.x, na.rm = TRUE))) %>% pivot_longer(cols = 2:5)) +geom_line(aes(x=LibSize, y=value, color=name))
```

We will now assess cross-map skill at smaller spatial scales, to see if
there are causal relationships in specific areas, strata, or regions
that are obscured when we average across all survey data.

## CCM by area - catch
```{r}
#CATCH

#Create the grid of all combinations we want to test
params_areas_ccm_combos <- expand.grid(jonah=c("jonah"), prey=c("scallop", "rock"), areaInput = areaList,  stringsAsFactors = FALSE) %>% 
    mutate(areaInput = as.integer(areaInput))

# By area, each individually - CATCH
RESULTS_ccm_by_area <- pmap_dfr(params_areas_ccm_combos,function(jonah,prey, areaInput){
  
  df_temp <- catchCCMdf %>% filter(area == areaInput)
  lib_vec <- paste(1, nrow(df_temp))

  rho_E_1<- EmbedDimension(dataFrame = df_temp, lib = lib_vec, pred = lib_vec, 
                           columns = jonah,target = prey, maxE = 7, showPlot = FALSE)
  E_out_1<-rho_E_1[which.max(rho_E_1$rho),"E"][1]
  out_1 <- CCM(dataFrame= df_temp, columns=jonah, target=prey, E = E_out_1, Tp=1,
               libSizes = paste(E_out_1+2, nrow(df_temp) - E_out_1, "1",sep=" "), sample=100, verbose=FALSE, showPlot = FALSE) %>%
    mutate(jonah=jonah,
           prey=prey,
           direction= paste("jonah","->","prey"),
           area = areaInput, 
           E = E_out_1)

  rho_E_2<- EmbedDimension(dataFrame = df_temp, lib = lib_vec, pred = lib_vec, 
                           columns = prey,target = jonah, maxE = 7, showPlot = FALSE)
  E_out_2<-rho_E_2[which.max(rho_E_1$rho),"E"][1]
  out_2 <- CCM(dataFrame= df_temp, columns=prey, target=jonah, E = E_out_2, Tp=1, 
               libSizes = paste(E_out_2+2, nrow(df_temp)-E_out_2, "1",sep=" "), sample=100, verbose=FALSE, showPlot = FALSE)  %>%
    mutate(jonah=jonah,
           prey=prey,
           direction= paste("prey","->","jonah"),
           area = areaInput,
           E = E_out_2)

  bind_rows(out_1,out_2)
  
}) %>% addDirection()

```

Let's try and visualize the results:

```{r}
ccm_col_order <- c("LibSize", "xmap", "area", "E", "jonah:scallop", "scallop:jonah","jonah:rock","rock:jonah")

p1<-ggplot(data = RESULTS_ccm_by_area %>% ungroup() %>% group_by(LibSize, area) %>% summarise(across(all_of(combos), ~mean(.x, na.rm = TRUE))) %>% pivot_longer(cols = 3:6) %>% mutate(Region = substr(area, 1, 1), Stratum = substr(area, 2, 2)) %>% arrange(area))+geom_line(aes(x=LibSize, y=value, color=name))+facet_grid(Region~Stratum)
p1
```

Just from a preliminary visual analysis, we can see clear spatial
heterogeneity in the cross-map results. For example, jonah:scallop shows
much higher predictive skill in 22 and 33 compared to the other areas.
As one might expect from the aggregate CCM, jonah:rock shows the
strongest relationship of the four cross-mappings for most areas, and
areas where jonah:scallop was not as strong (e.g., areas 11 and 42) did
not show strong relationships for any of the combinations we tested.
Stratum (column) 3 is also notable for the close clustering of all four
cross-map relationships around 0, while the four relationships generally
diverge from each other with increasing library size for areas in the
other three strata.

We'll repeat the analysis for the weight data.

## CCM by area - weight
```{r}
#WEIGHT

# By area, each individually - WEIGHT
RESULTS_ccm_wt_by_area <- pmap_dfr(params_areas_ccm_combos,function(jonah,prey, areaInput){
  
  df_temp <- wtCCMdf %>% filter(area == areaInput)
  lib_vec <- paste(1, nrow(df_temp))

  rho_E_1<- EmbedDimension(dataFrame = df_temp, lib = lib_vec, pred = lib_vec, 
                           columns = jonah,target = prey, maxE = 7, showPlot = FALSE)
  E_out_1<-rho_E_1[which.max(rho_E_1$rho),"E"][1]
  out_1 <- CCM(dataFrame= df_temp, columns=jonah, target=prey, E = E_out_1, Tp=1,
               libSizes = paste(E_out_1+2, nrow(df_temp) - E_out_1, "1",sep=" "), sample=100, verbose=FALSE, showPlot = FALSE) %>%
    mutate(jonah=jonah,
           prey=prey,
           direction= paste("jonah","->","prey"),
           area = areaInput, 
           E = E_out_1)

  rho_E_2<- EmbedDimension(dataFrame = df_temp, lib = lib_vec, pred = lib_vec, 
                           columns = prey,target = jonah, maxE = 7, showPlot = FALSE)
  E_out_2<-rho_E_2[which.max(rho_E_1$rho),"E"][1]
  out_2 <- CCM(dataFrame= df_temp, columns=prey, target=jonah, E = E_out_2, Tp=1, 
               libSizes = paste(E_out_2+2, nrow(df_temp)-E_out_2, "1",sep=" "), sample=100, verbose=FALSE, showPlot = FALSE)  %>%
    mutate(jonah=jonah,
           prey=prey,
           direction= paste("prey","->","jonah"),
           area = areaInput,
           E = E_out_2)

  bind_rows(out_1,out_2)
  
}) %>% addDirection()

```

Let's try and visualize the results:

```{r}

p2<-ggplot(data = RESULTS_ccm_wt_by_area %>% ungroup() %>% 
         group_by(LibSize, area) %>% 
         summarise(across(all_of(combos), ~mean(.x, na.rm = TRUE))) %>% 
         pivot_longer(cols = 3:6) %>% mutate(Region = substr(area, 1, 1), Stratum = substr(area, 2, 2))) +
         geom_line(aes(x=LibSize, y=value, color=name))+facet_grid(Region~Stratum)

(p1+ggtitle("catch"))+(p2+ggtitle("weight"))
```


Let's look at regions:

## CCM by region - catch
```{r}
params_regions_combos<- expand.grid(jonah=c("jonah"), prey=c("scallop", "rock"),
                                      region=c(1:5), stringsAsFactors = FALSE)

# By region, each individually - CATCH
catchCCMdf_reg <- catchCCMdf %>% group_by(Region, date) %>% summarise(scallop = mean(scallop), rock = mean(rock), jonah = mean(jonah))

RESULTS_ccm_by_reg <- pmap_dfr(params_regions_combos,function(jonah,prey, region){
  df_temp <- catchCCMdf_reg %>% filter(Region == region)
  lib_vec <- paste(1, nrow(df_temp))
  rho_E_1<- EmbedDimension(dataFrame = df_temp, lib = lib_vec, pred = lib_vec,
                           columns = jonah,target = prey, maxE = 7, showPlot = FALSE)
  E_out_1<-rho_E_1[which.max(rho_E_1$rho),"E"][1]
  out_1 <- CCM(dataFrame= df_temp, columns=jonah, target=prey, E = E_out_1, Tp=1,
               libSizes = paste(E_out_1+2, nrow(df_temp) - E_out_1, "1",sep=" "), sample=100, verbose=FALSE, showPlot = FALSE) %>%
    mutate(jonah=jonah,
           prey=prey,
           direction= paste("jonah","->","prey"),
           region=region,
           E = E_out_1)
  rho_E_2<- EmbedDimension(dataFrame = df_temp, lib = lib_vec, pred = lib_vec,
                           columns = prey,target = jonah, maxE = 7, showPlot = FALSE)
  E_out_2<-rho_E_2[which.max(rho_E_1$rho),"E"][1]
  out_2 <- CCM(dataFrame= df_temp, columns=prey, target=jonah, E = E_out_2, Tp=1,
               libSizes = paste(E_out_2+2, nrow(df_temp)-E_out_2, "1",sep=" "), sample=100, verbose=FALSE, showPlot = FALSE)  %>%
    mutate(jonah=jonah,
           prey=prey,
           direction= paste("prey","->","jonah"),
           region=region,
           E = E_out_2)

  bind_rows(out_1,out_2)
}) %>% addDirection()

p3<-ggplot(data = RESULTS_ccm_by_reg %>% ungroup() %>% 
         group_by(LibSize, region) %>% 
         summarise(across(all_of(combos), ~mean(.x, na.rm = TRUE))) %>% 
         pivot_longer(cols = 3:6)) +
         geom_line(aes(x=LibSize, y=value, color=name))+facet_wrap(~region)

#compare the 5 plots on the right (regions) to the 5 rows on the left
(p1+ggtitle("catch"))+p3
```

Again, we repeat for weight.

## CCM by region - weight
```{r}
# By region, each individually - WEIGHT
wtCCMdf_reg <- wtCCMdf %>% group_by(Region, date) %>% summarise(scallop = mean(scallop), rock = mean(rock), jonah = mean(jonah))

RESULTS_ccm_wt_by_reg <- pmap_dfr(params_regions_combos,function(jonah,prey, region){
  df_temp <- wtCCMdf_reg %>% filter(Region == region)
  lib_vec <- paste(1, nrow(df_temp))
  #jonah -> prey
  rho_E_1<- EmbedDimension(dataFrame = df_temp, lib = lib_vec, pred = lib_vec,
                           columns = jonah,target = prey, maxE = 7, showPlot = FALSE)
  E_out_1<-rho_E_1[which.max(rho_E_1$rho),"E"][1]
  out_1 <- CCM(dataFrame= df_temp, columns=jonah, target=prey, E = E_out_1, Tp=1,
               libSizes = paste(E_out_1+2, nrow(df_temp) - E_out_1, "1",sep=" "), sample=100, verbose=FALSE, showPlot = FALSE) %>%
    mutate(jonah=jonah,
           prey=prey,
           direction= paste("jonah","->","prey"),
           region=region,
           E = E_out_1)
  rho_E_2<- EmbedDimension(dataFrame = df_temp, lib = lib_vec, pred = lib_vec,
                           columns = prey,target = jonah, maxE = 7, showPlot = FALSE)
  E_out_2<-rho_E_2[which.max(rho_E_1$rho),"E"][1]
  out_2 <- CCM(dataFrame= df_temp, columns=prey, target=jonah, E = E_out_2, Tp=1,
               libSizes = paste(E_out_2+2, nrow(df_temp)-E_out_2, "1",sep=" "), sample=100, verbose=FALSE, showPlot = FALSE)  %>%
    mutate(jonah=jonah,
           prey=prey,
           direction= paste("prey","->","jonah"),
           region=region,
           E = E_out_2)

  bind_rows(out_1,out_2)
}) %>% addDirection()

p4<-ggplot(data = RESULTS_ccm_wt_by_reg %>% ungroup() %>% 
         group_by(LibSize, region) %>% 
         summarise(across(all_of(combos), ~mean(.x, na.rm = TRUE))) %>% 
         pivot_longer(cols = 3:6)) +
         geom_line(aes(x=LibSize, y=value, color=name))+facet_wrap(~region)

#compare the 5 plots on the right (regions) to the 5 rows on the left
(p2+ggtitle("wt"))+p4
(p3+ggtitle("catch"))+(p4+ggtitle("wt"))
```

Catch and weight are pretty similar for regions 3, 4, and 5, but
somewhat opposite for regions 1 and 2.

Now on to strata:

## CCM by stratum - catch
```{r}
params_strat_combos<- expand.grid(jonah=c("jonah"), prey=c("scallop", "rock"),
                                      stratum=c(1:4), stringsAsFactors = FALSE)

# By depth stratum, each individually - CATCH
catchCCMdf_strat <- catchCCMdf %>% group_by(Stratum, date) %>% summarise(scallop = mean(scallop), rock = mean(rock), jonah = mean(jonah))

RESULTS_ccm_by_strat <- pmap_dfr(params_strat_combos,function(jonah,prey, stratum){
  df_temp <- catchCCMdf_strat %>% filter(Stratum == stratum)
  lib_vec <- paste(1, nrow(df_temp))
  #jonah -> prey
  rho_E_1<- EmbedDimension(dataFrame = df_temp, lib = lib_vec, pred = lib_vec,
                           columns = jonah,target = prey, maxE = 7, showPlot = FALSE)
  E_out_1<-rho_E_1[which.max(rho_E_1$rho),"E"][1]
  out_1 <- CCM(dataFrame= df_temp, columns=jonah, target=prey, E = E_out_1, Tp=1,
               libSizes = paste(E_out_1+2, nrow(df_temp) - E_out_1, "1",sep=" "), sample=100, verbose=FALSE, showPlot = FALSE) %>%
    mutate(jonah=jonah,
           prey=prey,
           direction= paste("jonah","->","prey"),
           stratum=stratum,
           E = E_out_1)
  rho_E_2<- EmbedDimension(dataFrame = df_temp, lib = lib_vec, pred = lib_vec,
                           columns = prey,target = jonah, maxE = 7, showPlot = FALSE)
  E_out_2<-rho_E_2[which.max(rho_E_1$rho),"E"][1]
  out_2 <- CCM(dataFrame= df_temp, columns=prey, target=jonah, E = E_out_2, Tp=1,
               libSizes = paste(E_out_2+2, nrow(df_temp)-E_out_2, "1",sep=" "), sample=100, verbose=FALSE, showPlot = FALSE)  %>%
    mutate(jonah=jonah,
           prey=prey,
           direction= paste("prey","->","jonah"),
           stratum=stratum,
           E = E_out_2)

  bind_rows(out_1,out_2)
}) %>% addDirection()

p5<-ggplot(data = RESULTS_ccm_by_strat %>% ungroup() %>% 
         group_by(LibSize, stratum) %>% 
         summarise(across(all_of(combos), ~mean(.x, na.rm = TRUE))) %>% 
         pivot_longer(cols = 3:6)) +
         geom_line(aes(x=LibSize, y=value, color=name))+facet_wrap(~stratum)
p5
```

And for weight:

## CCM by stratum - weight
```{r}
# By depth stratum, each individually - WEIGHT
wtCCMdf_strat <- wtCCMdf %>% group_by(Stratum, date) %>% summarise(scallop = mean(scallop), rock = mean(rock), jonah = mean(jonah))

RESULTS_ccm_wt_by_strat <- pmap_dfr(params_strat_combos,function(jonah,prey, stratum){
  df_temp <- wtCCMdf_strat %>% filter(Stratum == stratum)
  lib_vec <- paste(1, nrow(df_temp))
  #jonah -> prey
  rho_E_1<- EmbedDimension(dataFrame = df_temp, lib = lib_vec, pred = lib_vec,
                           columns = jonah,target = prey, maxE = 7, showPlot = FALSE)
  E_out_1<-rho_E_1[which.max(rho_E_1$rho),"E"][1]
  out_1 <- CCM(dataFrame= df_temp, columns=jonah, target=prey, E = E_out_1, Tp=1,
               libSizes = paste(E_out_1+2, nrow(df_temp) - E_out_1, "1",sep=" "), sample=100, verbose=FALSE, showPlot = FALSE) %>%
    mutate(jonah=jonah,
           prey=prey,
           direction= paste("jonah","->","prey"),
           stratum=stratum,
           E = E_out_1)
  rho_E_2<- EmbedDimension(dataFrame = df_temp, lib = lib_vec, pred = lib_vec,
                           columns = prey,target = jonah, maxE = 7, showPlot = FALSE)
  E_out_2<-rho_E_2[which.max(rho_E_1$rho),"E"][1]
  out_2 <- CCM(dataFrame= df_temp, columns=prey, target=jonah, E = E_out_2, Tp=1,
               libSizes = paste(E_out_2+2, nrow(df_temp)-E_out_2, "1",sep=" "), sample=100, verbose=FALSE, showPlot = FALSE)  %>%
    mutate(jonah=jonah,
           prey=prey,
           direction= paste("prey","->","jonah"),
           stratum=stratum,
           E = E_out_2)

  bind_rows(out_1,out_2)
}) %>% addDirection()

p6<-ggplot(data = RESULTS_ccm_wt_by_strat %>% ungroup() %>% 
         group_by(LibSize, stratum) %>% 
         summarise(across(all_of(combos), ~mean(.x, na.rm = TRUE))) %>% 
         pivot_longer(cols = 3:6)) +
         geom_line(aes(x=LibSize, y=value, color=name))+facet_wrap(~stratum)
p6
```

At this point it is helpful to review where we've been and where we plan
to go. We've done cross-mapping at four different scales:

1.  aggregate (average everything) = run CCM on 1 37-point time series

2.  by area (no averaging except for tows in an area) = run CCM on 20
    different \~37-point time series INDIVIDUALLY

3.  by region (average across all areas in a region; i.e., across the
    different depth strata at a given latitude) to get 1 time series for
    each of the 5 regions = run CCM on 5 different \~37-point time
    series INDIVIDUALLY)

4.  by stratum (average across all areas in a depth stratum; i.e.,
    across the different latitude groupings at a given depth) to get 1
    time series for each of the 4 depth strata = run CCM on 4 different
    \~37-point time series INDIVIDUALLY)

We will now see which relationships demonstrate significant CCM, which is considered to occur
if: (1) Cross-map prediction skill at max lib size is statistically
significant based on a randomized null distribution: i.e., cross-map
skill is greater than 95% of cross-map values obtained from 100 random
rearrangements of the Jonah crab data while holding the order of the
other species constant, destroying any dynamic relationship between the
two time series. (2) The time series demonstrates convergence, meaning
prediction skill increases with increasing library size (e.g., more
information about the other species makes Jonah crab predictions better
because they are causally related), quantified as a delta rho (rho at
max lib size minus rho at min lib size) greater than zero.

Finally, we will make our investigation of inter-area dynamics more
comprehensive by using a geographic permutation approach where we
compute cross-map skill for a given area as the mean cross-map skill for
100 different "extended" areas created by adding a random sample of
neighboring tows to the tows within the area, so we get slightly
different results when averaging across tows for each area. This
approach has two benefits: first, it creates distributions of rho values
so we can perform statistical comparisons between areas more easily than
comparing single values. Second, averaging across more tows for each
area decreases the likelihood that all tows in an area will have caught
none of the target species. As having many zeros tends to decrease the
performance of Simplex models, this spatial extension approach may
improve model performance simply by lowering the number of zeros in a
given time series. We will also use the significance-testing functions
we develop to test significance for each of these 100 extensions, giving
us a number (out of 100) of significant permutations as an additional
method of comparison between areas.

# Significance testing

## Convergence
The first criteria for determining significant CCM is convergence, which we determine to be present when predictive skill at maximum library size is larger than predictive skill at minimum library size.
First we define a helper function:

```{r}
find_delta = function(x) {
  delta_rho = last(x)-first(x)
 return(delta_rho)
}
```

Let's revisit our aggregate CCM results:
```{r}
RESULTS_ccm_aggregate

conv_agg <- RESULTS_ccm_aggregate %>% group_by(xmap) %>% summarize(across(all_of(combos),find_delta))
conv_agg %>% pivot_longer(cols=2:5) %>% na.omit()

RESULTS_ccm_wt_aggregate %>% group_by(xmap) %>% summarize(across(all_of(combos),find_delta)) %>% pivot_longer(cols=2:5) %>% na.omit()
```

Okay, so we can see that the only direction that does not converge is rock:jonah, which decreases in predictive skill with increasing library size regardless of if the CCM is using the optimal E for rock or for Jonah crabs. Now let's revisit our CCM-by-area results

```{r}
min_max_lib <- RESULTS_ccm_by_area %>% group_by(area, xmap) %>% 
  summarise(across(.cols=all_of(combos), find_delta))

converges<- pivot_longer(min_max_lib, cols=3:6) %>% na.omit() %>% 
 # mutate(optE = substr(xmap, 1, 3)==substr(name, 1, 3)) %>% 
  group_by(area, name) %>% summarise(delta_rho = mean(value)) %>% filter(delta_rho >0)

table(converges$area, converges$name)
colSums(table(converges$area, converges$name))

#WEIGHT
min_max_lib_wt <- RESULTS_ccm_wt_by_area %>% group_by(area, xmap) %>% 
  summarise(across(.cols=all_of(combos), find_delta))

converges_wt<- pivot_longer(min_max_lib_wt, cols=3:6) %>% na.omit() %>% 
 # mutate(optE = substr(xmap, 1, 3)==substr(name, 1, 3)) %>% 
  group_by(area, name) %>% summarise(delta_rho = mean(value)) %>% filter(delta_rho >0)

table(converges_wt$area, converges_wt$name)
colSums(table(converges_wt$area, converges_wt$name))
```

We can see that jonah:rock and jonah:scallop converge for at least 18/20 areas. We will revisit the spatial distribution of the less-frequently converging relationships later. We will also add these results to our convergence table, assigning "TRUE" if a relationship converged for over half of the 20 areas.

Now let's revisit our region and stratum results. We will assign "TRUE" in our table if a relationship converges for 3 or more of the 5 regions or 4 strata.

```{r}
RESULTS_ccm_by_reg %>% group_by(region, xmap) %>% 
  summarise(across(.cols=all_of(combos), find_delta)) %>% 
  pivot_longer(cols=3:6) %>% na.omit() %>% 
  group_by(region, name) %>% summarise(delta_rho = mean(value)) #%>% filter(delta_rho >0)
```

We can see that jonah:rock and jonah:scallop converged for all 5 regions, while cross-maps in the opposite directions did not converge for any regions.

```{r}
RESULTS_ccm_wt_by_reg %>% group_by(region, xmap) %>% 
  summarise(across(.cols=all_of(combos), find_delta)) %>% 
  pivot_longer(cols=3:6) %>% na.omit() %>% 
  group_by(region, name) %>% summarise(delta_rho = mean(value)) %>% filter(delta_rho >0) %>% group_by(name) %>% summarise(num = length(name))
```

Again, only jonah:rock and jonah:scallop converge for over half of the regions.

We move on to strata:

```{r}
#CATCH
RESULTS_ccm_by_strat %>% group_by(stratum, xmap) %>% 
  summarise(across(.cols=all_of(combos), find_delta)) %>% 
  pivot_longer(cols=3:6) %>% na.omit() %>% 
  group_by(stratum, name) %>% summarise(delta_rho = mean(value)) %>% filter(delta_rho >0) %>% group_by(name) %>% summarise(num = length(name))

#WEIGHT
RESULTS_ccm_wt_by_strat %>% group_by(stratum, xmap) %>% 
  summarise(across(.cols=all_of(combos), find_delta)) %>% 
  pivot_longer(cols=3:6) %>% na.omit() %>% 
  group_by(stratum, name) %>% summarise(delta_rho = mean(value)) %>% filter(delta_rho >0) %>% group_by(name) %>% summarise(num = length(name))
```

For catch (unsurprisingly), jonah:rock and jonah:scallop converge for all four regions. Cross-mapping from rock:jonah converges for 3, while scallop:jonah converges for 2. However, for weight, only jonah:rock converges for more than 2/4 strata.

Let's create a new data frame to summarize the convergence results.

```{r}
#AGGREGATE
conv_res <- data.frame(xmap = combos, method="agg", metric="catch", conv=TRUE)
conv_res <- rbind(conv_res, conv_res %>% mutate(metric="wt"))
conv_res[conv_res$xmap=="rock:jonah","conv"]<- FALSE

#BY AREA - catch
conv_res <- rbind(conv_res, data.frame(xmap = combos, method="area", metric="catch", conv=TRUE))
conv_res[c(10, 12),"conv"]<- FALSE

#BY AREA - wt
conv_res <- rbind(conv_res, data.frame(xmap = combos, method="area", metric="wt", conv=TRUE))
conv_res[14,"conv"]<- FALSE

#BY REG - catch
conv_res <- rbind(conv_res, data.frame(xmap = combos, method="reg", metric="catch", conv=TRUE))
conv_res <- conv_res %>% mutate(conv = ifelse(
  metric =="catch" & method =="reg" & 
    (xmap %in% c("scallop:jonah", "rock:jonah")), FALSE, conv))

#BY REG - wt
conv_res <- rbind(conv_res, data.frame(xmap = combos, method="reg", metric="wt", conv=TRUE))
conv_res <- conv_res %>% mutate(conv = ifelse(
  metric =="wt" & method =="reg" & 
    (xmap %in% c("scallop:jonah", "rock:jonah")), FALSE, conv))

#BY STRAT - catch
conv_res <- rbind(conv_res, data.frame(xmap = combos, method="strat", metric="catch", conv=TRUE))
conv_res <- conv_res %>% mutate(conv = ifelse(
  metric =="catch" & method =="strat" & 
    (xmap %in% c("scallop:jonah")), FALSE, conv))

#BY STRAT - wt
conv_res <- rbind(conv_res, data.frame(xmap = combos, method="strat", metric="wt", conv=FALSE))
conv_res <- conv_res %>% mutate(conv = ifelse(
  metric =="wt" & method =="strat" & 
    (xmap %in% c("jonah:rock")), TRUE, conv))

conv_res

ggplot(data=conv_res %>% filter(conv==TRUE))+geom_bar(aes(x=xmap,fill=metric), position="dodge")+labs(x="Cross-map direction", y="Count")+scale_fill_manual(values=c("#D98880", "black"))+theme(axis.text=element_text(size=10),
        axis.title.x=element_text(size=12, margin = margin(t=10)), axis.title.y = element_text(margin = margin(r=10)))
#+ggtitle("Num CCM methods with \U0394\U03C1 >0 for over half of trials")

```