##k-means clustering----
#using all data
clustering_data_all <- scaled_full_data %>%
  select(-c(Country, Season, season_start, European_performance, lead_european_performance))
#considering optimal clusters
n <- 10
wss <- numeric(n)
set.seed(123)
for (i in 1:n) {
  km.out <- kmeans(clustering_data_all, centers = i, nstart = 20)
  wss[i] <- km.out$tot.withinss
}
wss_df <- tibble(clusters = 1:n, wss = wss)
ggplot(wss_df, aes(x = clusters, y = wss, group = 1)) +
  geom_point(size = 4)+
  geom_line() +
  scale_x_continuous(breaks = c(2, 4, 6, 8, 10)) +
  xlab('Number of clusters')
#perform clustering
set.seed(123)
kmeans_model_all <- kmeans(clustering_data_all, centers = 3, nstart = 50)
#joining to data and considering characteristics of clusters
clustered_data_all <- pivoted_scaled_full_data_comp_breakdown %>%
  mutate(cluster = as.factor(kmeans_model_all$cluster))
cluster_characteristics_all <- clustered_data_all %>%
  group_by(cluster) %>%
  summarise(across(where(is.numeric), mean, na.rm = TRUE)) %>%
  select(-season_start)


#using principal components
clustering_data_pc <- pca_results %>%
  select(contains("PC"))
#considering optimal clusters
n <- 10
wss <- numeric(n)
set.seed(123)
for (i in 1:n) {
  km.out <- kmeans(clustering_data_pc, centers = i, nstart = 20)
  wss[i] <- km.out$tot.withinss
}
wss_df <- tibble(clusters = 1:n, wss = wss)
ggplot(wss_df, aes(x = clusters, y = wss, group = 1)) +
  geom_point(size = 4)+
  geom_line() +
  scale_x_continuous(breaks = c(2, 4, 6, 8, 10)) +
  xlab('Number of clusters')
#perform clustering
set.seed(123)
kmeans_model_pc <- kmeans(clustering_data_pc, centers = 3, nstart = 50)
#joining to data and considering characteristics of clusters
clustered_data_pc <- pivoted_scaled_full_data_comp_breakdown %>%
  mutate(cluster = as.factor(kmeans_model_pc$cluster))
cluster_characteristics_pc <- clustered_data_pc %>%
  group_by(cluster) %>%
  summarise(across(where(is.numeric), mean, na.rm = TRUE)) %>%
  select(-season_start)


#using competitiveness indices
clustering_data_comp_index <- scaled_full_data_competitiveness_index %>%
  select(contains("index")) %>%
  select(1:3)
#considering optimal clusters
n <- 10
wss <- numeric(n)
set.seed(123)
for (i in 1:n) {
  km.out <- kmeans(clustering_data_comp_index, centers = i, nstart = 20)
  wss[i] <- km.out$tot.withinss
}
wss_df <- tibble(clusters = 1:n, wss = wss)
ggplot(wss_df, aes(x = clusters, y = wss, group = 1)) +
  geom_point(size = 4)+
  geom_line() +
  scale_x_continuous(breaks = c(2, 4, 6, 8, 10)) +
  xlab('Number of clusters')
#perform clustering
set.seed(123)
kmeans_model_comp_index <- kmeans(clustering_data_comp_index, centers = 3, nstart = 50)
#joining to data and considering characteristics of clusters
clustered_data_comp_index <- pivoted_scaled_full_data_comp_breakdown_competitiveness_index %>%
  mutate(cluster = as.factor(kmeans_model_comp_index$cluster))
cluster_characteristics_comp_index <- clustered_data_comp_index %>%
  group_by(cluster) %>%
  summarise(across(where(is.numeric), mean, na.rm = TRUE)) %>%
  select(-season_start)
plotting_cluster_data <- cluster_characteristics_comp_index %>% select(-11:-19)


##Evaluating clusters----
compare_clusters <- function(data, km_model, name) {
  dist_matrix <- dist(data)
  sil <- silhouette(km_model$cluster, dist_matrix)
  stats <- cluster.stats(dist_matrix, km_model$cluster)
  
  tibble(
    model = name,
    avg_silhouette = mean(sil[, 3]),
    dunn_index = stats$dunn,
    within_ss = km_model$tot.withinss,
    between_ss = km_model$betweenss,
    total_ss = km_model$totss,
    variance_explained = km_model$betweenss / km_model$totss
  )
}

comparison <- bind_rows(
  compare_clusters(clustering_data_all, kmeans_model_all, "All variables"),
  compare_clusters(clustering_data_pc, kmeans_model_pc, "Principal components"),
  compare_clusters(clustering_data_comp_index, kmeans_model_comp_index, "Theory-based indices")
) %>%
  arrange(desc(avg_silhouette), desc(dunn_index), desc(variance_explained))


##evaluating 'best' (indices) model
#full European performance
anova_model <- aov(European_performance ~ cluster, data = clustered_data_comp_index)
summary(anova_model)

kruskal.test(European_performance ~ cluster, data = clustered_data_comp_index)

panel_clusters <- plm(
  European_performance ~ cluster,
  data = clustered_data_comp_index,
  index = c("Country","season_start"),
  effect = "twoways",
  model = "within"
)
summary(panel_clusters)
coeftest(panel_clusters, vcov = vcovHC, type = "HC1")

#lead European performance
anova_model <- aov(lead_european_performance ~ cluster, data = clustered_data_comp_index)
summary(anova_model)

kruskal.test(lead_european_performance ~ cluster, data = clustered_data_comp_index)

panel_clusters_lead <- plm(
  lead_european_performance ~ cluster,
  data = clustered_data_comp_index,
  index = c("Country","season_start"),
  effect = "twoways",
  model = "within"
)
summary(panel_clusters_lead)
coeftest(panel_clusters_lead, vcov = vcovHC, type = "HC1")

#CL performance
anova_model <- aov(European_performance_CL ~ cluster, data = clustered_data_comp_index)
summary(anova_model)

kruskal.test(European_performance_CL ~ cluster, data = clustered_data_comp_index)

panel_clusters_CL <- plm(
  European_performance_CL ~ cluster,
  data = clustered_data_comp_index,
  index = c("Country","season_start"),
  effect = "twoways",
  model = "within"
)
summary(panel_clusters_CL)
coeftest(panel_clusters_CL, vcov = vcovHC, type = "HC1")

#EL performance
anova_model <- aov(European_performance_EL ~ cluster, data = clustered_data_comp_index)
summary(anova_model)

kruskal.test(European_performance_EL ~ cluster, data = clustered_data_comp_index)

panel_clusters_EL <- plm(
  European_performance_EL ~ cluster,
  data = clustered_data_comp_index,
  index = c("Country","season_start"),
  effect = "twoways",
  model = "within"
)
summary(panel_clusters_EL)
coeftest(panel_clusters_EL, vcov = vcovHC, type = "HC1")


##finding seasons closest to cluster centre
closest_n_per_cluster <- clustered_data_comp_index %>%
  group_by(cluster) %>%
  group_modify(~ {
    cluster_id <- as.character(.y$cluster)
    center_coords <- kmeans_model_comp_index$centers[cluster_id, , drop = FALSE]
    
    matrix_data <- .x %>% select(contains("index")) %>% as.matrix()
    
    center_matrix <- matrix(center_coords, nrow = nrow(matrix_data), ncol = ncol(matrix_data), byrow = TRUE)
    
    .x %>% mutate(
      Distance = sqrt(rowSums((matrix_data - center_matrix)^2))
    )
  }) %>%
  ungroup() %>%
  group_by(cluster) %>%
  slice_min(order_by = Distance, n = 10, with_ties = FALSE) %>%
  relocate(c(cluster, Season, Country, European_performance))