Agrupación de la información

Análisis de componentes principales

El Análisis de Componentes Principales (ACP) es una técnica estadística de síntesis de la información, o reducción de la dimensión (número de variables.

Dadas \(n\) observaciones de \(p\) variables, se buscan \(m<p\) variables que sean combinaciones lineales de las \(p\) originales y que estén incorreladas, recogiendo la mayor parte de la información o variabilidad de los datos. Hay que tener presente que, si las variables originales están incorreladas de partida, entonces no tiene sentido realizar un análisis de componentes principales.

Se considera una serie de variables \((x_1,x_2,...x_p)\) sobre un grupo de objetos o individuos y se trata de calcular, a partir de ellas, un nuevo conjunto de variables \((y_1,y_2,...y_p)\),incorreladas entre sí, cuyas varianzas vayan decreciendo progresivamente. Cada \(y_j\) (donde \(j=1,...,p\)) es una combinación lineal de las \(x_1,x_2,...x_p\) originales, es decir:

\(y_j=a_{j1}x_1+a_{j2}x_2+...+a_{jP}x_p=a'_jx\)

donde \(a'_j=(a_{j1},a_{j2},...,a_{jP})\) y

\[x= \begin{bmatrix} x_{1}\\ x_{2}\\ . . . x_{p} \end{bmatrix}\]

El objeto del análisis es maximizar la varianza del componente obtenido, pero manteniendo la ortogonalidad de la transformación, esto es \(a'_ja_j=\sum_{j=1}^p a_{kj}^2=1\).

El primer componente por tanto se obtendrá:

\(Var(y_1)=Var(a'_1 x)=a'_1\sigma a_1\)

El problema consiste en maximizar la función \(a'_1\sigma a_1\) donde \(\sigma\) es la matriz de covarianzas de las \(p\) variables,sujeta a la restricción \(a'_1a_1=1\), utilizandose para ello los multiplicadores de Lagrange:

\(L(a_1)=a'_1\sigma a_1-\lambda (a'_1a_1-1)\)

siendo ahora \(a_1\) la incognita, se busca el maximo:

\(\frac{L(\partial a_1)}{\partial a_1}=0 \longrightarrow (\sigma - \lambda I)a_1=0\)

La solución del sistema de ecuaciones requiere que el determinante \(\mid \sigma - \lambda I \mid\) sea no nulo, de manera que \(\lambda\) es un autovalor de \(\sigma\), que al ser definida positiva, tendrá \(p\) autovalores distintos.

Entonces, para maximizar la varianza de \(y_1\) se tiene que tomar el mayor autovalor, \(\lambda_1\) y el correspondiente autovector de \(a_1\).

El segundo componente principal,\(y_2=a'_2x\) ha de estar incorrelacionado con \(y_1\), es decir:

\(cov(y_1,y_2)=cov(a'_1x,a'_2x)=a'_1 \sigma a_2=0\).

Equivale a:

\(a'_1 \lambda a_2=0\)

Y a maximizar las \(Var(y_2)=a'_2\sigma a_2\) imponinedo ahora dos restricciones:

\(a'_2a_2=1\)

\(a'_2a_1=0\)

Se forma el lagrangiano \(L(a_2)=a'_12\sigma a_2-\lambda (a'_2a_2-1) - \delta (a'2a_1)\), se deriva respecto a \(a_2\) y se iguala a cero:

\(\frac{L(\partial a_2)}{\partial a_2}=0 \longrightarrow (\sigma - \lambda I)a_2=0\).

Elegimos \(\lambda_2\) como segundo autovalor de la matriz \(\sigma\) asociado al autovector asociado \(a_2\).

Entonces todos los p componentes de \(y\) se pueden expresar como el producto de una matriz formada por los autovectores, multiplicada por el vector \(x\) que contiene las variables originales:

\(y=Ax\)

donde la \(Var(y_j)=\lambda_j\), y la \(Cov(y_j,y_s)=0\).

Si sumamos todos los autovalores,tendremos entonces la varianza total de los componentes, que por las propiedades del operador traza, acaba siendo igual que la varianza de los variables originales:

\(\sum Var(y_j)= \sum (\lambda_j)=\sum Var(x_j)\)

El porcentaje de varianza total que recoge cada componentes principales, es el critierio utilizado para reducir las \(p\) variables de nuestro conjunto de datos, por \(m\) combinaciones lineales de dichas variables.

El ACP se emplea sobre todo en análisis exploratorio de datos (Analisis factorial) y para construir modelos predictivos cuando las variables explicativas tienen problemas de colinealidad.

El ACP comporta el cálculo de la descomposición en autovalores de la matriz de covarianza, normalmente tras centrar los datos en la media de cada atributo.

El análisis de componentes principales (PCA) puede realizarse con una función de la librería básica de «R»: prcomp o princomp.

require(graphics)
str(USArrests)
## 'data.frame':    50 obs. of  4 variables:
##  $ Murder  : num  13.2 10 8.1 8.8 9 7.9 3.3 5.9 15.4 17.4 ...
##  $ Assault : int  236 263 294 190 276 204 110 238 335 211 ...
##  $ UrbanPop: int  58 48 80 50 91 78 77 72 80 60 ...
##  $ Rape    : num  21.2 44.5 31 19.5 40.6 38.7 11.1 15.8 31.9 25.8 ...
prcomp(USArrests, scale = TRUE)
## Standard deviations (1, .., p=4):
## [1] 1.5748783 0.9948694 0.5971291 0.4164494
## 
## Rotation (n x k) = (4 x 4):
##                 PC1        PC2        PC3         PC4
## Murder   -0.5358995  0.4181809 -0.3412327  0.64922780
## Assault  -0.5831836  0.1879856 -0.2681484 -0.74340748
## UrbanPop -0.2781909 -0.8728062 -0.3780158  0.13387773
## Rape     -0.5434321 -0.1673186  0.8177779  0.08902432
prcomp(~ Murder + Assault + Rape, data = USArrests, scale = TRUE)
## Standard deviations (1, .., p=3):
## [1] 1.5357670 0.6767949 0.4282154
## 
## Rotation (n x k) = (3 x 3):
##                PC1        PC2        PC3
## Murder  -0.5826006  0.5339532 -0.6127565
## Assault -0.6079818  0.2140236  0.7645600
## Rape    -0.5393836 -0.8179779 -0.1999436
plot(prcomp(USArrests))

summary(prcomp(USArrests, scale = TRUE))
## Importance of components%s:
##                           PC1    PC2     PC3     PC4
## Standard deviation     1.5749 0.9949 0.59713 0.41645
## Proportion of Variance 0.6201 0.2474 0.08914 0.04336
## Cumulative Proportion  0.6201 0.8675 0.95664 1.00000
biplot(prcomp(USArrests, scale = TRUE))

Análisis factorial

El análisis factorial (AF) es tambien una técnica estadística de reducción de datos usada para explicar las correlaciones entre las variables observadas en términos de un número menor de variables no observadas llamadas factores.

El análisis factorial exploratorio, se utiliza para tratar de descubrir la estructura interna de un número relativamente grande de variables. La hipótesis a priori del investigador es que pueden existir una serie de factores asociados a grupos de variables. Las cargas de los distintos factores se utilizan para intuir la relación de éstos con las distintas variables. Es el tipo de análisis factorial más común.

El análisis factorial confirmatorio, AFC, trata de determinar si el número de factores obtenidos y sus cargas se corresponden con los que cabría esperar a la luz de una teoría previa acerca de los datos. La hipótesis a priori es que existen unos determinados factores preestablecidos y que cada uno de ellos está asociado con un determinado subconjunto de las variables. El análisis factorial confirmatorio entonces arroja un nivel de confianza para poder aceptar o rechazar dicha hipótesis.

El análisis de componentes principales es el método apropiado de extracción de factores, cuando el interés primordial se centra en la predicción o el número mínimo de factores necesarios para justificar la porción máxima de varianza representada en la serie de variables original, y cuando el conocimiento previo sugiere que la varianza específica y de error representan una porción relativamente pequeña de la varianza total. Por el contrario cuando se pretende identificar las dimensiones latentes o las construcciones representadas en las variables originales y se tiene poco conocimiento de la varianza específica y el error, lo más apropiado es utilizar el método factorial común. Si bien las complicaciones del análisis factorial común han contribuido al análisis generalizado de la técnica de componentes principales.

Para interpretar bien los factores se utiliza una rotación de ejes, ya que las soluciones factoriales no rotadas extraen factores según su orden de importancia. El primer factor tiende a ser un factor general por el que casi toda variable se ve afectada significativamente dando cuenta del mayor porcentaje de varianza. Los métodos de rotación ortogonales mas utilizados son VARIMAX, QUARTIMAX y EQUIMAX.

En la practica no se utiliza un criterio único de los factores a extraer, si bien se utiliza como primera aproximación el gráfico de autovalor para el criterio de contraste de caida, despues de interpretar los factores extraidos hay que valorar su caracter práctico.

Se leen la base de datos Seatbelts, que contiene datos mensuales de conductores de Gran Bretaña, muestos o heridos entre Enero de 1969 y Diciembre de 1983:

library(rela)
Belts <- Seatbelts[,1:7]
summary(Belts)
##  DriversKilled      drivers         front             rear      
##  Min.   : 60.0   Min.   :1057   Min.   : 426.0   Min.   :224.0  
##  1st Qu.:104.8   1st Qu.:1462   1st Qu.: 715.5   1st Qu.:344.8  
##  Median :118.5   Median :1631   Median : 828.5   Median :401.5  
##  Mean   :122.8   Mean   :1670   Mean   : 837.2   Mean   :401.2  
##  3rd Qu.:138.0   3rd Qu.:1851   3rd Qu.: 950.8   3rd Qu.:456.2  
##  Max.   :198.0   Max.   :2654   Max.   :1299.0   Max.   :646.0  
##       kms         PetrolPrice        VanKilled     
##  Min.   : 7685   Min.   :0.08118   Min.   : 2.000  
##  1st Qu.:12685   1st Qu.:0.09258   1st Qu.: 6.000  
##  Median :14987   Median :0.10448   Median : 8.000  
##  Mean   :14994   Mean   :0.10362   Mean   : 9.057  
##  3rd Qu.:17203   3rd Qu.:0.11406   3rd Qu.:12.000  
##  Max.   :21626   Max.   :0.13303   Max.   :17.000

Calculo la matriz de correlaciones:

correl=cor(Belts,use="pairwise.complete.obs") 
correl
##               DriversKilled    drivers      front       rear        kms
## DriversKilled     1.0000000  0.8888264  0.7067596  0.3533510 -0.3211016
## drivers           0.8888264  1.0000000  0.8084114  0.3436685 -0.4447631
## front             0.7067596  0.8084114  1.0000000  0.6202248 -0.3573823
## rear              0.3533510  0.3436685  0.6202248  1.0000000  0.3330069
## kms              -0.3211016 -0.4447631 -0.3573823  0.3330069  1.0000000
## PetrolPrice      -0.3866061 -0.4576675 -0.5392394 -0.1326272  0.3839004
## VanKilled         0.4070412  0.4853995  0.4724207  0.1217581 -0.4980356
##               PetrolPrice  VanKilled
## DriversKilled  -0.3866061  0.4070412
## drivers        -0.4576675  0.4853995
## front          -0.5392394  0.4724207
## rear           -0.1326272  0.1217581
## kms             0.3839004 -0.4980356
## PetrolPrice     1.0000000 -0.2885584
## VanKilled      -0.2885584  1.0000000

La prueba de esfericidad de Bartlett contrasta la hipótesis nula de que la matriz de correlaciones es una matriz identidad, en cuyo caso no existirían correlaciones significativas entre las variables y el modelo factorial no sería pertinente.

library(psych)
cortest.bartlett(Belts)
## R was not square, finding R from data
## $chisq
## [1] 968.7579
## 
## $p.value
## [1] 1.260101e-191
## 
## $df
## [1] 21

Para calcular los componentes principales basados en la matriz de correlaciones:

## Importance of components:
##                           Comp.1    Comp.2    Comp.3     Comp.4    Comp.5
## Standard deviation     1.9203515 1.2141729 0.8489513 0.78140438 0.5754445
## Proportion of Variance 0.5268214 0.2106022 0.1029598 0.08722754 0.0473052
## Cumulative Proportion  0.5268214 0.7374236 0.8403834 0.92761093 0.9749161
##                           Comp.6      Comp.7
## Standard deviation     0.3298839 0.258386585
## Proportion of Variance 0.0155462 0.009537661
## Cumulative Proportion  0.9904623 1.000000000
## 
## Loadings:
##               Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7
## DriversKilled -0.444         0.136  0.521 -0.360 -0.519  0.323
## drivers       -0.480         0.100  0.377         0.505 -0.596
## front         -0.475 -0.198        -0.103  0.417  0.414  0.612
## rear          -0.226 -0.686        -0.322  0.277 -0.406 -0.361
## kms            0.276 -0.627               -0.606  0.365  0.168
## PetrolPrice    0.327 -0.141  0.829  0.277  0.311              
## VanKilled     -0.334  0.259  0.526 -0.627 -0.387

A continuación se realiza un Analisis factorial con la libreria : psych, representando el gráfico de autovalor para el criterio de contraste de caida.

library(psych)
fa.parallel(Belts)
## The estimated weights for the factor scores are probably incorrect.  Try a different factor extraction method.
## Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate =
## rotate, : A loading greater than abs(1) was detected. Examine the loadings
## carefully.
## The estimated weights for the factor scores are probably incorrect.  Try a different factor extraction method.
## Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate =
## rotate, : An ultra-Heywood case was detected. Examine the results carefully

## Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate =
## rotate, : A loading greater than abs(1) was detected. Examine the loadings
## carefully.
## The estimated weights for the factor scores are probably incorrect.  Try a different factor extraction method.
## Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate =
## rotate, : An ultra-Heywood case was detected. Examine the results carefully

## Parallel analysis suggests that the number of factors =  3  and the number of components =  2
Belts.ps <- fa(Belts, fm="minres", nfactors=3, rotate="varimax")
## The estimated weights for the factor scores are probably incorrect.  Try a different factor extraction method.
## Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate =
## rotate, : An ultra-Heywood case was detected. Examine the results carefully
print(Belts.ps) # Factor loadings as $loadings
## Factor Analysis using method =  minres
## Call: fa(r = Belts, nfactors = 3, rotate = "varimax", fm = "minres")
## Standardized loadings (pattern matrix) based upon correlation matrix
##                 MR1   MR3   MR2   h2       u2 com
## DriversKilled  0.81  0.27  0.26 0.80  0.19785 1.4
## drivers        0.88  0.39  0.25 1.00  0.00454 1.6
## front          0.47  0.54  0.71 1.01 -0.01177 2.7
## rear           0.21 -0.12  0.82 0.72  0.27655 1.2
## kms           -0.19 -0.93  0.32 1.00  0.00012 1.3
## PetrolPrice   -0.27 -0.44 -0.20 0.30  0.69660 2.1
## VanKilled      0.30  0.49  0.11 0.34  0.65652 1.8
## 
##                        MR1  MR3  MR2
## SS loadings           1.89 1.83 1.45
## Proportion Var        0.27 0.26 0.21
## Cumulative Var        0.27 0.53 0.74
## Proportion Explained  0.37 0.35 0.28
## Cumulative Proportion 0.37 0.72 1.00
## 
## Mean item complexity =  1.7
## Test of the hypothesis that 3 factors are sufficient.
## 
## The degrees of freedom for the null model are  21  and the objective function was  5.16 with Chi Square of  968.76
## The degrees of freedom for the model are 3  and the objective function was  0.07 
## 
## The root mean square of the residuals (RMSR) is  0.02 
## The df corrected root mean square of the residuals is  0.04 
## 
## The harmonic number of observations is  192 with the empirical chi square  1.85  with prob <  0.6 
## The total number of observations was  192  with Likelihood Chi Square =  13.41  with prob <  0.0038 
## 
## Tucker Lewis Index of factoring reliability =  0.922
## RMSEA index =  0.137  and the 90 % confidence intervals are  0.067 0.212
## BIC =  -2.36
## Fit based upon off diagonal values = 1

Analsis factorial con factanal

Belts.fa<-factanal(Belts,factors=3, rotation="varimax",scores = "Bartlett")
Belts.fa
## 
## Call:
## factanal(x = Belts, factors = 3, scores = "Bartlett", rotation = "varimax")
## 
## Uniquenesses:
## DriversKilled       drivers         front          rear           kms 
##         0.197         0.005         0.005         0.227         0.051 
##   PetrolPrice     VanKilled 
##         0.660         0.645 
## 
## Loadings:
##               Factor1 Factor2 Factor3
## DriversKilled  0.820   0.264   0.244 
## drivers        0.883   0.391   0.249 
## front          0.485   0.531   0.691 
## rear           0.204  -0.121   0.846 
## kms           -0.191  -0.903   0.311 
## PetrolPrice   -0.241  -0.469  -0.249 
## VanKilled      0.301   0.507         
## 
##                Factor1 Factor2 Factor3
## SS loadings      1.915   1.813   1.482
## Proportion Var   0.274   0.259   0.212
## Cumulative Var   0.274   0.533   0.744
## 
## Test of the hypothesis that 3 factors are sufficient.
## The chi square statistic is 11.54 on 3 degrees of freedom.
## The p-value is 0.00914

Análisis factorial con libreria : rela

library(rela)
res <- paf(as.matrix(Belts))
summary(res) # Automatically calculate KMO with MSA, determine the number of factors,
## $KMO
## [1] 0.66899
## 
## $MSA
##                   MSA
## DriversKilled 0.75338
## drivers       0.72295
## front         0.65795
## rear          0.40214
## kms           0.53517
## PetrolPrice   0.83413
## VanKilled     0.88197
## 
## $Bartlett
## [1] 968.76
## 
## $Communalities
##               Initial Communalities Final Extraction
## DriversKilled               0.79824          0.66181
## drivers                     0.87486          0.85141
## front                       0.87225          0.91049
## rear                        0.76216          0.69016
## kms                         0.67184          0.89761
## PetrolPrice                 0.36120          0.29405
## VanKilled                   0.36294          0.34590
## 
## $Factor.Loadings
##                   [,1]      [,2]
## DriversKilled  0.80934 -0.082307
## drivers        0.92263 -0.012744
## front          0.92638 -0.228709
## rear           0.41852 -0.717642
## kms           -0.53335 -0.783039
## PetrolPrice   -0.53096 -0.110142
## VanKilled      0.55008  0.208105
## 
## $RMS
## [1] 0.050492
# calculate chi-square of Bartlett's sphericity test, communalities and
# factor loadings. Communalities are 1 minus uniquenesses.
barplot(res$Eigenvalues[,1]) # First column of eigenvalues.

resv <- varimax(res$Factor.Loadings) # Varimax rotation is possible later.
print(resv)
## $loadings
## 
## Loadings:
##               [,1]   [,2]  
## DriversKilled  0.773 -0.254
## drivers        0.898 -0.211
## front          0.856 -0.422
## rear           0.255 -0.791
## kms           -0.689 -0.650
## PetrolPrice   -0.542       
## VanKilled      0.582       
## 
##                 [,1]  [,2]
## SS loadings    3.309 1.343
## Proportion Var 0.473 0.192
## Cumulative Var 0.473 0.664
## 
## $rotmat
##         [,1]     [,2]
## [1,] 0.97667 -0.21474
## [2,] 0.21474  0.97667
barplot(sort(colSums(loadings(resv)^2),decreasing=TRUE)) # screeplot using rotated SS loadings.

scores <- as.matrix(Belts) %*% as.matrix(resv$loadings) # Get factor scores in a simple manner.

El índice KMO compara los valores de las correlaciones entre las variables y sus correlaciones parciales. Si el índice KMO está próximo a 1, el ACP se puede hacer. Si el índice es bajo (próximo a 0), el ACP no será relevante.

Análisis cluster

Introducción

El análisis cluster (AC) es un conjunto de técnicas multivariantes cuyo principal propósito es agrupar objetos basándose en las características que poseen. El AC clasifica los objetos en clases ó conglomerados de tal forma que cada objeto sea parecido a los que hay en el conjunto de su conglomerado. Los conglomerados resultantes tendrán que tener un alto grado de homogeneidad interna (dentro del conglomerado) y de heterogeneidad externa (entre conglomerados).

Para realizar un AC se necesita una medida de similitud, de correspondencia, o parecido entre objetos que van a ser analizados. Esta puede ser medida de varias formas, pero en el AC, tres son los métodos de simulitud que se utilizan: medidas de correlación (la matriz de correlación entre variables), medidas de distancia (la más utilizada es la distancia euclidea) y las medidas de asociación, cuando se tienen datos no métricos.

Una vez se tiene una matriz de simulitud calculada, hay que realizar el proceso de partición de los datos, para ello hay que decidir el algoritmo de aglomeración utilizado en la formación de conglomerados, para despues tomar una decisión acerca del número de conglomerados que se van a formar. Los algoritmos de obtención de conglomerados se subdividen en métodos jerárquicos (encadenamiento simple, encadenamiento completo, encadenamiento médio, método de Ward, método del centroide), de los no jerarquicos (umbral secuencial, umbral paralelo, optimización), o una combinación de los dos.

El dendrograma ó gráfico en forma de arbol, es una herramienta visual que ayudar a decidir el número de conglomerados que podrían representar mejor la estructura de los datos.

Para determinar el número final de conglomerados a formar o regla de parada, no hay un procedimiento determinado, ha de decidirlo el investigador en la fase de interpretación de los datos.

Cluster Jerarquico

R incluye la función hclust para la realización de clasificaciones jerárquicas.

Una clasificación precisa que el investigador le pase dos opciones básicas a R:

●Medida de similitud/disimilitud para calcular la matriz de distancias

●Estrategia de fusión

Posteriormente, en la interpretación de los datos el investigador deberá decidir el punto de corte y, con ello, el número de grupos.

Se puede usar la función nativa de R dist para calcular la matriz de distancias; por defecto el programa usa la distancia euclídea (euclidean), siendo posible elegir varias más con la opción method=“”: maximum, manhattan, canberra, binary, minkowski.

Las estrategias de fusión disponibles en hclust son las siguientes:

●“ward”: la distancia entre dos conglomerados es la suma de los cuadrados entre dos conglomerados sumados para todas las variables. En cada paso del procedimiento de aglomeración se minimiza la suma de cuadrados dentro del conglomerado para todas las particiones. Tiende a combinar los conglomerados con un número reducido de dimensiones.

●“single”: distancia mínima o vecino más próximo, encuentra dos objetos separados por la distancia más corta y los coloca en un primer conglomerado. A continuación se encuetra la distancia más corta, y se une a el un tercer objeto o se forma un nuevo conglomerado de dos miembros. El proceso continua hasta que todos los objetos están en un conglomerado.

●“complete”: Parecido al “single”(simple) pero utilizando como criterio de agregación la distancia máxima (aproximación al vecino más lejano).

●“average”: comienza igual que los métodos “single” y “complete”, pero el método de aproximación es la distancia media entre todos los individuos de un conglomerado y los de otro.

●“mcquitty”

●“median”

●“centroid”: la distancia entre dos conglomerados es la distancia (euclidia) entre centroides.

Realizamos un Cluster Jerarquico con las tasas de crimen en USA (USArrests), utilizando una medida de similaridad de distancias euclideas y el método de agrupación de encadenamiento medio:

require(graphics)
str(USArrests)
## 'data.frame':    50 obs. of  4 variables:
##  $ Murder  : num  13.2 10 8.1 8.8 9 7.9 3.3 5.9 15.4 17.4 ...
##  $ Assault : int  236 263 294 190 276 204 110 238 335 211 ...
##  $ UrbanPop: int  58 48 80 50 91 78 77 72 80 60 ...
##  $ Rape    : num  21.2 44.5 31 19.5 40.6 38.7 11.1 15.8 31.9 25.8 ...
hc <- hclust(dist(USArrests), "ave")
plot(hc)

plot(hc, hang = -1)

Se realiza el mismo cluster Jerarquico, se establece una parada de 10 cluster (cutree), y se realiza un nuevo cluster con los 10 conglomerados.

hc <- hclust(dist(USArrests)^2, "cen")
memb <- cutree(hc, k = 10)
cent <- NULL
for(k in 1:10){
  cent <- rbind(cent, colMeans(USArrests[memb == k, , drop = FALSE]))
}
str(cent)
##  num [1:10, 1:4] 11.47 13.5 9.95 11.5 5.59 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ : NULL
##   ..$ : chr [1:4] "Murder" "Assault" "UrbanPop" "Rape"
hc1 <- hclust(dist(cent)^2, method = "cen", members = table(memb))
opar <- par(mfrow = c(1, 2))
plot(hc,  labels = FALSE, hang = -1, main = "Original Tree")
plot(hc1, labels = FALSE, hang = -1, main = "Re-start from 10 clusters")

par(opar)

Determinar el número de cluster con k-means

Determinar el número de cluster o la regla de parada es una de las tareas más sensibles en el AC, una forma de incluir un test para decidir el numero de conglomerados en utilizar la función k-meas de R, y representar gráficamente el vector la suma de cuadrados entre conglomerados (withniss), para diferentes agrupamientos.

d=dist(USArrests)^2

# Determinar el numero de cluster

# Fija semilla inicial de numeros pseudoaleatorios
# para obtener misma serie aleatoria en ejecuciones.
set.seed(123)
 
# Crea vector "Errores", sin datos
# Crea variable "K_Max" con la cant. maxima de k a analizar
Errores <-NULL
K_Max   <-10
 

# Ejecuta kmeans con diferentes cluster, desde 1 hasta 10
# Luego guarda el error de cada ejecucion en el vector "Errores"
for (i in 1:K_Max)
  {
  Errores[i] <- sum(kmeans(d[-1], centers=i)$withinss)
  }
# Grafica el vector "Errores"
plot(1:K_Max, Errores, type="b", 
              xlab="Cantidad de Cluster", 
              ylab="Suma de error")

Cluster Jerárquico de series temporales con DTW

La agrupación en series temporales consiste, como se apreció en el apartardo anterior, en dividir el data-set de series temporales en grupos basados en similitudes o distancias, de modo que las series temporales de un mismo grupo sean similares.

Para la agrupación de series de tiempo, el primer paso es elaborar una métrica de distancia ó similitud apropiada, y luego, en el segundo paso, utilizar técnicas de agrupación existentes, como k-means, clustering jerárquico, clustering basado en densidad o clustering subespacial.

Para encontrar estructuras de agrupamiento, de series temporales una metrica adecuada es la distancia Dynamic Time Warping (DTW), que encuentra la alineación óptima entre dos series temporales.

El Alineamiento Temporal Dinámico (Dynamic Time Warping, DTW), es una técnica desarrollada en el ambito de reconocimiento de la voz, cuyo objetivo es comparar dos emisiones vocales aisladas a fin de determinar si corresponden o no a la misma palabra.

De manera resumida, el DTW se diseño para comprobar si existe una sincronización temporal (alineamiento temporal) entre dos grupos fónicos. La falta de alineamiento, por lo general, no obedece a una ley fija (p. e., un retardo constante), sino que se da de forma heterogénea, produciéndose así variaciones localizadas que aumentan o disminuyen la duración del tramo de análisis.

El problema de calcular una función de alineamiento se resuelve en el DTW mediante técnicas de programación dinámica, motivo por el que se conoce a esta técnica como alineamiento temporal dinámico.

En el ejemplo, se utiliza un conjunto de datos de la serie temporal de gráficos de control sintético, que contiene 600 ejemplos de gráficos de control. Cada gráfico de control es una serie de tiempo con 60 valores. Hay seis clases:

  1. 1-100 Normal,
  2. 101-200 Cíclico,
  3. 201-300 Tendencia creciente,
  4. 301-400 Tendencia decreciente,
  5. 401-500 Cambio hacia arriba y
  6. 501-600 Descenso cambio.

Los datos y el script se ha obtenido de :http://www.rdatamining.com/examples/time-series-clustering-classification.

sc <- read.table("http://kdd.ics.uci.edu/databases/synthetic_control/synthetic_control.data", header=F, sep="")
# randomly sampled n cases from each class, to make it easy for plotting
n <- 10
s <- sample(1:100, n)
idx <- c(s, 100+s, 200+s, 300+s, 400+s, 500+s)
sample2 <- sc[idx,]
observedLabels <- c(rep(1,n), rep(2,n), rep(3,n), rep(4,n), rep(5,n), rep(6,n))
# compute DTW distances
library(dtw)
## Loading required package: proxy
## 
## Attaching package: 'proxy'
## The following objects are masked from 'package:stats':
## 
##     as.dist, dist
## The following object is masked from 'package:base':
## 
##     as.matrix
## Loaded dtw v1.18-1. See ?dtw for help, citation("dtw") for use in publication.
distMatrix <- dist(sample2, method="DTW")
# hierarchical clustering
hc <- hclust(distMatrix, method="average")
plot(hc, labels=observedLabels, main="")

Cluster Jerárquico con span-tree

Si al clasificar un conjunto de datos no encontramos una decisión adecuada para establecer una regla de parada, y obtenemos como resultado un conjunto inclasificable, es util combinar la metodología de jerarquía y topología de los árboles de expansión mínima (MST, Minimal Spanning Tree).

ado un grafo, su árbol mínimo generador (o árbol de peso mínimo o árbol mínimo de expansión) es un árbol que pasa por todos los vértices y que la suma de sus aristas es la de menor peso.

La gráfica siguiente ilustra el proceso, se trata de un arbol de expansión de 7 nodos (acciones), con las medidas de distancia que hay entre ellos. Partiendo de la distancia menor, 29, entre el nodo 8 y 4, se inicia el arbol, que en el nodo 4, encuentra la segunda minima distancia en su enlace con el nodo 2, 31, y en el nodo 2, la menor distancia la encuentra con el nodo 5, operando de esta manera hasta tener enlazados todos los nodos. Como se puede observar en la figura, un nodo puede tener, más de un enlace con otros nodos, tal y como ocurre en el gráfico, nodo 3.

Curva ROC

Curva ROC

Basandonos en el trabajo pionero de Mantegna (1999) sobre tasas de retornos en los mercados financieros, realizamos un ejercicio clasificatorio con la función span-tree, de la librería vegan, que ordena y establece dependencias entre diversos activos financieros utilizando datos de rendimientos mensuales de activos procedentes de Berndt’s The Practice of Econometrics: http://web.stanford.edu/~clint/berndt/

# read prices from csv file
bolsa.df = load("berndt.RData")
colnames(berndt.df)
##  [1] "Año"    "MOBIL"  "TEXACO" "IBM"    "DEC"    "DATGEN" "CONED" 
##  [8] "PSNH"   "WEYER"  "BOISE"  "MOTOR"  "TANDY"  "PANAN"  "DELTA" 
## [15] "CONTIL" "CITGRP" "GERBER" "GENMIL" "MARKET" "RKFREE"
#[1] "Año"    "MOBIL"  "TEXACO" "IBM"    "DEC"    "DATGEN" "CONED"  "PSNH"   "WEYER"  "BOISE" 
#[11] "MOTOR"  "TANDY"  "PANAN"  "DELTA"  "CONTIL" "CITGRP" "GERBER" "GENMIL" "MARKET" "RKFREE

# create zooreg object - regularly spaced zoo object


library(zoo)
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
berndt.z = zooreg(berndt.df[,-1], start=c(1978, 1), end=c(1987,12),frequency=12)
index(berndt.z) = as.yearmon(index(berndt.z))

start(berndt.z)
## [1] "ene. 1978"
end(berndt.z)
## [1] "dic. 1987"
nrow(berndt.z)
## [1] 120
# note: coredata() function extracts data from zoo object
returns.mat = as.matrix(coredata(berndt.z))

# create the correlation coefficients

coef.corr <- cor(returns.mat)
coef.d <- (1-coef.corr^2) # compute distance (Mantegna, 1998)

# hierarchical cluster whir hclust

d <- as.dist(as.matrix(coef.d)) # find distance matrix 

#Function spantree finds a minimum spanning tree (MTA) connecting all points, but disregarding dissimilarities that are at or above the threshold or NA.


library(vegan)
## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.4-3
tr <- spantree(d)
## Add tree to a metric scaling 
plot(tr, cmdscale(d), type = "t")

## Find a configuration to display the tree neatly
plot(tr, type = "t")
## Initial stress        : 0.10907
## stress after   2 iters: 0.06448

## Depths of nodes
depths <- spandepth(tr)
plot(tr, type = "t", label = depths)
## Initial stress        : 0.10907
## stress after   2 iters: 0.06448

## Plot as a dendrogram
cl <- as.hclust(tr)
plot(cl)