Laboratorio 10.4 - Analisis de componentes principales
Cesar Tinoco - 13003387
Lab 1
A continuacion realizaremos PCA en el conjunto de datos de USArrests, que forma parte del paquete base r. Las filas del conjunto de datos contienen los 50 estados, en orden alfabetico.
states =row.names(USArrests )
states [1] "Alabama" "Alaska" "Arizona" "Arkansas" "California" "Colorado"
[7] "Connecticut" "Delaware" "Florida" "Georgia" "Hawaii" "Idaho"
[13] "Illinois" "Indiana" "Iowa" "Kansas" "Kentucky" "Louisiana"
[19] "Maine" "Maryland" "Massachusetts" "Michigan" "Minnesota" "Mississippi"
[25] "Missouri" "Montana" "Nebraska" "Nevada" "New Hampshire" "New Jersey"
[31] "New Mexico" "New York" "North Carolina" "North Dakota" "Ohio" "Oklahoma"
[37] "Oregon" "Pennsylvania" "Rhode Island" "South Carolina" "South Dakota" "Tennessee"
[43] "Texas" "Utah" "Vermont" "Virginia" "Washington" "West Virginia"
[49] "Wisconsin" "Wyoming" Las columnas del conjunto de datos contienen las cuatro variables.
names(USArrests )[1] "Murder" "Assault" "UrbanPop" "Rape" Primero examinamos brevemente los datos. Notamos que las variables tienen medios muy diferentes.
apply(USArrests , 2, mean) Murder Assault UrbanPop Rape
7.788 170.760 65.540 21.232 Tenga en cuenta que la funcion apply() permite aplicar la funcion mean(), a cada fila o columna del conjunto de datos. La segunda entrada indica si se desea calcular la media de las filas, 1, o las columnas, 2. Vemos que hay en promedio tres veces mas violaciones que asesinatos, y mas de ocho veces mas asaltos que violaciones. Tambien podemos examinar las variaciones de las cuatro variables usando la función apply().
apply(USArrests , 2, var) Murder Assault UrbanPop Rape
18.97047 6945.16571 209.51878 87.72916 No es sorprendente que las variables tambien tengan variaciones muy diferentes: la variable UrbanPop mide el porcentaje de la poblacion en cada estado que vive en un area urbana, que no es un numero comparable al numero de violaciones en cada estado por cada 100,000 personas. Si no logramos escalar las variables antes de realizar PCA, entonces la mayoria de los componentes principales que observamos estaria impulsada por la variable Assault, ya que tiene, por mucho, la media y la varianza mas grande. Por lo tanto, es importante estandarizar las variables para que tengan una media de cero y una desviacion estandar antes de realizar la PCA.
Ahora realizamos el analisis de componentes principales utilizando la funcion prcomp(), que es una de las varias funciones en R que realizan PCA.
pr.out =prcomp (USArrests , scale =TRUE)Por defecto, la funcion prcomp() centra las variables para que tengan media cero. Al usar la opcion scale = TRUE, escalamos las variables para tener una desviacion estandar de uno. La salida de prcomp() contiene una cantidad de cantidades utiles.
names(pr.out )[1] "sdev" "rotation" "center" "scale" "x" El centro y los componentes de escala corresponden a las medias y desviaciones estandar de las variables que se utilizaron para escalar antes de implementar PCA.
pr.out$center Murder Assault UrbanPop Rape
7.788 170.760 65.540 21.232 pr.out$scale Murder Assault UrbanPop Rape
4.355510 83.337661 14.474763 9.366385 La matriz de rotacion proporciona las cargas principales de los componentes; Cada columna de rotacion pr.out$ contiene el vector de carga del componente principal correspondiente.
pr.out$rotation 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.08902432Vemos que hay cuatro componentes principales distintos. Esto es para ser esperado porque en general hay min(n - 1, p) componentes principales informativos en un conjunto de datos con n observaciones y p variables. Usando la funcion prcomp(), no necesitamos multiplicar explicitamente los datos por el componente principal de carga de vectores con el fin de obtener el Componente principal en los vectores de puntuacion. Mas bien, la matriz x de 50×4 tiene como columnas los vectores de puntuacion del componente principal. Es decir, la columna kth es el kth vector de puntuacion del componente principal.
dim(pr.out$x )[1] 50 4Podemos trazar los dos primeros componentes principales de la siguiente manera:
biplot (pr.out , scale =0)
El argumento scale = 0 para biplot() asegura que las flechas se escalan para representar las cargas; Otros valores para la escala dan biplots ligeramente diferentes con diferentes interpretaciones. Observe que esta figura es una imagen reflejada de la Figura 10.1. Recordar que los componentes principales solo son unicos hasta un cambio de signo, por lo que podemos reproducir la Figura 10.1 realizando algunos pequeños cambios:
pr.out$rotation=-pr.out$rotation
pr.out$x=-pr.out$x
biplot (pr.out , scale =0)
La funcion prcomp() tambien genera la desviacion estandar de cada componente principal. Por ejemplo, en el conjunto de datos de USArrests, podemos acceder a estas desviaciones estandar de la siguiente manera:
pr.out$sdev[1] 1.5748783 0.9948694 0.5971291 0.4164494La varianza explicada por cada componente principal se obtiene al cuadrar estos:
pr.var =pr.out$sdev ^2
pr.var[1] 2.4802416 0.9897652 0.3565632 0.1734301Para calcular la proporcion de varianza explicada por cada componente principal, 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.62006039 0.24744129 0.08914080 0.04335752Vemos que el primer componente principal explica el 62.0% de la varianza en los datos, el siguiente componente principal explica el 24.7% de la varianza, y asi sucesivamente. Podemos trazar el PVE explicado por cada componente, asi como el PVE acumulado, de la siguiente manera:
plot(pve, xlab = "Principal Component", ylab = "Proportion of Variance Explained", ylim = c(0,1), type = 'b')
plot(cumsum(pve), xlab = "Principal component", ylab = "Cumulative Proportion of Variance Explained", ylim = c(0,1), type = 'b')
El resultado se muestra en la Figura 10.4. Tenga en cuenta que la función cumsum() calcula la suma acumulativa de los elementos de un vector numerico. Por ejemplo:
a=c(1,2,8,-3)
cumsum (a)[1] 1 3 11 8