ISLR - Lab 10.4: Principal Components Analysis

En este laboratorio usaremos el dataset USArrests que encontramos por defecto en R. En cada fila encontramos uno de los 50 estados ordenados alfabeticamente.

states =row.names(USArrests)
states

 [1] "Alabama"        "Alaska"         "Arizona"        "Arkansas"       "California"    
 [6] "Colorado"       "Connecticut"    "Delaware"       "Florida"        "Georgia"       
[11] "Hawaii"         "Idaho"          "Illinois"       "Indiana"        "Iowa"          
[16] "Kansas"         "Kentucky"       "Louisiana"      "Maine"          "Maryland"      
[21] "Massachusetts"  "Michigan"       "Minnesota"      "Mississippi"    "Missouri"      
[26] "Montana"        "Nebraska"       "Nevada"         "New Hampshire"  "New Jersey"    
[31] "New Mexico"     "New York"       "North Carolina" "North Dakota"   "Ohio"          
[36] "Oklahoma"       "Oregon"         "Pennsylvania"   "Rhode Island"   "South Carolina"
[41] "South Dakota"   "Tennessee"      "Texas"          "Utah"           "Vermont"       
[46] "Virginia"       "Washington"     "West Virginia"  "Wisconsin"      "Wyoming"       

names(USArrests)

[1] "Murder"   "Assault"  "UrbanPop" "Rape"    

cada una de las variables tiene una media muy diferente

apply(USArrests , 2, mean)

  Murder  Assault UrbanPop     Rape 
   7.788  170.760   65.540   21.232 

las violaciones son poco mas del doble que los asesinatos y asaltos es lo mas común, luego vemos las varianzas de cada variable

apply(USArrests , 2, var)

    Murder    Assault   UrbanPop       Rape 
  18.97047 6945.16571  209.51878   87.72916 

UrbanPop nos dice el porcentaje de las personas en cada estado que viven en areas urbanas, que no es comparable con el numero de violaciones por cada 100,000 habitantes. Por lo tanto debemos estandarizar nuestras variables para poder realizar el análisis PCA sino Assaults será la predominante porque tener la mayor media y varianza. Por lo tanto es importante hacer que las variables tenga media cero y desviacion estandar cero antes del análisis. El análisis lo haremos con la función prcomp() que es una de las muchas opciones disponibles en el mercado.

pr.out =prcomp (USArrests , scale =TRUE)

por defeto la función escala las variables a media cero y el parámetro scale=TRUE hace que la desviación estándar sea uno

names(pr.out )

[1] "sdev"     "rotation" "center"   "scale"    "x"       

el componente center son las medias y scale son las desviaciones estandar de las variables que se utilizaron antes del 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 rotation contiene en cada columna el vector PCA corresponiente

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

no necestimas multiplicar los vectores por nuestras variables originales para obtener los vectores de score pca, sino que ya se encuentran en la matriz x donde cada columna es un vector que se encuentran en el mismo orden del dataset que estamos utilizando

dim(pr.out$x )

[1] 50  4

podemos graficar los primeros dos componentes

biplot (pr.out , scale =0)

el parámetro scale=0 nos asegura que las flechas esten en la escala para representar correctamente las cargas

pr.out$rotation=-pr.out$rotation
pr.out$x=-pr.out$x
biplot (pr.out , scale =0)

también podemos extraer las desviaciones estandar de cada componente

pr.out$sdev

[1] 1.5748783 0.9948694 0.5971291 0.4164494

la varianza de cada componente

pr.var =pr.out$sdev ^2
pr.var

[1] 2.4802416 0.9897652 0.3565632 0.1734301

la proporción de la varianza que explicada cada uno de los vectores principales

pve=pr.var/sum(pr.var )
pve

[1] 0.62006039 0.24744129 0.08914080 0.04335752

el primer vector explica el 62% de la varianza, el segudno explica el 24.7% y así sucesivamente. Estos datos tambien los podemos examinar de forma gráfica.

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'
)

la función cumsum() da la suma acumulada de un vector numerico

a=c(1,2,8,-3)
cumsum (a)

[1]  1  3 11  8