---
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)
# }
```