# Create correlation plot with consistent styling
plot_correlation_matrix <- function(cor_matrix, title) {
corrplot(
cor_matrix,
method = "color",
type = "lower",
order = "hclust",
tl.col = "black",
tl.srt = 45,
tl.cex = 0.9,
addCoef.col = "black",
number.cex = 0.7,
title = title,
mar = c(0, 0, 2, 0)
)
}
# Calculate correlation matrix from wide data
calculate_correlation <- function(df_wide, exclude_cols = "year") {
mat <- df_wide |>
select(-all_of(exclude_cols)) |>
as.matrix()
cor(mat, use = "pairwise.complete.obs")
}
# Summarize correlation statistics
summarize_correlation <- function(cor_matrix) {
lower_tri <- cor_matrix[lower.tri(cor_matrix)]
tibble(
mean = mean(lower_tri, na.rm = TRUE),
min = min(lower_tri, na.rm = TRUE),
max = max(lower_tri, na.rm = TRUE)
)
}
# Perform hierarchical clustering and return results
perform_clustering <- function(df_wide,
province_col = "province",
dist_method = "euclidean",
hclust_method = "ward.D2") {
mat <- df_wide |>
select(-all_of(province_col)) |>
as.matrix()
rownames(mat) <- df_wide[[province_col]]
dist_mat <- dist(mat, method = dist_method)
hc <- hclust(dist_mat, method = hclust_method)
list(
matrix = mat,
distance = dist_mat,
hclust = hc
)
}
# Plot dendrogram with optional cluster rectangles
plot_dendrogram <- function(hc, main_title, k = NULL, colors = NULL) {
plot(
hc,
main = main_title,
xlab = "Province",
ylab = "Height",
sub = ""
)
if (!is_null(k)) {
rect.hclust(hc, k = k, border = colors)
}
}
# Create cluster assignments for multiple k values
create_cluster_assignments <- function(hc, k_values) {
cluster_list <- map(k_values, ~cutree(hc, k = .x))
names(cluster_list) <- str_c("cluster_", k_values)
tibble(province = names(cluster_list[[1]])) |>
bind_cols(map_dfc(cluster_list, as.factor))
}
# Plot time series with cluster assignments
plot_timeseries_with_clusters <- function(df_asi, df_clusters,
cluster_col = "cluster",
title, subtitle) {
n_clusters <- n_distinct(df_clusters[[cluster_col]])
linetypes <- c("solid", "dashed", "dotted", "dotdash", "longdash")[1:n_clusters]
df_asi |>
left_join(df_clusters, by = "province") |>
ggplot(aes(x = year, y = data, color = province, linetype = .data[[cluster_col]])) +
geom_line(linewidth = 1) +
geom_point(size = 2) +
scale_linetype_manual(values = set_names(linetypes, as.character(1:n_clusters))) +
labs(
title = title,
subtitle = subtitle,
x = "Year",
y = "ASI",
color = "Province",
linetype = "Cluster"
) +
theme(legend.position = "bottom")
}
# Summarize clusters
summarize_clusters <- function(df_asi, df_clusters, cluster_col = "cluster") {
df_asi |>
left_join(df_clusters, by = "province") |>
group_by(.data[[cluster_col]]) |>
summarise(
provinces = paste(unique(province), collapse = ", "),
n_provinces = n_distinct(province),
mean_asi = round(mean(data, na.rm = TRUE), 1),
sd_asi = round(sd(data, na.rm = TRUE), 1),
.groups = "drop"
)
}
# Create threshold flags for multiple RP values
create_threshold_flags <- function(df_rp, rp_thresholds = c(3, 4, 5)) {
for (rp in rp_thresholds) {
flag_col <- str_c("flag_rp", rp)
df_rp <- df_rp |>
mutate({{flag_col}} := as.integer(rp >= !!rp))
}
df_rp
}
# Pivot threshold flags to wide format for all RP values
pivot_threshold_flags_wide <- function(df_flags, rp_thresholds = c(3, 4, 5)) {
rp_thresholds |>
map(~{
flag_col <- str_c("flag_rp", .x)
df_flags |>
select(year, province, all_of(flag_col)) |>
pivot_wider(names_from = province, values_from = all_of(flag_col)) |>
arrange(year)
}) |>
set_names(str_c("rp", rp_thresholds))
}
# Plot threshold heatmap
plot_threshold_heatmap <- function(df_flags, title_suffix = "") {
df_flags |>
select(year, province, starts_with("flag_")) |>
pivot_longer(
cols = starts_with("flag"),
names_to = "threshold",
values_to = "exceeded"
) |>
mutate(
threshold = case_when(
threshold == "flag_rp3" ~ "3-Year RP",
threshold == "flag_rp4" ~ "4-Year RP",
threshold == "flag_rp5" ~ "5-Year RP",
TRUE ~ threshold
),
threshold = factor(threshold, levels = c("3-Year RP", "4-Year RP", "5-Year RP"))
) |>
ggplot(aes(x = year, y = province, fill = factor(exceeded))) +
geom_tile(color = "white", linewidth = 0.5) +
scale_x_continuous(breaks=1984:2025)+
facet_wrap(~threshold, ncol = 1) +
scale_fill_manual(
values = c("0" = "gray90", "1" = "tomato"),
labels = c("0" = "No", "1" = "Yes"),
name = "Threshold\nExceeded"
) +
labs(
title = str_c("Return Period Threshold Exceedances", title_suffix),
subtitle = "End-of-Season ASI (May Dekad 3)",
x = "Year",
y = "Province"
) +
theme(
axis.text.x = element_text(angle = 45, hjust = 1),
panel.grid = element_blank()
)
}
# Print correlation summary for multiple thresholds
print_correlation_summary <- function(cor_matrices, title_suffix = "") {
cat(str_c("Correlation Summary by Threshold", title_suffix, ":\n\n"))
walk2(cor_matrices, names(cor_matrices), ~{
stats <- summarize_correlation(.x)
threshold_name <- str_replace(.y, "rp", "") |>
str_c("-Year RP Threshold")
cat(threshold_name, ":\n", sep = "")
cat(" Mean correlation:", round(stats$mean, 3), "\n")
cat(" Range:", round(stats$min, 3), "to", round(stats$max, 3), "\n\n")
})
}
