###############################################################################
## 
## COMISIÓN DE REGULACION DE COMUNICACIONES - CRC
##
## ANÁLISIS DE COMPONENTES PRINCIPALES Y DE CLÚSTER DE MUNICIPIOS PARA SERVICIOS FIJOS
## Elaboró: Laura Hernandez - Isabella Russi
## Fecha de Revisión : 28/02/2022
## 
################################################################################

rm(list=ls())       

# Librerias ----

library(tidyverse)  
library(caret)      
library(zeallot)    
library(tidyr)
library(dplyr)
library(readr)
library(readxl)
library(stats)      
library(FactoMineR) 
library(factoextra)
library(corrplot)
library(ggplot2)
library(gridExtra)
library(cluster)
library(clValid)
library(NbClust)
library(fpc)
library(clustertend)
library(dbscan)
library(Rcpp)
library(corrr)
library(skimr)
library(huxtable)
library(psych)  
library(GPArotation)
library(gg)


#####################################
#     Carga de datos                #
#####################################

Municipal <- read_delim("DATOS/Municipal_data.txt","|",
                        escape_double = FALSE, trim_ws = TRUE)

datos <- Municipal %>%
  select(-c("DP","DPNOM","MPIO"))%>%
  column_to_rownames(var="id_municipio")

# Descriptivas 

skim(datos)

# Correlaciones por pares

correlaciones <- datos %>% select_if(is.numeric) %>%
  correlate(method = "pearson") %>% 
  stretch(remove.dups = TRUE,na.rm = T)
corr_full <- correlaciones%>% mutate(r_abs = abs(r)) %>% arrange(desc(r_abs))
corr_full

# Prueba Barlett ----
matriz.cor<-cor(datos)
print(cortest.bartlett(matriz.cor, n=nrow(datos))) 

 # Este test Examina si la matriz de correlaciones es significativamente diferente
 # de una matriz de identidad (identity matrix). Si es significativo (p<0.05) = 
 # las correlaciones entre las variables son (en general) significativamente 
 # diferentes de cero. Eso pasa en estos datos, sin embargo el indicador se puede
 # haber visto afectado por el alto número de observaciones.

# Prueba KMO ----

      KMO(matriz.cor)
  
  # El test KMO (Kaiser, Meyer y Olkin) relaciona los coeficientes de correlación, 
  # rjh, observados entre las variables Xj y Xh, y ajh son los coeficientes de 
  # correlación parcial entre las variables Xj y Xh. Cuanto más cerca de 1 tenga 
  # el valor obtenido del test KMO, implica que la relación entres las variables 
  # es alta. 
 
## Ajuste de la base

datos <- Municipal %>%
  select(-c("DP","DPNOM","MPIO","Vel_sub_prom","desvest_s",
            "penetracion_telefonia_2020","Acc_17_20","desvest_B",
            "o_acto_terror_2019"))%>%
  column_to_rownames(var="id_municipio") 

matriz.cor<-cor(datos)
KMO(matriz.cor)


###########################################
# Analisis de componentes principales ----
###########################################

# Se calculan los CP. 
datos <- datos %>%
  rename("Oper_Tv"="op_tv_17_20",
         "Oper_Telef"="op_tel_17_20",
         "Penetr_Tv"="penetracion_tv_2019",
         "Empaquet"="porc_suscriptores_empaquetados",
         "HHI_Internet"="hhi_accesos",
         "HHI_TV"="hhi_susc_tv",
         "HHI_Telef"="hhi_lineas",
         "Suscrip_TV"="TV_17_20",
         "Suscrip_Telef"="Lineas_17_20",
         "Ruralidad"="Part_pobrural_18_21",
         "Dens_pob"="Dens_pob_21",
         "Valor_agr"="Valor_agr_percap_ctes2018_2016-2019",
         "Ing_mpal"="Ing_totpercap_imput",
         "IPM"="IPM",
         "Dist_capital"="DIST_CAPITAL_imput",
         "Amenaza_acc"="Proporcion_amenaza",
         "Penetr_Internet"="Acc_hogares",
         "Tec_inal"="Tec_inal_imput",
         "Velbajada"="Vel_baja_prom",
         "Oper_Internet"="op_int_17_20")


acp <- prcomp(datos, center = TRUE, scale = TRUE)


dim(acp$rotation) 
summary(acp)



################################################
# Eigenvalores y varianzas ----
################################################

# Los eigenvalores miden la varianza explicada por cada componente principal


get_eigenvalue(acp)


fviz_eig(acp, addlabels = T, main="",
         xlab="Dimensiones",ylab="Porcentaje explicado")

#####################################################
# Resultados para variables ----
#####################################################

var<-get_pca_var(acp)

# 1. Calculo de correlacion ----

# Las coordenadas del calculo son las correlaciones entre las variables y el CP

var$coord
fviz_pca_var(acp,col.var="cos2", repel=T,axes = c(1,2))

# 2. Contribucion a las variables ----

var$contrib
grid.arrange(fviz_contrib(acp,choice = "var", axes=1, top = 10, title="Dimensión 1"),
             fviz_contrib(acp,choice = "var", axes=2, top = 10, title="Dimensión 2"),
             ncol = 2)
grid.arrange(fviz_contrib(acp,choice = "var", axes=3, top = 10, title="Dimensión 3"),
             fviz_contrib(acp,choice = "var", axes=4, top = 10, title="Dimensión 4"),
             ncol = 2)
grid.arrange(fviz_contrib(acp,choice = "var", axes=5, top = 10, title="Dimensión 5"),
             ncol = 2)
             

fviz_contrib(acp, choice = "var", axes = c(1:5))


# 3. Calidad y contribucion ----

fviz_pca_ind(acp,col.var="cos2", repel=T,axes = c(1,2))

fviz_contrib(acp, choice = "ind", axes = 1, repel=T)

########################################################################
## Rotación PCA ----
########################################################################


acprot2 <- principal(datos, nfactors = 5 ,rotate = "varimax") 
acprot2$loadings


datos_scale<-scale(datos)
datos_rot <- as.data.frame((datos_scale) %*% (acprot2$loadings))

# Correlaciones por pares
correlaciones <- datos_rot %>% select_if(is.numeric) %>%
  correlate(method = "pearson") %>% 
  stretch(remove.dups = TRUE,na.rm = T)
corr_full <- correlaciones%>% mutate(r_abs = abs(r)) %>% arrange(desc(r_abs))
corr_full

# Prueba Barlett ----
matriz.cor<-cor(datos_rot)
print(cortest.bartlett(matriz.cor, n=nrow(datos_rot))) 

# Prueba KMO ----

KMO(matriz.cor)


########################################################################
## Clusters ----
##########################################################################

cp <- datos_rot[,1:3]

# Plot de correlaciones de los CP 

plot(datos_rot, pch=16, col=rgb(0,0,0,0.5)) 


# 1. Prueba de Hopkins: validación de tendencia de clustering ----
set.seed(123)
hopkins(cp,nrow(cp)-1)

 # Valores cercanos a cero (y alejados de 0.5) muestran evidencias de que 
 # las observaciones del set de datos no siguen una distribución espacial 
 # uniforme, su estructura, contiene algún tipo de agrupación.Por lo tanto se    
 # justifica hacer clusters.

# 2. Eleccion algoritmo ----
 set.seed(123) 
  comparacion1 <- clValid(
    obj        = cp,
    nClust     = 3:10,
    clMethods  = c("pam","kmeans", "hierarchical"),
    validation = c("stability", "internal")
  )
  y
  summary(comparacion1)
  

  
  ######################################################
  #  Cluster por kmeans ----
  set.seed(123)
  fviz_nbclust(x = cp, FUNcluster = kmeans, method = "wss", k.max = 15,
               diss = dist(cp, method = "euclidean")) 
  fviz_nbclust(x = cp, FUNcluster = kmeans, method = "silhouette", k.max = 15,
               diss = dist(cp, method = "euclidean"))  
  fviz_nbclust(x = cp, FUNcluster = kmeans, method = "gap_stat", k.max = 7,
               diss = dist(cp, method = "euclidean",nboot = 50))  
  
  
  # Se estima modelos con 5 clusters
  
  set.seed(123)
  
  km_clusters_5 <- kmeans(cp, centers = 5, nstart = 50)


  # Pruebas de validación
  
  # k=5
  km_clusters_5_p <- eclust(cp, FUNcluster = "kmeans", k = 5, seed = 123,
                            nstart = 50, graph = T)
  fviz_silhouette(sil.obj = km_clusters_5_p, print.summary = TRUE, palette = "jco",
                  ggtheme = theme_classic())
  km_clusters_5_p$silinfo$clus.avg.widths
  
  km_indices5 <- cluster.stats(d = dist(cp, method = "manhattan"), 
                               clustering = km_clusters_5_p$cluster)
  km_indices5$dunn

  
  # Exportamos resultados
  
  results <- as.data.frame(cbind(results,km5=km_clusters_5_p$cluster))
  results$id_municipio <- rownames(results)
  quick_xlsx(results,file="results.xlsx")