James
1422 lines
55 KiB
Plaintext
1422 lines
55 KiB
Plaintext
---
|
|
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)
|
|
# }
|
|
```
|
|
|
|
|