Analisis Principal de Componentes

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:

LS0tCnRpdGxlOiAiRmlhYmlsaWRhZCAtIExhYm9yYXRvcmlvIDMiCm91dHB1dDogaHRtbF9ub3RlYm9vawotLS0KCiMjIEFuYWxpc2lzIFByaW5jaXBhbCBkZSBDb21wb25lbnRlcwoKRW4gZXN0ZSBsYWJvcmF0b3JpbywgaGFyZW1vcyBQQ0EgKEFuYWxpc2lzIFByaW5jaXBhbCBkZSBDb21wb25lbnRlcykgc29icmUgZWwgZGF0YXNldCBgVVNBcnJlc3RzYCwgZWwgY3VhbCB2aWVuZSBpbnRlZ3JhZG8gZW4gUiwgZXN0YXMgY29udGllbmVuIGxvcyA1MCBlc3RhZG9zIG9yZGVuYWRvcyBhbGZhYmV0aWNhbWVudGUgY29tbyB2ZXJlbW9zIGEgY29udGludWFjaW9uOgoKYGBge3J9CnN0YXRlcyA9IHJvdy5uYW1lcyhVU0FycmVzdHMpCnN0YXRlcwpgYGAKCkVsIGRhdGFzZXQgY29udGllbmUgNCBjb2x1bW5hcyBxdWUgc29uOgoKYGBge3J9Cm5hbWVzKFVTQXJyZXN0cykKYGBgCgpFeGFtaW5lbW9zIHVuIHBvY28gbGEgZGF0YToKCmBgYHtyfQphcHBseShVU0FycmVzdHMsIDIsIG1lYW4pCmBgYAoKbGEgZnVuY2lvbiBgYXBwbHlgIG5vcyBwZXJtaXRlIGl0ZXJhciBzb2JyZSBjYWRhIGZpbGEgbyBjb2x1bW5hIGRlbCBkYXRhc2V0LCBlbCBwYXJhbWV0cm8gMiBpbmRpY2EgZXhwbGljaXRhbWVudGUgcXVlIGl0ZXJhcmVtb3MgcG9yIGNvbHVtbmEsIHNpIHNlIGxlIGVudmlhIDEgZW50b25jZXMgc2UgaXRlcmEgcG9yIGZpbGEsIHkgZmluYWxtZW50ZSBlbnZpYW1vcyBsYSBmdW5jaW9uIGBtZWFuYCBwYXJhIGdlbmVyYXIgZWwgcHJvbWVkaW8gZGUgY2FkYSBjb2x1bW5hLiBBbCBhbmFsaXphciBsYSBkYXRhIHZlbW9zIHF1ZSBoYXkgdW5hIHZhcmlhbnphIHJlbGF0aXZhbWVudGUgZ3JhbmRlLCBhaG9yYSB2ZWFtb3MgbG8gcXVlIHN1Y2VkZSBjdWFuZG8gaGFjZW1vcyB1biBhbmFsaXNpcyBkZSBjb21wb25lbnRlczoKCmBgYHtyfQpwci5vdXQgPC0gcHJjb21wKFVTQXJyZXN0cywgc2NhbGU9VFJVRSkKYGBgCgpQb3IgZGVmYXVsdCBsYSBmdW5jaW9uIGBwcmNvbXBgIGNlbnRyYSBsYXMgdmFyaWFibGVzIHBhcmEgcXVlIGVzdGFuIHRlbmdhbiB1bmEgbWVkaWEgZGUgY2VybywgYWRlbWFzIGFsIHV0aWxpemFyIGVsIHBhcmFtZXRybyBgc2NhbGVgIHBhcmEgdGVuZXIgdW5hIGRlc3ZpYWNpb24gZXN0YW5kYXIgZGUgdW5vLiAKCmBgYHtyfQpuYW1lcyhwci5vdXQpCmBgYAoKVmVtb3MgcXVlIGVuIGVsIHNldCBkZSBkYXRvcyBsYXMgdmFyaWFibGVzIGBjZW50ZXJgIHkgYHNjYWxlYCBjb250aWVuZW4gbGEgbWVkaWEgeSBsYSBkZXN2aWFjaW9uIGVzdGFuZGFyIAoKYGBge3J9CnByLm91dCRjZW50ZXIKYGBgCmBgYHtyfQpwci5vdXQkc2NhbGUKYGBgCgpBaG9yYSB2ZWFtb3MgbGEgbWF0cmlzIGRlIHJvdGFjaW9uOgoKYGBge3J9CnByLm91dCRyb3RhdGlvbgpgYGAKClZlYW1vcyBjb21vIHF1ZWRhcmlhIGVsIGNhbGN1bG8gZGUgbGFzIGRpbWVuc2lvbmVzIGRlbCB2ZWN0b3IgYHhgCgpgYGB7cn0KZGltKHByLm91dCR4KQpgYGAKClZlYW1vcyBjb21vIHF1ZWRhIGFsIGdyYWZpY2FybG86CgpgYGB7cn0KYmlwbG90KHByLm91dCwgc2NhbGU9MCkKYGBgCgpMYSBpZGVhIGRlIGNvbG9jYXIgZWwgcGFyYW1ldHJvIGBzY2FsZWAgY29uIHZhbG9yIGNlcm8gbm9zIGFzZWd1cmEgcXVlIGxhcyBmbGVjaGFzIHJvamFzIHNlIGVuY3VlbnRyZW4gZW4gbGEgbWlzbWEgZXNjYWxhLiBBaG9yYSBhcnJlZ2xlbW9zIHVuIHBvY28gZXN0YSBncmFmaWNhCgpgYGB7cn0KcHIub3V0JHJvdGF0aW9uIDwtIC1wci5vdXQkcm90YXRpb24KcHIub3V0JHggPSAtcHIub3V0JHgKYmlwbG90KHByLm91dCwgc2NhbGU9MCkKYGBgCgpWZW1vcyBxdWUgaGVtb3Mgcm90YWRvIGxvcyBlamVzIGRlIGxhcyBmbGVjaGFzIHBhcmEgZGFybGUgdW4gcG9jbyBtYXMgZGUgc2VudGlkbyBhIGxhIGdyYWZpY2EuIExhIGZ1bmNpb24gYHByY29tcGAgdGFtYmllbiBub3MgcHJvdmVlIGxhcyBkZXN2aWFjaW9lbnMgZXN0YW5kYXIsIHZlYW1vczoKCmBgYHtyfQpwci5vdXQkc2RldgpgYGAKCkFob3JhIGNhbGN1bGVtb3MgbGEgdmFyaWFuemE6CgpgYGB7cn0KcHIudmFyIDwtIHByLm91dCRzZGV2XjIKcHIudmFyCmBgYAoKQWhvcmEgY2FsY3VsZW1vcyBlbCBQVkU6CgpgYGB7cn0KcHZlIDwtIHByLnZhci9zdW0ocHIudmFyKQpwdmUKYGBgCkNvbW8gdmVtb3MgbGEgcHJpbWVyYSBQVkUgZXhwbGljYSBlbCBtb2RlbCBlbiB1biA2MiUsIGxhIHNlZ3VuZGEgZW4gdW4gMjQlIHkgbGFzIHVsdGltYXMgZG9zIGVuIG1lbm9zIGRlbCAxMCUuIEFob3JhIHZlYW1vcyB1bmEgZ3JhZmljYSBkZWwgUFZFOgoKYGBge3J9CnBsb3QocHZlLCB4bGF2PSJQcmluY2lwYWwgQ29tcG9uZW50IiwgeWxhYj0iUHJvcG9ydGlvbiBvZiBWYXJpYW5jZSBFeHBsYWluZWQiLCB5bGltPWMoMCwxKSwgdHlwZT0iYiIpCmBgYAoKYGBge3J9CnBsb3QoY3Vtc3VtKHB2ZSksIHhsYXY9IlByaW5jaXBhbCBDb21wb25lbnQiLCB5bGFiPSJDdW11bGF0aXZlIFByb3BvcnRpb24gb2YgVmFyaWFuY2UgRXhwbGFpbmVkIiwgeWxpbT1jKDAsMSksIHR5cGU9ImIiKQpgYGAKCg==