diff --git a/covid19/covid_PCA.Rmd b/covid19/covid_PCA.Rmd
new file mode 100644
index 0000000..a58c0d8
--- /dev/null
+++ b/covid19/covid_PCA.Rmd
@@ -0,0 +1,1421 @@
+---
+title: "Covid_fin_EWS"
+output: html_document
+---
+
+```{r setup}
+library(tidyverse)
+
+# NOTE: In these new ones, vars with > 10 missings are dropped, and missing rows
+# omitted completely
+
+until_date <- "2022-01-18"
+
+country_in_question <- "Netherlands"
+
+# ### FULL DATA
+# figure_identifier <- paste0(country_in_question, 
+                            # " Full data until ", 
+                            # until_date)
+# mattidata_orig <- readr::read_csv(paste0(
+#   "../shield-complexity/data/",
+#   country_in_question, "_worldsurvey_nonmissing_since_2021-06-08_to_",
+#   until_date, ".csv"))
+
+# ### FULL DATA OMITTING SYMPTOMS
+# figure_identifier <- paste0(country_in_question,
+#                             " Full data omitting symptoms until ", 
+#                             until_date)
+# 
+# mattidata_orig <- readr::read_csv(
+#   paste0("../shield-complexity/data/",
+#          country_in_question,
+#          "_worldsurvey_nonmissing_since_2021-06-08_to_",
+#          until_date, ".csv")) %>%
+#   dplyr::select(-contains("symp"),
+#                 -contains("covid"),
+#                 -contains("cli"),
+#                 -contains("_anos"))
+
+### COEFFICIENT OF VARIATION DATA
+figure_identifier <- paste0(country_in_question,
+                            " Coefficient of variation until ", 
+                            until_date)
+mattidata_orig <- readr::read_csv(
+  paste0("../shield-complexity/data/",
+  country_in_question,
+  "_worldsurvey_nonmissing_c_of_v_since_2021-06-08_to_",
+  until_date, ".csv"))
+
+# ### COEFFICIENT OF VARIATION DATA OMITTING SYMPTOMS
+# figure_identifier <-
+#   paste0(country_in_question,
+#          " Coefficient of variation, omitting symptoms until ",
+#          until_date)
+# mattidata_orig <- readr::read_csv(
+#   paste0("../shield-complexity/data/",
+#          country_in_question,
+#          "_worldsurvey_nonmissing_c_of_v_since_2021-06-08_to_",
+#          until_date, ".csv")) %>%
+#   dplyr::select(-contains("symp"),
+#                 -contains("covid"),
+#                 -contains("cli"),
+#                 -contains("_anos"))
+
+# ### AAFT SURROTAGE COEFFICIENT OF VARIATION DATA OMITTING SYMPTOMS
+# figure_identifier <-
+  # paste0(country_in_question,
+  #        " AAFT coefficient of variation, omitting symptoms until ",
+  #        until_date)
+# mattidata_orig_presurrogate <- readr::read_csv(
+#   paste0("../shield-complexity/data/",
+#          country_in_question,
+#          "_worldsurvey_nonmissing_c_of_v_since_2021-06-08_to_",
+#          until_date, ".csv")) %>%
+#   dplyr::select(-contains("symp"),
+#                 -contains("covid"),
+#                 -contains("cli"),
+#                 -contains("_anos"))
+# set.seed(10)
+# mattidata_orig <- mattidata_orig_presurrogate %>%
+#   dplyr::mutate(across(-date, ~tseries::surrogate(.,
+#                                                   ns = 1,
+#                                                   fft = TRUE,
+#                                                   amplitude = TRUE)))
+
+# ### AAFT SURROGATE DATA OMITTING SYMPTOMS
+# figure_identifier <- paste0(country_in_question,
+                            # " AAFT surrogates omitting symptoms until ",
+#                             until_date)
+# mattidata_orig_presurrogate <- readr::read_csv(
+#   paste0("../shield-complexity/data/",
+#          country_in_question,
+#          "_worldsurvey_nonmissing_since_2021-06-08_to_",
+#          until_date, ".csv")) %>%
+#   dplyr::select(-contains("symp"),
+#                 -contains("covid"),
+#                 -contains("cli"),
+#                 -contains("_anos"))
+# set.seed(10)
+# mattidata_orig <- mattidata_orig_presurrogate %>%
+#   dplyr::mutate(across(-date, ~tseries::surrogate(.,
+#                                                   ns = 1,
+#                                                   fft = TRUE,
+#                                                   amplitude = TRUE)))
+
+# ### RANDOM SURROGATE DATA OMITTING SYMPTOMS
+# figure_identifier <- paste0(country_in_question,
+                            # " Random surrogates omitting symptoms until ",
+#                             until_date)
+# mattidata_orig_presurrogate <- readr::read_csv(
+#   paste0("../shield-complexity/data/",
+#          country_in_question,
+#          "_worldsurvey_nonmissing_since_2021-06-08_to_",
+#          until_date, ".csv")) %>%
+#   dplyr::select(-contains("symp"),
+#                 -contains("covid"),
+#                 -contains("cli"),
+#                 -contains("_anos"))
+# set.seed(10)
+# mattidata_orig <- mattidata_orig_presurrogate %>%
+#   # Shuffle each column
+#   dplyr::mutate(across(-date, ~sample(., replace = FALSE)))
+
+
+# ### OLD: EPIDEMIC INDICATORS FOR FINLAND
+# epidemic_indicators <-  readr::read_csv("../shield-complexity/data/solanpaa_epidemic_indicators.csv") %>% 
+#   dplyr::filter(date >= min(mattidata_orig$date, na.rm = TRUE) &
+#                   date <= max(mattidata_orig$date, na.rm = TRUE)) %>% 
+#   dplyr::select(date, new_cases, Rt)
+
+dir.create(paste0("ews-output/", figure_identifier))
+```
+
+# Analysis on shuffled data
+
+```{r}
+number_of_surrogates <- 19
+
+set.seed(10)
+mattidata_withdate_surrogates <- mattidata_orig %>% 
+  dplyr::mutate(id = "real_data") %>% 
+  dplyr::group_by(id) %>% 
+  tidyr::nest() %>% 
+  # Duplicate real data; it will be shuffled later
+  slice(rep(1:n(), each = number_of_surrogates + 1)) %>% 
+  dplyr::group_by(id) %>% 
+  # Rename duplicate rows to contain "surrogate_xx"
+  dplyr::mutate(id = dplyr::case_when(
+    dplyr::row_number() != 1 ~ paste0("surrogate_data_",
+                                      dplyr::row_number()-1),
+    TRUE ~ id)) %>% 
+  # For the rows where id contains "surrogate", randomise columns 
+  dplyr::mutate(data = dplyr::case_when(
+    stringr::str_detect(string = id,
+                        pattern = "surrogate") ~ 
+      purrr::map(.x = data,
+                 .f = ~.x %>% 
+                   # # Shuffle each column except date; KEEPING daily configuration
+                   dplyr::slice(sample(1:n())) %>% 
+                   dplyr::mutate(date = sort(date))),
+                   # # Shuffle each column except date; every day RANDOM values
+                   # dplyr::mutate(across(-date, ~sample(., replace = FALSE)))),
+    # The non-surrogate stays intact
+    TRUE ~ data))
+
+max_varpercs_store <- list()
+max_eigvals_store <- list()
+min_varpercs_store <- list()
+all_varpercs_store <- list()
+list_of_varpercs <- list()
+
+for(analysis_number in 1:(number_of_surrogates + 1)){ 
+  
+  omitted_from_pca <- c("date")
+  
+  mattidata_withdate_shuffle <- 
+    mattidata_withdate_surrogates$data[[analysis_number]]
+  
+  mattidata_shuffle <- mattidata_withdate_shuffle %>% 
+    dplyr::select(-all_of(omitted_from_pca))
+  
+  window_lengths <- seq(from = 7, to = 28, by = 7)
+  
+  loadings_list <- list()
+  eigval_list <- list()
+  varperc_list <- list()
+  mattidata_eigval <- list()
+  max_eigvals <- list()
+  max_varpercs <- list()
+  min_varpercs <- list()
+  
+  for(k in 1:length(window_lengths)){
+    
+    eig.val <- NA
+    varperc <- NA
+    loadings <- 
+      matrix(NA, 
+             nrow = 1, 
+             ncol = ncol(mattidata_withdate_shuffle %>% 
+                           dplyr::select(-all_of(omitted_from_pca)))) %>% 
+      data.frame()
+    
+    colnames(loadings) <- names(mattidata_shuffle)
+    
+    win <- window_lengths[[k]]
+    
+    for(i in (win:nrow(mattidata_withdate_shuffle))){
+      res.pca <- prcomp(
+        mattidata_withdate_shuffle %>% 
+          # For e.g. window of 28, keep observations 1-28, 2-29, 3-30 etc.
+          dplyr::filter(row_number() >= i-win+1 &
+                          row_number() <= i) %>% 
+          dplyr::select(-all_of(omitted_from_pca)),
+        scale = TRUE)
+      # Pull loadings of variables on the 1st principal component
+      loadings[(i-win+1), ] <- res.pca$rotation[ , 1]
+      # Pull eigenvalue of the 1st principal component
+      eig.val[i-win+1] <- factoextra::get_eigenvalue(res.pca)$eigenvalue[1] %>% 
+        sum()
+      # Pull % variance explained by first dimension
+      varperc[i-win+1] <- factoextra::get_eigenvalue(res.pca)$variance.percent[1] %>%
+        sum()
+      
+      loadings_list[[k]] <- loadings
+      eigval_list[[k]] <-  eig.val 
+      varperc_list[[k]] <- varperc
+    }
+    
+    mattidata_eigval[[k]] <- dplyr::bind_cols(mattidata_withdate_shuffle, 
+                                              eigval = c(rep(NA, win - 1), 
+                                                         eig.val),
+                                              varperc = c(rep(NA, win - 1), 
+                                                          varperc))
+    
+    max_eigvals[[k]] <- 
+      mattidata_eigval[[k]] %>% 
+      dplyr::filter(eigval == max(eigval, na.rm = TRUE)) %>% 
+      dplyr::pull(eigval)
+    
+    max_varpercs[[k]] <- 
+      mattidata_eigval[[k]] %>%
+      dplyr::filter(varperc == max(varperc, na.rm = TRUE)) %>% 
+      dplyr::pull(varperc)
+    
+    min_varpercs[[k]] <- 
+      mattidata_eigval[[k]] %>% 
+      dplyr::filter(varperc == min(varperc, na.rm = TRUE)) %>% 
+      dplyr::pull(varperc)
+  }
+  
+  
+  list_of_varpercs[[analysis_number]] <- purrr::map(.x = mattidata_eigval,
+                                 .f = ~.x %>% dplyr::select(date, varperc))
+  
+  max_eigvals_store[[analysis_number]] <-
+    purrr::set_names(max_eigvals, paste0("win_", window_lengths)) %>% 
+    tibble::as_tibble()
+  
+  min_varpercs_store[[analysis_number]] <- 
+    purrr::set_names(min_varpercs, paste0("win_", window_lengths)) %>% 
+    tibble::as_tibble()
+  
+  max_varpercs_store[[analysis_number]] <- 
+    purrr::set_names(max_varpercs, paste0("win_", window_lengths)) %>% 
+    tibble::as_tibble()
+  
+  all_varpercs_store[[analysis_number]] <- 
+    purrr::set_names(list_of_varpercs[[analysis_number]], paste0("win_", window_lengths))
+}
+
+### Plot result
+
+real_percentile <- dplyr::bind_rows(max_varpercs_store) %>% 
+  dplyr::mutate(rownum = dplyr::row_number()) %>% 
+  tidyr::pivot_longer(cols = -rownum) %>% 
+  dplyr::group_by(name) %>% 
+  dplyr::mutate(percentile = 1 - dplyr::percent_rank(value)) %>% 
+  dplyr::filter(rownum == 1) %>% 
+  dplyr::mutate(percentile = round(percentile, digits = 2))
+
+real_varperc <- dplyr::bind_rows(max_varpercs_store) %>% 
+  slice(1) %>% 
+  tidyr::pivot_longer(cols = everything()) %>% 
+  dplyr::mutate(
+    name = factor(name, 
+                  levels = paste0("win_", window_lengths),
+                  labels = paste0("Window length: ", window_lengths,
+                                  " (p = ", real_percentile$percentile, ")")))
+
+randomised_varperc <- dplyr::bind_rows(max_varpercs_store) %>% 
+  slice(-1) %>% 
+  tidyr::pivot_longer(cols = everything()) %>% 
+  dplyr::mutate(
+    name = factor(name, 
+                  levels = paste0("win_", window_lengths),
+                  labels = paste0("Window length: ", window_lengths,
+                                  " (p = ", real_percentile$percentile, ")")))
+
+randomised_varpercs <- all_varpercs_store[-1] %>% # Omit the real data
+  dplyr::bind_rows() %>% 
+  purrr::map(as.vector)
+
+randomised_varpercs_daily <- randomised_varpercs %>% 
+  purrr::map(
+    .f = ~.x %>% 
+      na.omit() %>% 
+      dplyr::group_by(date) %>% 
+      dplyr::summarise(min_varperc_for_day = min(varperc, na.rm = TRUE),
+                       max_varperc_for_day = max(varperc, na.rm = TRUE)))
+
+# varperc_minmax <- dplyr::bind_rows(max_varpercs_store) %>% 
+#   slice(-1) %>% 
+#   tidyr::pivot_longer(cols = everything()) %>% 
+#   dplyr::mutate(
+#     name = factor(name, 
+#                   levels = paste0("win_", window_lengths),
+#                   labels = paste0("Window length: ", window_lengths,
+#                                   " (p = ", real_percentile$percentile, ")")))
+
+# Displaying rectangle around real data doesn't work well unless this is done:
+source("statbin2.R")
+
+randomised_varperc %>% 
+  ggplot(aes(x = value,
+             fill = name)) +
+  geom_histogram(binwidth = 1) +
+  scale_y_continuous(breaks = seq(from = 0, to = 1000, by = 5)) +
+  scale_x_continuous(breaks = seq(from = 0, to = 1000, by = 5)) +
+  scale_fill_viridis_d(end = 0.8, option = "inferno") +
+  # geom_vline(data = real_varperc,
+  #            aes(xintercept = value),
+  #            linetype = "dashed") +
+  geom_histogram(data = real_varperc,
+                 stat = StatBin2,
+                 binwidth = 1,
+                 colour = "red",
+                 linetype = "dashed") +
+  labs(title = paste0("Maximum % variance explained in time series (",
+                      number_of_surrogates,
+                      " surrogates)"),
+       x = "Value of highest peak",
+       y = "Count",
+       caption = 
+         paste0("Dashed red square for real data; others are randomised copies.
+                P-value indicates probability of getting equally or more extreme result as the real data, if order was random.")) +
+  theme_bw() +
+  theme(legend.position = "none") +
+  facet_wrap(~name, scales = "free_x")
+
+ggsave(filename = paste0("ews-output/", 
+                         figure_identifier, "/", 
+                         "random_largest_eigenvalues_",
+                         figure_identifier, ".png"),
+       height = 7, width = 10)
+
+# To be used in the next plot: save min and max % variance explained
+ 
+varpercs_minmax <- dplyr::bind_rows(
+  dplyr::bind_rows(min_varpercs_store) %>% slice(-1), 
+  dplyr::bind_rows(max_varpercs_store) %>% slice(-1)) %>% 
+  tidyr::pivot_longer(cols = everything()) %>% 
+  dplyr::group_by(name) %>% 
+  dplyr::summarise(minvalue = min(value),
+                   maxvalue = max(value)) %>% 
+  dplyr::mutate(
+    name = factor(name, 
+                  levels = paste0("win_", window_lengths), ordered = TRUE)) %>% 
+  dplyr::arrange(name)
+
+
+```
+
+# Eigenvalue of the largest eigenvector
+
+```{r}
+#############ANALYSIS ############################
+### Eigenvalue of the largest eigenvector
+#using pca
+
+mattidata_withdate <- mattidata_orig 
+  # dplyr::mutate(across(is.numeric, 
+  #                      ~casnet::ts_standardise(., type = "median.mad")))
+
+omitted_from_pca <- c("date")
+
+mattidata <- mattidata_withdate %>% 
+  dplyr::select(-all_of(omitted_from_pca))
+
+date_of_max_eigval_1st <- c()
+window_lengths <- seq(from = 7, to = 4*7, by = 7)
+
+loadings_list <- list()
+eigval_list <- list()
+varperc_list <- list()
+mattidata_eigval <- list()
+date_of_max_eigval_1st <- list()
+date_of_max_eigval_2nd <- list()
+date_of_max_eigval_3rd <- list()
+mattidata_eigval_plots <- list()
+  
+for(k in 1:length(window_lengths)){
+
+  eig.val <- NA
+  varperc <- NA
+  loadings <- 
+    matrix(NA, 
+           nrow = 1, 
+           ncol = ncol(mattidata_withdate %>% 
+                         dplyr::select(-all_of(omitted_from_pca)))) %>% 
+    data.frame()
+  
+colnames(loadings) <- names(mattidata)
+
+  win <- window_lengths[[k]]
+
+for(i in (win:nrow(mattidata_withdate))){
+  res.pca <- prcomp(
+    mattidata_withdate %>% 
+      # For e.g. window of 28, keep observations 1-28, 2-29, 3-30 etc.
+      dplyr::filter(row_number() >= i-win+1 &
+                      row_number() <= i) %>% 
+      dplyr::select(-all_of(omitted_from_pca)),
+    scale = TRUE)
+  # Pull loadings of variables on the 1st principal component
+  loadings[(i-win+1), ] <- res.pca$rotation[ , 1]
+  # Pull eigenvalue of the 1st principal component
+  eig.val[i-win+1] <- factoextra::get_eigenvalue(res.pca)$eigenvalue[1] %>% 
+    sum()
+  # Pull % variance explained by first dimension
+  varperc[i-win+1] <- factoextra::get_eigenvalue(res.pca)$variance.percent[1] %>%
+    sum()
+
+  loadings_list[[k]] <- loadings
+  eigval_list[[k]] <-  eig.val 
+  varperc_list[[k]] <- varperc
+}
+
+mattidata_eigval[[k]] <- dplyr::bind_cols(mattidata_withdate, 
+                                          eigval = c(rep(NA, win - 1), 
+                                                     eig.val),
+                                          varperc = c(rep(NA, win - 1), 
+                                                     varperc))
+
+# Pull date of maximum eigenvalue in the first 1/3rd of data
+# How many days is a third of the data? 
+one_third_of_data <- # win + # window size to date
+  # Now take the remaining data and figure out, how many days is 
+  # a third of it. The 1st peak is grabbed from the 1st third of data.
+  (((nrow(mattidata_eigval[[k]]) - win) / 3) %>% floor(.))
+
+date_of_max_eigval_1st[[k]] <- 
+  mattidata_eigval[[k]] %>% 
+  # Only include the first third of data
+  dplyr::filter(
+    date < (min(mattidata_eigval[[k]]$date) + # first date available 
+            win + one_third_of_data)) %>% # add one third of remaining days
+  # Pull date of maximum eigenvalue in the first third of data
+  dplyr::filter(eigval == max(eigval, na.rm = TRUE)) %>% 
+  dplyr::pull(date)
+
+# For visualisation: pull date when the window, which leads to the peak,
+#  begins. Hence, peak contains info between these two dates.
+date_of_max_eigval_1st[[k]][[2]] <- date_of_max_eigval_1st[[k]] - win
+
+# Pull date of maximum eigenvalue in the second 1/3rd of data 
+date_of_max_eigval_2nd[[k]] <- 
+  mattidata_eigval[[k]] %>% 
+      dplyr::filter(
+        date >= (min(mattidata_eigval[[k]]$date)) + win + 
+          one_third_of_data &
+          date < min(mattidata_eigval[[k]]$date) + win +
+          2 * one_third_of_data) %>% 
+          dplyr::filter(eigval == max(eigval, na.rm = TRUE)) %>% 
+          dplyr::pull(date)
+        
+# For visualisation: grab date when the window, which leads to the peak,
+#  begins. Hence, peak contains info between these two dates.
+date_of_max_eigval_2nd[[k]][[2]] <- date_of_max_eigval_2nd[[k]] - win
+
+# Pull date of maximum eigenvalue in the last 1/3rd of data 
+date_of_max_eigval_3rd[[k]] <- 
+  mattidata_eigval[[k]] %>% 
+      dplyr::filter(date >= min(mattidata_eigval[[k]]$date) + win +
+                              2 * one_third_of_data) %>% 
+      dplyr::filter(eigval == max(eigval, na.rm = TRUE)) %>% 
+      dplyr::pull(date)
+    
+# For visualisation: grab date when the window, which leads to the peak,
+#  begins. Hence, peak contains info between these two dates.
+date_of_max_eigval_3rd[[k]][[2]] <- date_of_max_eigval_3rd[[k]] - win
+
+mattidata_eigval_plots[[k]] <- mattidata_eigval[[k]] %>% 
+  # dplyr::full_join(., epidemic_indicators) %>% 
+  # dplyr::mutate(across(c(eigval, new_cases),
+  #                      ~./max(., na.rm = TRUE))) %>% 
+  # dplyr::rename(`Eigenvalue\n(unit scale)` = eigval, 
+  #               `New cases\n(unit scale)` = new_cases) %>% 
+  # tidyr::pivot_longer(cols = c(`Eigenvalue\n(unit scale)`, 
+  #                              Rt, 
+  #                              `New cases\n(unit scale)`)) %>% 
+  # dplyr::mutate(across(c(eigval),
+  #                      ~./max(., na.rm = TRUE)),
+  #               varperc = varperc / 100) %>% 
+  dplyr::rename(`Eigenvalue` = eigval,
+                `Percentage of variance explained` = varperc) %>%
+  dplyr::select(-Eigenvalue) %>% 
+  tidyr::pivot_longer(cols = c(#`Eigenvalue`,
+                               `Percentage of variance explained`)) %>%
+  tidyr::drop_na(value) %>% 
+  dplyr::group_by(name) %>% 
+  dplyr::mutate(label = dplyr::case_when(date == max(date) ~ name,
+                                          TRUE ~ as.character(NA))) %>% 
+  dplyr::ungroup() %>% 
+  ggplot(aes(x = date,
+             y = value,
+             colour = name)) +
+  # ### same max and min shading for the whole plot
+  # geom_ribbon(aes_string(ymin = varpercs_minmax$minvalue[[k]],
+  #                 ymax = varpercs_minmax$maxvalue[[k]],
+  #                 alpha = 0.2,
+  #                 colour = NA),
+  #             fill = "lightgrey") +
+  ### max and min shading based on max and min of the day
+  geom_ribbon(aes_string(ymin = randomised_varpercs_daily[[k]]$min_varperc_for_day,
+                         ymax = randomised_varpercs_daily[[k]]$max_varperc_for_day,
+                         alpha = 0.2,
+                         colour = NA),
+              fill = "lightgrey") +
+  geom_line() +
+  geom_vline(xintercept = date_of_max_eigval_1st[[k]][[2]],
+             color = "grey",
+             linetype = "dotted") +
+  geom_vline(xintercept = date_of_max_eigval_1st[[k]][[1]],
+             color = "red",
+             linetype = "dashed") +
+  geom_vline(xintercept = date_of_max_eigval_2nd[[k]][[2]],
+             color = "grey",
+             linetype = "dotted") +
+  geom_vline(xintercept = date_of_max_eigval_2nd[[k]][[1]],
+             color = "red",
+             linetype = "dashed") +
+  geom_vline(xintercept = date_of_max_eigval_3rd[[k]][[2]],
+             color = "grey",
+             linetype = "dotted") +
+  geom_vline(xintercept = date_of_max_eigval_3rd[[k]][[1]],
+             color = "red",
+             linetype = "dashed") +
+  labs(title = paste0("Window size: ", win),
+       y = NULL) +
+  scale_x_date(name = NULL,
+               date_breaks = "1 week", 
+               date_labels = "%F",
+               # expand = expansion(mult = c(0.1, 0.20))
+               ) +
+  # scale_y_continuous(breaks = seq(from = 0, to = 100, by = 5),
+  #                    labels = seq(from = 0, to = 100, by = 5)) +
+  scale_colour_viridis_d(option = "viridis", end = 0.8) +
+  # ggrepel::geom_text_repel(aes(x = date,
+  #                              y = value,
+  #                              label = label),
+  #                          # xlim = c(max(mattidata_eigval[[k]]$date) + 2, 
+  #                          #          max(mattidata_eigval[[k]]$date) + 3),
+  #                          box.padding = 0.5,
+  #                          nudge_x = 2,
+  #                          force_pull = 0,
+  #                          direction = "y",
+  #                          hjust = "left",
+  #                          min.segment.length = 100,
+  #                          segment.linetype = "dotted",
+  #                          na.rm = TRUE,
+  #                          max.overlaps = Inf,
+  #                          size = 3) +
+  theme_bw() +
+  theme(axis.text.x = element_text(angle = 45, 
+                                   vjust = 1, 
+                                   hjust = 1),
+        panel.grid.minor.x = element_blank(),
+        legend.position = "none",
+        legend.title = element_blank()) +
+  facet_wrap(~name, 
+             scales = "free_y",
+             ncol = 1)
+
+}
+
+mattidata_eigval_plots
+
+purrr::pmap(list(..1 = mattidata_eigval_plots,
+                 ..2 = window_lengths),
+            .f = ~ggsave(filename = 
+                           paste0("ews-output/", 
+                                  figure_identifier, "/", 
+                                  "loadings_1st_peak_win_",
+                                  ..2, figure_identifier, ".png"),
+                         plot = ..1,
+                         height = 7, width = 10))
+
+```
+
+### First peak loadings
+
+```{r}
+
+position_of_1st_peak <- list()
+peak_plots_1st <- list()
+
+peak_loadings_1st_temp <- list()
+peak_loadings_1st <- list()
+peak_viz <- list()
+top_loadings_at_1st_peak <- list()
+alphavalues <- list()
+maxvalue <- list()
+
+# How many days around the peak to include in visualisation?
+inspection_span <- 20
+
+# For each window length, pull row number of the date on which the peak 
+# (i.e. maximum eigenvalue) resides
+for(i in 1:length(window_lengths)){
+position_of_1st_peak[[i]] <- mattidata_eigval[[i]] %>% 
+  # Need to substract window length to make row number compatible with 
+  # loadings_list and eigval_list
+  dplyr::mutate(rownum = dplyr::row_number() - window_lengths[[i]] + 1) %>% 
+  dplyr::filter(date %in% date_of_max_eigval_1st[[i]][[1]]) %>% 
+  pull(rownum)
+}
+
+# For every window length: Pull the loadings for the peak day, 
+# as well as the surrounding inspection_span days
+for(i in 1:length(window_lengths)){
+  for(k in (-1 * inspection_span):inspection_span){ 
+    # When k == -20, placeholder in list equals 1
+    peak_loadings_1st_temp[[k + (inspection_span + 1)]] <- 
+      loadings_list[[i]] %>% 
+      # if_else assures row number to be pulled is not outside the bounds
+      # of the data frame
+      dplyr::filter(row_number() == dplyr::if_else(
+        condition = 
+          (position_of_1st_peak[[i]] + k < 0) |
+            (position_of_1st_peak[[i]] + k > nrow(loadings_list[[i]])),
+        true = 0,
+        false = position_of_1st_peak[[i]] + k,
+        missing = NULL)) %>% 
+      tidyr::pivot_longer(cols = everything()) %>% 
+      dplyr::arrange(desc(abs(value))) %>% 
+      dplyr::mutate(days_since_peak = k,
+                    date = date_of_max_eigval_1st[[i]][[1]] + k,
+                    window_size = window_lengths[[i]])
+    
+    # If number to be pulled is outside data, don't do it
+    if((position_of_1st_peak[[i]] + k < 0) |
+            (position_of_1st_peak[[i]] + k > nrow(loadings_list[[i]]))){
+    print("Out of bounds")
+    # Otherwise, add in the corresponding eigenvalue
+    } else {
+    peak_loadings_1st_temp[[k + (inspection_span + 1)]] <- 
+      peak_loadings_1st_temp[[k + (inspection_span + 1)]] %>% 
+      dplyr::mutate(eigval = eigval_list[[i]][position_of_1st_peak[[i]] + k])
+    }
+    
+  }
+  peak_loadings_1st[[i]] <- purrr::reduce(.x = peak_loadings_1st_temp,
+                                             .f = full_join)
+}
+
+top_loadings_at_1st_peak <- purrr::map(.x = peak_loadings_1st,
+                                   .f = ~.x %>% 
+                                     dplyr::filter(days_since_peak == 0) %>% 
+                                     dplyr::arrange(desc(abs(value))) %>% 
+                                     dplyr::mutate(value = row_number()) %>% 
+                                     head(15) %>% 
+                                     dplyr::select(name, colourvalue = value))
+ 
+peak_viz <- purrr::pmap(list(..1 = peak_loadings_1st,
+                             ..2 = top_loadings_at_1st_peak),
+                        .f = ~dplyr::full_join(..1, ..2) %>% 
+  dplyr::mutate(absvalue = abs(value)) %>% 
+  dplyr::group_by(name) %>% 
+  dplyr::mutate(label = dplyr::case_when(absvalue == max(absvalue) 
+                                         ~ name,
+                                         TRUE ~ as.character(NA))) %>% 
+  dplyr::ungroup() %>% 
+  dplyr::mutate(sign = dplyr::case_when(sign(value) >= 0 ~ "positive",
+                                        sign(value) < 0 ~ "negative"),
+                label2 = dplyr::case_when(!is.na(colourvalue) & 
+                                               date == max(date) ~ name,
+                                          TRUE ~ as.character(NA)),
+                peak_date = dplyr::case_when(days_since_peak == 0 ~ date,
+                                             TRUE ~ as.Date(NA)),
+                eigval = scales::rescale(eigval, to = c(0, 0.35))))
+
+alphavalues <- purrr::map(.x = peak_viz,
+                          .f = ~ifelse(!is.na(.x$colourvalue),
+                      1, 0.75))
+maxvalue <- purrr::map(.x = peak_viz,
+                       .f = ~max(.x$eigval))
+
+peak_plots_1st <- purrr::pmap(list(..1 = peak_viz,
+                                   ..2 = alphavalues,
+                                   ..3 = maxvalue,
+                                   ..4 = window_lengths),
+                              .f = ~..1 %>%  
+  ggplot(aes(x = date,
+             y = absvalue,
+             colour = colourvalue,
+             group = name)) +
+  geom_vline(aes(xintercept = peak_date),
+             linetype = "dashed",
+             colour = "red") + 
+  geom_line(aes(x = date,
+                y = eigval),
+            size = 2,
+            colour = "grey",
+            linetype = "dotted") +
+  geom_line(aes(alpha = ..2)) +
+  scale_x_date(name = NULL,
+               date_breaks = "2 days", 
+               date_labels = "%F",
+               expand = expansion(mult = c(0.1, 0.55))) +
+  scale_colour_viridis_c(option = "inferno",
+                         end = 0.8,
+                         direction = -1) +
+  scale_alpha_continuous(range = c(0.25, 1)) +
+  ggrepel::geom_text_repel(aes(x = date,
+                               y = absvalue,
+                               label = label2),
+                           xlim = c(max(..1$date) + 3, 
+                                    max(..1$date) + 30),
+                           box.padding = 0.5,
+                           nudge_x = 2,
+                           force_pull = 0,
+                           direction = "y",
+                           hjust = "left",
+                           min.segment.length = 0,
+                           segment.linetype = "dotted",
+                           na.rm = TRUE,
+                           max.overlaps = Inf,
+                           size = 3) +
+  coord_cartesian(ylim = c(0, ..3)) +
+  labs(x = NULL,
+       y = "Absolute value of loading",
+       title = paste("Loadings on the first principal component around 1st peak | Window size:",
+                     ..4)) +
+  theme_bw() +
+  theme(legend.position = "none",
+        axis.text.x = element_text(angle = 30, 
+                                   vjust = 1, 
+                                   hjust = 1),
+        panel.grid.minor.x = element_blank()))
+
+peak_plots_1st
+
+purrr::pmap(list(..1 = peak_plots_1st,
+                 ..2 = window_lengths),
+            .f = ~ggsave(filename = paste0("ews-output/", figure_identifier, "/", "loadings_1st_peak_win_", 
+                                           ..2, figure_identifier, ".png"),
+                         plot = ..1,
+                         height = 7, width = 10))
+
+```
+
+### Second peak loadings
+
+```{r}
+
+position_of_2nd_peak <- list()
+peak_plots_2nd <- list()
+
+peak_loadings_2nd_temp <- list()
+peak_loadings_2nd <- list()
+peak_viz <- list()
+top_loadings_at_2nd_peak <- list()
+alphavalues <- list()
+maxvalue <- list()
+
+# How many days around the peak to include in visualisation?
+inspection_span <- 20
+
+# For each window length, pull row number of the date on which the peak 
+# (i.e. maximum eigenvalue) resides
+for(i in 1:length(window_lengths)){
+position_of_2nd_peak[[i]] <- mattidata_eigval[[i]] %>% 
+  # Need to substract window length to make row number compatible with 
+  # loadings_list and eigval_list
+  dplyr::mutate(rownum = dplyr::row_number() - window_lengths[[i]] + 1) %>% 
+  dplyr::filter(date %in% date_of_max_eigval_2nd[[i]][[1]]) %>% 
+  pull(rownum)
+}
+
+# For every window length: Pull the loadings for the peak day, 
+# as well as the surrounding inspection_span days
+for(i in 1:length(window_lengths)){
+  for(k in (-1 * inspection_span):inspection_span){ 
+    # When k == -20, placeholder in list equals 1
+    peak_loadings_2nd_temp[[k + (inspection_span + 1)]] <- 
+      loadings_list[[i]] %>% 
+      # if_else assures row number to be pulled is not outside the bounds
+      # of the data frame
+      dplyr::filter(row_number() == dplyr::if_else(
+        condition = 
+          (position_of_2nd_peak[[i]] + k < 0) |
+            (position_of_2nd_peak[[i]] + k > nrow(loadings_list[[i]])),
+        true = 0,
+        false = position_of_2nd_peak[[i]] + k,
+        missing = NULL)) %>% 
+      tidyr::pivot_longer(cols = everything()) %>% 
+      dplyr::arrange(desc(abs(value))) %>% 
+      dplyr::mutate(days_since_peak = k,
+                    date = date_of_max_eigval_2nd[[i]][[1]] + k,
+                    window_size = window_lengths[[i]])
+    
+    # If number to be pulled is outside data, don't do it
+    if((position_of_2nd_peak[[i]] + k < 0) |
+            (position_of_2nd_peak[[i]] + k > nrow(loadings_list[[i]]))){
+    print("Out of bounds")
+    # Otherwise, add in the corresponding eigenvalue
+    } else {
+    peak_loadings_2nd_temp[[k + (inspection_span + 1)]] <- 
+      peak_loadings_2nd_temp[[k + (inspection_span + 1)]] %>% 
+      dplyr::mutate(eigval = eigval_list[[i]][position_of_2nd_peak[[i]] + k])
+    }
+    
+  }
+  peak_loadings_2nd[[i]] <- purrr::reduce(.x = peak_loadings_2nd_temp,
+                                             .f = full_join)
+}
+
+top_loadings_at_2nd_peak <- purrr::map(.x = peak_loadings_2nd,
+                                   .f = ~.x %>% 
+                                     dplyr::filter(days_since_peak == 0) %>% 
+                                     dplyr::arrange(desc(abs(value))) %>% 
+                                     dplyr::mutate(value = row_number()) %>% 
+                                     head(15) %>% 
+                                     dplyr::select(name, colourvalue = value))
+ 
+peak_viz <- purrr::pmap(list(..1 = peak_loadings_2nd,
+                             ..2 = top_loadings_at_2nd_peak),
+                        .f = ~dplyr::full_join(..1, ..2) %>% 
+  dplyr::mutate(absvalue = abs(value)) %>% 
+  dplyr::group_by(name) %>% 
+  dplyr::mutate(label = dplyr::case_when(absvalue == max(absvalue) 
+                                         ~ name,
+                                         TRUE ~ as.character(NA))) %>% 
+  dplyr::ungroup() %>% 
+  dplyr::mutate(sign = dplyr::case_when(sign(value) >= 0 ~ "positive",
+                                        sign(value) < 0 ~ "negative"),
+                label2 = dplyr::case_when(!is.na(colourvalue) & 
+                                               date == max(date) ~ name,
+                                          TRUE ~ as.character(NA)),
+                peak_date = dplyr::case_when(days_since_peak == 0 ~ date,
+                                             TRUE ~ as.Date(NA)),
+                eigval = scales::rescale(eigval, to = c(0, 0.35))))
+
+alphavalues <- purrr::map(.x = peak_viz,
+                          .f = ~ifelse(!is.na(.x$colourvalue),
+                      1, 0.75))
+maxvalue <- purrr::map(.x = peak_viz,
+                       .f = ~max(.x$eigval))
+
+peak_plots_2nd <- purrr::pmap(list(..1 = peak_viz,
+                                   ..2 = alphavalues,
+                                   ..3 = maxvalue,
+                                   ..4 = window_lengths),
+                              .f = ~..1 %>%  
+  ggplot(aes(x = date,
+             y = absvalue,
+             colour = colourvalue,
+             group = name)) +
+  geom_vline(aes(xintercept = peak_date),
+             linetype = "dashed",
+             colour = "red") + 
+  geom_line(aes(x = date,
+                y = eigval),
+            size = 2,
+            colour = "grey",
+            linetype = "dotted") +
+  geom_line(aes(alpha = ..2)) +
+  scale_x_date(name = NULL,
+               date_breaks = "2 days", 
+               date_labels = "%F",
+               expand = expansion(mult = c(0.1, 0.55))) +
+  scale_colour_viridis_c(option = "inferno",
+                         end = 0.8,
+                         direction = -1) +
+  scale_alpha_continuous(range = c(0.25, 1)) +
+  ggrepel::geom_text_repel(aes(x = date,
+                               y = absvalue,
+                               label = label2),
+                           xlim = c(max(..1$date) + 2, 
+                                    max(..1$date) + 30),
+                           box.padding = 0.5,
+                           nudge_x = 2,
+                           force_pull = 0,
+                           direction = "y",
+                           hjust = "left",
+                           min.segment.length = 0,
+                           segment.linetype = "dotted",
+                           na.rm = TRUE,
+                           max.overlaps = Inf,
+                           size = 3) +
+  coord_cartesian(ylim = c(0, ..3)) +
+  labs(x = NULL,
+       y = "Absolute value of loading",
+       title = paste("Loadings on the first principal component around 2nd peak | Window size:",
+                     ..4)) +
+  theme_bw() +
+  theme(legend.position = "none",
+        axis.text.x = element_text(angle = 30, 
+                                   vjust = 1, 
+                                   hjust = 1),
+        panel.grid.minor.x = element_blank()))
+
+peak_plots_2nd
+
+purrr::pmap(list(..1 = peak_plots_2nd,
+                 ..2 = window_lengths),
+            .f = ~ggsave(filename = paste0("ews-output/", figure_identifier, "/", "loadings_2nd_peak_win_", 
+                                           ..2, figure_identifier, ".png"),
+                         plot = ..1,
+                         height = 7, width = 10))
+
+```
+
+### Third peak loadings
+
+```{r}
+
+position_of_3rd_peak <- list()
+peak_plots_3rd <- list()
+
+peak_loadings_3rd_temp <- list()
+peak_loadings_3rd <- list()
+peak_viz <- list()
+top_loadings_at_3rd_peak <- list()
+alphavalues <- list()
+maxvalue <- list()
+
+# How many days around the peak to include in visualisation?
+inspection_span <- 20
+
+# For each window length, pull row number of the date on which the peak 
+# (i.e. maximum eigenvalue) resides
+for(i in 1:length(window_lengths)){
+position_of_3rd_peak[[i]] <- mattidata_eigval[[i]] %>% 
+  # Need to substract window length to make row number compatible with 
+  # loadings_list and eigval_list
+  dplyr::mutate(rownum = dplyr::row_number() - window_lengths[[i]] + 1) %>% 
+  dplyr::filter(date %in% date_of_max_eigval_3rd[[i]][[1]]) %>% 
+  pull(rownum)
+}
+
+# For every window length: Pull the loadings for the peak day, 
+# as well as the surrounding inspection_span days
+for(i in 1:length(window_lengths)){
+  for(k in (-1 * inspection_span):inspection_span){ 
+    # When k == -20, placeholder in list equals 1
+    peak_loadings_3rd_temp[[k + (inspection_span + 1)]] <- 
+      loadings_list[[i]] %>% 
+      # if_else assures row number to be pulled is not outside the bounds
+      # of the data frame
+      dplyr::filter(row_number() == dplyr::if_else(
+        condition = 
+          (position_of_3rd_peak[[i]] + k < 0) |
+            (position_of_3rd_peak[[i]] + k > nrow(loadings_list[[i]])),
+        true = 0,
+        false = position_of_3rd_peak[[i]] + k,
+        missing = NULL)) %>% 
+      tidyr::pivot_longer(cols = everything()) %>% 
+      dplyr::arrange(desc(abs(value))) %>% 
+      dplyr::mutate(days_since_peak = k,
+                    date = date_of_max_eigval_3rd[[i]][[1]] + k,
+                    window_size = window_lengths[[i]])
+    
+    # If number to be pulled is outside data, don't do it
+    if((position_of_3rd_peak[[i]] + k < 0) |
+            (position_of_3rd_peak[[i]] + k > nrow(loadings_list[[i]]))){
+    print("Out of bounds")
+    # Otherwise, add in the corresponding eigenvalue
+    } else {
+    peak_loadings_3rd_temp[[k + (inspection_span + 1)]] <- 
+      peak_loadings_3rd_temp[[k + (inspection_span + 1)]] %>% 
+      dplyr::mutate(eigval = eigval_list[[i]][position_of_3rd_peak[[i]] + k])
+    }
+    
+  }
+  peak_loadings_3rd[[i]] <- purrr::reduce(.x = peak_loadings_3rd_temp,
+                                             .f = full_join)
+}
+
+top_loadings_at_3rd_peak <- purrr::map(.x = peak_loadings_3rd,
+                                   .f = ~.x %>% 
+                                     dplyr::filter(days_since_peak == 0) %>% 
+                                     dplyr::arrange(desc(abs(value))) %>% 
+                                     dplyr::mutate(value = row_number()) %>% 
+                                     head(15) %>% 
+                                     dplyr::select(name, colourvalue = value))
+ 
+peak_viz <- purrr::pmap(list(..1 = peak_loadings_3rd,
+                             ..2 = top_loadings_at_3rd_peak),
+                        .f = ~dplyr::full_join(..1, ..2) %>% 
+  dplyr::mutate(absvalue = abs(value)) %>% 
+  dplyr::group_by(name) %>% 
+  dplyr::mutate(label = dplyr::case_when(absvalue == max(absvalue) 
+                                         ~ name,
+                                         TRUE ~ as.character(NA))) %>% 
+  dplyr::ungroup() %>% 
+  dplyr::mutate(sign = dplyr::case_when(sign(value) >= 0 ~ "positive",
+                                        sign(value) < 0 ~ "negative"),
+                label2 = dplyr::case_when(!is.na(colourvalue) & 
+                                               date == max(date) ~ name,
+                                          TRUE ~ as.character(NA)),
+                peak_date = dplyr::case_when(days_since_peak == 0 ~ date,
+                                             TRUE ~ as.Date(NA)),
+                eigval = scales::rescale(eigval, to = c(0, 0.35))))
+
+alphavalues <- purrr::map(.x = peak_viz,
+                          .f = ~ifelse(!is.na(.x$colourvalue),
+                      1, 0.75))
+maxvalue <- purrr::map(.x = peak_viz,
+                       .f = ~max(.x$eigval))
+
+peak_plots_3rd <- purrr::pmap(list(..1 = peak_viz,
+                                   ..2 = alphavalues,
+                                   ..3 = maxvalue,
+                                   ..4 = window_lengths),
+                              .f = ~..1 %>%  
+  ggplot(aes(x = date,
+             y = absvalue,
+             colour = colourvalue,
+             group = name)) +
+  geom_vline(aes(xintercept = peak_date),
+             linetype = "dashed",
+             colour = "red") + 
+  geom_line(aes(x = date,
+                y = eigval),
+            size = 2,
+            colour = "grey",
+            linetype = "dotted") +
+  geom_line(aes(alpha = ..2)) +
+  scale_x_date(name = NULL,
+               date_breaks = "2 days", 
+               date_labels = "%F",
+               expand = expansion(mult = c(0.1, 0.55))) +
+  scale_colour_viridis_c(option = "inferno",
+                         end = 0.8,
+                         direction = -1) +
+  scale_alpha_continuous(range = c(0.25, 1)) +
+  ggrepel::geom_text_repel(aes(x = date,
+                               y = absvalue,
+                               label = label2),
+                           xlim = c(max(..1$date) + 2, 
+                                    max(..1$date) + 30),
+                           box.padding = 0.5,
+                           nudge_x = 2,
+                           force_pull = 0,
+                           direction = "y",
+                           hjust = "left",
+                           min.segment.length = 0,
+                           segment.linetype = "dotted",
+                           na.rm = TRUE,
+                           max.overlaps = Inf,
+                           size = 3) +
+  coord_cartesian(ylim = c(0, ..3)) +
+  labs(x = NULL,
+       y = "Absolute value of loading",
+       title = paste("Loadings on the first principal component around 3rd peak | Window size:",
+                     ..4)) +
+  theme_bw() +
+  theme(legend.position = "none",
+        axis.text.x = element_text(angle = 30, 
+                                   vjust = 1, 
+                                   hjust = 1),
+        panel.grid.minor.x = element_blank()))
+
+peak_plots_3rd
+
+purrr::pmap(list(..1 = peak_plots_3rd,
+                 ..2 = window_lengths),
+            .f = ~ggsave(filename = paste0("ews-output/", figure_identifier, "/", "loadings_3rd_peak_win_", 
+                                           ..2, figure_identifier, ".png"),
+                         plot = ..1,
+                         height = 7, width = 10))
+
+```
+
+### Combine figures
+
+```{r}
+
+combined_pca_plots <- purrr::pmap(
+  list(..1 = mattidata_eigval_plots,
+       ..2 = peak_plots_1st,
+       ..3 = peak_plots_2nd,
+       ..4 = peak_plots_3rd,
+       ..5 = as.list(window_lengths)),
+  .f = ~egg::ggarrange(
+    plots = list(..1 + ggtitle("Eigenvalue of the largest eigenvector"),
+                 ..2 + ggtitle("Loadings around the first peak"),
+                 ..3 + ggtitle("Loadings around the second peak"),
+                 ..4 + ggtitle("Loadings around the third peak")), 
+    top = paste("Window length:", ..5),
+    ncol = 1))
+           
+for (i in 1:length(window_lengths)){
+ggsave(plot = combined_pca_plots[[i]], 
+       filename = paste0("ews-output/", figure_identifier, 
+                         "/", "combined_figure_win_", 
+                         window_lengths[[i]], figure_identifier, 
+                         ".png"), 
+       height = 15)
+}
+
+```
+
+
+# Inspect time series
+
+## Extract change points 
+
+```{r}
+# Which variables are in any window size's top-15 at the 3rd peak? 
+top_loaders <- top_loadings_at_3rd_peak %>% 
+  dplyr::bind_rows() %>%
+  dplyr::pull(name) %>% 
+  unique()
+  
+data_for_cp <- mattidata_withdate %>% 
+  dplyr::select(date, 
+                # Select just the top loading variables
+                all_of(top_loaders))
+
+### Data frame for change points
+
+changepoints <- data_for_cp %>%
+  tidyr::pivot_longer(-date) %>% 
+  dplyr::group_by(name) %>% 
+  tidyr::nest() %>% 
+  dplyr::mutate(
+    # change_points_tsmcp = purrr::map(
+    # .x = data,
+    # .f = ~TSMCP::tsmcplm(
+    #   Y = .x[["value"]],
+    #   X = NULL,
+    #   method = "adapt",
+    #   c = 1
+    # )),
+    # tsmcp_dates = purrr::pmap(list(..1 = name,
+    #                                ..2 = data,
+    #                                ..3 = change_points_tsmcp),
+    #                           .f = ~..2 %>% dplyr::slice(..3) %>%
+    #                             dplyr::mutate(name = ..1)),
+    change_points_ecp = purrr::map(
+      .x = data,
+      .f = ~ecp::e.divisive(data.matrix(.x$value), 
+                            R = 1000,
+                            min.size = 7,
+                            alpha = 1)),
+    ecp_dates = purrr::pmap(list(..1 = name,
+                                   ..2 = data,
+                                   ..3 = change_points_ecp),  
+                              .f = ~..2 %>% dplyr::slice(..3$estimates) %>% 
+                                dplyr::mutate(name = ..1)))
+
+# tsmcp_changedates <- changepoints$tsmcp_dates %>% 
+#   dplyr::bind_rows()
+
+ecp_changedates <- changepoints$ecp_dates %>% 
+  dplyr::bind_rows()
+
+```
+
+
+## Build figure
+
+```{r}
+sample_sizes <- readr::read_csv(paste0("../shield-complexity/data/",
+                                       country_in_question,
+                                       "_sample_sizes_since_2021-06-08.csv"))
+
+std_errors <- readr::read_csv(paste0("../shield-complexity/data/",
+                                       country_in_question,
+                                       "_std_errors_since_2021-06-08.csv"))
+
+for (i in 1:length(window_lengths)){
+### Data frame for data in question for this window length
+data_for_viz <- mattidata_withdate %>% 
+  dplyr::select(date, 
+                # Pull names of top loadings and select just those variables
+                top_loadings_at_2nd_peak[[i]] %>% dplyr::pull(name)) %>% 
+  tidyr::pivot_longer(-date) %>% 
+  dplyr::mutate(name = forcats::fct_relevel(name, 
+                                            top_loadings_at_2nd_peak[[i]] %>%
+                                              dplyr::pull(name))) 
+
+### Data frame for visualisible sample sizes
+minvalues <- data_for_viz %>% 
+  dplyr::group_by(name) %>% 
+  dplyr::summarise(minvalue = min(value),
+                   maxvalue = max(value)) %>% 
+  dplyr::mutate(diff = maxvalue - minvalue)
+
+sample_sizes_for_viz <- sample_sizes %>% 
+  dplyr::filter(name %in% top_loadings_at_2nd_peak[[i]]$name) %>% 
+  dplyr::full_join(., minvalues) %>% 
+  dplyr::group_by(name) %>% 
+  dplyr::mutate(minvalue_for_viz = min(minvalue) - min(diff) * 0.5,
+                value = scales::rescale(x = sample_size,
+                                        to = c(min(minvalue_for_viz), 
+                                               min(minvalue))),
+                date = as.Date(date)) %>% 
+  dplyr::ungroup() %>% 
+  dplyr::mutate(name = forcats::fct_relevel(name, 
+                                            top_loadings_at_2nd_peak[[i]] %>%
+                                              dplyr::pull(name)))
+
+std_errors_for_viz <- std_errors %>% 
+  dplyr::filter(name %in% top_loadings_at_2nd_peak[[i]]$name) %>% 
+  dplyr::full_join(., minvalues) %>% 
+  dplyr::group_by(name) %>% 
+  dplyr::mutate(minvalue_for_viz = min(minvalue) - min(diff) * 0.5,
+                value = scales::rescale(x = std_error,
+                                        to = c(min(minvalue_for_viz), 
+                                               min(minvalue))),
+                date = as.Date(date)) %>% 
+  dplyr::ungroup() %>% 
+  dplyr::mutate(name = forcats::fct_relevel(name, 
+                                            top_loadings_at_2nd_peak[[i]] %>%
+                                              dplyr::pull(name)))
+
+# tsmcp_changedates_viz <- tsmcp_changedates %>% 
+#   dplyr::filter(name %in% top_loadings_at_2nd_peak[[i]]$name) %>% 
+#   dplyr::mutate(name = forcats::fct_relevel(name, 
+#                                             top_loadings_at_2nd_peak[[i]] %>%
+#                                               dplyr::pull(name)))
+
+ecp_changedates_viz <- ecp_changedates %>% 
+  dplyr::filter(name %in% top_loadings_at_2nd_peak[[i]]$name) %>% 
+  dplyr::mutate(name = forcats::fct_relevel(name, 
+                                            top_loadings_at_2nd_peak[[i]] %>%
+                                              dplyr::pull(name)))
+
+### Draw plot
+data_for_viz %>% 
+  ggplot(aes(x = date,
+             y = value,
+             colour = name)) +
+  # geom_smooth(aes(x = date,
+  #                 y = value,
+  #                 colour = name),
+  #             se = FALSE,
+  #             span = 0.15) +
+  geom_line(aes(x = date,
+                y = value,
+                colour = name)) +
+  annotate(geom = "rect",
+           xmin = date_of_max_eigval_2nd[[i]][[1]],
+           xmax = date_of_max_eigval_2nd[[i]][[2]],
+           ymin = -Inf,
+           ymax = Inf,
+           alpha = 0.1,
+           fill = "red") + 
+  geom_vline(xintercept = date_of_max_eigval_2nd[[i]],
+             linetype = c("dashed"),
+             colour = "red",
+             alpha = 0.5) +  
+  scale_colour_viridis_d(option = "inferno",
+                         end = 0.8,
+                         direction = -1) +
+  scale_x_date(date_labels = "%y-%m-%d",
+               date_breaks = "1 month") +
+  labs(x = NULL,
+       y = NULL,
+       title = paste0(figure_identifier, 
+                      " | Variables with highest loadings at 2nd peak (window size: ",
+                      window_lengths[[i]], ")"),
+       caption = "Grey bars: rescaled sample size | Grey points: change points | Dashed lines: the window at 2nd peak
+       Darker colours in solid lines indicate lower loadings (black = variable with 15th highest loading)") +
+  theme_bw() +
+  theme(legend.position = "none",
+        axis.text.x = element_text(angle = 30, hjust = 1)) +
+  geom_segment(data = sample_sizes_for_viz, #std_errors_for_viz,
+               aes(xend = date,
+                   y = -Inf, #minvalue_for_viz,
+                   yend = value),
+            colour = "grey") +
+  geom_point(data = ecp_changedates_viz,
+             aes(x = date,
+                 y = value),
+             shape = 18,
+             size = 2,
+             colour = "darkgrey") +
+  # geom_point(data = tsmcp_changedates_viz,
+  #            aes(x = date,
+  #                y = value),
+  #            shape = 4,
+  #            size = 2,
+  #            colour = "darkgrey") +
+  facet_wrap(~name, scales = "free_y", nrow = 3)
+
+ggsave(paste0("ews-output/", 
+              figure_identifier, "/", 
+              "values_of_top_vars_at_2nd_peak", "_win_", 
+              window_lengths[[i]], 
+              figure_identifier,
+              ".png"),
+       width = 12)
+}
+
+data_omni <- dplyr::full_join(data_for_viz,
+                 std_errors_for_viz,
+                 by = c("name", "date")) %>% 
+  dplyr::full_join(.,
+                   sample_sizes_for_viz,
+                 by = c("name", "date")) %>% 
+  dplyr::group_by(name) %>% 
+  tidyr::nest() %>% 
+  dplyr::mutate(
+    meanvalue = purrr::map_dbl(
+      .x = data,
+      .f = ~mean(.x[["value.x"]], na.rm = TRUE)),
+    correlation_with_std_err = purrr::map_dbl(
+      .x = data,
+      .f = ~cor(x = .x[["value.x"]], 
+                y = .x[["std_error"]],
+                use = "pairwise.complete.obs")[[1]]
+    )) 
+  
+cor(data_omni$meanvalue, data_omni$correlation_with_std_err,
+      use = "everything")
+
+data_omni %>% 
+  dplyr::select(name, meanvalue, correlation_with_std_err) %>% 
+  dplyr::ungroup() %>% 
+  ggplot(aes(x = meanvalue, y = correlation_with_std_err)) +
+  geom_point() + 
+  ggrepel::geom_text_repel(aes(label = name),
+                           size = 1.4) +
+  # geom_label(aes(label =))
+  theme_bw() 
+
+```
+
+# To export data
+
+```{r}
+
+dates_of_max_eigvals <- list(date_of_max_eigval_1st, 
+                             date_of_max_eigval_2nd, 
+                             date_of_max_eigval_3rd)  %>% 
+  purrr::map(
+    .f = ~purrr::set_names(x = .x,
+                           nm = paste0("win_", window_lengths)))
+
+saveRDS(object = dates_of_max_eigvals,
+        paste0("ews-output/", 
+              figure_identifier, "/", 
+              "dates_of_max_eigvals", 
+              figure_identifier,
+              ".RDS"))
+
+```
+
+
+# Old plots w/o sample sizes
+
+```{r}
+
+# for (i in 1:length(window_lengths)){
+# mattidata_withdate %>% 
+#   dplyr::select(date, 
+#                 # Pull names of top loadings and select just those variables
+#                 top_loadings_at_3rd_peak[[i]] %>% dplyr::pull(name)) %>% 
+#   tidyr::pivot_longer(-date) %>% 
+#   dplyr::mutate(name = forcats::fct_relevel(name, 
+#                                             top_loadings_at_3rd_peak[[i]] %>%
+#                                               dplyr::pull(name))) %>% 
+#   ggplot(aes(x = date,
+#              y = value,
+#              colour = name)) +
+#   # geom_smooth(aes(x = date,
+#   #                 y = value,
+#   #                 colour = name),
+#   #             se = FALSE,
+#   #             span = 0.15) +
+#   geom_line(aes(x = date,
+#                 y = value,
+#                 colour = name)) +
+#   annotate(geom = "rect",
+#            xmin = date_of_max_eigval_3rd[[i]][[1]],
+#            xmax = date_of_max_eigval_3rd[[i]][[2]],
+#            ymin = -Inf,
+#            ymax = Inf,
+#            alpha = 0.1,
+#            fill = "red") +
+#   geom_vline(xintercept = date_of_max_eigval_3rd[[i]],
+#              linetype = c("dashed"),
+#              colour = "red",
+#              alpha = 0.5) +  
+#   scale_colour_viridis_d(option = "inferno",
+#                          end = 0.8,
+#                          direction = -1) +
+#   scale_x_date(date_labels = "%y-%m-%d") +
+#   labs(x = NULL,
+#        y = NULL,
+#        title = paste0(figure_identifier, 
+#                       " | Variables with highest loadings at 3rd peak (window size: ", window_lengths[[i]], ")"),
+#        caption = "Dashed lines indicate the window at 3rd peak. 
+#        Darker colours indicate lower loadings (black = variable with 15th highest loading)") +
+#   theme_bw() +
+#   theme(legend.position = "none",
+#         axis.text.x = element_text(angle = 30, hjust = 1)) +
+#   facet_wrap(~name, scales = "free_y", nrow = 3)
+# 
+# ggsave(paste0("ews-output/", figure_identifier, "/", "values_of_top_vars_at_3rd_peak", "_win_", window_lengths[[i]], 
+#               figure_identifier,
+#               ".png"),
+#        width = 12)
+# }
+```
+
+