#** 1.- DATOS**
library(ISLR)
names(NCI60)
## [1] "data" "labs"
Hacemos un str para comprobar la estructura de los datos del Data Frame.
# Dimensión de los datos
#-------------------------------
datos.nci <- NCI60$data
dim(datos.nci)
## [1] 64 6830
Tenemos un total de 64 variables con 6.830 datos en total.
# Estructura de los datos.
#------------------------
str(datos.nci)
## num [1:64, 1:6830] 0.3 0.68 0.94 0.28 0.485 ...
## - attr(*, "dimnames")=List of 2
## ..$ : chr [1:64] "V1" "V2" "V3" "V4" ...
## ..$ : chr [1:6830] "1" "2" "3" "4" ...
# Headde los datos.
#------------------------
head(datos.nci)[, 1:6]
## 1 2 3 4 5 6
## V1 0.300000 1.180000 0.550000 1.140000 -0.265000 -7.000000e-02
## V2 0.679961 1.289961 0.169961 0.379961 0.464961 5.799610e-01
## V3 0.940000 -0.040000 -0.170000 -0.040000 -0.605000 0.000000e+00
## V4 0.280000 -0.310000 0.680000 -0.810000 0.625000 -1.387779e-17
## V5 0.485000 -0.465000 0.395000 0.905000 0.200000 -5.000000e-03
## V6 0.310000 -0.030000 -0.100000 -0.460000 -0.205000 -5.400000e-01
#Estandarizamos los datos (normalizamos, esto es obligatorio para el análisis cluster)
#------------------------------------------------------
datos.nci <- scale(datos.nci, center = TRUE, scale = TRUE)
head(datos.nci)[, 1:6]
## 1 2 3 4 5 6
## V1 0.7229554 1.594614647 1.3152906 1.3450554 -0.6001006 -0.21892339
## V2 1.5838967 1.739790603 0.4382214 0.6489885 0.9047460 1.63581692
## V3 2.1731106 -0.016089747 -0.3463542 0.2643754 -1.3010255 -0.01917014
## V4 0.6776381 -0.372557113 1.6153098 -0.4408142 1.2346734 -0.01917014
## V5 1.1421409 -0.577195786 0.9575754 1.1298352 0.3585172 -0.03343823
## V6 0.7456141 -0.002887252 -0.1848054 -0.1202735 -0.4764080 -1.56012378
# Hacemos el summary, pues si es correcto, la media de las variables deve ser cero y la desviación típica 1
#----------------------------------------------------------
summary(datos.nci)[, 1:6]
## 1 2 3 4
## Min. :-2.3586 Min. :-2.85463 Min. :-3.90043 Min. :-2.0893
## 1st Qu.:-0.8008 1st Qu.:-0.49796 1st Qu.:-0.39825 1st Qu.:-0.9102
## Median : 0.0432 Median : 0.03672 Median : 0.04598 Median : 0.3010
## Mean : 0.0000 Mean : 0.00000 Mean : 0.00000 Mean : 0.0000
## 3rd Qu.: 0.7456 3rd Qu.: 0.50211 3rd Qu.: 0.42098 3rd Qu.: 0.9352
## Max. : 2.1731 Max. : 2.99408 Max. : 2.69999 Max. : 1.6748
## 5 6
## Min. :-1.75456 Min. :-2.01670
## 1st Qu.:-0.51764 1st Qu.:-0.46505
## Median :-0.05379 Median :-0.01917
## Mean : 0.00000 Mean : 0.00000
## 3rd Qu.: 0.37913 3rd Qu.: 0.50874
## Max. : 3.48176 Max. : 3.29103
# Tipo s de cáncer distintos en el set de datos
#---------------------------------------------
unique(NCI60$labs)
## [1] "CNS" "RENAL" "BREAST" "NSCLC" "UNKNOWN"
## [6] "OVARIAN" "MELANOMA" "PROSTATE" "LEUKEMIA" "K562B-repro"
## [11] "K562A-repro" "COLON" "MCF7A-repro" "MCF7D-repro"
# Número de muestras por tipo de cáncer
#-------------------------------------------
table(NCI60$labs)
##
## BREAST CNS COLON K562A-repro K562B-repro LEUKEMIA
## 7 5 7 1 1 6
## MCF7A-repro MCF7D-repro MELANOMA NSCLC OVARIAN PROSTATE
## 1 1 8 9 6 2
## RENAL UNKNOWN
## 9 1
#2.- CLUSTER JERÁRQUICO A continuación, se ,uestra un ejemplo usando complete, single y average linkage, escogiendo la distancia euclídea como medida de similitud.
# Matriz de distancia euclídea entre observaciones
#--------------------------------------------------
datos.nci.euc <- dist(datos.nci, method = "euclidean")
# Representación de dendogramas
#--------------------------------------------
par(mfrow = c(3,1)) # Matriz de gráficos (3x1)
#Plots
#--------------------------------
plot(hclust(datos.nci.euc, method = "complete"),
labels = NCI60$labs,
main = "Complete linkage",
xlab = "",
ylab = "",
cex = 0.3,
sub = "")
plot(hclust(datos.nci.euc, method = "average"),
labels = NCI60$labs,
main = "Average linkage",
xlab = "",
ylab = "",
cex = 0.3,
sub = "")
plot(hclust(datos.nci.euc, method = "single"),
labels = NCI60$labs,
main = "Single linkage",
xlab = "",
ylab = "",
cex = 0.3,
sub = "")
Aplicando cada uno de los métodos de enlace, se obtienen soluciones de
agrupamiento diferentes. En este caso, usaremos el método complete
linkage, ya que es el método donde la solución parece más robusta (menor
número de divisiones por nivel y no encontrarmos subdivisiones aisladas
como en las otras dos)
# Retornamos el mpar a su estado inicial
#-------------------------------------------
par(mfrow = c(1,1)) # Matriz de gráficos (1x1)
#Complete linkage
plot(hclust(datos.nci.euc, method = "complete"),
labels = NCI60$labs,
main = "Complete linkage",
xlab = "",
ylab = "",
cex = 0.3,
sub = "")
# Clustering jerárquico con complete linkage
#-------------------------------------------------
clust.comp <- hclust(dist(datos.nci))
clust.comp
##
## Call:
## hclust(d = dist(datos.nci))
##
## Cluster method : complete
## Distance : euclidean
## Number of objects: 64
# Corte del dendograma resultante en 4 clústeres. Almacena el clúster asignado a cada observación
#---------------------------------------------------------------------
clusteres.hc <- cutree(tree = clust.comp, k = 4)
clusteres.hc
## V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 V12 V13 V14 V15 V16 V17 V18 V19 V20
## 1 1 1 1 2 2 2 2 1 1 1 1 1 1 1 1 1 2 2 2
## V21 V22 V23 V24 V25 V26 V27 V28 V29 V30 V31 V32 V33 V34 V35 V36 V37 V38 V39 V40
## 1 1 1 1 1 1 1 1 1 1 1 1 1 3 3 3 3 3 3 3
## V41 V42 V43 V44 V45 V46 V47 V48 V49 V50 V51 V52 V53 V54 V55 V56 V57 V58 V59 V60
## 3 1 4 1 4 4 4 4 4 4 4 4 1 1 1 1 1 1 1 1
## V61 V62 V63 V64
## 1 1 1 1
# Representación del dendograma con el corte en 4 clusters
#---------------------------------------------------------------------------
plot(clust.comp, labels = NCI60$labs, cex = 0.4)
abline(h = 139, col = "red")
# Resultados del análisis clúster
#--------------------------------------------
table(clusteres.hc, NCI60$labs)
##
## clusteres.hc BREAST CNS COLON K562A-repro K562B-repro LEUKEMIA MCF7A-repro
## 1 2 3 2 0 0 0 0
## 2 3 2 0 0 0 0 0
## 3 0 0 0 1 1 6 0
## 4 2 0 5 0 0 0 1
##
## clusteres.hc MCF7D-repro MELANOMA NSCLC OVARIAN PROSTATE RENAL UNKNOWN
## 1 0 8 8 6 2 8 1
## 2 0 0 1 0 0 1 0
## 3 0 0 0 0 0 0 0
## 4 1 0 0 0 0 0 0
Todas las líneas celulares de leucemia se agrupan en el clúster 3, mientras que las líneas de cáncer de mama se reparten en tres clústeres diferentes.
Otra opción sería llevar a cabo el clustering jerárquico sobre los primeros componentes principales, de la siguiente manera:
# Cálculo de componentes principales
#---------------------------------------
pca.nci <- prcomp(datos.nci, scale = TRUE)
# Cinco primeros vectores de scores
#---------------------------------------
head(pca.nci$x)[,1:5]
## PC1 PC2 PC3 PC4 PC5
## V1 -19.68245 3.527748 -9.7354382 0.8177816 -12.511081
## V2 -22.90812 6.390938 -13.3725378 -5.5911088 -7.972471
## V3 -27.24077 2.445809 -3.5053437 1.3311502 -12.466296
## V4 -42.48098 -9.691742 -0.8830921 -3.4180227 -41.938370
## V5 -54.98387 -5.158121 -20.9291076 -15.7253986 -10.361364
## V6 -26.96488 6.727122 -21.6422924 -13.7323153 7.934827
# Clustering jerárquico sobre los primeros 5 componentes principales
#---------------------------------------------
clust.jer <- hclust(dist(pca.nci$x[,1:5]))
plot(clust.jer, labels = NCI60$labs, cex = 0.4,
main = "clust jerárquico sobre componentes principales")
# Resultados obtenidos para un cut de 4 clusters
#--------------------------------------------------
table(cutree(clust.jer, 4), NCI60$labs)
##
## BREAST CNS COLON K562A-repro K562B-repro LEUKEMIA MCF7A-repro MCF7D-repro
## 1 0 2 7 0 0 2 0 0
## 2 5 3 0 0 0 0 0 0
## 3 0 0 0 1 1 4 0 0
## 4 2 0 0 0 0 0 1 1
##
## MELANOMA NSCLC OVARIAN PROSTATE RENAL UNKNOWN
## 1 1 8 5 2 7 0
## 2 7 1 1 0 2 1
## 3 0 0 0 0 0 0
## 4 0 0 0 0 0 0
Los resultados difieren a los obtenidos utilizando el set de datos completo. Podríamos entender la aplicación de PCA como un proceso de eliminación de ruido en los datos.
#3.- K-Means
El resultado del clustering jerárquico (con el corte en el dendograma para obtener 4 clústeres) y el resultado de la aplicación de k-means clustering con K = 4 para el mismo problema pueden diferir. Veamos qué ocurriría en este caso, empleando el set de datos completo (aunque podría llevarse a cabo también sobre los primeros componentes principales):
# Asignación de semilla para resultados reproducibles
#------------------------------------------------
set.seed(12345)
# K-Means clustering con K = 4 y 20 asignaciones aleatorias de clústeres iniciales
#-----------------------------------------------------------
k.means <- kmeans(x = datos.nci, centers = 4, nstart = 20)
# Suma de cuadrados intra-cluster total
#-----------------------------------
k.means$tot.withinss
## [1] 344566.9
# Asignación de clústeres de cada observación
#-----------------------------------------
clusteres.km <- k.means$cluster
# Comparación clústeres de k-means y clustering jerarquico
#----------------------------------------------
table(clusteres.km, clusteres.hc)
## clusteres.hc
## clusteres.km 1 2 3 4
## 1 20 7 0 0
## 2 0 0 8 0
## 3 9 0 0 0
## 4 11 0 0 9
#BIPLOTS DE LOS CLUSTERS
#------------------------------------------------
library(FactoMineR)
# Cálculo componentes principals con la función PCA()
#-------------------------------------------------
pca.nci <- PCA(X = datos.nci, scale.unit = TRUE, ncp = 64, graph = TRUE)
library(factoextra)
## Loading required package: ggplot2
## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa
fviz_pca_ind(pca.nci, geom.ind = "point",
col.ind = as.factor(k.means$cluster),
axes = c(1, 2),
pointsize = 1.5)