library(readxl)
library(factoextra)## Loading required package: ggplot2## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBalibrary(cluster)
#Establacer el espacio de trabajo
dir = "/home/tati/Documents/PhD/2022-I/Metodos_Multivariados/"
setwd(dir)
dfa = read_excel("chaa.xlsx")
#Generación de la matrix de distancia usando el método de distancias Euclidiana
dist_mat = dist(dfa[,-4], method = "euclidean")
dim(as.matrix(dist_mat))## [1] 90 90#Generación del heatmap
heatmap(as.matrix(dist_mat))
# Ahora evaluamos las correlaciones entre los datos
# la correlación es otra forma de distancia
#Recuerda quitar la última etiqueta que es está en la última columna.
R= cor(dfa[,-4])
R## humedad almidon proteina
## humedad 1.0000000 0.8697394 -0.6366744
## almidon 0.8697394 1.0000000 -0.7928747
## proteina -0.6366744 -0.7928747 1.0000000heatmap(R)
library(corrplot)## corrplot 0.92 loaded#Haciendo uso de la libreria corrplot se explora de una manera visual la matrix de correlación.
#En este caso se hace uso del método "number" para que los números de los coeficientes presenten diferente color.
corrplot::corrplot(R, method = "number")
#En cuanto a la interpretación de este gráfico, los valores negativos indican una correlacion inversa, por ejemplo la relación entre proteína y húmedad (-0.64).Mientras que los valores positivos indican una correlación directa entre las variables.Existen otras librerias como PerformanceAnalytics que permiten la visualización de la Matriz de Correlación y el p-valor asociado. En la gráfica se observan: 1. En la diagonal se indica la distribución de cada variables 2. En la parte inferior se observan los diagramas de dispersión bivariados con una línea ajustada. 3. En la parte superior se observa el valor de la correlación junto con el nivel de significancia p-values, el cual es presentado con estrellas
library(PerformanceAnalytics)## Loading required package: xts## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric##
## Attaching package: 'PerformanceAnalytics'## The following object is masked from 'package:graphics':
##
## legendchart.Correlation(dfa[,-4])
#En cuanato a la interpretación, se destacan 3 grupos entre las variables de almidón y proteina Clustering jerárquico #Agrupación jerárquica aglomerativa. #A continuación se trabajó con agrupamiento jerárquico, especificamente con cuatro de los métodos implementados en el paquete hclust(Hierarchical Clustering), estos son: “single” - Single linkage method “average” Average linkage method “complete” - Complete linkage method “ward.D” - Variant of Ward’s method
En cuanto al método de enlace único “single” y el método de enlace completo “complete” en este se calculan la distancia entre el Cluster 1 y el Cluster 2 considerando, respectivamente, las distancias mínima y máxima entre las unidades asignadas a los dos conglomerados.
#A continuación se presenta un ejemplousando el método de "Single linkage method" o el vecino más cercano para formar el cluster.
hierar_cl = hclust(dist_mat, method = "single")
hierar_cl##
## Call:
## hclust(d = dist_mat, method = "single")
##
## Cluster method : single
## Distance : euclidean
## Number of objects: 90plot(hierar_cl)
abline(h=25, col="red")
# Ahora cambiamos el método a "complete" - Complete linkage o vecino más alejado
hierar_cl = hclust(dist_mat, method = "complete")
hierar_cl##
## Call:
## hclust(d = dist_mat, method = "complete")
##
## Cluster method : complete
## Distance : euclidean
## Number of objects: 90plot(hierar_cl)
abline(h=15, col = "red")
El siguiente método es “average”. En este se considera el valor promedio de las distancias entre las unidades asignadas a los dos clusters evaluados.
hierar_cl = hclust(dist_mat, method = "average")
hierar_cl##
## Call:
## hclust(d = dist_mat, method = "average")
##
## Cluster method : average
## Distance : euclidean
## Number of objects: 90plot(hierar_cl)
abline(h=10, col="red")
El último método a consider es ward.D. o suma de cuadrados. En este se mide la distancia entre el Cluster 1 y Cluster 2 en términos de la distancia entre las medias de conglomerados de los dos cluster. Se debe tener en cuenta que el método de Ward debe aplicarse a variables cuantitativas.
hierar_cl = hclust(dist_mat, method = "ward.D")
hierar_cl##
## Call:
## hclust(d = dist_mat, method = "ward.D")
##
## Cluster method : ward.D
## Distance : euclidean
## Number of objects: 90plot(hierar_cl)
#Luego de la evaluación de los datos con hclust, se empleó la función cutree especificando el número de grupos deseado, el cual es simbolizado con la letra k. Para el presente ejemplo k=3.
mod1= cutree(hierar_cl, k=3)
mod1## [1] 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 2 2 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2
## [39] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 2 3 3 3 3 3 3 3 3 3 3 3
## [77] 3 3 3 3 3 3 3 3 3 3 3 3 3 3rect.hclust(hierar_cl, k = 3, border = "green")
Matriz de confusión (o matriz de error) Esta matriz es una forma de resumir el rendimiento de un clasificador para tareas de clasificación.
##matriz de confusion
table(dfa$harina, mod1)## mod1
## 1 2 3
## ama 0 30 0
## arr 0 1 29
## cha 27 3 0#100*(27+29+30)/90library(caret)## Loading required package: latticelibrary(tibble)
library(cvms)
# Visualización de la matriz de confusión
df_cm = tibble(prediction=factor(mod1, c(2,1,3), c("cha","ama","arr"), ordered = T), target=factor(dfa$harina, ordered = T))
cm =as_tibble(table(df_cm))
plot_confusion_matrix(cm,target_col = 'target',
prediction_col = "prediction",
counts_col = 'n') ## Warning in plot_confusion_matrix(cm, target_col = "target", prediction_col =
## "prediction", : 'ggimage' is missing. Will not plot arrows and zero-shading.## Warning in plot_confusion_matrix(cm, target_col = "target", prediction_col =
## "prediction", : 'rsvg' is missing. Will not plot arrows and zero-shading.
Con el objetivo de comparar los cuatro métodos de hclust realizamos el análisis con dos datos sin estandarizar y luego estandarizados.
dfb=dfa[,1:3]
#Calculamos el número de cluster con los datos sin estandarizar
gap_stat = clusGap(dfb, FUN=hcut, nstart=25, K.max = 10, B =50)
#gap_stat
#Generar plot
fviz_gap_stat(gap_stat)
#Nota al calcular #Evaluamos todos los métodos sin estandarizar la información
dfb=dfa[,1:3]
dfb## # A tibble: 90 × 3
## humedad almidon proteina
## <dbl> <dbl> <dbl>
## 1 7.72 45.9 10.1
## 2 7.65 46.0 10.6
## 3 9.12 49.3 11.9
## 4 10.3 49.3 10.9
## 5 9.73 53.7 11.8
## 6 8.03 44.9 10.6
## 7 7.67 44.4 10.3
## 8 8.30 47.3 11.0
## 9 7.08 42.4 10.1
## 10 8.37 46.4 10.9
## # … with 80 more rowsm= c("average", "single", "complete", "ward")
names(m) = c("average", "single","complete", "ward")
ac = function(x){
agnes(dfb, method =x )$ac}
sapply(m, ac)## average single complete ward
## 0.9585437 0.8229332 0.9749480 0.9926597Ahora con los datos estandarizados empleando la función Scale() la cual centra y/o escala las columnas de una matriz numérica.
# Estandarizarción de los datos
dfb=scale(dfa[,1:3])
m= c("average", "single", "complete", "ward")
names(m) = c("average", "single","complete", "ward")
ac = function(x){
agnes(dfb, method =x )$ac
}
sapply(m, ac)## average single complete ward
## 0.9183692 0.8615611 0.9460516 0.9868489#Calculamos el número de cluster con los datos estandarizados
gap_stat = clusGap(dfb, FUN=hcut, nstart=25, K.max = 10, B =50)
#Generar plot
fviz_gap_stat(gap_stat)
#Usando la función summarise se pueden obtener algunos estadisticos con respecto a los clustes generados
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:xts':
##
## first, last## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, uniondfa$clust = mod1
b=dfa%>%
group_by(clust)%>%
summarise(mh=mean(humedad),
ma=mean(almidon),
mp= mean(proteina))
b## # A tibble: 3 × 4
## clust mh ma mp
## <int> <dbl> <dbl> <dbl>
## 1 1 7.97 45.5 10.6
## 2 2 10.6 57.7 8.66
## 3 3 11.2 68.7 7.48y=cbind(dfa$humedad, dfa$almidon, dfa$proteina); y ## [,1] [,2] [,3]
## [1,] 7.721578 45.90128 10.139913
## [2,] 7.652724 46.04431 10.598239
## [3,] 9.123838 49.27112 11.897356
## [4,] 10.268620 49.29338 10.925669
## [5,] 9.725898 53.65593 11.762685
## [6,] 8.032232 44.87319 10.588989
## [7,] 7.668702 44.41084 10.269907
## [8,] 8.297506 47.25879 11.049940
## [9,] 7.079020 42.37882 10.128377
## [10,] 8.372860 46.35557 10.890065
## [11,] 9.114401 48.73582 11.170282
## [12,] 7.624699 44.52850 10.061102
## [13,] 7.320832 43.71180 10.838241
## [14,] 7.918740 47.71021 10.994711
## [15,] 7.511419 47.13314 10.973178
## [16,] 7.432836 42.15872 10.452795
## [17,] 7.048295 42.45157 10.653445
## [18,] 9.495883 54.26860 12.213783
## [19,] 9.566397 52.33493 11.850099
## [20,] 8.428448 46.49959 10.900030
## [21,] 7.137224 41.30813 9.902352
## [22,] 7.996107 43.88332 10.601019
## [23,] 7.857408 46.03987 10.083055
## [24,] 8.475420 46.82790 10.963417
## [25,] 7.770122 44.57957 10.869224
## [26,] 8.337881 46.23931 10.272636
## [27,] 7.322514 43.44046 10.863603
## [28,] 7.120598 44.54915 10.050192
## [29,] 8.086297 47.17898 10.392180
## [30,] 8.476982 46.36449 10.825985
## [31,] 11.406200 61.52768 8.854925
## [32,] 11.865100 62.42881 9.582678
## [33,] 10.994501 56.17152 8.378672
## [34,] 8.547076 54.75710 7.944448
## [35,] 10.139319 58.35008 8.047774
## [36,] 11.251462 59.55176 8.516003
## [37,] 10.721676 56.58344 8.620696
## [38,] 11.169609 55.60496 8.312420
## [39,] 11.709281 58.70984 8.901657
## [40,] 10.631209 59.47148 8.164276
## [41,] 11.837275 61.97836 8.696064
## [42,] 10.899388 60.19528 8.676086
## [43,] 11.034178 59.11214 8.234401
## [44,] 10.364361 56.64930 8.074072
## [45,] 9.453628 55.83020 8.448407
## [46,] 11.304953 59.98648 8.208168
## [47,] 11.817252 58.19939 8.569444
## [48,] 8.854873 53.88221 7.792924
## [49,] 11.705227 56.49285 9.178956
## [50,] 9.864662 57.04852 7.720985
## [51,] 10.976432 58.43758 8.014868
## [52,] 9.758869 53.97065 7.716782
## [53,] 11.104813 54.70654 8.511586
## [54,] 10.752708 62.26898 8.319577
## [55,] 10.770620 56.78398 8.050882
## [56,] 10.874112 59.70302 9.000629
## [57,] 10.591282 60.77411 9.171865
## [58,] 11.443537 60.78443 9.100978
## [59,] 10.878284 57.91546 8.019242
## [60,] 7.694603 52.56606 7.072489
## [61,] 11.639415 69.06842 7.613599
## [62,] 10.560772 66.15191 7.179469
## [63,] 12.900343 70.87005 8.780711
## [64,] 11.536620 73.11278 8.425234
## [65,] 10.154163 61.85292 6.785652
## [66,] 11.339808 68.91843 7.454255
## [67,] 10.168894 66.09209 6.453378
## [68,] 12.010602 70.16086 8.223884
## [69,] 11.194967 72.04267 7.661631
## [70,] 10.411941 66.11417 6.885927
## [71,] 10.869293 65.10160 7.136257
## [72,] 10.965178 71.73988 7.507435
## [73,] 10.592406 66.33288 7.683299
## [74,] 10.677790 66.10951 6.824915
## [75,] 11.148594 70.48496 7.224137
## [76,] 11.253639 66.65250 7.103698
## [77,] 12.338184 71.15704 8.292266
## [78,] 11.690715 69.72544 7.711396
## [79,] 9.977745 68.58767 7.420820
## [80,] 10.762089 66.60586 7.493518
## [81,] 11.902305 70.32080 7.424477
## [82,] 11.073847 66.09679 6.724237
## [83,] 10.378991 69.68403 7.365961
## [84,] 11.920298 67.07839 7.663896
## [85,] 10.146433 71.52956 7.488857
## [86,] 12.372345 71.70629 7.511267
## [87,] 10.926828 65.26694 6.679849
## [88,] 11.028630 65.52440 7.168652
## [89,] 11.240799 68.61351 7.801128
## [90,] 12.772005 71.14063 7.872976mod2 = manova(y~dfa$harina)
summary(mod2)## Df Pillai approx F num Df den Df Pr(>F)
## dfa$harina 2 1.4658 78.65 6 172 < 2.2e-16 ***
## Residuals 87
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Creamos el objeto Y con la información de la variable húmedad y almidón.Con este realizamos un manova, la cual analiza simultaneamente la relación entre varias variables de respuesta y un conjunto común de predictores.
#solo con dos variables
y=cbind(dfa$humedad, dfa$almidon); y ## [,1] [,2]
## [1,] 7.721578 45.90128
## [2,] 7.652724 46.04431
## [3,] 9.123838 49.27112
## [4,] 10.268620 49.29338
## [5,] 9.725898 53.65593
## [6,] 8.032232 44.87319
## [7,] 7.668702 44.41084
## [8,] 8.297506 47.25879
## [9,] 7.079020 42.37882
## [10,] 8.372860 46.35557
## [11,] 9.114401 48.73582
## [12,] 7.624699 44.52850
## [13,] 7.320832 43.71180
## [14,] 7.918740 47.71021
## [15,] 7.511419 47.13314
## [16,] 7.432836 42.15872
## [17,] 7.048295 42.45157
## [18,] 9.495883 54.26860
## [19,] 9.566397 52.33493
## [20,] 8.428448 46.49959
## [21,] 7.137224 41.30813
## [22,] 7.996107 43.88332
## [23,] 7.857408 46.03987
## [24,] 8.475420 46.82790
## [25,] 7.770122 44.57957
## [26,] 8.337881 46.23931
## [27,] 7.322514 43.44046
## [28,] 7.120598 44.54915
## [29,] 8.086297 47.17898
## [30,] 8.476982 46.36449
## [31,] 11.406200 61.52768
## [32,] 11.865100 62.42881
## [33,] 10.994501 56.17152
## [34,] 8.547076 54.75710
## [35,] 10.139319 58.35008
## [36,] 11.251462 59.55176
## [37,] 10.721676 56.58344
## [38,] 11.169609 55.60496
## [39,] 11.709281 58.70984
## [40,] 10.631209 59.47148
## [41,] 11.837275 61.97836
## [42,] 10.899388 60.19528
## [43,] 11.034178 59.11214
## [44,] 10.364361 56.64930
## [45,] 9.453628 55.83020
## [46,] 11.304953 59.98648
## [47,] 11.817252 58.19939
## [48,] 8.854873 53.88221
## [49,] 11.705227 56.49285
## [50,] 9.864662 57.04852
## [51,] 10.976432 58.43758
## [52,] 9.758869 53.97065
## [53,] 11.104813 54.70654
## [54,] 10.752708 62.26898
## [55,] 10.770620 56.78398
## [56,] 10.874112 59.70302
## [57,] 10.591282 60.77411
## [58,] 11.443537 60.78443
## [59,] 10.878284 57.91546
## [60,] 7.694603 52.56606
## [61,] 11.639415 69.06842
## [62,] 10.560772 66.15191
## [63,] 12.900343 70.87005
## [64,] 11.536620 73.11278
## [65,] 10.154163 61.85292
## [66,] 11.339808 68.91843
## [67,] 10.168894 66.09209
## [68,] 12.010602 70.16086
## [69,] 11.194967 72.04267
## [70,] 10.411941 66.11417
## [71,] 10.869293 65.10160
## [72,] 10.965178 71.73988
## [73,] 10.592406 66.33288
## [74,] 10.677790 66.10951
## [75,] 11.148594 70.48496
## [76,] 11.253639 66.65250
## [77,] 12.338184 71.15704
## [78,] 11.690715 69.72544
## [79,] 9.977745 68.58767
## [80,] 10.762089 66.60586
## [81,] 11.902305 70.32080
## [82,] 11.073847 66.09679
## [83,] 10.378991 69.68403
## [84,] 11.920298 67.07839
## [85,] 10.146433 71.52956
## [86,] 12.372345 71.70629
## [87,] 10.926828 65.26694
## [88,] 11.028630 65.52440
## [89,] 11.240799 68.61351
## [90,] 12.772005 71.14063mod2 = manova(y~dfa$harina)
summary(mod2)## Df Pillai approx F num Df den Df Pr(>F)
## dfa$harina 2 1.2255 68.834 4 174 < 2.2e-16 ***
## Residuals 87
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(dfa$humedad, dfa$almidon, col=c("red","blue","orange") )
Una alternativa para establecer el número de cluster es empleando la función fviz_nbclust, la cual determina y realiza la visualización del número óptimo de clústers. También empleamos el algoritmo de cluster K-means, el cual agrupa objetos en k grupos basándose en sus características. El agrupamiento se realiza minimizando la suma de distancias entre cada objeto y el centroide de su grupo o cluster. la función fviz_nbclust permite evaluar el número de clústerse a través de tres métodos: “silhouette” (para el ancho promedio), “wss” (suma de cuadrados) y “gap_stat” (estadisticas).
### suma de cuadrados
fviz_nbclust(dfb, kmeans,method = "wss" )
#set.seed(1)
km = kmeans(dfb, centers = 3, nstart = 25)
km## K-means clustering with 3 clusters of sizes 30, 30, 30
##
## Cluster means:
## humedad almidon proteina
## 1 0.7438022 1.14051632 -0.9437156
## 2 0.4212798 0.04391756 -0.3166988
## 3 -1.1650820 -1.18443389 1.2604144
##
## Clustering vector:
## [1] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 2 2 2 2 2 2 2 2
## [39] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [77] 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##
## Within cluster sum of squares by cluster:
## [1] 12.92903 17.16446 15.70853
## (between_SS / total_SS = 82.8 %)
##
## Available components:
##
## [1] "cluster" "centers" "totss" "withinss" "tot.withinss"
## [6] "betweenss" "size" "iter" "ifault"fviz_cluster(km,data = dfb)
#ahora con 2 cluster
set.seed(1)
km = kmeans(dfb, centers = 2, nstart = 25)
km## K-means clustering with 2 clusters of sizes 60, 30
##
## Cluster means:
## humedad almidon proteina
## 1 0.582541 0.5922169 -0.6302072
## 2 -1.165082 -1.1844339 1.2604144
##
## Clustering vector:
## [1] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 1
## [39] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [77] 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##
## Within cluster sum of squares by cluster:
## [1] 55.58899 15.70853
## (between_SS / total_SS = 73.3 %)
##
## Available components:
##
## [1] "cluster" "centers" "totss" "withinss" "tot.withinss"
## [6] "betweenss" "size" "iter" "ifault"fviz_cluster(km,data = dfb)
# ahora 4 cluster
set.seed(1)
km = kmeans(dfb, centers = 4, nstart = 25)
km## K-means clustering with 4 clusters of sizes 24, 6, 30, 30
##
## Cluster means:
## humedad almidon proteina
## 1 -1.3855555 -1.31425807 1.1170586
## 2 -0.2831879 -0.66513715 1.8338374
## 3 0.4212798 0.04391756 -0.3166988
## 4 0.7438022 1.14051632 -0.9437156
##
## Clustering vector:
## [1] 1 1 2 2 2 1 1 1 1 1 2 1 1 1 1 1 1 2 2 1 1 1 1 1 1 1 1 1 1 1 3 3 3 3 3 3 3 3
## [39] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
## [77] 4 4 4 4 4 4 4 4 4 4 4 4 4 4
##
## Within cluster sum of squares by cluster:
## [1] 4.175005 1.211874 17.164460 12.929032
## (between_SS / total_SS = 86.7 %)
##
## Available components:
##
## [1] "cluster" "centers" "totss" "withinss" "tot.withinss"
## [6] "betweenss" "size" "iter" "ifault"fviz_cluster(km,data = dfb)
# Finalmente podemos encontrar la media en cada cluster
aggregate(dfb, by=list(cluster=km$cluster), mean)## cluster humedad almidon proteina
## 1 1 -1.3855555 -1.31425807 1.1170586
## 2 2 -0.2831879 -0.66513715 1.8338374
## 3 3 0.4212798 0.04391756 -0.3166988
## 4 4 0.7438022 1.14051632 -0.9437156