#** 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)