Reduciendo la dimensión: Componentes Principales
Gabriel Martos
Técnicas de Investigación - Septiembre, 2014
En esta clase retomaremos el repaso de R, y nos adentraremos en el Análisis de Componentes Principales (PCA por sus siglas en Inglés).
Analisis Multivariante en R
Vamos a trabajar con el conjunto de datos iris. Comenzaremos a analizar la estructura de los datos con los comandos: dim() y names().
data(iris)
dim(iris) # Número de filas y de columnas del conjunto de datos.## [1] 150 5names(iris) # Nombre de las columnas## [1] "Sepal.Length" "Sepal.Width" "Petal.Length" "Petal.Width"
## [5] "Species"head(iris, 3) # Para visualizar las primeras 3 observaciones de cada columna.## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 5.1 3.5 1.4 0.2 setosa
## 2 4.9 3.0 1.4 0.2 setosa
## 3 4.7 3.2 1.3 0.2 setosastr(iris) # Otro comando útil que nos permite ver los primeros valores de cada variable y sus atributos## 'data.frame': 150 obs. of 5 variables:
## $ Sepal.Length: num 5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
## $ Sepal.Width : num 3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
## $ Petal.Length: num 1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
## $ Petal.Width : num 0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
## $ Species : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...Explorando las variables de manera individual
Con el comando table() generamos tablas de frecuencias absolutas:
table(iris$Species)##
## setosa versicolor virginica
## 50 50 50Con la función summary() podemos obtener los estadísticos descriptivos básicos para todas las variables (columnas) de nuestra matriz de datos.
summary(iris)## Sepal.Length Sepal.Width Petal.Length Petal.Width
## Min. :4.30 Min. :2.00 Min. :1.00 Min. :0.1
## 1st Qu.:5.10 1st Qu.:2.80 1st Qu.:1.60 1st Qu.:0.3
## Median :5.80 Median :3.00 Median :4.35 Median :1.3
## Mean :5.84 Mean :3.06 Mean :3.76 Mean :1.2
## 3rd Qu.:6.40 3rd Qu.:3.30 3rd Qu.:5.10 3rd Qu.:1.8
## Max. :7.90 Max. :4.40 Max. :6.90 Max. :2.5
## Species
## setosa :50
## versicolor:50
## virginica :50
##
##
## La media, mediana y el rango pueden obtenerse igualmente con las funciones: mean(), median() y range() respectivamente. Los cuartiles y percentiles con la función quantile(). Las medidas de dispersión (varianzas y desviación estandard) se obtienen así:
attach(iris) # Para trabajar con los nombres de las columnas del data set "iris"
var(iris$Sepal.Length) # varianza de la variable Sepal.Length## [1] 0.6857sd(iris$Sepal.Length) # desviación estandard de la variable Sepal.Length## [1] 0.8281Para conocer la distribución de cada una de las variables podemos utilizar hitogramas (variables numéricas continuas), gráficos de barras (variables numéricas discretas y varaibles cualitativas) y gráficos de pastel:
hist(iris$Sepal.Length) # Histograma
plot(density(iris$Sepal.Length), main= 'Estimación de la Densidad') # Para estimar la densidad (una versión "suave"" del histograma)
barplot(table(iris$Species)) # Gráfico de barras
pie(table(iris$Species)) # Gráfico de torta
Mientras que las medidas de relación (covarianza y correlación) se obtienen con las funciones cov() y cor() respectivamente:
cov(iris$Sepal.Length, iris$Sepal.Width) # co-varianza entre las variables Sepal.Length y Sepal.Width## [1] -0.04243cor(iris$Sepal.Length, iris$Sepal.Width) # co-relación entre las variables Sepal.Length y Sepal.Width## [1] -0.1176Si quisieramos computar la matriz de co-varianzas y de correlación:
cov(iris[,1:4]) # Matriz de covarianzas para las columnas 1 a 4## Sepal.Length Sepal.Width Petal.Length Petal.Width
## Sepal.Length 0.68569 -0.04243 1.2743 0.5163
## Sepal.Width -0.04243 0.18998 -0.3297 -0.1216
## Petal.Length 1.27432 -0.32966 3.1163 1.2956
## Petal.Width 0.51627 -0.12164 1.2956 0.5810cor(iris[,1:4]) # Matriz de correlaciones para las columnas 1 a 4## Sepal.Length Sepal.Width Petal.Length Petal.Width
## Sepal.Length 1.0000 -0.1176 0.8718 0.8179
## Sepal.Width -0.1176 1.0000 -0.4284 -0.3661
## Petal.Length 0.8718 -0.4284 1.0000 0.9629
## Petal.Width 0.8179 -0.3661 0.9629 1.0000Una manera gráfica de abordar las relaciones entre variables es utilizando una matriz de scatterplots:
plot(iris[,1:4], pch='*',col=c("red", "green", "blue")[unclass(iris[,5])])
pairs(iris[1:4], main = "Datos Iris", pch = 21, bg = c("red", "green3", "blue")[unclass(iris$Species)], lower.panel=NULL, labels=c("SL","SW","PL","PW"), font.labels=2, cex.labels=4.5) 
library(car) # Una alternativa un tanto mas vistosa
scatterplotMatrix(iris[,1:4], groups=iris$Species)
Si queremos comparar la distribución de una variable entre grupos, podemos utilizar un boxplot:
boxplot(Sepal.Length ~ Species, col = 'red', main= 'Sepal Lenght por Especie')
Otras herramientas de exploración de datos
Un scatterplot en 3 dimensiones y un “mapa de calor” de los puntos (para ello utilizamos una función que calcula la distancia o “discimilaridad” entre las observaciones)
library(scatterplot3d)
scatterplot3d(iris$Petal.Width, iris$Sepal.Length, iris$Sepal.Width, color=c("red", "green", "blue")[unclass(iris[,5])])
distMatrix <- as.matrix(dist(iris[,1:4]))
heatmap(distMatrix)
Por útlimo, la librería lattice nos permite obtener gráficos de paralelos, donde podemos observar la distribución de las variables numéricas en función de algún atributo (en este caso las especies)
library(lattice)
parallelplot(~iris[1:4] | Species, data=iris)
Finalmente, si queremos guardar los gráficos en formato PDF o PS, simplemente indicamos antes de cada gráfica el formato específico con el que salvaremos el gráfico y la ruta y nombre del fichero:
# Guardar figura como un fichero PDF
pdf("C:/.../myPlot.pdf")
plot(...)
graphics.off()
# Guardar figura como un fichero postscript
postscript("C:/.../myPlot2.ps")
plot(...)
graphics.off()Ejercicio en clase: utiliza el conjunto de datos _Europa_ (ver detalle de la lectura a continuación) y empleando las herramientas gráficas antes expuestas analiza la distribución de las variables involucradas en este conjunto de datos. Europa=read.csv("http://www.instantr.com/wp-content/uploads/2013/01 /europe.csv", header=TRUE)Analisis de Componentes Principales
Con el desarrollo de la tecnología ha resultado muy sensillo para las empresas, los bancos, el sector público, la extracción y almacenamiento de cantidades significativas de datos. Estos datos se organizan habitualmente en una matriz \(X_{n \times p}\), donde n es el número de filas (el número de observaciones, por ejemplo el número de compradores, el número de meses, dias u horas en que observamos una variable, etc) y p representa el número de columnas (el número de variables o atributos que observamos). El objetivo del PCA es sintetizar la información contenida en las p columnas de la matriz de datos \(X\) (es decir, las variables o atributos en cuestión) tal que la dimensión del problema se reduce.
El PCA se basa en un procedimiento matemático que utiliza una transformación ortogonal (en terminos estadísticos: ortogonalidad = independencia) tal que transformamos el conjunto de variales originales (probalemente correlacionadas) \(x_1 , x_2 , \dots , x_p\) (las columnas de la matriz \(X\)) en un conjunto de nuevas variables que llamaremos componentes principales, y cuya principal característica es que estan incorreladas. El número de componentes principales, que denotaremos por \(k\), generalmente es menor o igual que \(p\) (habitualmente \(k\ll p\)).
Esta transformación se lleva a cabo de manera tal que la primera componente principal (que es una combinación lineal de las \(p\) varaibles originales) tiene la máxima varianza entre todas las componentes. Y cada componente adicional presenta una variabilidad que resultará menor o igual a la componente precedente, dada la restricción de que cada componente adicional resulta ortogonal (independiente) a las componentes precedentes.
El PCA es una técnica habitual para analizar información de alta dimensión, y se lleva a cabo a través de la descomposición en valores singulares de la matriz de varianzas y covarianzas (una matriz de \(p\times p\)) o a través de la descomposición en valores singulares de la matriz de correlacion (es decir, despúes de centrar y normalizar cada columna de la matriz de varianzas y covarianzas).
En estas notas de clase omitiremos los detalles matemáticos de la mencionada descomposición, y solo nos enfocaremos en resolver el problema desde una perspectiva práctica.
Comandos útiles para llevar adelante el PCA en R
La función de R que nos permite realizar el PCA se llama princomp:
princomp(dataset)
Es importante que el conjunto de datos solo contenga variables numéricas, puesto que para el cálculo de las componentes necesitamos estimar una matriz de correlación o de covarianzas. Si hubiera variables NO numéricas en el conjunto de datos, debemos excluirlas con el comando [ , ].
La función princomp nos mostrará las varianzas de cada una de las componentes. Sin embargo, existen más resultados que involucran al PCA y no son visualizados por defecto, entre ellos las cargas y los factores. Para organizar correctamente el output de este análisis podemos guadar todas los elementos del output en un objeto como se muestra a continuación:
Mi.modelo.PCA<-princomp(dataset)
Una vez guardado los resultados en un objeto, poderemos recuperar todos los outputs de manera conveniente. Por ejemplo, puedes utilizar la función summary() para visualizar las varianzas de cada componente, las cargas y los factores:
summary(Mi.modelo.PCA)
Para visualizar solo las cargas: > Mi.modelo.PCA$loadings
Para visualizar los scores: > Mi.modelo.PCA$scores
Ejemplo I: PCA datos Económicos de Europa
Europa=read.csv("http://www.instantr.com/wp-content/uploads/2013/01/europe.csv", header=TRUE)
head(Europa, 3)## Country Area GDP Inflation Life.expect Military Pop.growth
## 1 Austria 83871 41600 3.5 79.91 0.8 0.03
## 2 Belgium 30528 37800 3.5 79.65 1.3 0.06
## 3 Bulgaria 110879 13800 4.2 73.84 2.6 -0.80
## Unemployment
## 1 4.2
## 2 7.2
## 3 9.6Como podemos ver, este conjunto de datos tiene 1 variable categórica (Country) y las restantes numéricas. Por ello deberíamos quitar la primera columna para poder llevar adelante el análisis:
Europa2=Europa[,2:8]
head(Europa2, 3)## Area GDP Inflation Life.expect Military Pop.growth Unemployment
## 1 83871 41600 3.5 79.91 0.8 0.03 4.2
## 2 30528 37800 3.5 79.65 1.3 0.06 7.2
## 3 110879 13800 4.2 73.84 2.6 -0.80 9.6A continuación utilizamos la función prcomp() para obtener las componentes principales:
pca.Europa2 <- prcomp(Europa2,scale=FALSE) # Scale = False -> utilizamos la matriz de COVARIANZA para obtener las componentes! Veamos que ocurre:
summary(pca.Europa2)## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 1.66e+05 1.44e+04 4.02 2.32 0.882 0.747 0.241
## Proportion of Variance 9.93e-01 7.47e-03 0.00 0.00 0.000 0.000 0.000
## Cumulative Proportion 9.93e-01 1.00e+00 1.00 1.00 1.000 1.000 1.000plot(pca.Europa2)
La primera componente captura ‘casi toda’ la variabilidad de los datos! Esto quiere decir que podríamos reducir las \(7\) variables originales a una sola variable (componente principal) manteniendo (prácticamente) constante la cantidad de información disponible con respecto al conjunto de datos originales.
Es esto normal? Porque creen que ocurre esto?
Lo que ocurre es que estamos utilizando variables medidas en escalas diferentes (area en \(m^2\), GDP en millones, y el resto de variables en porcentajes!). Cuando ocurre esto, un subconjunto de variables (en este caso GDP y Area) tienden a ‘monopolizar’ la variabilidad en los datos! Pero esto es debido exclusivamente al efecto de ‘escala’ con el que hemos medido los datos. Por ello, en vez de utilizar la matriz de covarianzas para hacer PCA, se utiliza la matriz de correlación:
pca.Europa2 <- prcomp(Europa2,scale=T) # Scale = True -> utilizamos la matriz de CORRELACIÓN para obtener las componentes!
#Veamos que ocurre:
summary(pca.Europa2)## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 1.796 1.090 1.031 0.878 0.6766 0.4107 0.354
## Proportion of Variance 0.461 0.170 0.152 0.110 0.0654 0.0241 0.018
## Cumulative Proportion 0.461 0.631 0.782 0.893 0.9579 0.9820 1.000plot(pca.Europa2)
Al homogeneizar la ‘escala’ en la que hemos medido las variables, la distribución de la variabilidad entre las com ponentes parece más racional. Podemos ver más elementos del output del PCA:
pca.Europa2$sdev # Varianza de cada componente.## [1] 1.7964 1.0896 1.0311 0.8777 0.6766 0.4107 0.3545pca.Europa2$rotation # cargas de cada componente.## PC1 PC2 PC3 PC4 PC5 PC6
## Area 0.1249 0.17287 -0.8982967 0.04485 -0.32402 0.1901
## GDP -0.5005 0.13014 -0.0839558 -0.08426 0.39063 0.6387
## Inflation 0.4065 0.36966 -0.1981947 0.16469 0.68950 -0.3239
## Life.expect -0.4829 -0.26525 -0.2460825 0.02677 -0.10179 -0.6064
## Military 0.1881 -0.65827 -0.2436794 -0.56237 0.36815 0.0356
## Pop.growth -0.4757 -0.08262 -0.1636972 0.39246 0.34787 -0.1209
## Unemployment 0.2717 -0.55320 -0.0005001 0.70197 0.01016 0.2597
## PC7
## Area 0.06664
## GDP -0.39741
## Inflation -0.22670
## Life.expect -0.50703
## Military 0.13731
## Pop.growth 0.67115
## Unemployment -0.24466head(pca.Europa2$x, 3) # Matriz de datos (solo primeras 3 filas) con las componentes (en columnas las nuevas variables).## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## [1,] -1.0623 1.2472 0.5055 -0.4064 0.001229 -0.374667 -0.56518
## [2,] -0.6688 0.4085 0.6749 -0.2783 0.268590 -0.372298 -0.47230
## [3,] 2.5629 -0.2648 0.6118 -1.3079 0.001911 0.003763 -0.02468Como interpretar, por ejemplo, la primera componente?\[y_{1,i} = 0.1248739 Area_i - 0.5005059GDP_i + ... + 0.2716558 Unemployment_i \quad i=1,...,n\]
Algunas preguntas para reflexionar en clase: Con que criterio seleccionar el número de componentes que vamos a utilizar finalmente?
Ejemplo II: Evolución mensual de la tasa de paro en USA
Aquí la matriz de datos \(X\) tiene \(n = 50\) filas (cada fila representa uno de los \(50\) estados de américa del norte) y \(p = 416\) columnas (observaciones mensuales sobre la tasa de paro). El objetivo es reducir la información de los \(416\) meses a un conjunto más reducido de variables (reduciendo así la dimensión de los datos).
De esta manera, podremos visualizar los datos (es decir que utilizaremos a lo sumo 2 o 3 dimensiones) , de manera tal que podamos descubrir si existen ‘clusters’ de estados (grupos de estados que tengan comportamientos similares respecto de la evolución de la tasa de paro).
library(cluster) ## Vamos a utilizar este paquete para hacer 'clusters', un tema que estudiaremos en detalle mas adelante.
states=c("AL","AK","AZ","AR","CA","CO","CT","DE","FL","GA","HI","ID", "IL","IN","IA","KS","KY","LA","ME","MD","MA","MI","MN","MS","MO", "MT","NE","NV","NH","NJ","NM","NY","NC","ND","OH","OK","OR","PA", "RI","SC","SD","TN","TX","UT","VT","VA","WA","WV","WI","WY")
states## [1] "AL" "AK" "AZ" "AR" "CA" "CO" "CT" "DE" "FL" "GA" "HI" "ID" "IL" "IN"
## [15] "IA" "KS" "KY" "LA" "ME" "MD" "MA" "MI" "MN" "MS" "MO" "MT" "NE" "NV"
## [29] "NH" "NJ" "NM" "NY" "NC" "ND" "OH" "OK" "OR" "PA" "RI" "SC" "SD" "TN"
## [43] "TX" "UT" "VT" "VA" "WA" "WV" "WI" "WY"raw <- read.csv("http://www.biz.uiowa.edu/faculty/jledolter/DataMining/unempstates.csv")
raw[1:3,1:5] # Y las filas continuan hasta la dimensión 416## AL AK AZ AR CA
## 1 6.4 7.1 10.5 7.3 9.3
## 2 6.3 7.0 10.3 7.2 9.1
## 3 6.1 7.0 10.0 7.1 9.0## transponemos la matriz 50 filas (cada uno de los estados) y 416 columnas (cada tasa de paro)
rawt=matrix(nrow=50,ncol=416)
rawt=t(raw)
rawt[1:3,1:5] # Y las columnas continuan hasta la dimensión 416## [,1] [,2] [,3] [,4] [,5]
## AL 6.4 6.3 6.1 6.0 6.0
## AK 7.1 7.0 7.0 7.0 7.0
## AZ 10.5 10.3 10.0 9.8 9.6pcaunemp <- prcomp(rawt,scale=FALSE) # computamos las componentes ( porque en este caso scaling = False ???)
head(pcaunemp$sdev) # Aqui vemos las varianzas de las componentes## [1] 24.130 13.213 9.521 6.487 5.567 4.880plot(pcaunemp, main="Varianza de las componentes")
mtext(side=1,"Unemployment: 50 states",line=1,font=2)
pcaunemp$rotation[1:10,1] ## Cargas de la primera componentes. Solo se visualizan las primeras 10 variables (recordemos que aquí tendremos una ecuación con 416 cargas para cada una de las variables en columna de la matriz de datos original)## [1] -0.03261 -0.03209 -0.03204 -0.03153 -0.03165 -0.03201 -0.03212
## [8] -0.03212 -0.03233 -0.03248### Asegurate de entender la manera en la que se contruye la primera componente. Toma un papel y un lápiz e intenta escribir la ecuación de la primera componente.
# Calculamos la tasa media de paro para todos los estos para cada uno de los meses en el análisis.
ave=dim(416)
for (j in 1:416) {
ave[j]=mean(rawt[,j])
} ### El vector 'ave' contiene 416 promedios.
par(mfrow=c(1,2))
## Gráfica de los valores (negativos) de las cargas de la primera componente.
plot(-pcaunemp$rotation[,1], main='Cargas de la primera componente')
## Gráfica de los valores medios de paro para los estados
plot(ave,type="l",ylim=c(3,10),xlab="month",ylab="Evol. de la tasa de paro promedio")
# Pregunta: Que conclusiones obtienes de este análisis????
# Calulemos la correlación entre los factores de la primera componente y las tasas promedio de paro:
abs(cor(ave,pcaunemp$rotation[,1]))## [1] 0.8257unemppc <- predict(pcaunemp)
# A continuación construímos una gráfica de los estados utilizando solo la información de las primeras 2 componentes principales. Llevamos a cabo un análisis de clusters y pintamos los estados que pertenecen a cada cluster de diferentes colores:
set.seed(123456789)
grpunemp3 <- kmeans(rawt,centers=3,nstart=10)
par(mfrow=c(1,1))
plot(unemppc[,1:2],type="n")
text(x=unemppc[,1],y=unemppc[,2],labels=states,col=rainbow(7)
[grpunemp3$cluster])
Que tipo de conclusiones podríamos sacar de este tipo de análisis? Piensa en como aplicar esta técnica con datos de paro (o de cualquier otra serie económica a España).
Ejercicio en clase: con los datos del Heptatlon: 'data("heptathlon", package = "HSAUR")', lleva adelante el análisis de componentes principales (OMITE LA VARIABLE SCORE, que resume la performance global de cada participante).
1) Es necesario normalizar los datos?
2) Cuantas componentes principales debemos utilizar?
3) Se puede interpretar el significado de las primeras componentes?
4) Analiza la correlación entre la primera componente y la variable score (computa el coeficiente de correlación y dibuja ambas variables en un scatterplot). QUe conclusiones obtienies?Ejemplo III: PCA datos decathlon
En este ejemplo vamos a utilizar (solo a modo informativo) otras librerias que complementan el analisis PCA con los Biplots. El conjunto de datos hace referencia a los registros obtenidos por diferentes deportistas (33 atletas) en los juegos olímpicos para diferentes disciplinas. El conjunto de datos pertenece a la librería FactoMineR.
library(FactoMineR)
data(decathlon)
head(decathlon, 3)## 100m Long.jump Shot.put High.jump 400m 110m.hurdle Discus
## SEBRLE 11.04 7.58 14.83 2.07 49.81 14.69 43.75
## CLAY 10.76 7.40 14.26 1.86 49.37 14.05 50.72
## KARPOV 11.02 7.30 14.77 2.04 48.37 14.09 48.95
## Pole.vault Javeline 1500m Rank Points Competition
## SEBRLE 5.02 63.19 291.7 1 8217 Decastar
## CLAY 4.92 60.15 301.5 2 8122 Decastar
## KARPOV 4.92 50.31 300.2 3 8099 Decastarres <- PCA(decathlon,quanti.sup=11:12,quali.sup=13)
plot(res,invisible="quali")
plot(res,choix="var",invisible="quanti.sup")
plot(res,habillage=13)
aa=cbind.data.frame(decathlon[,13],res$ind$coord)
bb=coord.ellipse(aa,bary=TRUE)
plot.PCA(res,habillage=13,ellipse=bb)