############################################################################################################# ############################################################################################################# ############################################################################################################# ############################################################################################################# ############################# ############################# ############################# Comenzamos ############################# ############################# ############################# ############################################################################################################# ############################################################################################################# ############################################################################################################# ############################################################################################################# ############################################################################################################# # Mundo No Supervisado versus Mundo Supervisado. # Llamamos #install.packages("readxl") # Instalamos paquete library(readxl) # Corremos libreria Master <- read_excel("E:/PROYECTOS/DOCENCIA_/INAP/Versión3_/Clustering_/Data_set_a_Utilizar_/Indicadores_Econo_y_Salud_.xlsx", sheet = 1) # La funcion o comando header= encabezado, al colocarle TRUE, significa que la data lleva una etiqueta o nombre de las variables. dim(Master) # Revisamos Dimensiones str(Master) # Revisamos Clase de las Variables variable.names(Master) CP=143*13 CP # Respaldamos, Convertimos en tipo datframe, Indexamos, Convertimos todo en clase numeric y Eliminamos Vectores. Master_respaldo1 <- Master Master <- as.data.frame(Master) # Convierto Matriz a tipo dataframe dado que si continuo con tipo tabble, podr?an presentarse problemas. str(Master) # Consulto conversi?n y clase. View(Master) rownames(Master) <- Master$ID #Indexamos View(Master) # Consultamos Index str(Master) # Revisamos Clases # Eliminamos Variables View(Master) # Consultamos Index variable.names(Master) Master <- Master[ ,!colnames(Master)=="Pais"] Master <- Master[ ,!colnames(Master)=="ID"] variable.names(Master) # Revisamos str(Master) # Revisamos Clases View(Master) # Consultamos Index #View(Master$`Tasa_Mortalidad_Muertes/1000_habitantes`) ################# ################# ################# Escalamos ################# ################# ################# Master_respaldo2 <- Master # Respaldamos # Master <- Master_respaldo2 variable.names(Master) str(Master) Master <- as.data.frame(scale(Master)) # Escalamos. Explicar el escalamiento (funci?n scale y por que no otra) View(Master) # Revisamos str(Master) ####################################################################################################################### ####################################################################################################################### ############################## ##################################### ############################## Estudiamos las Correlaciones ##################################### ############################## ##################################### ####################################################################################################################### ####################################################################################################################### # Corresponder?a un estudio de correlaciones si mi objetivo es determinar cuales son las variables que est?n explicando la # variable dependiente. Lo que no es objeto de este estudio. Ya que es un proceso No Supervisado. # No obstante si quiero clusterizar sobrevalorando el peso de la informaci?n de ciertas variables para un modelo no supervisado, dado que # la realidad que envuelve a un fen?meno determinado presenta variables que entregan la misma informaci?n (la repite), puedo evadir # el estudio de correlaciones y dejarlas, argumentando que el fen?meno en estudio se relaciona con estas variables y son imprescindibles # conservarlas por motivo de negocio. # Para argumentar lo anterior ser?a relevante volver a la etapa de an?lisis exploratorio y volver al equipo de analista y al equipo # de negocio y/o ?pice estrat?gico para discutir sobre la ganancia y perdida de la inclusi?nd e estas variables, # seleccionando as? estas variables (y sin duda, tener claro tambi?n por que se eliminaron otras). Esta discusi?n entre el analista de datos # y el conocedor del negocio # es fundamental. (explicar con el ejemplo de cuando jug? Chile versus la contaminaci?n ambiental, aqu?l d?a se dispararon todos los # ?ndices de contaminaci?n, influyendo en el hacinamiento de centros de urgencia) ya que muchas veces la realidad, s? # le da importancia a una determinada variable repitiendo incluso la informaci?n de otra de ellas. # Por el contrario, si el objetivo de mi estudio es "depurar", "limpiar", o buscar la parsimonia, para as? determinar cuales variables que # est?n influyendo de mayor manera en la variable dependiente o simplemente "como est?n determinando o sesgando el modelo de datos" para # fines exploratorios, descriotivos 8no supervisados), elimino aquellas variables que est?n correlacionadas # un buen criterio es correlaci?n fuerte (igual o superior a 0.8), dejando la que tiene mayor poder # predictivo, y en segundo termino, la que tiene mayor informaci?n. Para esto la funci?n "cor" y la Matriz "CorrMat" es la que trabajaremos. # El poder predictivo se puede calcular con un information value de la libreria woe en r. # En definitiva, es relevante presentar dos estudios; uno manteniendo las variables correlacionadas y un segundo, sin variables correlacionadas. # Para este caso, haremos un estudio de correlaci?n y eliminar? aquellas que se encuentran correlacionadas por criterio de correlaci?n fuerte, # mayor a 80%. # Estudiamos Correlacion. Preparamos. #Redundancia (no se ve correlacion con las variables como estan) #install.packages("corrplot") library(corrplot) str(Master) #install.packages("dplyr") library(dplyr) sum(is.na(Master)) # Reviso nuevamente que no hayan NA # Corremos funci?n de correlaci?n variable.names(Master) str(Master) CorrMat <- cor(Master[,1:11]) #Correlacion View(CorrMat) corrplot.mixed(CorrMat) # Matriz de Correlación str(CorrMat) #Matriz de Correlaciones Matriz_de_Correlaciones_ <- data.frame(CorrMat) Matriz_de_Correlaciones_respaldo <- Matriz_de_Correlaciones_ # Respaldo Matriz_de_Correlaciones_ <- abs(Matriz_de_Correlaciones_) # Dejamos en Valores Absolutos str(Matriz_de_Correlaciones_) View(Matriz_de_Correlaciones_) # Crearemos un Dataframe a partir de la matriz de correlaci?n CorrMat para identificar rapidamente # las variables correlacionadas y su grado de correlaci?n y "con cuantas variables se correlaciona". # Reemplazamos Valores (todo lo que sea menor a 0.8 no me sirve; Ya que necesito identificar todas aquellas # variables que tienen correlación fuerte. Esas son las que necesito seleccionar para revisar cual elimino. Matriz_de_Correlaciones_[Matriz_de_Correlaciones_ < 0.8] <- NA str(Matriz_de_Correlaciones_) View(Matriz_de_Correlaciones_) # Aplicamos un ciclo FOR para eliminar la diagonal de la matriz, es decir las variables correlacionadas # consigo mismas (valor 1) # for = Esto significa Iteración o ciclo # i significa "1". Es un vector. Es un vector que comienza en el 1. Y "in 1:NROW", # me está diciendo que desde la fila 1 (ya que NROW significa eso) pase por toda la # fila de 128 filas de la Master Matriz de Correlaciones. El "[i,i]" es una cordenada # matricial. Significa lo que solemos entender como el "campo 1", es decir, el primer ciclo, # for, o iteración, comienza en la cordenada "1,1" (posición) es decir, que tome el valor #del campo intersectado por la fila 1 y la columna 1. Luego, la iteracion 2 siguira a la # cordenada "2,2". Luego, la iteración 3, sigue la cordenada "3,3", y asi sucecivamente. # Ojo, que todo es posible dado que sabemos que estámos en presenaic a de una matriz cuadrada (la # misma cantidad de filas es igual a la misma cantidad de columnas) ya que es una matriz de # correlación. Finalmente la salida es "N" que eso ya lo sabemos inicio <- Sys.time() # Con esto le indico que mida el inicio en tiempo del proceso for (i in 1:NROW(Matriz_de_Correlaciones_)) { Matriz_de_Correlaciones_[i,i] <- NA } fin <- Sys.time() # Con esto le indico que mida el fin en tiempo del proceso View(Matriz_de_Correlaciones_) # Contamos Frecuencia # Contar veces que un registro o etiqueta presenta informacion en sus siguientes variables Matriz_de_Correlaciones_$Fcorrelacion <- rowSums(!is.na(Matriz_de_Correlaciones_)) View(Matriz_de_Correlaciones_) variable.names(Matriz_de_Correlaciones_) Matriz_de_Correlaciones_prueba <- Matriz_de_Correlaciones_ Matriz_de_Correlaciones_prueba$Variables <- rownames(Matriz_de_Correlaciones_) # Des Indexar RUN ####### Hasta aquí llegamos el día 15 de diciembre de 2022, con los alumnos ##### # Eliminamos Variables Correlacionadas. Si no es por iv, ni informaci?n, aplicar criterio de negocio versus analistas # del preprocesamiento. Master_respaldo3 <- Master # Respaldamos # Master <- Master_respaldo3 Master <- Master[ ,!colnames(Master)=="Personas_vacunadas"] Master <- Master[ ,!colnames(Master)=="Muertes"] Master <- Master[ ,!colnames(Master)=="Tasa_mortalidad_en_porcentajes"] # Eliminamos esta variable a modo de simular una soplicitud # de negocio. Es tarea para la pr?xima clases, los alumnos investigar que tipo de correlaci?n utiliza la funci?n cor de R. str(Master) variable.names(Master) # Volvemos a Estudiar la Correlaci?n. sum(is.na(Master)) # Reviso nuevamente que no hayan NA # Corremos funci?nd de correlaci?n variable.names(Master) CorrMat2 <- cor(Master[,1:8]) #Correlación corrplot.mixed(CorrMat2) # Matriz de Correlación #Matriz de Correlaciones Matriz_de_Correlaciones_2 <- data.frame(CorrMat2) Matriz_de_Correlaciones_respaldo2 <- Matriz_de_Correlaciones_2 # Respaldo Matriz_de_Correlaciones_2 <- abs(Matriz_de_Correlaciones_2) # Dejamos en Valores Absolutos str(Matriz_de_Correlaciones_2) View(Matriz_de_Correlaciones_2) # Crearemos un Dataframe a partir de la matriz de correlaci?n CorrMat para identificar r?pidamente # las variables correlacionadas y su grado de correlaci?n y "con cuantas variables se correlaciona". # Reemplazamos Valores (todo lo que sea menor a 0.8 no me sirve; Ya que necesito identificar todas aquellas # variables que tienen correlaci?n fuerte. Esas son las que necesito seleccionar para revisar cual elimino. Matriz_de_Correlaciones_2[Matriz_de_Correlaciones_2 < 0.8] <- NA str(Matriz_de_Correlaciones_2) View(Matriz_de_Correlaciones_2) # Ya no existen correlaciones. # Aplicamos un ciclo FOR para eliminar la diagonal de la matriz, es decir las variables correlacionadas # consigo mismas (valor 1) # for = Esto significa Iteración o ciclo # i significa "1". Es un vector. Es un vector que comienza en el 1. Y "in 1:NROW", # me esta diciendo que desde la fila 1 (ya que NROW significa eso) pase por toda la # fila de 128 filas de la Master Matriz de Correlaciones. El "[i,i]" es una coordenada # matricial. Significa lo que solemos entender como el "campo 1", es decir, el primer ciclo, # for, o iteracion, comienza en la coordenada "1,1" (posicion) es decir, que tome el valor #del campo interceptando por la fila 1 y la columna 1. Luego, la iteracion 2 seguira a la # coordenada "2,2". Luego, la iteracion 3, sigue la coordenada "3,3", y asi sucecivamente. # Ojo, que todo es posible dado que sabemos que estamos en presencia de una matriz cuadrada (la # misma cantidad de filas es igual a la misma cantidad de columnas) ya que es una matriz de # correlacion. Finalmente la salida es "N" que eso ya lo sabemos inicio <- Sys.time() # Con esto le indico que mida el inicio en tiempo del proceso for (i in 1:NROW(Matriz_de_Correlaciones_2)) { Matriz_de_Correlaciones_2[i,i] <- NA } fin <- Sys.time() # Con esto le indico que mida el fin en tiempo del proceso View(Matriz_de_Correlaciones_2) # No existen correlaciones superiores a 0.8 Master_respaldo4 <- Master # Respaldamos variable.names(Master) View(Master) ########## Otros Plot ######### # instalar paquetes: PerformanceAnalytics, inspectdf, dplyr, corrplot, GGally, car, leaps, MPV, boot install.packages("PerformanceAnalytics") install.packages("inspectdf") install.packages("dplyr") install.packages("corrplot") install.packages("GGally") install.packages("car") install.packages("leaps") install.packages("MPV") install.packages("boot") library(inspectdf) library(dplyr) # calcula correlaciones seleccionando solo las variables numericas cor(crimeUS %>% select_if(is.numeric)) # entrega las correlaciones y su significancia crimeUS %>% inspect_cor() # cargamos en memoria el paquete corrplot library(corrplot) corrplot(cor(crimeUS),type="upper") # cargamos en memoria el paquete PerformanceAnalytics library(PerformanceAnalytics) chart.Correlation(crimeUS) # Otra nomesclatura cor(prueba1, method='pearson') cor(prueba1, method='kendall') cor(prueba1, method='spearman') ################# ################# PRIMER ACERCAMIENTO ################# ################# Plotearemos a trav?s de PCA el Modelo de Datos de Manera cuasi Natural ################# ################# antes de Probar cu?l es el mejor algoritmo de Clustering ################# ################# que Utilizaremos ################# ################# ################# # El objetivo de esto es encontrar una dispersi?n natural del Modelo de DAtos en el espacio cartesiano. # Con el objeto de Mirar si existen clasterS CUASI-naturales. # Reducimos en dos Variables a trav?s de la t?cnica de Componentes Principales. # Componentes Principales PCA = prcomp(Master) PCA # Reviso la Proporci?n Acumulada sobre la Varianza Total por PCA. Es decir, cuanto representa la aportaci?n de la varianza (acumulada) # con n PCA sobre la Varianza Total del Modelo de Datos. summary(PCA) # Tendremos la contribuci?n o participacion de la varianza de los pa?ses para la formaci?n de la componente # 1 y componente 2. # Esto lo haremos para ver si aparecen cuasi Clusters Naturales o No como ya se?alamos. # Traemos la contribuci?n u observaciones para la construcci?n del PCA por individuo. View(Master) l <- PCA$x View(l) # Seleccionamos las dos primeras PCA, al igual que lo anterior, pero por individuo y convertimos en Dataframe. variable.names(l) l <- as.data.frame(l[,-3:-8]) str(l) View(l) # Si nos fijamos vemos la aportaci?n de de los pa?ses, por cada pa?s, en la Varianza Total del Modelo de Datos. l$Estados <- rownames(l) # Desendixamos str(l) View(l) # Instalamos paquetes y corremos librer?as #install.packages("ggplot2") library(ggplot2) str(l) l$Estados <- as.factor(l$Estados) # Cambiamos a Tipo factor. # con colores. str(l) View(l) # Ploteamos con colores las dos primeras PCA. plot(x = l[,"PC1"], y = l[,"PC2"], col = l[,"Estados"], xlab = "PC1", ylab = "PCA2") rownames(Master_respaldo1) <- Master_respaldo1$Pais #Indexamos labels <- rownames(Master_respaldo1) text(l$PC1, l$PC2, labels, cex = 0.5, pos = 1) # Ploteamos sin colores las dos primeras PCA. # Requisto: Indexar la Variable de interes plot(x = l[,"PC1"], y = l[,"PC2"], xlab = "PC1", ylab = "PCA2") rownames(l) <- l$Estados #Indexamos labels <- rownames(l) text(l$PC1, l$PC2, labels, cex = 0.5, pos = 1) # En ID ingresamos el numero de interes para saber que país es. print((paisDeInteres <- subset(Master_respaldo1, ID==67))[1,1]) ################# ################# ################# Ploteamos 3d ################# ################# ################# # Traemos la contribuci?n u observaciones para la construcci?n del PCA por individuo. View(Master) l <- PCA$x View(l) # Seleccionamos las dos primeras PCA, al igual que lo anterior, pero por individuo y convertimos en Dataframe. variable.names(l) l <- as.data.frame(l[,-4:-8]) str(l) View(l) # Si nos fijamos vemos la aportaci?n de de los pa?ses, por cada pa?s, en la Varianza Total del Modelo de Datos. l$Estados <- rownames(l) # Desendixamos str(l) View(l) install.packages("rgl") library(rgl) #l$Estados <- rownames(l) # Des Indexar View(l) mycolors <- c('oldlace') l$Estados <- as.character(l$Estados) variable.names(l) str(l$Estados) dim(l) l <- cbind(l, color = rep("oldlace",143)) #Agrego filas str(l$color) View(l) # Plot 3D uno. plot3d( x=l$PC1, y=l$PC2, z=l$PC3, col = l$color, type = 's', text, radius = 0.2, #aspect = !add, # Encenderlo y se verá sin rectas. xlab="PCA1", ylab="PCA2", zlab="PCA3") # Plot 3D dos. install.packages("plotly") library(plotly) fig <- l %>% plot_ly(x = l$PC1, y = l$PC2, z = l$PC3, type = "scatter3d", mode = 'makers', marker = list(size = 6)) fig <- fig %>% layout(tittle = "3D Scatter plo") fig # Buscaremos el x o cordenada en el cubo dinámico del punto outliers hipotetico de interes, en el objeto l # para identificar su ID. View(l) # En ID ingresamos el numero de interes para saber que país es. print((paisDeInteres2 <- subset(Master_respaldo1, ID==67))[1,1]) # Guardamos Medio ambiente a Modo de Respaldo RData. #save.image("E:/PROYECTOS/DOCENCIA_/INAP/Versión4_/Clustering_/Data_set_a_Utilizar_/Environment_1.RData") # Respaldamos Medio Ambiente Rdata. Guardar en formato RDS. #load("E:/PROYECTOS/DOCENCIA_/INAP/Versión3_/Clustering_/Data_set_a_Utilizar_/Environment_1.RData") # Para cargar # Recordar aislar los outliers. Y luego plotear nuevamente. # Es recomedable presentar dos estudios; con los outliers aislados y sin los outliers aislados. ####################################################################################################################### ####################################################################################################################### ############################## ##################################### ############################## Evaluamos el mejor Algoritmo de Clustering ##################################### ############################## y tambien el optimo de K ##################################### ############################## ##################################### ####################################################################################################################### ####################################################################################################################### install.packages("clValid") library(clValid) dim(Master) # Revisamos Dimensiones de la Matriz variable.names(Master) str(Master) View(Master) Master_respaldo5 <- Master #Master <- Master_respaldo5 ############################# ############################# ############################# Probaremos el mejor algoritmod e Clustering ############################# ############################# para mi set de datos y el ?ptimo de K ############################# ############################# ############################# # Ahora inicio la prueba de 2 a 30 Clusters (K). Con esto ver? cu?l es el algoritmo mas apropiado. inicio <- Sys.time() # Con esto le indico que mida el inicio en tiempo del proceso comparacion <- clValid(obj = Master, nClust = 2:30, # Numero de Clusters a Probar clMethods = c("hierarchical", "kmeans", "pam"), validation = c("stability", "internal")) #y# Estos es por si nos pide el alto rendimiento o recurso que ocupar? de nuestro hardware o m?quina. Mejor entrenarlo en un servidor. summary(comparacion) fin <- Sys.time() # Con esto le indico que mida el fin en tiempo del proceso para evaluar cuanto demora sino # Clustering Methods: # hierarchical kmeans pam # # Cluster sizes: # 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 # # Validation Measures: # 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 # # hierarchical APN 0.0029 0.0096 0.0095 0.0121 0.0180 0.0349 0.0575 0.0600 0.0599 0.0639 0.1141 0.1874 0.1746 0.2959 0.3321 0.3333 0.2075 0.2396 0.2551 0.2604 0.2490 0.3121 0.2769 0.2835 0.2536 0.2683 0.2900 0.2962 0.2815 # AD 3.1809 3.1306 3.0167 2.9093 2.8189 2.7872 2.7106 2.6585 2.5757 2.5224 2.4949 2.4522 2.3649 2.3461 2.3007 2.2736 2.0487 2.0044 1.9380 1.9129 1.8773 1.8397 1.7657 1.7425 1.6556 1.6238 1.6008 1.5754 1.5452 # ADM 0.0738 0.1367 0.0999 0.1491 0.1403 0.2722 0.2572 0.2709 0.3091 0.2800 0.3676 0.5106 0.5221 0.7445 0.8360 0.8262 0.6713 0.6542 0.6850 0.7035 0.6816 0.7081 0.6787 0.6718 0.6832 0.6857 0.6925 0.7402 0.7242 # FOM 0.9787 0.9653 0.9622 0.9520 0.9465 0.9493 0.9412 0.9245 0.9266 0.8759 0.8778 0.8796 0.8815 0.8806 0.8808 0.8807 0.8724 0.8716 0.8727 0.8705 0.8687 0.8616 0.8043 0.8038 0.8052 0.8067 0.8053 0.8048 0.8077 # Connectivity 5.2869 6.7869 9.7159 14.0500 16.9790 21.1980 23.1980 28.2226 36.5579 39.4869 39.4869 45.0948 54.6813 55.1813 65.6734 67.6734 87.1310 87.1310 101.3798 104.6254 105.1183 105.3683 109.5357 110.5357 123.8694 125.3552 127.6052 136.2956 139.9813 # Dunn 0.4871 0.4677 0.6278 0.3490 0.4142 0.3748 0.3748 0.3748 0.2581 0.2581 0.2581 0.2581 0.1630 0.1630 0.1665 0.1665 0.1624 0.1624 0.1932 0.1932 0.1932 0.1932 0.1932 0.1932 0.1992 0.1992 0.1992 0.1992 0.1992 # Silhouette 0.6851 0.6665 0.6363 0.5900 0.5722 0.4763 0.4714 0.3260 0.3048 0.2437 0.2135 0.1617 0.1252 0.1183 0.1041 0.1016 0.1819 0.1665 0.1647 0.1531 0.1379 0.1303 0.1462 0.1397 0.1722 0.1660 0.1596 0.1631 0.1649 # kmeans APN 0.0114 0.0110 0.0131 0.0354 0.0626 0.1372 0.4609 0.3295 0.2860 0.3431 0.4193 0.4019 0.3766 0.3683 0.3587 0.3990 0.3987 0.3504 0.3832 0.3700 0.3496 0.3273 0.3263 0.3195 0.3201 0.2868 0.2724 0.3161 0.3234 # AD 3.1703 3.0936 2.9521 2.8787 2.7634 2.7308 2.7221 2.4122 2.2819 2.2485 2.2169 2.0929 2.0219 1.9506 1.9120 1.8786 1.8474 1.7829 1.7808 1.7396 1.6984 1.6418 1.6107 1.5709 1.5466 1.4964 1.4502 1.4533 1.4283 # ADM 0.0905 0.2206 0.1563 0.2073 0.2465 0.4215 1.0087 0.7706 0.7658 0.8477 0.9784 0.9032 0.8788 0.8424 0.8484 0.8470 0.8728 0.7881 0.8609 0.8555 0.7971 0.7824 0.7500 0.7133 0.7242 0.6865 0.6408 0.7010 0.6932 # FOM 0.9558 0.9741 0.9627 0.9479 0.9421 0.9449 0.9426 0.9396 0.9237 0.9097 0.9091 0.9124 0.9139 0.8641 0.8688 0.8679 0.8622 0.8639 0.8637 0.8600 0.8557 0.8564 0.8581 0.8532 0.7990 0.8011 0.8010 0.8040 0.8064 # Connectivity 5.2869 10.5976 13.5266 15.0266 24.6270 24.7270 26.7270 64.7607 90.8242 90.9429 99.9845 100.8345 114.2254 127.0992 122.5325 124.5325 118.6730 154.3385 146.1111 154.9587 173.9095 161.9198 163.6790 163.8079 154.0885 162.8139 165.5702 164.3702 165.1496 # Dunn 0.4871 0.1980 0.2681 0.2967 0.1563 0.1773 0.1773 0.0984 0.0976 0.1113 0.1149 0.1128 0.1248 0.1258 0.1643 0.1643 0.1550 0.1503 0.1412 0.1503 0.1232 0.1482 0.1482 0.1482 0.2322 0.1848 0.1876 0.1876 0.1835 # Silhouette 0.6851 0.5404 0.5453 0.5458 0.4298 0.3719 0.3670 0.2328 0.2119 0.2157 0.2042 0.1993 0.1929 0.1954 0.1956 0.1908 0.1968 0.1786 0.1838 0.1907 0.1751 0.1984 0.1999 0.2009 0.2050 0.2096 0.2064 0.2062 0.1933 # pam APN 0.1560 0.1473 0.2211 0.3688 0.3825 0.2662 0.3072 0.2658 0.3415 0.2957 0.2801 0.3379 0.3264 0.3526 0.3145 0.3664 0.3556 0.3657 0.3806 0.3271 0.3529 0.3788 0.3801 0.3919 0.3734 0.3831 0.3818 0.3597 0.3654 # AD 3.1448 3.0042 2.7910 2.7689 2.5959 2.4448 2.4369 2.2728 2.2321 2.1528 2.0671 1.9958 1.9187 1.9189 1.8528 1.8106 1.7719 1.7477 1.7343 1.6605 1.6551 1.6392 1.6168 1.5911 1.5529 1.5274 1.4987 1.4673 1.4465 # ADM 0.5002 0.5486 0.6220 0.9825 0.9231 0.7479 0.9731 0.7705 0.8600 0.7956 0.7198 0.7392 0.6983 0.8081 0.7316 0.7549 0.7176 0.7292 0.7603 0.6635 0.7054 0.7549 0.7570 0.7594 0.7261 0.7265 0.7135 0.6904 0.7025 # FOM 0.9984 0.9989 0.9605 0.9614 0.9448 0.9420 0.9424 0.9415 0.9427 0.9442 0.9444 0.9355 0.9339 0.9370 0.9398 0.9322 0.9301 0.9336 0.9355 0.9246 0.9252 0.8618 0.8651 0.8591 0.8592 0.8544 0.8546 0.8536 0.8535 # Connectivity 60.7095 84.0857 85.2468 87.2750 87.1750 83.0361 101.9333 110.0536 110.2536 113.1183 114.6183 115.4266 121.6313 134.3056 143.8952 145.8952 155.8357 162.0171 163.4933 156.0873 158.0873 167.3925 168.6643 170.2321 170.0262 172.2484 174.3317 174.8528 174.5456 # Dunn 0.0359 0.0359 0.0468 0.0449 0.0469 0.0374 0.0469 0.0469 0.0469 0.0469 0.0656 0.0976 0.0976 0.0976 0.1014 0.1284 0.1323 0.1323 0.1433 0.1309 0.1309 0.1429 0.1429 0.1429 0.1377 0.1377 0.1377 0.1470 0.1470 # Silhouette 0.1915 0.1254 0.1523 0.1257 0.1460 0.1760 0.1549 0.1732 0.1820 0.1852 0.1859 0.1986 0.1978 0.1894 0.1906 0.1861 0.1816 0.1808 0.1750 0.1726 0.1702 0.1760 0.1814 0.1835 0.1706 0.1704 0.1703 0.1699 0.1744 # # Optimal Scores: # # Score Method Clusters # APN 0.0029 hierarchical 2 # AD 1.4283 kmeans 30 # ADM 0.0738 hierarchical 2 # FOM 0.7990 kmeans 26 # Connectivity 5.2869 hierarchical 2 # Dunn 0.6278 hierarchical 4 # Silhouette 0.6851 hierarchical 2 # El APN debe ser lo mas cercano a 0. # El AD debe ser lo mas cercano a 0. # El ADM debe ser lo mas cercano a 0. # El FOM debe ser lo mas cercano a 0. # El CONNECTIVITY debe ser lo mas cercano a 0. # El DUNN debe ser lo mas alto posible. # El SILHOUETTE debe ser lo mas alto posible. # Preeliminarmente me quedo con 4 hierarchical (y para un segundo estudio, 26 Kmeans ya que sabemos # que es sensible a los outliers). # Hay que recordar que nunca se debe elegir como optimo el minimo de clusteres ni el # máximo de clústeres testeado o arrojado por la técnica de testeo del optimo. Dado que es lógico obtener # los mejores resultados desde el punto de vista de la maximización de la inercia entre-clústers con el # menor número de clusters. De la misma manera, es lógico también obtener la mayor minimización de la # inercia intra-clústers al aumentar al máximo el número de clústers ya que las # observaciones dentro de un clústers comenzaran a compactarse. Entonces los # parametros que miden por separado cada distancia siempre sugeriran en primer lugar, # con minimo y maximo de clusters. # Probamos con 20. # Ahora inicio la prueba de 2 a 20 Clusters (K). Con esto ver? cu?l es el algoritmo mas apropiado. inicio <- Sys.time() # Con esto le indico que mida el inicio en tiempo del proceso comparacion <- clValid(obj = Master, nClust = 2:20, # Numero de Clusters a Probar clMethods = c("hierarchical", "kmeans", "pam"), validation = c("stability", "internal")) #y# Estos es por si nos pide el alto rendimiento o recurso que ocupar? de nuestro hardware o m?quina. Mejor entrenarlo en un servidor. summary(comparacion) fin <- Sys.time() # Con esto le indico que mida el fin en tiempo del proceso para evaluar cuanto demora sino # Clustering Methods: # hierarchical kmeans pam # # Cluster sizes: # 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 # # Validation Measures: # 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 # # hierarchical APN 0.0029 0.0096 0.0095 0.0121 0.0180 0.0349 0.0575 0.0600 0.0599 0.0639 0.1141 0.1874 0.1746 0.2959 0.3321 0.3333 0.2075 0.2396 0.2551 # AD 3.1809 3.1306 3.0167 2.9093 2.8189 2.7872 2.7106 2.6585 2.5757 2.5224 2.4949 2.4522 2.3649 2.3461 2.3007 2.2736 2.0487 2.0044 1.9380 # ADM 0.0738 0.1367 0.0999 0.1491 0.1403 0.2722 0.2572 0.2709 0.3091 0.2800 0.3676 0.5106 0.5221 0.7445 0.8360 0.8262 0.6713 0.6542 0.6850 # FOM 0.9787 0.9653 0.9622 0.9520 0.9465 0.9493 0.9412 0.9245 0.9266 0.8759 0.8778 0.8796 0.8815 0.8806 0.8808 0.8807 0.8724 0.8716 0.8727 # Connectivity 5.2869 6.7869 9.7159 14.0500 16.9790 21.1980 23.1980 28.2226 36.5579 39.4869 39.4869 45.0948 54.6813 55.1813 65.6734 67.6734 87.1310 87.1310 101.3798 # Dunn 0.4871 0.4677 0.6278 0.3490 0.4142 0.3748 0.3748 0.3748 0.2581 0.2581 0.2581 0.2581 0.1630 0.1630 0.1665 0.1665 0.1624 0.1624 0.1932 # Silhouette 0.6851 0.6665 0.6363 0.5900 0.5722 0.4763 0.4714 0.3260 0.3048 0.2437 0.2135 0.1617 0.1252 0.1183 0.1041 0.1016 0.1819 0.1665 0.1647 # kmeans APN 0.0114 0.0110 0.0131 0.0354 0.0626 0.1372 0.4609 0.3295 0.2860 0.3431 0.4193 0.4019 0.3766 0.3683 0.3587 0.3990 0.3987 0.3504 0.3832 # AD 3.1703 3.0936 2.9521 2.8787 2.7634 2.7308 2.7221 2.4122 2.2819 2.2485 2.2169 2.0929 2.0219 1.9506 1.9120 1.8786 1.8474 1.7829 1.7808 # ADM 0.0905 0.2206 0.1563 0.2073 0.2465 0.4215 1.0087 0.7706 0.7658 0.8477 0.9784 0.9032 0.8788 0.8424 0.8484 0.8470 0.8728 0.7881 0.8609 # FOM 0.9558 0.9741 0.9627 0.9479 0.9421 0.9449 0.9426 0.9396 0.9237 0.9097 0.9091 0.9124 0.9139 0.8641 0.8688 0.8679 0.8622 0.8639 0.8637 # Connectivity 5.2869 10.5976 13.5266 15.0266 24.6270 24.7270 26.7270 64.7607 90.8242 90.9429 99.9845 100.8345 114.2254 127.0992 122.5325 124.5325 118.6730 154.3385 146.1111 # Dunn 0.4871 0.1980 0.2681 0.2967 0.1563 0.1773 0.1773 0.0984 0.0976 0.1113 0.1149 0.1128 0.1248 0.1258 0.1643 0.1643 0.1550 0.1503 0.1412 # Silhouette 0.6851 0.5404 0.5453 0.5458 0.4298 0.3719 0.3670 0.2328 0.2119 0.2157 0.2042 0.1993 0.1929 0.1954 0.1956 0.1908 0.1968 0.1786 0.1838 # pam APN 0.1560 0.1473 0.2211 0.3688 0.3825 0.2662 0.3072 0.2658 0.3415 0.2957 0.2801 0.3379 0.3264 0.3526 0.3145 0.3664 0.3556 0.3657 0.3806 # AD 3.1448 3.0042 2.7910 2.7689 2.5959 2.4448 2.4369 2.2728 2.2321 2.1528 2.0671 1.9958 1.9187 1.9189 1.8528 1.8106 1.7719 1.7477 1.7343 # ADM 0.5002 0.5486 0.6220 0.9825 0.9231 0.7479 0.9731 0.7705 0.8600 0.7956 0.7198 0.7392 0.6983 0.8081 0.7316 0.7549 0.7176 0.7292 0.7603 # FOM 0.9984 0.9989 0.9605 0.9614 0.9448 0.9420 0.9424 0.9415 0.9427 0.9442 0.9444 0.9355 0.9339 0.9370 0.9398 0.9322 0.9301 0.9336 0.9355 # Connectivity 60.7095 84.0857 85.2468 87.2750 87.1750 83.0361 101.9333 110.0536 110.2536 113.1183 114.6183 115.4266 121.6313 134.3056 143.8952 145.8952 155.8357 162.0171 163.4933 # Dunn 0.0359 0.0359 0.0468 0.0449 0.0469 0.0374 0.0469 0.0469 0.0469 0.0469 0.0656 0.0976 0.0976 0.0976 0.1014 0.1284 0.1323 0.1323 0.1433 # Silhouette 0.1915 0.1254 0.1523 0.1257 0.1460 0.1760 0.1549 0.1732 0.1820 0.1852 0.1859 0.1986 0.1978 0.1894 0.1906 0.1861 0.1816 0.1808 0.1750 # # Optimal Scores: # # Score Method Clusters # APN 0.0029 hierarchical 2 # AD 1.7343 pam 20 # ADM 0.0738 hierarchical 2 # FOM 0.8622 kmeans 18 # Connectivity 5.2869 hierarchical 2 # Dunn 0.6278 hierarchical 4 # Silhouette 0.6851 hierarchical 2 # Preeliminarmente, sigo seleccionando 4 hierarchical (y para un segundo estudio, no 26 sino 18 Kmeans. # Es interesante observar como a medido que aumento el numero de K Kmeans selecciona como optimo # un numero mayor. Esto se explica por la sensibilidad de la medida de distancia a los outliers). # Probaremos por tercera vez, pero con 7. # Ahora inicio la prueba de 2 a 7 Clusters (K). Con esto ver? cu?l es el algoritmo mas apropiado. inicio <- Sys.time() # Con esto le indico que mida el inicio en tiempo del proceso comparacion <- clValid(obj = Master, nClust = 2:7, # Numero de Clusters a Probar clMethods = c("hierarchical", "kmeans", "pam"), validation = c("stability", "internal")) #y# Estos es por si nos pide el alto rendimiento o recurso que ocupar? de nuestro hardware o m?quina. Mejor entrenarlo en un servidor. summary(comparacion) fin <- Sys.time() # Con esto le indico que mida el fin en tiempo del proceso para evaluar cuanto demora sino # Clustering Methods: # hierarchical kmeans pam # # Cluster sizes: # 2 3 4 5 6 7 # # Validation Measures: # 2 3 4 5 6 7 # # hierarchical APN 0.0029 0.0096 0.0095 0.0121 0.0180 0.0349 # AD 3.1809 3.1306 3.0167 2.9093 2.8189 2.7872 # ADM 0.0738 0.1367 0.0999 0.1491 0.1403 0.2722 # FOM 0.9787 0.9653 0.9622 0.9520 0.9465 0.9493 # Connectivity 5.2869 6.7869 9.7159 14.0500 16.9790 21.1980 # Dunn 0.4871 0.4677 0.6278 0.3490 0.4142 0.3748 # Silhouette 0.6851 0.6665 0.6363 0.5900 0.5722 0.4763 # kmeans APN 0.0114 0.0110 0.0131 0.0354 0.0626 0.1372 # AD 3.1703 3.0936 2.9521 2.8787 2.7634 2.7308 # ADM 0.0905 0.2206 0.1563 0.2073 0.2465 0.4215 # FOM 0.9558 0.9741 0.9627 0.9479 0.9421 0.9449 # Connectivity 5.2869 10.5976 13.5266 15.0266 24.6270 24.7270 # Dunn 0.4871 0.1980 0.2681 0.2967 0.1563 0.1773 # Silhouette 0.6851 0.5404 0.5453 0.5458 0.4298 0.3719 # pam APN 0.1560 0.1473 0.2211 0.3688 0.3825 0.2662 # AD 3.1448 3.0042 2.7910 2.7689 2.5959 2.4448 # ADM 0.5002 0.5486 0.6220 0.9825 0.9231 0.7479 # FOM 0.9984 0.9989 0.9605 0.9614 0.9448 0.9420 # Connectivity 60.7095 84.0857 85.2468 87.2750 87.1750 83.0361 # Dunn 0.0359 0.0359 0.0468 0.0449 0.0469 0.0374 # Silhouette 0.1915 0.1254 0.1523 0.1257 0.1460 0.1760 # # Optimal Scores: # # Score Method Clusters # APN 0.0029 hierarchical 2 # AD 2.4448 pam 7 # ADM 0.0738 hierarchical 2 # FOM 0.9420 pam 7 # Connectivity 5.2869 hierarchical 2 # Dunn 0.6278 hierarchical 4 # Silhouette 0.6851 hierarchical 2 # Preeliminarmente, continúo con mi selección, 4 hierarchical. # Y por el estudio anterior, mantengo en mi memoria 18 Kmeans o 26 (zona) # a medida que aumento el K. Esto es sintoma de sensibilidad a los outliers. # Ahora probaremos otras tecnicas complementarias para tomar una decision fundada o robusta. ################################################################################################################### ################################################################################################################### ############################## ##################################### ############################## Utilizaremos Otras Medidas de ##################################### ############################## elecci?n del ?ptimo de K ##################################### ############################## "Optimizando la suma de errores al cuadrado" ##################################### ############################## ##################################### ################################################################################################################### ################################################################################################################### # "Metodo" del Codo o Elbow Method a fin de determinar dicho valor. Basicamente, este m?todo busca seleccionar la # cantidad ideal de grupos a partir de la optimizacion de la WCSS (Within Clusters Summed Squares -optimizando la # suma errores al cuadrado). # La base del algoritmo k-means, buscar optimizar la suma errores al cuadrado. La suma de errores al cuadrado es # igual a la distancia euclidea entre las observaciones. # La idea del algoritmo es agrupar cada observacion (fila de una tabla) en k grupos, calcular el numero de clusters # optimo. Para ello, el algoritmo sigue los siguientes pasos: # 1. Se inicializan k puntos aleatorios, llamados centroides. # 2. Para cada observacion se calcula la suma de errores al cuadrado de esa observacion respecto a cada uno de los k # centroides. # 3. A cada observacion se le asigna el centroide que menos error tenga. # Mediremos el optimo de K con sslo el algoritmo K-Means, a partir del cual K la ganancia en la reduccion # en la suma total de varianza intra-cluster deja de ser sustancial. Buscamos que esta varianza sea la mas pequeña. install.packages("factoextra") library(factoextra) # Optimo de K bajo el Metodo Wss y algoritmo kmeans. # Mediremos la ganancia al aumentar en un numero el K o Numero de Clusters con 30 K. fviz_nbclust(x = Master, FUNcluster = kmeans, method = "wss", k.max = 30) + labs(title = "Numero optimo de clusters") # Se puede apreciar que la maxima ganancia en la disminucion de la varianza intra clusters es # cuando saltamos de 2 a 3 o 4 Clusters con kmeans. Este es un problema de pendiente. La pendiente # mas pronunciada se obtiene en el salto de 2 a 3 o 4. Esto habla con las tecnicas anteriores al obtener # 4 k con medoidso PAM. La segunda gran disminucion es cuando saltamos de 8 a 9 K. Luego pareciese observar # una ganancia marginal e inestable. # Guardamos Medio ambiente a Modo de Respaldo RData. #save.image("E:/PROYECTOS/DOCENCIA_/INAP/Versión4_/Clustering_/Data_set_a_Utilizar_/Environment_2.RData") # Respaldamos Medio Ambiente Rdata. Guardar en formato RDS. #load("E:/PROYECTOS/DOCENCIA_/INAP/Versión4_/Clustering_/Data_set_a_Utilizar_/Environment_2.RData") # Para cargar # Optimo de K bajo el Metodo Silhouette y algoritmo kmeans. # Probaremos otra medida de optimo con 30 K, el silhouette que mide el grado de confianza de la asigna-cion de una observacion, # con solo K means. Su valor puede estar entre -1 y 1, siendo valores altos un indicativo de que la observacion se ha # asignado al cluster correcto. fviz_nbclust(x = Master, FUNcluster = kmeans, method = "silhouette", k.max = 30) + labs(title = "Numero optimo de clusters") # Se puede apreciar que el grado de confianza da un salto sustantivo a partir de 4 a 7 Clusters. # El punto mas alto se obtiene en 17 como una de la stecnicas anteriores vistas. # Optimo de K bajo el Metodo wss y algoritmo pam (manhattan) # Probaremos otra medida del optimo de K pero con el algoritmo PAM, bajo el parametro wss # Una forma sencilla de estimar el numero K optimo de clusters cuando no se dispone de informacion adicional en la que # basarse es aplicar el algoritmo para un rango de valores de K, identificando aquel a partir del cual la # reduccion en la suma total de varianza intra-cluster deja de ser sustancial. # La funcion fviz_nbclust() automatiza este proceso. En este caso, dado que se sospecha de la presencia de outliers, # se emplea la distancia de Manhattan como medida de similitud. # Probaremos con 30 fviz_nbclust(x = Master, FUNcluster = pam, k.max = 30, method = "wss", diss = dist(Master, method = "manhattan")) # Utilizamos distancia Manhatan en este argumento. # Se aprecia una ganancia significativa cuando saltamos de 3 a 4 Clusters como lo habiamos detectado en las tecnicas # anteriores de manera numerica. Podteriormente pareciese estabilizarse marginalmente la ganancia. # Recomendaciones: 4 estudios # E1: 4K Medois con oUtliers # E2: 4K Medois sin oUtliers # E3: 17K Kmeans con outliers # E4: 17K Kmeans sin outliers # Aqui debe acompañarse las primeras medidas comerciales. # Guardamos a Modo de Respaldo. #save.image("E:/PROYECTOS/DOCENCIA_/INAP/Versión4_/Clustering_/Data_set_a_Utilizar_/Environment_4.RData") # Respaldamos Medio Ambiente Rdata. Guardar en formato RDS. #load("E:/PROYECTOS/DOCENCIA_/INAP/Versión4_/Clustering_/Data_set_a_Utilizar_/Environment_4.RData") # Para cargar ####################################################################################################################### ####################################################################################################################### ############################## ##################################### ############################## Generacion de Clusters 4 M y 17 K ##################################### ############################## ##################################### ####################################################################################################################### ####################################################################################################################### # Generaremos ############## ############## ############## Comenzaremos con 4 K, algoritmo PAM ############## ############## ############## variable.names(Master) str(Master) # Corro el Clusters set.seed(123) # Fijamos la Semillas. Esto es un hiperparámetro. Clusters_PAM_4K_ <- pam(x = Master, k = 4, metric = "manhattan", nstart = 1000) Clusters_PAM_4K_$medoids # Revisamos los objetos que se contienen dentro class(Clusters_PAM_4K_$cluster) View(Clusters_PAM_4K_$cluster) str(Clusters_PAM_4K_) Clusters_PAM_4K_$cluster # Ploteamos las dos primeras PCA. La funcion fviz_cluster realiza la reduccion PCA automaticamente, # plotea las dos primeras PCA y utiliza como etiqueta los clusters ya realizados. fviz_cluster(object = Clusters_PAM_4K_, data = Master, geom = "point", ellipse = FALSE, show.clust.cent = FALSE, pallete = "jco") + theme_bw() + theme(legend.position = "none") # Traemos los Clusters a un Datframe y luego plotearemos en 3D Clusters_PAM_4K_ <- as.data.frame(Clusters_PAM_4K_$clustering) # Creamos un Dat.frame con los Clustering (objeto) de los SS.OO. View(Clusters_PAM_4K_) Clusters_PAM_4K_$id_paises <- rownames(Clusters_PAM_4K_) # Des Indexo ID paises variable.names(Clusters_PAM_4K_) colnames(Clusters_PAM_4K_) <- c("Clusters_4K_PAM","id_paises") # Renombramos View(Clusters_PAM_4K_) variable.names(Clusters_PAM_4K_) str(Clusters_PAM_4K_) #### #### #### Plotearemos en 3D #### #### #### # Traemos la contribucion u observaciones para la construccion de las 3 primeras PCAs por individuo. View(Master) l <- PCA$x View(l) # Seleccionamos las dos primeras PCA, al igual que lo anterior, pero por individuo y convertimos en Dataframe. variable.names(l) l <- as.data.frame(l[,-4:-8]) str(l) View(l) # Si nos fijamos vemos la aportaci?n de de los pa?ses, por cada pa?s, en la Varianza Total del Modelo de Datos. l$Estados <- rownames(l) # Desendixamos str(l) View(l) # Ploteamos install.packages("rgl") library(rgl) #l$Estados <- rownames(l) # Des Indexar str(l) View(l) mycolors <- c('oldlace') s <- Clusters_PAM_4K_ s$Clusters_4K_PAM <- as.character(s$Clusters_4K_PAM) str(s$Clusters_4K_PAM) dim(s) library(dplyr) str(s) View(Juguete1) s <- rename(s,c(Estados="id_paises")) # funcion rename, cambia el nombre de la variable str(s) l <- merge(l, s, by="Estados", all.x = TRUE) View(l) l <- cbind(l, color = rep("oldlace",143)) # Creamos vector de colores str(l$color) # Reemplazamos l$color <- ifelse(l$Clusters_4K_PAM == 1, "yellow", l$color) l$color <- ifelse(l$Clusters_4K_PAM == 2, "red", l$color) l$color <- ifelse(l$Clusters_4K_PAM == 3, "lightblue", l$color) l$color <- ifelse(l$Clusters_4K_PAM == 4, "green", l$color) View(l) # Plot 3D uno. plot3d( x=l$PC1, y=l$PC2, z=l$PC3, col = l$color, type = 's', text, radius = 0.2, #aspect = !add, # Encenderlo y se verá sin rectas. xlab="PCA1", ylab="PCA2", zlab="PCA3") # Plot 3D dos con cordenadas. install.packages("plotly") library(plotly) fig <- l %>% plot_ly(x = l$PC1, y = l$PC2, z = l$PC3, type = "scatter3d", mode = 'makers', marker = list(size = 6), color=l$color) fig <- fig %>% layout(tittle = "3D Scatter plo") fig # Buscaremos el x o cordenada en el cubo dinámico del punto outliers hipotetico de interes, en el objeto l # para identificar su ID. View(l) # En ID ingresamos el numero de interes para saber que país es. print((paisDeInteres2 <- subset(Master_respaldo1, ID==67))[1,1]) # Plot con el data set original pero escalado, escogiendo dos variables. # Se seleccioanaron dos variables de manera aleatoria, para transmitir a los alumnos la conveniencia de proyectar el la dispersi?n de los datos # con todas las variables o de a pares (como es este caso) para enriquecer el an?lisis exploratorio. As? tambi?n, es conveniente proyectar # un plot tipo cubo o scatterplot en 3D. Pueden realizarlo siguiendo el Script del link que comparto a continuaci?n. Recomiendo # generar el S3D con las 3 primeras PCA para seber si efectivamente el traslape que pod?amos observar entre dos clusters al proyectarlos con # las dos primeras PCA, podr?a no observarse con la misma intensidad. O tal vez con m?s. De cualquier manera, enriquecer? el an?lisis # exploratorio. # http://www.sthda.com/english/wiki/scatterplot3d-3d-graphics-r-software-and-data-visualization # Ejemplo 1 2D: str(Master) plot(x = Master[,"Tasa_Mortalidad_Muertes/1000_habitantes"], y = Master[,"Deficit_en_Millones_de_Euros"], col = Clusters_PAM_4K_[,"Clusters_4K_PAM"], xlab = "Tasa_Mortalidad_Muertes/1000_habitantes", ylab = "Deficit_en_Millones_de_Euros") # Ejemplo 2 3D: # Plot 3D dos con cordenadas. #install.packages("plotly") library(plotly) fig <- l %>% plot_ly(x = Master$`Tasa_Mortalidad_Muertes/1000_habitantes`, y = Master$IDH, z = Master$Numero_de_Homicidios, type = "scatter3d", mode = 'makers', marker = list(size = 6), color=l$color) fig <- fig %>% layout(tittle = "3D Scatter plo") fig # Buscaremos el x o cordenada en el cubo dinámico del punto outliers hipotetico de interes, en el objeto l # para identificar su ID. View(l) # En ID ingresamos el numero de interes para saber que país es. print((paisDeInteres2 <- subset(Master_respaldo1, ID==67))[1,1]) ############## ############## ############## ############## ############## Seguimos con 17 K, algoritmo Kmeans ############## ############## ############## ############## ############## # Comenzaremos con 3 K, algoritmo Kmeans. str(Master) # Revisamos set.seed(88) # Fijamos semilla. Hiperparámetro. Clusters_Kmeans_17K_ <- kmeans(x = Master, centers = 17, nstart = 1000) # El argumento nstart corresponde al n?mero de iteraciones # # La opci?n nstart # intenta m?ltiples configuraciones iniciales e informa sobre la mejor. Por ejemplo, agregar nstart = 25 # generar? 25 centroides aleatorios iniciales y elegir? el mejor para el algoritmo. Tambi?n podr?amos indicar el # n?mero de iteraciones con el argumento "iter.max = N" # La base del algoritmo k-means: optimizando la suma errores al cuadrado. La suma de errores al cuadrado es igual a la distancia # eucl?dea entre las observaciones. # La idea del algoritmo es agrupar cada observaci?n (fila de una tabla) en k grupos, calcular el n?mero de clusters ?ptimo. # Para ello, el algoritmo sigue los siguientes pasos: # Se inicializan k puntos aleatorios, llamados centroides. # Para cada observaci?n se calcula la suma de errores al cuadrado de esa observaci?n respecto a cada uno de los k centroides. # A cada observaci?n se le asigna el centroide que menos error tenga. Clusters_Kmeans_17K_ # Revisamos los objetos que se contienen dentro class(Clusters_Kmeans_17K_$cluster) View(Clusters_Kmeans_17K_$cluster) str(Clusters_Kmeans_17K_) # Ploteamos las dos primeras PCA. La funcion fviz_cluster realiza la reduccion PCA automaticamente, # plotea las dos primeras PCA y utiliza como etiqueta los clusters ya realizados. fviz_cluster(object = Clusters_Kmeans_17K_, data = Master, geom = "point", ellipse = FALSE, show.clust.cent = FALSE, pallete = "jco") + theme_bw() + theme(legend.position = "none") # Traemos los Clusters a un Datframe Clusters_Kmeans_17K_ <- as.data.frame(Clusters_Kmeans_17K_$cluster) # Creamos un Dat.frame con los Clustering (objeto) de los SS.OO. View(Clusters_Kmeans_17K_) Clusters_Kmeans_17K_$id_paises <- rownames(Clusters_Kmeans_17K_) # Des Indexar View(Clusters_Kmeans_17K_) colnames(Clusters_Kmeans_17K_) <- c("Clusters_Kmeans_17K_","id_paises") # Renombramos View(Clusters_Kmeans_17K_) str(Clusters_Kmeans_17K_) # Puede plotear también en 3D lo anterior, pero no será eficiente dicha visualización al existir 17 colores, # es decir, 17 clusters o agrupaciones. ############## ############## ############## ############## ############## Seguimos con 7K, hierarchica ############## ############## (Jerarquizado) ############## ############## ############## # Plotearemos 7K Jeraquizado solo a modo de juguete para enseñar como funciona la sintaxis en R. # Finalmente seguimos con hierarchical, bajo el criterio Dendograma con tres enlaces o Linkage; complete, single y average. # Con 2 K para la generaci?n de Clusters. # Ahora, visualmente determinaremos si ocuparemos Linkage complete, Linkage single o Linkage average, para hierarchical. # Al aplicar un hierarchical clustering se tiene que escoger una medida de distancia (1-similitud. Es decir, euclidiana # o Manhattan para el caso de este curso) y un tipo de linkage determinado que tambi?n debemos ocupar. # En este caso, emplearemos la funci?n hclust() indicando la # distancia eucl?dea como medida de similitud y se comparan los linkages complete, single y average. # Posteriormente evaluaremos que linkage es el m?s apropiado. # Generamos primero la Matriz de distancia. # Para esto el algoritmo necesita calcular una matriz de distancia antes, en este caso euclidiana. matriz_distancias <- dist(x = Master, method = "euclidean") View(matriz_distancias) # Generamos los 3 Linkage o Enlaces. set.seed(56) # Fijamos la semilla para hacer el experimento reproducible. h_cluster_completo <- hclust(d = matriz_distancias, method = "complete") h_cluster_single <- hclust(d = matriz_distancias, method = "single") h_cluster_average <- hclust(d = matriz_distancias, method = "average") # Validaci?n del M?todo para distinguir si utilizar "complete", "single" o "average". # Una vez creado el dendrograma, hay que evaluar hasta qu? punto su estructura refleja las distancias originales # entre observaciones. # Una forma de hacerlo es empleando el coeficiente de correlaci?n entre las distancias cophenetic del dendrograma # (altura de los nodos) y la matriz de distancias original. Cuanto m?s cercano es el valor a 1, # mejor refleja el dendrograma la verdadera similitud entre las observaciones. # Valores superiores a 0.75 suelen considerarse como buenos. Esta medida puede emplearse como # criterio de ayuda para escoger entre los distintos m?todos de linkage. # En R, la funci?n cophenetic() (Correlaci?n cofen?tica) calcula las distancias cophenetic # de un hierarchical clustering. cor(x = matriz_distancias, cophenetic(h_cluster_completo)) # Probamos "complete" # 0.8113958 cor(x = matriz_distancias, cophenetic(h_cluster_single)) # Probamos "single" # 0.9020753 cor(x = matriz_distancias, cophenetic(h_cluster_average)) # Probamos "average" # 0.9395717 # Escogemos el Enlace promedio o Linkage Average. ################ ################ ################ Plotearemos lo anterior ################ ################ y Seleccionamos Visualmente ################ ################ donde Cortar el ?ptimo de K ################ ################ ################ #par(mfrow = c(3,1)) plot(x = h_cluster_completo, cex = 0.6, xlab = "", ylab = "", sub = "", main = "Linkage complete") plot(x = h_cluster_single, cex = 0.6, xlab = "", ylab = "", sub = "", main = "Linkage single") plot(x = h_cluster_average, cex = 0.6, xlab = "", ylab = "", sub = "", main = "Linkage average") # Los argumentos "cex = " nos ayudar?n a modificar la escala del plot para indicar la l?nea donde # deseo cortar. # Por Criterio Visual Corto en 53 para Obtener 3 Clusters. plot(x = h_cluster_average, cex = 0.7, xlab = "", ylab = "", sub = "", main = "Linkage average") # 4 Clusters abline(h = 5.5, lty = 2) # El argumento lty indica el tipo de l?nea gr?fica. El argumento h indica a la # altura que cortar? en el eje Y para seleccionar los Clusters. Miro el dendograna y por sobre todo # la matriz o dataframe. ################ ################ ################ ################ ################ Generamos 7 Clusters ################ ################ con hierarchical ################ ################ ################ #set.seed(65) # Fijamos la Semillas Clusters_hierarchical_7K_ <- as.data.frame(cutree(h_cluster_average, h = 5.5)) # Cutree es la funci?n para recibir los par?metros del dendograma # en 8.3 logro cortar el enlace para la generaci?n de 3 clusters. View(Clusters_hierarchical_7K_) Clusters_hierarchical_7K_$id_paises <- rownames(Clusters_hierarchical_7K_) # Des Indexar colnames(Clusters_hierarchical_7K_) <- c("Clusters_hierarchical_7K_","id_paises") # Renombro View(Clusters_hierarchical_7K_) str(Clusters_hierarchical_7K_) # Lo mismo; podría traer a un plot 3d como en el ejercicio anterior, pero la eficiencia de visualizar 7 agrupaciones # tan pocas observaciones (143) no sería eficiente para efectos de visualizacion. # Traemos los clusters que generamos a mi set de datos Master para unificar la informacion. Master$id_paises <- rownames(Master) # Des Indexar Master <- merge(Master, Clusters_PAM_4K_,"id_paises",all.x = TRUE) #Antes identificar en una columna o un vector en cada data.frame Master <- merge(Master, Clusters_Kmeans_17K_,"id_paises",all.x = TRUE) #Antes identificar en una columna o un vector en cada data.frame Master <- merge(Master, Clusters_hierarchical_7K_,"id_paises",all.x = TRUE) #Antes identificar en una columna o un vector en cada data.frame library(dplyr) str(Master) str(Master_respaldo1) Master_respaldo1 <- rename(Master_respaldo1, c(id_paises="ID")) # funcion rename, cambia el nombre de la variable variable.names(Master_respaldo1) str(Master) Master <- merge(Master, Master_respaldo1[,-3:-13],"id_paises",all.x = TRUE) #Antes identificar en una columna o un vector en cada data.frame View(Master) #save.image("E:/PROYECTOS/DOCENCIA_/INAP/Versión4_/Clustering_/Data_set_a_Utilizar_/Environment_5.RData") # Respaldamos Medio Ambiente Rdata. Guardar en formato RDS. #load("E:/PROYECTOS/DOCENCIA_/INAP/Versión4_/Clustering_/Data_set_a_Utilizar_/Environment_5.RData") # Para cargar ######################## ######################## ######################## ######################## ######################## Inspeccion de Clusters ######################## ######################## ######################## ######################## ######################## # Preprocesamos para consultar variable.names(Master_respaldo1) variable.names(Master) Master_estudios <- merge(Master_respaldo1[,-1], Master[,-2:-9],"id_paises",all.x = TRUE) #Antes identificar en una columna o un vector en cada data.frame View(Master_estudios) str(Master_estudios) variable.names(Master_estudios) Master_estudios[,1:15] <- lapply(Master_estudios[,1:15], as.numeric) # Convertimos todas variables en numeric salvo pais. str(Master_estudios) variable.names(Master_estudios) library(dplyr) library(inspectdf) Master_estudios %>% inspect_na dim(Master_estudios) # Asignamos contador Master_estudios <- cbind(Master_estudios, Freq = rep(1,143)) #Agrego vector poblado de 1. esto ya se enseñó. # Inspeccioanremos el Clusters 4 bajo la tecnica PAM. Consulta <- subset(Master_estudios, Clusters_4K_PAM==4) # Buscamos en que Clusters estan dim(Consulta) variable.names(Consulta) str(Consulta$`Tasa_Mortalidad_Muertes/1000_habitantes`) Tas_Mor_1000_H_ <- (format(summary(Consulta$`Tasa_Mortalidad_Muertes/1000_habitantes`), scientific = FALSE, big.mark = ",")) Personas_vacunadas <- (format(summary(Consulta$Personas_vacunadas), scientific = FALSE, big.mark = ",")) Def_Mill_Euros_ <- (format(summary(Consulta$Deficit_en_Millones_de_Euros), scientific = FALSE, big.mark = ",")) Def_PIB_ <- (format(summary(Consulta$Deficit_en_porc_del_PIB), scientific = FALSE, big.mark = ",")) IDH_ <- (format(summary(Consulta$IDH), scientific = FALSE, big.mark = ",")) F_Muertes_ <- (format(summary(Consulta$Muertes), scientific = FALSE, big.mark = ",")) Tasa_Morta_ <- (format(summary(Consulta$Tasa_mortalidad_en_porcentajes), scientific = FALSE, big.mark = ",")) F_Suicidios_ <- (format(summary(Consulta$Suicidios), scientific = FALSE, big.mark = ",")) VarPor_Suicidios_ <- (format(summary(Consulta$Varia_porc_Suicidios), scientific = FALSE, big.mark = ",")) F_Homicidios_ <- (format(summary(Consulta$Numero_de_Homicidios), scientific = FALSE, big.mark = ",")) VarPor_Homicidios <- (format(summary(Consulta$Varia_porc_Homicidios), scientific = FALSE, big.mark = ",")) print(cbind <- cbind(Tas_Mor_1000_H_, Personas_vacunadas, Def_Mill_Euros_, Def_PIB_, IDH_, F_Muertes_, Tasa_Morta_, F_Suicidios_, VarPor_Suicidios_, F_Homicidios_, VarPor_Homicidios)) print(paste("La composicion de este Clusters se caracteriza por poseer:", sep=" ", (paste("", sep="", sum(Consulta$Freq))),"Países ;", (paste("", sep="", (c <- (format((c <- sum(Consulta$Personas_vacunadas)), scientific = FALSE, big.mark = ","))))), "personas vacunadas ;", (paste("", sep="", (d <- (format((d <- sum(Consulta$Muertes)), scientific = FALSE, big.mark = ","))))), "muertes ; Y ", (paste("", sep="", (b <- (format((b <- sum(Consulta$Suicidios)), scientific = FALSE, big.mark = ","))))), "suicidios.")) #print(cbind) # Con big.mark = "," coloco la coma en en vez de punto en los miles. ####################################################################################################################### ####################################################################################################################### ############################## ##################################### ############################## DBSCAN ##################################### ############################## ##################################### ####################################################################################################################### ####################################################################################################################### # Idea intuitiva # Density-based spatial clustering of applications with noise (DBSCAN) fue presentado en 1996 por Ester et al. como una forma # de identificar clusters siguiendo el modo intuitivo en el que lo hace el cerebro humano, identificando regiones con alta # densidad de observaciones separadas por regiones de baja densidad. # Usamos esta t?cnica cuando esperamos que la distribuci?n espacial de los grupos no sea esf?rica, install.packages("fpc") install.packages("dbscan") install.packages("factoextra") library(fpc) library(dbscan) library(factoextra) variable.names(Master) # El algoritmo DBSCAN necesita dos parámetros: # Epsilon (e): Radio que define la region vecina a una observacion, tambien llamada epsilon-neighborhood. # Minimum points (minPts): Numero minimo de observaciones dentro de la region epsilon. # Empleando estos dos parámetros, cada observación del set de datos se puede clasificar en una de las siguientes tres categorías: # Core point: observación que tiene en su 𝜖-neighborhood un numero de observaciones vecinas igual o mayor a minPts. # Border point: observacion no satisface el mínimo de observaciones vecinas para ser core point pero que pertenece al # 𝜖-neighborhood de otra observación que sí es core point. # Noise u outlier: observación que no es core point ni border point. # Por ultimo, empleando las tres categorías anteriores se pueden definir tres niveles de conectividad entre observaciones: # 1. Directamente alcanzable (direct density reachable): una observacion 𝐴 es directamente alcanzable desde otra observación # 𝐵 si 𝐴 forma parte del 𝜖-neighborhood de 𝐵 y 𝐵 es un core point. Por definición, las observaciones solo pueden ser # directamente alcanzables desde un core point. # 2. Alcanzable (density reachable): una observación 𝐴 es alcanzable desde otra observación 𝐵 si existe una secuencia de # core points que van desde 𝐵 a 𝐴. # 3. Densamente conectadas (density conected): dos observaciones 𝐴 y 𝐵 están densamente conectadas si existe una observación core point 𝐶 tal que 𝐴 y 𝐵 son alcanzables desde 𝐶. # Seleccion de parametros. # Como ocurre en muchas otras tecnicas estadisticas, en DBSCAN no existe una forma unica y exacta de encontrar el # valor adecuado de epsilon epsilon (e) y minpts. A modo orientativo se pueden seguir las siguientes premisas: # minPts : cuanto mayor sea el tamaño del set de datos, mayor debe ser el valor minimo de observaciones vecinas. # En el libro Practical Guide to Cluster Analysis in R recomiendan no bajar nunca de 3. Si los datos contienen niveles # altos de ruido (outliers por ejemplo), aumentar minpts favorecera la creacion de clusters significativos menos # influidos por outliers. # epsilon: una buena forma de escoger el valor epsilon es estudiando las distancias promedio entre # las k = minPts observaciones mas proximas. # Al representar estas distancias en funcion de epsilon, el punto de inflexion de la curva suele ser un valor optimo. # Si el valor de epsilojn (e) escogido es muy pequeño, una proporcion alta de las observaciones no se asignaran a ningun cluster. # Por el contrario, si el valor es demasiado grande, la mayor?a de observaciones se agruparan en un unico cluster. # Seleccion del valor ?ptimo de epsilon. Para esto probaremos con los K de las recomendaciones y de # nuestros an?lisis. 3, 4 y 7 # Comenzaremos con 4 K # Selecci?n del epsilon. variable.names(Master) dbscan::kNNdistplot(Master[,2:9], k = 4) # Se selecciona K # La curva tiene el punto de inflexion en torno a 3 por lo que se escoge este valor como epsilon para DBSCAN. # Generamos el Cluster variable.names(Master) set.seed(321) # DBSCAN con epsilon = 3 y minPts = 4 dbscan_cluster <- fpc::dbscan(data = Master[,2:9], eps = 2.8, MinPts = 4) # dado que indica la distancia. # Resultados de la asignaci?n head(dbscan_cluster$cluster) str(dbscan_cluster) # Visualizaci?n de los clusters fviz_cluster(object = dbscan_cluster, data = Master[,2:9], stand = FALSE, geom = "point", ellipse = FALSE, show.clust.cent = FALSE, pallete = "jco") + theme_bw() + theme(legend.position = "bottom") # Seguiremos con 7 K para determinar si cambia el epsilon. Si no cambia, de juguete fijaremos 0.8 dbscan::kNNdistplot(Master[,2:9], k = 7) # La curva tiene el punto de inflexi?n en torno a 3. No cambi? en # relaci?n a la prueba anterior. Por ende es probable que se mantenga la selecci?n no esf?rica. # Mientras menor sea el epsilon mayor ser? el traslape o detecci?n. # Generamos el Cluster set.seed(321) # DBSCAN con epsilon = 3 y minPts = 2 dbscan_cluster <- fpc::dbscan(data = Master[,2:9], eps = 0.8, MinPts = 7) # Resultados de la asignaci?n head(dbscan_cluster$cluster) str(dbscan_cluster) # Visualizaci?n de los clusters fviz_cluster(object = dbscan_cluster, data = Master[,2:9], stand = FALSE, geom = "point", ellipse = FALSE, show.clust.cent = FALSE, pallete = "jco") + theme_bw() + theme(legend.position = "bottom") # Lo interesante es que a diferencia de kmeans, detecta los outliers. # Guardamos a Modo de Respaldo. #save.image("E:/PROYECTOS/DOCENCIA_/INAP/Versión4_/Clustering_/Data_set_a_Utilizar_/Environment_6.RData") # Respaldamos Medio Ambiente Rdata. Guardar en formato RDS. #load("E:/PROYECTOS/DOCENCIA_/INAP/Versión4_/Clustering_/Data_set_a_Utilizar_/Environment_6.RData") # Para cargar ####################################################################################################################### ####################################################################################################################### ############################## #################################### ############################## Ploteo Sofisticad deClusters #################################### ############################## librer?a gganimate #################################### ############################## #################################### ####################################################################################################################### ####################################################################################################################### install.packages("plotly") install.packages("gganimate") library(plotly) library(gganimate) # Llamamos un dataframe en formato rds Clusters_para_Plotly <- readRDS(file = "E:/PROYECTOS/DOCENCIA_/INAP/Versión1_/Adrian_Modulos_/Clustering_/Data_set_a_Utilizar_/Clusters_Años_.rds") # Data set de juguete. str(Clusters_para_Plotly) variable.names(Clusters_para_Plotly) View(Clusters_para_Plotly) plot1 <- ggplot(data=Clusters_para_Plotly, aes(x = Saldo_de_Credito, y = Avaluo_de_Propiedades, size = Clusters, color = Clusters, frame = Año)) + geom_point(aes(frame=Año,ids=Run_Falso)) ggplotly(plot1) # Un ploteo utilizando la renderizaci?n a trav?s de la funci?n transition_states ggplot(data=Clusters_para_Plotly, aes(x = Saldo_de_Credito, y = Avaluo_de_Propiedades, size = Clusters, color = Clusters, frame = Año)) + geom_point(aes(frame=Año,ids=Run_Falso)) + transition_states(Clusters, transition_length = 2, state_length = 1) # Otro tipo de Animaci?n. ggplot(Clusters_para_Plotly, aes(x = Saldo_de_Credito, y = Avaluo_de_Propiedades)) + geom_point(aes(colour = Clusters, group = 1L)) + transition_states(Clusters, transition_length = 2, state_length = 1) # Para Mayor informaci?n # https://translate.google.com/translate?hl=es&sl=en&u=https://gganimate.com/articles/gganimate.html&prev=search&pto=aue # https://community.rstudio.com/t/warning-message-file-renderer-failed-to-copy-frames-to-the-destination-directory/45261/4