Preparación

En primer lugar, cargamos las librerías necesarias.

if (!require(readxl)) install.packages('readxl', dependencies = T)
## Loading required package: readxl
library(readxl)

Ejercicio 1: PCA

Cargamos la base de datos que necesitaremos para el primer ejercicio.

PROSTATA2 <- read_excel("PROSTATA2.xlsx")

Análisis exploratorio de los datos

El análisi exploratorio lo comenzaremos con la instrucción dim() para conocer el número de variables (5) y de muestras (97) de la base de datos. A continucación, la variable summary() nos ofrece información sobre las medidas de cada variable.

dim(PROSTATA2)
## [1] 97  5
summary(PROSTATA2)
##      lcavol           lweight           lbph              lcp         
##  Min.   :-1.3471   Min.   :2.375   Min.   :-1.3863   Min.   :-1.3863  
##  1st Qu.: 0.5128   1st Qu.:3.376   1st Qu.:-1.3863   1st Qu.:-1.3863  
##  Median : 1.4469   Median :3.623   Median : 0.3001   Median :-0.7985  
##  Mean   : 1.3500   Mean   :3.629   Mean   : 0.1004   Mean   :-0.1794  
##  3rd Qu.: 2.1270   3rd Qu.:3.876   3rd Qu.: 1.5581   3rd Qu.: 1.1787  
##  Max.   : 3.8210   Max.   :4.780   Max.   : 2.3263   Max.   : 2.9042  
##       lpsa        
##  Min.   :-0.4308  
##  1st Qu.: 1.7317  
##  Median : 2.5915  
##  Mean   : 2.4784  
##  3rd Qu.: 3.0564  
##  Max.   : 5.5829

Para comparar las medias y las variaznas de las variables usamos la función apply que nos permite aplicar otras funciones, en este caso mean y var a cada variable.

apply(PROSTATA2, 2 , mean)
##     lcavol    lweight       lbph        lcp       lpsa 
##  1.3500096  3.6289427  0.1003556 -0.1793656  2.4783869

Respecto a la media, podemos observar que hay mucha diferencia entre unos datos y otros. Así, por ejemplo, observamos que la variables lweight es puntúa 36 veces más alto que lbph.

apply(PROSTATA2, 2 , var)
##    lcavol   lweight      lbph       lcp      lpsa 
## 1.3891566 0.1835362 2.1048399 1.9551020 1.3324756

También se observan diferencias muy notables entre las varianzas de las diferentes variables.

Por estos motivos resulta importante estandarizar las variables para que tengan media cero y desviación típica uno antes de realizar el Análisis de Componentens Principales.

Análisis de Componentes Principales

Tal como hemos visto en el apartado anterior, comenzamos estandarizando las variables. Por defecto, con la función prcomp() se centran las variables para que tengan media cero. En la formulación, usando la opción scale = TRUE, escalamos las variables para que tengan una desviación estándar de uno.

pr.out <- prcomp(PROSTATA2 , scale = TRUE)
names(pr.out)
## [1] "sdev"     "rotation" "center"   "scale"    "x"

Matriz de rotación

Observamos la matriz de rotación que nos proporciona las cargas de los componentes principales. Cada columna de pr.out$rotation contiene el vector de carga del componente principal correspondiente.

pr.out$rotation
##               PC1        PC2         PC3        PC4         PC5
## lcavol  0.5473154 -0.2587215 -0.01397918 -0.2490854  0.75582412
## lweight 0.3603542  0.5456492  0.62576090  0.4185211  0.07533353
## lbph    0.1770236  0.7197733 -0.65703021 -0.1204011  0.06636225
## lcp     0.4787795 -0.3406613 -0.39098637  0.6619246 -0.25239953
## lpsa    0.5567976 -0.0347343  0.15384681 -0.5569164 -0.59577285

Representación de los dos primeros componentes principales

A través de la función biplo() graficamos los resultados.

biplot(pr.out , scale = 0)

Se puede observar que los pacientes con una alta puntuación en el Componente 1, como el 97 o el 89, tienen un alto nivel de antígeno prostático específico, mientras que pacientes como el 1 o el 3 que puntuaron muy bajo en el Componente 1 tienen un bajo nivel de antígeno prostático específico.

Desviación típica, varianza y PVE explicadas por cada componente principal

Una vez más usamos la función summary para que nos de los diferentes datos del Análisis de Componentes Principales que hemos realizado. Así podemos estudiar la desviación estandar de cada componente.

summary(pr.out)
## Importance of components:
##                           PC1    PC2    PC3     PC4     PC5
## Standard deviation     1.5911 1.1489 0.7301 0.62405 0.47531
## Proportion of Variance 0.5063 0.2640 0.1066 0.07789 0.04518
## Cumulative Proportion  0.5063 0.7703 0.8769 0.95482 1.00000

Sin embargo,. de esta manera redondéa los resultados. Así pues, para obtener la numeración completa podemos llamar a dicho dato.

pr.out$sdev
## [1] 1.5910606 1.1489468 0.7301351 0.6240472 0.4753054

En cualquier caso, tal como podemos obervar, la SD va disminuyendo con cada componente

La varianza explicada por cada componente principal se obtiene mediante sus cuadrados

pr.var <- pr.out$sdev^2
pr.var
## [1] 2.5314738 1.3200787 0.5330972 0.3894350 0.2259153

Al igual que la SD, la vairanza explicada por cada componente principal va disminuyendo

Si bien la función summary() ya nos ofreció proporción de la varianza explicada por cada componente principal, esta estaba redondeada. Para obtener el dato completo simplemente dividimos la varianza explicada por cada componente principal por la varianza total explicada por los cuatro componentes principales

pve <- pr.var / sum(pr.var) 
pve
## [1] 0.50629476 0.26401575 0.10661944 0.07788699 0.04518305

Tal como se puede observar, el primer componente ofrece el 50.63% de la varianza de los datos, el siguiente el 26.40%, el siguiente el 10.66%, el cuarto el 7.78% y el último el 4.52%

Gráfica del PVE explicado por cada componente y el PVE acumulado

Para trazar el PVE explicado por cada componente, así como el PVE acumulado, lo realizamos de la siguiente manera

par(mfrow = c(1, 2))
plot(pve , xlab = "Componente Principal", ylab = "Proporción de Varianza Explicada", ylim = c(0, 1), type = "b")
plot(cumsum(pve), xlab = "Componente Principal", ylab = "Proporción de Varianza Explicada Acumulada", ylim = c(0, 1), type = "b")

Ejercicio 2: Cluster

Cargamos la base de datos que necesitaremos para el segundo ejercicio.

EUROEMPLEO <- read_excel("EUROEMPLEO.xlsx")

Para poder realizar el resto del ercicio observamos la base de datos

dim(EUROEMPLEO)
## [1] 25 10

Arreglamos la base de datos para que el nombre del pais sea el nombre de la fila y no solo una variable

EUROEMPLEO1 <- EUROEMPLEO
row.names (EUROEMPLEO1) <- c("Belgium", "Denmark", "France", "Germany", "Ireland", 
                             "Italy", "Luxembourg","Netherlands", "United Kingdom",
                             "Austria","Finland","Greece","Norway","Portugal","Spain",
                             "Sweden","Switzerland","Turkey","Bulgaria","Czechia",
                             "Hungary","Poland","Romania","Russia", "Serbia")
## Warning: Setting row names on a tibble is deprecated.
EUROEMPLEO1 <- EUROEMPLEO1 [ , -1]

Al igual que en el ejercicio anterior, vamos a comparar la mediana y la varianza

apply(EUROEMPLEO1, 2 , mean)
##    Agr    Min    Man     PS    Con     SI    Fin    SPS     TC 
## 19.728  1.188 26.440  0.892  8.188 13.028  4.112 19.940  6.472
apply(EUROEMPLEO1, 2 , var)
##         Agr         Min         Man          PS         Con          SI 
## 242.1062667   0.8627667  42.4266667   0.1407667   2.8069333  21.6712667 
##         Fin         SPS          TC 
##   7.8652667  48.3991667   1.8679333

Observamos que hay mucha diferencia entre las diferentes variables. Por tanto normalizamos esos datos con la función scale()

EUROEMPLEO2 <- scale(EUROEMPLEO1)

Si bien en el ejercicio ya nos marca cuantos cluster debemos hacer, para ampliar conocimientos analizaré según diversos estadísticos cuantos clusters sería lo más adecuado. Para ello instalamos el paquete “NbClust”

if (!require(NbClust)) install.packages('NbClust', dependencies = T)
## Loading required package: NbClust
library(NbClust) 
resnumclust<-NbClust(EUROEMPLEO2, distance = "euclidean", min.nc=2, max.nc=10, method = "kmeans", index = "alllong")
## Warning in pf(beale, pp, df2): NaNs produced

## *** : The Hubert index is a graphical method of determining the number of clusters.
##                 In the plot of Hubert index, we seek a significant knee that corresponds to a 
##                 significant increase of the value of the measure i.e the significant peak in Hubert
##                 index second differences plot. 
## 

## *** : The D index is a graphical method of determining the number of clusters. 
##                 In the plot of D index, we seek a significant knee (the significant peak in Dindex
##                 second differences plot) that corresponds to a significant increase of the value of
##                 the measure. 
##  
## ******************************************************************* 
## * Among all indices:                                                
## * 9 proposed 2 as the best number of clusters 
## * 12 proposed 3 as the best number of clusters 
## * 1 proposed 5 as the best number of clusters 
## * 1 proposed 6 as the best number of clusters 
## * 2 proposed 8 as the best number of clusters 
## * 1 proposed 9 as the best number of clusters 
## * 2 proposed 10 as the best number of clusters 
## 
##                    ***** Conclusion *****                            
##  
## * According to the majority rule, the best number of clusters is  3 
##  
##  
## *******************************************************************

Tal como podemos obervar, lo más adecuado swería hacer un cluster de 3 grupos ya que hay 12 esradísticos que así lo indican. EN cualquier caso, otros 9 estadísticos consideran que 2 cluster es adecuado

Clustering K-means con 2 grupos

Con la función kmenas realizamos el clúster

set.seed (2)
km.out <- kmeans(EUROEMPLEO2, 2, nstart = 20)

La asignaciones de cluster están contenidas km.out$cluster.

km.out$cluster
##  [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 2

A continuación podemos ver graficamente la dsitribución

plot(EUROEMPLEO2, col = (km.out$cluster+1),
main = "Resultados de la agrupación K-Means con K = 2",
xlab = "", ylab = "", pch = 20, cex = 2)

Para mejorar la represención gráfica, usaremos la biblioteca “factoextra” y la instrucción fviz_cluster

if (!require(factoextra)) install.packages('factoextra', dependencies = T)
## Loading required package: factoextra
## Loading required package: ggplot2
## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa
library(factoextra) 
fviz_cluster(km.out, data =EUROEMPLEO2,repel = TRUE,star.plot = TRUE, palette = "Set2", main = "Cluster con 2 grupos")

Clustering K-means con 3 grupos

Con la función kmenas realizamos el clúster

set.seed (3)
km.out <- kmeans(EUROEMPLEO2, 3, nstart = 20)

La asignaciones de cluster están contenidas km.out$cluster.

km.out$cluster
##  [1] 1 1 1 1 1 1 1 1 1 1 1 3 1 1 1 1 1 2 3 3 3 3 3 3 2

A continuación podemos ver graficamente la dsitribución.

plot(EUROEMPLEO2, col = (km.out$cluster + 1),
main = "Resultados de la agrupación K-Means con K = 3",
xlab = "", ylab = "", pch = 20, cex = 2)

Al igual que antes, graficaremos con fviz_cluster.

fviz_cluster(km.out, data =EUROEMPLEO2,repel = TRUE,star.plot = TRUE, palette = "Set2", main = "Cluster con 3 grupos")

Análisis comparativos con diferentes nstart

Para realizar la comparativa, realizaremos los diferentes cluster.

set.seed (3)
km.out_1 <- kmeans(EUROEMPLEO2, 2, nstart = 1)
km.out_1$tot.withinss
## [1] 155.1147
km.out_20 <- kmeans(EUROEMPLEO2, 2, nstart = 20)
km.out_20$tot.withinss
## [1] 155.1147
km.out_50 <- kmeans(EUROEMPLEO2, 2, nstart = 50)
km.out_50$tot.withinss
## [1] 155.1147

En este caso observamos que no las diferentes nstart realizadas siempre nos ofrece el mismo valor.

Agrupación jerarquica

Para realizar la agrupación jerárquica usamos la función hclust(). Para prácticar, usaremos las tres propuestas que ofrece el manual.

set.seed (3)
hc.complete <- hclust(dist(EUROEMPLEO2), method = "complete")
hc.average <- hclust(dist(EUROEMPLEO2), method = "average")
hc.single <- hclust(dist(EUROEMPLEO2), method = "single")

Reepresentación gráfica del dendograma

A continuación podemos graficar el dendograma

par(mfrow = c(1, 3))
plot(hc.complete , main = "Complete Linkage",
     xlab = "", sub = "", cex = .9)
plot(hc.average , main = "Average Linkage",
     xlab = "", sub = "", cex = .9)
plot(hc.single, main = "Single Linkage",
     xlab = "", sub = "", cex = .9)

Detemrinación del número de etiquetas

Para determinar las etiquetas de los clústers para cada observación asociada con un corte dado del dendrograma utilizaremos la función cutree() con k=2 que es el número de clusters realizados.

set.seed (3)
cutree(hc.complete , 2)
##  [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 2
cutree(hc.average , 2)
##  [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 2
cutree(hc.single, 2)
##  [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1

Finalmente, graficamos todo.

plot(hclust(dist(EUROEMPLEO2), method = "complete"),
       main = "Agrupación jerárquica")

Al igual que se presenta en el manual, la distancia basada en la correlación puede calcularse utilizando la función as.dist(), que convierte una matriz simétrica cuadrada arbitraria en una forma que la función hclust() reconoce como una matriz de distancia.Sería de la siguiente manera

set.seed (3)
dd <- as.dist (1 - cor(t(EUROEMPLEO2)))
plot(hclust(dd, method = "complete"),
       main = "Complete Linkage with Correlation -Based Distance", xlab = "", sub = "")