Markdown Document Overview

This markdown doc is set up for reviewing each analysis step I’ve done in parallel for each of the following three sections of the continuous plankton recorder dataset:

All available years (1961 - 2017)
The early period (1961 - 2003)
The Most recent and warmest period (2003 - 2017)

The first sections will mirror the steps taken in Pershing et al., 2005 in setting up the PCA and determining the weights and timelines for each time-period as shown in the following figure.

Data and Packages

####  Packages  ####
library(cluster)
library(factoextra)
library(ggbiplot)
library(tidyverse)
library(here)



####  Functions  ####
source(here::here("R", "cpr_helper_funs.R"))

#Set ggplot theme
theme_set(theme_classic())

####  Load Data  ####

# CPR species Data with the corrresponding SST and lagged SST
species_periods_long <- read_csv(str_c(cpr_boxpath, "data", "processed_data", "cpr_with_SSTlags.csv", sep = "/"))  %>% 
  mutate(period = case_when(
    period == "annual" ~ "Annual",
    period == "jf" ~ "January - February",
    period == "ma" ~ "March - April",
    period == "mj" ~ "May - June",
    period == "ja" ~ "July - August",
    period == "so" ~ "September - October",
    period == "nd" ~ "November - December"),
period = factor(period, 
                levels = c("Annual", "January - February", "March - April", "May - June", "July - August", "September - October", "November - December")))

Analysis Code

#Split into analysis units
full_ts <- species_periods_long
pre_split <- species_periods_long %>% filter(year <= 2003)
post_split <- species_periods_long %>% filter(year > 2003)

#Put the full data sets into a nested list
analysis_sets <- list("Full Timeseries" = list(data = full_ts),
                      "1961 - 2003" = list(data = pre_split),
                      "2004 - 2017" = list(data = post_split))

#Function for extracting %variance explained from PCA object
pull_deviance <- function(pca_sdev) {
  
  eigs <- pca_sdev ^ 2
  
  deviance_df <- rbind(
    SD = sqrt(eigs),
    Proportion = eigs/sum(eigs),
    Cumulative = cumsum(eigs)/sum(eigs))
  
  pca_dev_out <- data.frame(
    "PC1" = str_c(as.character(round(deviance_df[2,1] * 100, 2)), "% of Variance"),
    "PC2" = str_c(as.character(round(deviance_df[2,2] * 100, 2)), "% of Variance"))
  
  return(pca_dev_out)
  
}

# Perform all analyses on this list for consistency and for organization

####  Analysis Loop  ####
for (i in 1:length(analysis_sets)) {
  
  #Keep years with the data going into distance matrices for PCA and clustering
  analysis_sets[[i]]$pca_data <- analysis_sets[[i]]$data %>% 
    filter(period == "Annual") %>% 
    pivot_wider(names_from = species, values_from = anomaly) %>% 
    mutate(decade = factor(floor_decade(year))) %>% 
    select(year, decade, calanus:para_pseudocalanus) %>% 
    drop_na()
  
  #Pull out just the values used for PCA and clustering
  analysis_sets[[i]]$anom_mat <- analysis_sets[[i]]$pca_data %>% select(-year, -decade)
  
  
  #PCA Object
  analysis_sets[[i]]$pca_obj <- prcomp(analysis_sets[[i]]$anom_mat, center = F, scale. = F)
  
  #PCA Leading Modes
  analysis_sets[[i]]$leading_modes <- rownames_to_column(as.data.frame(analysis_sets[[i]]$pca_obj$rotation)) %>%
    dplyr::select(species = rowname, PC1, PC2)
  
  #Percent Deviance Explained (To slide into plots)
  analysis_sets[[i]]$deviance_explained <- pull_deviance(analysis_sets[[i]]$pca_obj$sdev)
  
  #Figure 2a
  analysis_sets[[i]]$fig2a <- analysis_sets[[i]]$leading_modes %>%
    gather(key = "PC", value =  "Principal Component Weight", PC1, PC2) %>%
    mutate(species = factor(species,
                            levels = c("calanus", "centropages", "oithona","para_pseudocalanus",
                                       "metridia", "calanus1to4", "euphausiacea")),
           PC = if_else(PC == "PC1", 
                       as.character(analysis_sets[[i]]$deviance_explained$PC1),
                       as.character(analysis_sets[[i]]$deviance_explained$PC2))
           ) %>%
    ggplot(aes(species, `Principal Component Weight` * -1, fill = PC)) +
      geom_col(position  = "dodge") +
      scale_fill_gmri(palette = "mixed") +
      labs(x = "") +
      theme(#legend.position = c(0.85, 0.85),
            axis.text.x = element_text(angle = 45, hjust = 1))
  
}

PCA Analysis

PCA Weights

Full TS

analysis_sets$`Full Timeseries`$fig2a

Early Period

analysis_sets$`1961 - 2003`$fig2a

Last 14 Years

analysis_sets$`2004 - 2017`$fig2a

PCA Timelines

Taking the PCA weights from the 1961-2003 PCA (mimicking the 2005 paper) we can recreate the original timeline and extend those weights out through 2017 to create the following figures:

Original Time-period

Extended Through 2017

Applying 1961-2005 Weights to Bi-Monthly Periods

When applying 1961-2005 PCA weights to the bi-monthly period anomalies and extend them out through 2017 you get this figure.

Abundance Anomaly Timelines

The following plots show the timelines of abundance anomalies for each of the seven taxa. There are timelines for each bi-monthly period and each quarterly period for each taxa.

Bi-monthly Periods

Calanus V+

Calanus I - IV

Centropages

Euphausiacea

Metridia

Oithona

Para Pseudocalanus

Quarterly Periods

Calanus V+

Calanus I - IV

Centropages

Euphausiacea

Metridia

Oithona

Para Pseudocalanus

Cluster Analysis

for (i in 1:length(analysis_sets)) {
  analysis_sets[[i]]$distance <- get_dist(analysis_sets[[i]]$pca_data)
}

Optimal Clusters

Full TS

Elbow Method

fviz_nbclust(analysis_sets$`Full Timeseries`$anom_mat, kmeans, method = "wss")

Average Silhouette

fviz_nbclust(analysis_sets$`Full Timeseries`$anom_mat, kmeans, method = "silhouette")

Gap Statistic

gap_stat <- clusGap(analysis_sets$`Full Timeseries`$anom_mat, FUN = kmeans, nstart = 25,
                    K.max = 10, B = 50)

fviz_gap_stat(gap_stat)

Early Period

Elbow Method

fviz_nbclust(analysis_sets$`1961 - 2003`$anom_mat, kmeans, method = "wss")

Average Silhouette

fviz_nbclust(analysis_sets$`1961 - 2003`$anom_mat, kmeans, method = "silhouette")

Gap Statistic

gap_stat <- clusGap(analysis_sets$`1961 - 2003`$anom_mat, FUN = kmeans, nstart = 25,
                    K.max = 10, B = 50)

fviz_gap_stat(gap_stat)

Last 14 Years

Elbow Method

fviz_nbclust(analysis_sets$`2004 - 2017`$anom_mat, kmeans, method = "wss")

Average Silhouette

fviz_nbclust(analysis_sets$`2004 - 2017`$anom_mat, kmeans, method = "silhouette")

Gap Statistic

gap_stat <- clusGap(analysis_sets$`2004 - 2017`$anom_mat, FUN = kmeans, nstart = 25,
                    K.max = 10, B = 50)

fviz_gap_stat(gap_stat)

Inter-Taxa Relationships Over Time

Using the Optimal number of clusters identified for each period, we can plot which years cluster together on a population anomaly plot with an axis for each taxa. The units of each axis is the deviance from annual mean abundance. To aid in visuallization a grid has been overlayed which splits the plot into four quadrants. Years falling in the Mutual + quadrant had positive abundance anomalies for both species, with years in the Mutual - quadrant exhibiting lower than average abundance for both species.

K-means Code

#Using 2 Clusters for All Periods
for (i in 1:length(analysis_sets)) {
  
  #Get the cluster identification for each period
  analysis_sets[[i]]$kmeans <- kmeans(analysis_sets[[i]]$anom_mat, centers = 2, nstart = 25)
  
  #Add the cluster ID to the dataset that has year information
  analysis_sets[[i]]$pca_data <- analysis_sets[[i]]$pca_data %>% 
  mutate(cluster = factor(analysis_sets[[i]]$kmeans$cluster))
  
  
  
  ####__####
  ####  Inter-taxa relationships  #####
  ####__####
  
  #Centropages and Calanus
  analysis_sets[[i]]$centro_x_cal <- analysis_sets[[i]]$pca_data %>% 
  ggplot(aes(calanus, centropages, color = cluster)) +
    geom_hline(yintercept = 0, color = "gray", linetype = 2, alpha = 0.5) +
    geom_vline(xintercept = 0, color = "gray", linetype = 2, alpha = 0.5) +
    geom_text(aes(label = year)) +
    ggforce::geom_mark_ellipse() +
    labs(x = "Calanus V+", y = "Centropages") +
    xlim(-2, 2) +
    ylim(-2, 2)
  
  #Table of Quadrat Tally
  analysis_sets[[i]]$centro_x_cal_t <- analysis_sets[[i]]$pca_data %>% 
    mutate(
      `Plot Quadrant` = if_else(calanus > 0 & centropages > 0, "Mutual +", "Mutual -"),
      `Plot Quadrant` = if_else(calanus > 0 & centropages < 0, "Calanus Favored", `Plot Quadrant`),
      `Plot Quadrant` = if_else(calanus < 0 & centropages > 0, "Competitor Favored", `Plot Quadrant`),
      `Plot Quadrant` = factor(`Plot Quadrant`,
                               levels = c("Mutual +", "Calanus Favored", "Competitor Favored", "Mutual -"))) %>%
    group_by(decade) %>%
    count(`Plot Quadrant`) %>%
    tidyr::complete(decade, `Plot Quadrant`, fill = list(n = 0)) %>%
    ggplot(aes(decade, n, fill = `Plot Quadrant`)) +
      geom_col(position = "fill") +
      labs(x = NULL, y = "Proportion") + 
      scale_fill_gmri(palette = "mixed") +
      theme(legend.title = element_blank())
  
  #Oithona and Calanus
  analysis_sets[[i]]$oithona_x_cal <- analysis_sets[[i]]$pca_data %>% 
  ggplot(aes(calanus, oithona, color = cluster)) +
    geom_hline(yintercept = 0, color = "gray", linetype = 2, alpha = 0.5) +
    geom_vline(xintercept = 0, color = "gray", linetype = 2, alpha = 0.5) +
    geom_text(aes(label = year)) +
    ggforce::geom_mark_ellipse() +
    labs(x = "Calanus V+", y = "Oithona") +
    xlim(-2, 2) +
    ylim(-2, 2)
  
  #Table of Quadrat Tally
  analysis_sets[[i]]$oithona_x_cal_t <- analysis_sets[[i]]$pca_data %>% 
    mutate(
      `Plot Quadrant` = if_else(calanus > 0 & oithona > 0, "Mutual +", "Mutual -"),
      `Plot Quadrant` = if_else(calanus > 0 & oithona < 0, "Calanus Favored", `Plot Quadrant`),
      `Plot Quadrant` = if_else(calanus < 0 & oithona > 0, "Competitor Favored", `Plot Quadrant`),
      `Plot Quadrant` = factor(`Plot Quadrant`,
                               levels = c("Mutual +", "Calanus Favored", "Competitor Favored", "Mutual -"))) %>%
    group_by(decade) %>%
    count(`Plot Quadrant`) %>%
    tidyr::complete(decade, `Plot Quadrant`, fill = list(n = 0)) %>%
    ggplot(aes(decade, n, fill = `Plot Quadrant`)) +
      geom_col(position = "fill") +
      labs(x = NULL, y = "Proportion") + 
      scale_fill_gmri(palette = "mixed") +
      theme(legend.title = element_blank())
  
  #Pseudocalanus and Calanus
  analysis_sets[[i]]$pseudo_x_cal <- analysis_sets[[i]]$pca_data %>% 
  ggplot(aes(calanus, para_pseudocalanus, color = cluster)) +
    geom_hline(yintercept = 0, color = "gray", linetype = 2, alpha = 0.5) +
    geom_vline(xintercept = 0, color = "gray", linetype = 2, alpha = 0.5) +
    geom_text(aes(label = year)) +
    ggforce::geom_mark_ellipse() +
    labs(x = "Calanus V+", y = "Pseudocalanus") +
    xlim(-2, 2) +
    ylim(-2, 2)
  
  #Does it match the corellations?
  analysis_sets[[1]]$pca_data %>% 
  ggplot(aes(calanus, para_pseudocalanus)) +
    geom_hline(yintercept = 0, color = "gray", linetype = 2, alpha = 0.5) +
    geom_vline(xintercept = 0, color = "gray", linetype = 2, alpha = 0.5) +
    geom_text(aes(color = cluster, label = year)) +
    geom_smooth(method = "lm")+
    labs(x = "Calanus V+", y = "Pseudocalanus") +
    xlim(-2, 2) +
    ylim(-2, 2)
  
  #Table of Quadrat Tally
  analysis_sets[[i]]$pseudo_x_cal_t <- analysis_sets[[i]]$pca_data %>% 
    mutate(
      `Plot Quadrant` = if_else(calanus > 0 & para_pseudocalanus > 0, "Mutual +", "Mutual -"),
      `Plot Quadrant` = if_else(calanus > 0 & para_pseudocalanus < 0, "Calanus Favored", `Plot Quadrant`),
      `Plot Quadrant` = if_else(calanus < 0 & para_pseudocalanus > 0, "Competitor Favored", `Plot Quadrant`),
      `Plot Quadrant` = factor(`Plot Quadrant`,
                               levels = c("Mutual +", "Calanus Favored", "Competitor Favored", "Mutual -"))) %>%
    group_by(decade) %>%
    count(`Plot Quadrant`) %>%
    tidyr::complete(decade, `Plot Quadrant`, fill = list(n = 0)) %>%
    ggplot(aes(decade, n, fill = `Plot Quadrant`)) +
      geom_col(position = "fill") +
      labs(x = NULL, y = "Proportion") + 
      scale_fill_gmri(palette = "mixed") +
      theme(legend.title = element_blank())
  
  #Calanus I -IV and Calanus
  analysis_sets[[i]]$cal2_x_cal <- analysis_sets[[i]]$pca_data %>% 
  ggplot(aes(calanus, calanus1to4, color = cluster)) +
    geom_hline(yintercept = 0, color = "gray", linetype = 2, alpha = 0.5) +
    geom_vline(xintercept = 0, color = "gray", linetype = 2, alpha = 0.5) +
    geom_text(aes(label = year)) +
    ggforce::geom_mark_ellipse() + 
    labs(x = "Calanus V+", y = "Calanus I - IV") +
    xlim(-2, 2) +
    ylim(-2, 2)
  
  analysis_sets[[i]]$cal2_x_cal_t <- analysis_sets[[i]]$pca_data %>% 
    mutate(
      `Plot Quadrant` = if_else(calanus > 0 & calanus1to4 > 0, "Mutual +", "Mutual -"),
      `Plot Quadrant` = if_else(calanus > 0 & calanus1to4 < 0, "Calanus Favored", `Plot Quadrant`),
      `Plot Quadrant` = if_else(calanus < 0 & calanus1to4 > 0, "Competitor Favored", `Plot Quadrant`),
      `Plot Quadrant` = factor(`Plot Quadrant`,
                               levels = c("Mutual +", "Calanus Favored", "Competitor Favored", "Mutual -"))) %>%
    group_by(decade) %>%
    count(`Plot Quadrant`) %>%
    tidyr::complete(decade, `Plot Quadrant`, fill = list(n = 0)) %>%
    ggplot(aes(decade, n, fill = `Plot Quadrant`)) +
      geom_col(position = "fill") +
      labs(x = NULL, y = "Proportion") + 
      scale_fill_gmri(palette = "mixed") +
      theme(legend.title = element_blank())
  
  
}

Full TS

Calanus x Centropages

analysis_sets$`Full Timeseries`$centro_x_cal

analysis_sets$`Full Timeseries`$centro_x_cal_t

Calanus x Oithona

analysis_sets$`Full Timeseries`$oithona_x_cal

analysis_sets$`Full Timeseries`$oithona_x_cal_t

Calanus x Pseudocalanus

analysis_sets$`Full Timeseries`$pseudo_x_cal

analysis_sets$`Full Timeseries`$pseudo_x_cal_t

Calanus x Calanus I-IV

analysis_sets$`Full Timeseries`$cal2_x_cal

analysis_sets$`Full Timeseries`$cal2_x_cal_t

Early Period

Calanus x Centropages

analysis_sets$`1961 - 2003`$centro_x_cal

analysis_sets$`1961 - 2003`$centro_x_cal_t

Calanus x Oithona

analysis_sets$`1961 - 2003`$oithona_x_cal

analysis_sets$`1961 - 2003`$oithona_x_cal_t

Calanus x Pseudocalanus

analysis_sets$`1961 - 2003`$pseudo_x_cal

analysis_sets$`1961 - 2003`$pseudo_x_cal_t

Calanus x Calanus I-IV

analysis_sets$`1961 - 2003`$cal2_x_cal

analysis_sets$`1961 - 2003`$cal2_x_cal_t

Last 14 Years

Calanus x Centropages

analysis_sets$`2004 - 2017`$centro_x_cal

analysis_sets$`2004 - 2017`$centro_x_cal_t

Calanus x Oithona

analysis_sets$`2004 - 2017`$oithona_x_cal

analysis_sets$`2004 - 2017`$oithona_x_cal_t

Calanus x Pseudocalanus

analysis_sets$`2004 - 2017`$pseudo_x_cal

analysis_sets$`2004 - 2017`$pseudo_x_cal_t

Calanus x Calanus I-IV

analysis_sets$`2004 - 2017`$cal2_x_cal

analysis_sets$`2004 - 2017`$cal2_x_cal_t

Comparing PCA modes and k-means cluster

Original Period, PCA Mode Timeline

Kmeans Cluster Year Assignments

cluster_1 <- analysis_sets$`1961 - 2003`$pca_data %>% 
  filter(cluster == 1) %>% 
  select(year) %>%   
  mutate(`Cluster 1` = "+")

cluster_2 <- analysis_sets$`1961 - 2003`$pca_data %>% 
  filter(cluster == 2) %>% 
  select(year) %>% 
  mutate(`Cluster 2` = "+")

#Plot them as a timeline
analysis_sets$`1961 - 2003`$pca_data %>% 
  select(year) %>% 
  left_join(cluster_1, by = "year") %>% 
  left_join(cluster_2, by = "year") %>% 
  mutate(`Cluster 1` = if_else(is.na(`Cluster 1`) == T, "-", `Cluster 1`),
         `Cluster 2` = if_else(is.na(`Cluster 2`) == T, "-", `Cluster 2`)) %>% 
  gather(key = "Cluster ID", value = "Selected +/-", `Cluster 1`, `Cluster 2`) %>% 
  ggplot(aes(x = year, 
             y = fct_rev(`Cluster ID`), 
             label = `Selected +/-`,
             color = `Selected +/-`)
         ) +
    geom_text(size = 5, show.legend = FALSE) +
    scale_color_gmri(palette = "mixed", reverse = T) +
    labs(x = NULL, y = NULL) + theme_minimal()

Quarterly SST Regressions

Data Setup for Quarterly Analysis

#Target species
species_05 <- c("calanus", "centropages", "oithona", "para_pseudocalanus",
                "metridia", "calanus1to4", "euphausiacea")

#Set up a list for the quarterly groups
quarterly_long <- read_csv(str_c(cpr_boxpath,"data", "processed_data", "cpr_allspecies_long_quarters.csv", sep = "/")) %>% 
  mutate(
    period = case_when(
      period == "annual" ~"Annual",
      period == "q1" ~"Q1",
      period == "q2" ~"Q2",
      period == "q3" ~"Q3",
      period == "q4" ~"Q4",
      TRUE ~ "Missed One"
    )
  ) %>% 
  filter(species %in% species_05)

#Split into analysis units
full_ts_q    <- quarterly_long
pre_split_q  <- quarterly_long %>% filter(year <= 2003)
post_split_q <- quarterly_long %>% filter(year > 2003)

#Put the full data sets into a nested list
quarterly_sets <- list("Full Timeseries" = list(Annual = list(data = filter(full_ts_q, period == "Annual")),
                                                Q1     = list(data = filter(full_ts_q, period == "Q1")),
                                                Q2     = list(data = filter(full_ts_q, period == "Q2")),
                                                Q3     = list(data = filter(full_ts_q, period == "Q3")),
                                                Q4     = list(data = filter(full_ts_q, period == "Q4"))),
                       "1961 - 2003"     = list(Annual = list(data = filter(pre_split_q, period == "Annual")),
                                                Q1     = list(data = filter(pre_split_q, period == "Q1")),
                                                Q2     = list(data = filter(pre_split_q, period == "Q2")),
                                                Q3     = list(data = filter(pre_split_q, period == "Q3")),
                                                Q4     = list(data = filter(pre_split_q, period == "Q4"))),
                       "2004 - 2017"     = list(Annual = list(data = filter(post_split_q, period == "Annual")),
                                                Q1     = list(data = filter(post_split_q, period == "Q1")),
                                                Q2     = list(data = filter(post_split_q, period == "Q2")),
                                                Q3     = list(data = filter(post_split_q, period == "Q3")),
                                                Q4     = list(data = filter(post_split_q, period == "Q4")))
)

Analysis Loop for Quarterly Time Periods

####  Analysis Loop  ####
# i for period
# q for quarter
for (i in 1:length(quarterly_sets)) {
  for (q in 1:5) {
  
  #Keep years with the data going into distance matrices for PCA and clustering
  quarterly_sets[[i]][[q]]$pca_data <- quarterly_sets[[i]][[q]]$data %>% 
    pivot_wider(names_from = species, values_from = anomaly) %>% 
    mutate(decade = factor(floor_decade(year))) %>% 
    dplyr::select(year, decade, calanus:para_pseudocalanus) %>% 
    drop_na()
  
  #Pull out just the values used for PCA and clustering
  quarterly_sets[[i]][[q]]$anom_mat <- quarterly_sets[[i]][[q]]$pca_data %>% select(-year, -decade)
  
  
  #PCA Object
  quarterly_sets[[i]][[q]]$pca_obj <- prcomp(quarterly_sets[[i]][[q]]$anom_mat, center = F, scale. = F)

  #PCA Leading Modes
  quarterly_sets[[i]][[q]]$leading_modes <- rownames_to_column(as.data.frame(quarterly_sets[[i]][[q]]$pca_obj$rotation)) %>%
    dplyr::select(species = rowname, PC1, PC2)

  #Percent Deviance Explained (To slide into plots)
  quarterly_sets[[i]][[q]]$deviance_explained <- pull_deviance(quarterly_sets[[i]][[q]]$pca_obj$sdev)

  #Figure 2a
  quarterly_sets[[i]][[q]]$fig2a <- quarterly_sets[[i]][[q]]$leading_modes %>%
    gather(key = "PC", value =  "Principal Component Weight", PC1, PC2) %>%
    mutate(species = factor(species,
                            levels = c("calanus", "centropages", "oithona","para_pseudocalanus",
                                       "metridia", "calanus1to4", "euphausiacea")),
           PC = if_else(PC == "PC1",
                       as.character(quarterly_sets[[i]][[q]]$deviance_explained$PC1),
                       as.character(quarterly_sets[[i]][[q]]$deviance_explained$PC2))
           ) %>%
    ggplot(aes(species, `Principal Component Weight` * -1, fill = PC)) +
      geom_col(position  = "dodge") +
      scale_fill_gmri(palette = "mixed") +
      labs(x = "") +
      theme(#legend.position = c(0.85, 0.85),
            axis.text.x = element_text(angle = 45, hjust = 1))
  }
  
}


col <- colorRampPalette(c("#BB4444", "#EE9988", "#FFFFFF", "#77AADD", "#4477AA"))

Quarterly Abundance and SST Regressions

The following plots show the results of simple linear regressions between plankton population anomalies and the corresponding SST anomalies for each quarterly period (left) and the same regressions using temperature anomalies from the previous quarter (right). For plankton anomalies in Q1, temperature anomalies from Q4 of the previous year are used.

1981-2017

Calanus V +

Calanus I - IV

Centropages

Euphausiacea

Metridia

Oithona

Para Pseudocalannus

2004 - 2017

Calanus V +

Calanus I - IV

Centropages

Euphausiacea

Metridia

Oithona

Para Pseudocalannus

Quarterly PCA Weights

Full TS

Annual

quarterly_sets$`Full Timeseries`$Annual$fig2a + labs(caption = "Annual")

quarterly_sets$`Full Timeseries`$Annual$anom_mat %>% cor %>% 
  corrplot::corrplot(method="color", 
                     col=col(200),  
                     diag=TRUE, 
                     type = "lower", 
                     addCoef.col = "black", # Add coefficient of correlation
                     mar=c(0,0,2,2) # http://stackoverflow.com/a/14754408/54964
)

Q1

quarterly_sets$`Full Timeseries`$Q1$fig2a + labs(caption = "Quarter 1")

quarterly_sets$`Full Timeseries`$Q1$anom_mat %>% cor %>% 
  corrplot::corrplot(method="color", 
                     col=col(200),  
                     diag=TRUE, 
                     type = "lower", 
                     addCoef.col = "black", # Add coefficient of correlation
                     mar=c(0,0,2,2) # http://stackoverflow.com/a/14754408/54964
)

Q2

quarterly_sets$`Full Timeseries`$Q2$fig2a + labs(caption = "Quarter 2")

quarterly_sets$`Full Timeseries`$Q2$anom_mat %>% cor %>% 
  corrplot::corrplot(method="color", 
                     col=col(200),  
                     diag=TRUE, 
                     type = "lower", 
                     addCoef.col = "black", # Add coefficient of correlation
                     mar=c(0,0,2,2) # http://stackoverflow.com/a/14754408/54964
)

Q3

quarterly_sets$`Full Timeseries`$Q3$fig2a + labs(caption = "Quarter 3")

quarterly_sets$`Full Timeseries`$Q3$anom_mat %>% cor %>% 
  corrplot::corrplot(method="color", 
                     col=col(200),  
                     diag=TRUE, 
                     type = "lower", 
                     addCoef.col = "black", # Add coefficient of correlation
                     mar=c(0,0,2,2) # http://stackoverflow.com/a/14754408/54964
)

Q4

quarterly_sets$`Full Timeseries`$Q4$fig2a + labs(caption = "Quarter 4")

quarterly_sets$`Full Timeseries`$Q4$anom_mat %>% cor %>% 
  corrplot::corrplot(method="color", 
                     col=col(200),  
                     diag=TRUE, 
                     type = "lower", 
                     addCoef.col = "black", # Add coefficient of correlation
                     mar=c(0,0,2,2) # http://stackoverflow.com/a/14754408/54964
)

Early Period

Annual

quarterly_sets$`1961 - 2003`$Annual$fig2a + labs(caption = "Annual")

quarterly_sets$`1961 - 2003`$Annual$anom_mat %>% cor %>% 
  corrplot::corrplot(method="color", 
                     col=col(200),  
                     diag=TRUE, 
                     type = "lower", 
                     addCoef.col = "black", # Add coefficient of correlation
                     mar=c(0,0,2,2) # http://stackoverflow.com/a/14754408/54964
)

Q1

quarterly_sets$`1961 - 2003`$Q1$fig2a + labs(caption = "Quarter 1")

quarterly_sets$`1961 - 2003`$Q1$anom_mat %>% cor %>% 
  corrplot::corrplot(method="color", 
                     col=col(200),  
                     diag=TRUE, 
                     type = "lower", 
                     addCoef.col = "black", # Add coefficient of correlation
                     mar=c(0,0,2,2) # http://stackoverflow.com/a/14754408/54964
)

Q2

quarterly_sets$`1961 - 2003`$Q2$fig2a + labs(caption = "Quarter 2")

quarterly_sets$`1961 - 2003`$Q2$anom_mat %>% cor %>% 
  corrplot::corrplot(method="color", 
                     col=col(200),  
                     diag=TRUE, 
                     type = "lower", 
                     addCoef.col = "black", # Add coefficient of correlation
                     mar=c(0,0,2,2) # http://stackoverflow.com/a/14754408/54964
)

Q3

quarterly_sets$`1961 - 2003`$Q3$fig2a + labs(caption = "Quarter 3")

quarterly_sets$`1961 - 2003`$Q3$anom_mat %>% cor %>% 
  corrplot::corrplot(method="color", 
                     col=col(200),  
                     diag=TRUE, 
                     type = "lower", 
                     addCoef.col = "black", # Add coefficient of correlation
                     mar=c(0,0,2,2) # http://stackoverflow.com/a/14754408/54964
)

Q4

quarterly_sets$`1961 - 2003`$Q4$fig2a + labs(caption = "Quarter 4")

quarterly_sets$`1961 - 2003`$Q4$anom_mat %>% cor %>% 
  corrplot::corrplot(method="color", 
                     col=col(200),  
                     diag=TRUE, 
                     type = "lower", 
                     addCoef.col = "black", # Add coefficient of correlation
                     mar=c(0,0,2,2) # http://stackoverflow.com/a/14754408/54964
)

Last 14 Years

Annual

quarterly_sets$`2004 - 2017`$Annual$fig2a + labs(caption = "Annual")

quarterly_sets$`2004 - 2017`$Annual$anom_mat %>% cor %>% 
  corrplot::corrplot(method="color", 
                     col=col(200),  
                     diag=TRUE, 
                     type = "lower", 
                     addCoef.col = "black", # Add coefficient of correlation
                     mar=c(0,0,2,2) # http://stackoverflow.com/a/14754408/54964
)

Q1

quarterly_sets$`2004 - 2017`$Q1$fig2a + labs(caption = "Quarter 1")

quarterly_sets$`2004 - 2017`$Q1$anom_mat %>% cor %>% 
  corrplot::corrplot(method="color", 
                     col=col(200),  
                     diag=TRUE, 
                     type = "lower", 
                     addCoef.col = "black", # Add coefficient of correlation
                     mar=c(0,0,2,2) # http://stackoverflow.com/a/14754408/54964
)

Q2

quarterly_sets$`2004 - 2017`$Q2$fig2a + labs(caption = "Quarter 2")

quarterly_sets$`2004 - 2017`$Q2$anom_mat %>% cor %>% 
  corrplot::corrplot(method="color", 
                     col=col(200),  
                     diag=TRUE, 
                     type = "lower", 
                     addCoef.col = "black", # Add coefficient of correlation
                     mar=c(0,0,2,2) # http://stackoverflow.com/a/14754408/54964
)

Q3

quarterly_sets$`2004 - 2017`$Q3$fig2a + labs(caption = "Quarter 3")

quarterly_sets$`2004 - 2017`$Q3$anom_mat %>% cor %>% 
  corrplot::corrplot(method="color", 
                     col=col(200),  
                     diag=TRUE, 
                     type = "lower", 
                     addCoef.col = "black", # Add coefficient of correlation
                     mar=c(0,0,2,2) # http://stackoverflow.com/a/14754408/54964
)

Q4

quarterly_sets$`2004 - 2017`$Q4$fig2a + labs(caption = "Quarter 4")

quarterly_sets$`2004 - 2017`$Q4$anom_mat %>% cor %>% 
  corrplot::corrplot(method="color", 
                     col=col(200),  
                     diag=TRUE, 
                     type = "lower", 
                     addCoef.col = "black", # Add coefficient of correlation
                     mar=c(0,0,2,2) # http://stackoverflow.com/a/14754408/54964
)

CPR Clustering Analysis

Markdown Document Overview

PCA Analysis

PCA Weights

Full TS

Early Period

Last 14 Years

PCA Timelines

Original Time-period

Extended Through 2017

Applying 1961-2005 Weights to Bi-Monthly Periods

Abundance Anomaly Timelines

Bi-monthly Periods

Calanus V+

Calanus I - IV

Centropages

Euphausiacea

Metridia

Oithona

Para Pseudocalanus

Quarterly Periods

Calanus V+

Calanus I - IV

Centropages

Euphausiacea

Metridia

Oithona

Para Pseudocalanus

Cluster Analysis

Optimal Clusters

Full TS

Elbow Method

Average Silhouette

Gap Statistic

Early Period

Elbow Method

Average Silhouette

Gap Statistic

Last 14 Years

Elbow Method

Average Silhouette

Gap Statistic

Inter-Taxa Relationships Over Time

Full TS

Calanus x Centropages

Calanus x Oithona

Calanus x Pseudocalanus

Calanus x Calanus I-IV

Early Period

Calanus x Centropages

Calanus x Oithona

Calanus x Pseudocalanus

Calanus x Calanus I-IV

Last 14 Years

Calanus x Centropages

Calanus x Oithona

Calanus x Pseudocalanus

Calanus x Calanus I-IV

Comparing PCA modes and k-means cluster

Original Period, PCA Mode Timeline

Kmeans Cluster Year Assignments

Quarterly SST Regressions

1981-2017

Calanus V +

Calanus I - IV

Centropages

Euphausiacea

Metridia

Oithona

Para Pseudocalannus

2004 - 2017

Calanus V +

Calanus I - IV

Centropages

Euphausiacea

Metridia

Oithona

Para Pseudocalannus

Quarterly PCA Weights

Full TS