En este laboratorio, haremos PCA (Analisis Principal de Componentes) sobre el dataset USArrests, el cual viene integrado en R, estas contienen los 50 estados ordenados alfabeticamente como veremos a continuacion:
states = row.names(USArrests)
states
[1] "Alabama" "Alaska" "Arizona" "Arkansas" "California" "Colorado" "Connecticut"
[8] "Delaware" "Florida" "Georgia" "Hawaii" "Idaho" "Illinois" "Indiana"
[15] "Iowa" "Kansas" "Kentucky" "Louisiana" "Maine" "Maryland" "Massachusetts"
[22] "Michigan" "Minnesota" "Mississippi" "Missouri" "Montana" "Nebraska" "Nevada"
[29] "New Hampshire" "New Jersey" "New Mexico" "New York" "North Carolina" "North Dakota" "Ohio"
[36] "Oklahoma" "Oregon" "Pennsylvania" "Rhode Island" "South Carolina" "South Dakota" "Tennessee"
[43] "Texas" "Utah" "Vermont" "Virginia" "Washington" "West Virginia" "Wisconsin"
[50] "Wyoming"
El dataset contiene 4 columnas que son:
names(USArrests)
[1] "Murder" "Assault" "UrbanPop" "Rape"
Examinemos un poco la data:
apply(USArrests, 2, mean)
Murder Assault UrbanPop Rape
7.788 170.760 65.540 21.232
la funcion apply nos permite iterar sobre cada fila o columna del dataset, el parametro 2 indica explicitamente que iteraremos por columna, si se le envia 1 entonces se itera por fila, y finalmente enviamos la funcion mean para generar el promedio de cada columna. Al analizar la data vemos que hay una varianza relativamente grande, ahora veamos lo que sucede cuando hacemos un analisis de componentes:
Por default la funcion prcomp centra las variables para que estan tengan una media de cero, ademas al utilizar el parametro scale para tener una desviacion estandar de uno.
names(pr.out)
[1] "sdev" "rotation" "center" "scale" "x"
Vemos que en el set de datos las variables center y scale contienen la media y la desviacion estandar
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
Ahora veamos la matris de rotacion:
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.08902432
Veamos como quedaria el calculo de las dimensiones del vector x
dim(pr.out$x)
[1] 50 4
Veamos como queda al graficarlo:
La idea de colocar el parametro scale con valor cero nos asegura que las flechas rojas se encuentren en la misma escala. Ahora arreglemos un poco esta grafica
pr.out$rotation <- -pr.out$rotation
pr.out$x = -pr.out$x
biplot(pr.out, scale=0)
Vemos que hemos rotado los ejes de las flechas para darle un poco mas de sentido a la grafica. La funcion prcomp tambien nos provee las desviacioens estandar, veamos:
pr.out$sdev
[1] 1.5748783 0.9948694 0.5971291 0.4164494
Ahora calculemos la varianza:
pr.var <- pr.out$sdev^2
pr.var
[1] 2.4802416 0.9897652 0.3565632 0.1734301
Ahora calculemos el PVE:
pve <- pr.var/sum(pr.var)
pve
[1] 0.62006039 0.24744129 0.08914080 0.04335752
Como vemos la primera PVE explica el model en un 62%, la segunda en un 24% y las ultimas dos en menos del 10%. Ahora veamos una grafica del PVE: