#####################################################-
### A3. Cluster-PCA management policies ###

#Prepare PCA and cluster over UE managament policies
#Outputs: 

#library
library(ggplot2); library(dplyr); library(tidyverse); library(ecospat); library(ade4); library(ggdendro)
library(factoextra); library(corrplot); library(cluster); library(ggrepel); library(car)


#Dades
manage_data <- read.csv("Brian data_scripts/EU_managementscore_NRP.csv",stringsAsFactor=F)
colnames(manage_data)[1] = "Country"

#Paisos
alldat_count_list = c("Austria","Belgium","Bulgaria","Cyprus","Czechia","Denmark","Estonia","Finland","France","Germany","Greece","Hungary","Ireland","Italy","Latvia","Lithuania","Netherlands","Poland","Portugal","Romania","Slovakia","Slovenia","Spain","Sweden","United Kingdom")


#### 0. Data clean ####
#Clean data.frame
pos=which(names(manage_data) %in% c("X4.2.11.Declassification.system","X4.3.10.Criterion.B10","X4.7.5.Criterion.F.5"))
manage_data=manage_data[,-pos]

colnames(manage_data)[3:14] = c("X4.2.1.IAS.definition.A1","X4.2.2.Scope.coverage.A2","X4.3.3..RA.procedure.A3","X4.2.4.RA.During.emergencies.A4","X4.2.5.Listing.IAS.A5",
                                "X4.2.6.Transition.measures.A6","X4.2.7.Horizon.Scanning.tool.A7","X4.2.8.Alert.list.A8","X4.2.9.Pathway.RA.A9","X4.2.10.Identify.pathways.A10",
                                "X4.2.11.Priorisation.tool.A11","X4.2.12.Information.system.A12")

manage_data=manage_data[manage_data$Country %in% alldat_count_list,]

#Normalitzar valors per columnes
manage_data=manage_data[order(manage_data$Country),]
pos=ncol(manage_data)

manage2=apply(manage_data[,-c(1:2,(pos-2):pos)],2,scale)
rownames(manage2) = manage_data$Country


#Per type
manage2.A = manage2[,c(1:12)]
manage2.B = manage2[,c(13:23)]
manage2.C = manage2[,c(24:29)]
manage2.D = manage2[,c(30:36)]
manage2.G = manage2[,c(45:49)]

manage2.list = list(A=manage2.A,B=manage2.B,C=manage2.C,D=manage2.D,G=manage2.G)



#### 1. PCA management ####

###### PCA All measures ####
#PCA Brian style
pc=prcomp(manage2,scale=T)
summary(pc)$importance

# Percentatge de contribució de cada variable a cada component
var_contrib_percent <- 100 * (pc$sdev^2) / sum(pc$sdev^2)
# Contribució acumulada
var_contrib_acum <- cumsum(var_contrib_percent)

#Eigen
eigen_data=as.data.frame(pc$x[,1:6])
eigen_data$Country = rownames(eigen_data)
eigen_data$PCsum = rowSums(pc$x[,1:12])

eigen_loading=lapply(1:6,function(i){
  x=pc$rotation[,i]
  return(x[rev(order(abs(x)))][1:5])
})

#Save
saveRDS(pc,"dades/Management/PCA_Europe.rds")
saveRDS(eigen_loading,"dades/Management/PCA_EigenLoading_Europe.rds")
#write.csv(eigen_data, "dades/Management/PCA_Countries_Europe.csv")


#Plots PCA
# Plots PCA
fviz_pca_ind(pc)
fviz_pca_var(pc)
fviz_screeplot(pc)
fviz_contrib(pc, choice="var")
fviz_contrib(pc, choice="ind")



###### PCA by criterions ####

for (i in 1:5){
  manage3 = manage2.list[[i]]
  
  pc.type=prcomp(manage3,scale=T)
  summary(pc.type)
  
  num_pc <- which(cumsum(pc.type$sdev^2) / sum(pc.type$sdev^2) > 0.8)[1]
  
  eigen_data.type=as.data.frame(pc.type$x[,1:num_pc])
  eigen_data.type$Country = rownames(eigen_data.type)
  eigen_data.type$PCsum = rowSums(pc.type$x[,1:num_pc])
  colnames(eigen_data.type)[1:num_pc] = paste0(colnames(eigen_data.type)[1:num_pc],"_type",names(manage2.list[i]))
  names(eigen_data.type)[num_pc+2] = paste0("PCsum_type",names(manage2.list[i]))
  
  #Save
  eigen_data = merge(eigen_data, eigen_data.type, by="Country")
  
  #Plots PCA
  # Plots PCA
  fviz_pca_ind(pc.type)+ ggtitle(paste("Individuals - pc.typeA: Management type", names(manage2.list[i])))
  fviz_pca_var(pc.type)+ ggtitle(paste("Variables PCA: Management type", names(manage2.list[i])))
  fviz_screeplot(pc.type)+ ggtitle(paste("Scree plot: Management type", names(manage2.list[i])))
  fviz_contrib(pc.type, choice="var")+ ggtitle(paste("Variables contribution: Management type", names(manage2.list[i])))
  fviz_contrib(pc.type, choice="ind")+ ggtitle(paste("Individuals contribution: Management type", names(manage2.list[i])))
}

#Save
write.csv(eigen_data, "dades/Management/PCA_Countries_Europe.csv")



#### 2. Cluster management ####

manage2.BCD = cbind(manage2.B, manage2.C, manage2.D)

###### Cluster All measures ####
#Correlation
res <- round(cor(manage2.BCD),2)
#Correlacions altes: valor absolut, -55 per les autocorrelacions i /2 ja que la taula està duplicada
(sum(res > abs(0.8))-55)/2
corrplot(res, method = "color")

manage2.BCD <- manage2.BCD[, !colnames(manage2.BCD) %in% c("X4.5.2.Criterion.D2")] #Selected from VIF and cor analysis


#Matriu distancies
distance = dist(manage2.BCD,method = "euclidean")

#Nombre de clusters
#Avarage silhouette, sum of squares and gap_stat
fviz_nbclust(manage2.BCD, kmeans, method = "silhouette")
fviz_nbclust(manage2.BCD, kmeans, method = "wss")
fviz_nbclust(manage2.BCD, kmeans, method = "gap_stat")
#2 groups as optimal
k=2

#Cluster
#k-means clustering
set.seed(123)
kc = kmeans(manage2.BCD,k); kc
#cluster jerarquic
hc1 <- hclust(distance, method = "complete" )
hc2<-hclust(distance, method="average") 
#cluster Partitioning Around Medoids
pamvshortset <- pam(distance,k, diss = FALSE)

#We chose PAM or hc1??
#Cluster groups
member = cutree(hc1,k)
table(member)
a = aggregate(manage2.BCD,list(member),mean); a

table(pamvshortset$clustering)
b = aggregate(manage2.BCD,list(pamvshortset$clustering),mean); b

#save
cluster.df = as.data.frame(pamvshortset$clustering)
cluster.df$Country = rownames(cluster.df)
colnames(cluster.df)[1] = "ClusterPAM_all"

kk=as.data.frame(cutree(hc1,k))
kk$Country = rownames(kk)
colnames(kk)[1] = "ClusterHC1_all"
cluster.df = merge(cluster.df, kk, by="Country")


#Plots
#hc1 - cluster
plot(hc1,hang=-1,main='Default from hclust: complete')
rect.hclust(hc1, k = k, border = 2)

#Clusplot
clusplot(pamvshortset, shade = FALSE,labels=2,col.clus="purple",col.p="darkorange",span=FALSE,main="Cluster Mapping",cex=1.2)

#Valor de les variables normalitzades segon el grup
barplot(t(a[-1]), beside = TRUE, col = rainbow(length(unique(member))))

kk= as.data.frame(t(b[-1]));kk$RowNames <- row.names(kk)
ggplot(kk %>% pivot_longer(cols = -RowNames, names_to = "Columns", values_to = "Values"), 
       aes(x = RowNames, y = Values, fill = Columns)) +
  geom_bar(stat = "identity", position = "dodge") +
  labs(x = "Management", y = "Value") +
  theme_minimal()+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

#Valor de les variables reals segon el grup
cluster_means <- left_join(manage_data, cluster.df, by="Country") %>%
  select(-c(Reactive.Score, Proactive.Score,TOTAL))%>%
  pivot_longer(cols = -c(Country, ClusterPAM_all, ClusterHC1_all, Member.State), names_to = "Variable", values_to = "Valor")%>%
  group_by(ClusterPAM_all, Variable) %>% # o ClusterHC1_all
  summarise(Media = mean(Valor))
ggplot(cluster_means, aes(x = Variable, y = Media, fill = factor(ClusterPAM_all))) + # o ClusterHC1_all
  geom_bar(stat = "identity", position = "dodge", width = 0.7) +
  labs(x = "Management", y = "Mean", title = "Management mean for each cluster group") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
  scale_fill_discrete(name = "Grup de Clúster")

# Puntuació de les variables en els clústers
var_score <- apply(manage2.BCD, 2, function(x) { tapply(x, cutree(hc1,k), mean) })
barplot(var_score, beside = TRUE, col = rainbow(k), main = "Variable Scores by Cluster", las = 2)

var_score <- apply(manage2.BCD, 2, function(x) { tapply(x, pamvshortset$clustering, mean) })
barplot(var_score, beside = TRUE, col = rainbow(k), main = "Variable Scores by Cluster", las = 2)







###### Cluster by criterions ####

k.vec = c(2,2,2,2,2)
manage_type=c("A. Definition of IAS","B. Prevention","C. Early warning","D. Control and restoration","G. Capacity building")

for (j in 1:5){
  
  manage3 = manage2.list[[j]]
  k = k.vec[j]
  
  #Correlation
  # res <- round(cor(manage3),2)
  # corrplot(res, method = "color")
  # #VIF
  # vif_model <- lm(as.formula(paste(colnames(manage3)[1]," ~ ", paste(colnames(manage3)[-1], collapse = " + "))),
  #                 data = as.data.frame(manage3))
  # vif(vif_model)
  # #VIF matrix
  # vif_values <- sapply(colnames(manage3), function(var) {
  #   formula <- as.formula(paste(var, "~ ."))
  #   model <- try(lm(formula, data = as.data.frame(manage3)), silent = TRUE)
  #   if(inherits(model, "try-error")) {
  #      return(NA)
  #    } else {
  #      return(vif(model))
  #    }
  #   })
  # rowMeans(vif_values, na.rm = TRUE)
  
  #Clean manage3
  manage3 <- manage3[, !colnames(manage3) %in% c("X4.5.2.Criterion.D2", "X4.2.12.Information.system.A12")] #Selected from VIF and cor analysis
  
  #Cluster building
  distance = dist(manage3,method = "euclidean")
  
  #Avarage silhouette
  # fviz_nbclust(manage3, kmeans, method = "silhouette")+ ggtitle(paste("Optimal number of clusters: Management type", names(manage2.list[j])))
  
  #Cluster
  #k-means clustering
  set.seed(123)
  # kc = kmeans(manage3,k)
  # hc1 <- hclust(distance, method = "complete" )
  pamvshortset <- pam(distance,k, diss = T) #diss=T since distance is a matrix of distances; if manage3, then diss=F
  
  #Score
  cluster_data <- as.data.frame(cbind(manage3, Cluster = pamvshortset$clustering))
  average_score <- aggregate(cluster_data[, -ncol(cluster_data)], by = list(Cluster = cluster_data$Cluster), FUN = mean)
  
  #PLOT PER MANUSCRIPT
  # Convertim la primera columna en rownames
  rownames(average_score) <- average_score[,1]
  average_score <- average_score[,-1]
  # Transformem el dataframe en un format adequat per a ggplot
  df <- as.data.frame(t(average_score))
  df$Cluster <- rownames(df)
  df_long <- tidyr::pivot_longer(df, cols = -Cluster, names_to = "Variable", values_to = "Average_Score")
  
  # Plot amb ggplot2
  ggplot(df_long, aes(x = Cluster, y = Average_Score, fill = Variable)) +
    geom_bar(stat = "identity", position = "dodge") +
    labs(title = manage_type[j], x = "", y = "Average Score") +
    theme_minimal() +
    theme(legend.position = "top", axis.text.x = element_text(angle = 45, hjust = 1))
    #theme(axis.text.x = element_blank(), legend.position = "none")
  
  #Plots
  # #hc1 - cluster
  # plot(hc1,hang=-1,main=paste('Default from hclust: complete: Type',names(manage2.list[j])))
  # rect.hclust(hc1, k = k, border = 2)
  #PAM
  #clusplot(pamvshortset, shade = FALSE,labels=2,col.clus="purple",col.p="darkorange",span=FALSE,main=paste("Cluster Mapping PAM: Type",names(manage2.list[j])),cex=1.2)
  # fviz_cluster(pamvshortset, data = manage3, geom = "point", show.clust.cent = FALSE, palette = c("coral3", "aquamarine4")) +
  #   labs(title = manage_type[j]) +
  #   geom_text_repel(aes(label = rownames(manage3)), size = 5, color="grey20") +
  #   theme_minimal(base_size = 16) +
  #   theme(legend.position = "none")
  #  
  #Dendrograma
  hc <- hclust(as.dist(distance), method = "ward.D2")
  hc$labels <- c("CY", "GR", "CZ", "FI", "LT", "PL", "SK", "EE", "LV", "DE", "IT", "SI", "AT", "BG", 
                 "RO", "NL", "PT", "IE", "FR", "DK", "SE", "HU", "ES", "BE", "GB")
  dendro <- as.dendrogram(hc)
  # Plotar el dendrograma
  plot(dendro, main = manage_type[j], xlab = "", ylab = "")
  
  
  # #kmeans
  # clusplot(manage3, kc$cluster, color = TRUE, shade = TRUE, labels = 2, col.p = "darkblue",main = paste("Cluster Plot for k-means: Type",names(manage2.list[j])))
   
  #save PAM
  cluster.df.type = as.data.frame(pamvshortset$clustering)
  cluster.df.type$Country = rownames(cluster.df.type)
  colnames(cluster.df.type)[1] = paste0("Cluster_",names(manage2.list[j]))
  
  cluster.df = merge(cluster.df, cluster.df.type, by="Country")
  
  
}


#Change to letters
cluster.df2 <- cluster.df %>%
  mutate_if(is.numeric, 
            ~ case_when(
              . == 1 ~ "B",
              . == 2 ~ "A",
              . == 3 ~ "C",
              . == 4 ~ "D",
              . == 5 ~ "E",
              TRUE ~ as.character(.)
            )
  )
cluster.df2$Cluster_C = ifelse(cluster.df2$Cluster_C=="B","A","B")

#Save
write.csv(cluster.df2, "dades/Management/Cluster_Countries_Europe.csv")




#### 3. Years of policies (SONIGO) ####
library(readxl)

file_path="C:/Users/Quim/OneDrive - Universitat de Girona/AlienScenarios/Brian paper/Info seleccio/Data/BIO_ISPolicies_Table_16092011.xls"
sheet_names <- excel_sheets("C:/Users/Quim/OneDrive - Universitat de Girona/AlienScenarios/Brian paper/Info seleccio/Data/BIO_ISPolicies_Table_16092011.xls")
# Inicialitzar una llista per guardar els resultats
data <- list()

# Processar cada pestanya (excepte les últimes 7)
for (sheet in sheet_names[1:(length(sheet_names) - 7)]) {
  # Llegir la pestanya
  sheet_data <- read_excel(file_path, sheet = sheet)
  
  # Extreure la fila 2 (en realitat és la 3a en R)
  if (nrow(sheet_data) >= 2) {
    row_1 <- sheet_data[1, ]
    # Cercar anys dins del text
    years <- unlist(regmatches(as.character(row_1), gregexpr("\\d{4}", as.character(row_1))))
    # Afegir els resultats a la llista
    if (length(years) > 0) {
      for (year in years) {
        data <- rbind(data, data.frame(Country = sheet, Year = year))
      }
    }
  }
}

# Convertir a DataFrame
result_df <- as.data.frame(data)
result_df$Year <- as.numeric(result_df$Year)
result_df2 <- result_df[!is.na(result_df$Year) & result_df$Year > 1950, ]
result_df2 <- result_df2[!is.na(result_df2$Year) & result_df2$Year < 2017, ]

hist(result_df2$Year)
summary(result_df2$Year)