##correlation coefficient----
cor_matrix <- season_league_full_dataset %>%
  select(-c(Country, Season, season_start)) %>%
  cor(use = "pairwise.complete.obs")

plot_ly(
  x = colnames(cor_matrix),
  y = rownames(cor_matrix),
  z = cor_matrix,
  type = "heatmap",
  colorscale = "RdBu",
  zmin = -1,
  zmax = 1,
  showscale = TRUE
) %>%
  add_trace(
    x = rep(colnames(cor_matrix), each = nrow(cor_matrix)),
    y = rep(rownames(cor_matrix), times = ncol(cor_matrix)),
    text = sprintf("%.2f", as.vector(cor_matrix)),
    type = "scatter",
    mode = "text"      ,
    textfont = list(color = "black"),
    hoverinfo = "none",
    showlegend = FALSE
  ) %>%
  layout(
    xaxis = list(tickangle = 45),
    yaxis = list(autorange = "reversed")
  )

#drop multicollinearity variables
season_league_full_dataset_no_multicollinearity <- season_league_full_dataset %>%
  select(-c(Winners_GD, First_and_relegated))


##competitiveness index----
scaled_full_data_competitiveness_index <- scaled_full_data %>%
  mutate(competitive_title_index =
           rowMeans(cbind(-First_and_second, -First_and_CL)),
         elite_dominance_index =
           rowMeans(cbind(CL_and_relegated, Top_4_total_GD, Gini_coefficient,
                          Average_goals_by_top_3_players)),
         competitive_balance_index =
           rowMeans(cbind(-First_and_second, -First_and_CL, -CL_and_relegated,
                          -Top_4_total_GD, -Gini_coefficient)),
         player_rotation_index = Average_number_of_players,
         multilevel_competitiveness_index =
           rowMeans(cbind(-First_and_second, -First_and_CL, -CL_and_relegated))
  )
#by competition breakdown
scaled_full_data_comp_breakdown_competitiveness_index <- scaled_full_data_comp_breakdown %>%
  mutate(competitive_title_index =
           rowMeans(cbind(-First_and_second, -First_and_CL)),
         elite_dominance_index =
           rowMeans(cbind(CL_and_relegated, Top_4_total_GD, Gini_coefficient,
                          Average_goals_by_top_3_players)),
         competitive_balance_index =
           rowMeans(cbind(-First_and_second, -First_and_CL, -CL_and_relegated,
                          -Top_4_total_GD, -Gini_coefficient)),
         player_rotation_index = Average_number_of_players,
         multilevel_competitiveness_index =
           rowMeans(cbind(-First_and_second, -First_and_CL, -CL_and_relegated))
  )
#for clustering
pivoted_scaled_full_data_comp_breakdown_competitiveness_index <- scaled_full_data_comp_breakdown_competitiveness_index %>%
  pivot_wider(names_from = competition, values_from = c(European_performance, lead_european_performance)) %>%
  mutate(European_performance = 2*European_performance_CL + European_performance_EL,
         lead_european_performance = 2*lead_european_performance_CL + lead_european_performance_EL)

##principal component analysis----
pca_data <- season_league_full_dataset %>%
  select(-c(Season, Country, season_start, European_performance, lead_european_performance))
pca <- prcomp(pca_data, center = TRUE, scale. = TRUE)
#weights of each variable in the principal components
loadings <- as.data.frame(pca$rotation) %>%
  rownames_to_column("Variable")
#calculating proportion of variance explained by each component
pca.var <- pca$sdev^2
propve <- pca.var / sum(pca.var)
plot(cumsum(propve), xlab = "Principal Component",
     ylab = "Cumulative Proportion of Variance Explained",
     ylim = c(0, 1), type = "b")
abline(h=0.9, lty = "dashed")
which(cumsum(propve) >= 0.9)[1] #number of principal components that explain at least 90% of variance
#values of the observations on each principal component
pc_scores <- as.data.frame(pca$x[, 1:5])
pc_scores["PC1"] <- pc_scores["PC1"] * -1
pca_results <- bind_cols(
  scaled_full_data,
  pc_scores
)
#by competition breakdown
comp_breakdown_pc_scores <- pc_scores %>%
  slice(rep(1:n(), each = 2))
pca_results_comp_breakdown <- bind_cols(
  scaled_full_data_comp_breakdown,
  comp_breakdown_pc_scores
)
#reformatting data for plotting composition of components
loadings_plotting <- loadings %>%
  select(Variable, PC1:PC5) %>%
  pivot_longer(
    cols = starts_with("PC"),
    names_to = "PC",
    values_to = "Loading"
  ) %>%
  mutate(Loading = ifelse(PC == "PC1", Loading*-1, Loading))
#plotting
ggplot(loadings_plotting %>% filter(PC == "PC4"),
       aes(x = reorder(Variable, abs(Loading)),
           y = Loading,
           fill = Loading > 0)) +
  geom_col() +
  coord_flip() +
  scale_fill_manual(values = c("TRUE" = "steelblue",
                               "FALSE" = "firebrick")) +
  labs(
    x = NULL,
    y = "Loading",
    title = "PCA Loadings for Principal Components explaining 90% of variance"
  ) +
  theme_minimal() +
  theme(
    legend.position = "none",
    strip.text = element_text(face = "bold")
  )