Lab 10.4: Principal components Analysis

Fiabilidad

Se usará PCA en el dataset de USArrests, el cual contiene \(50\) estados ordenados alfabeticamente.

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"       

Las variables del dataset son las siguientes:

names(USArrests)

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

Se usará la función apply() que permite aplicar una función a cada columna o fila del dataset. En este caso usaremos la función mean(). El segundo parámetro de la función apply() denota si queremos obtener la media de las filas o de las columnas, para las filas se utiliza un 1 y para las columnas un 2.

apply(USArrests, 2, var)

    Murder    Assault   UrbanPop       Rape 
  18.97047 6945.16571  209.51878   87.72916 

Las variables también tienen diferente varianza. La variable UrbanPop mide el porcentaje de población de cada estado que vive en un área urbana, que no se puede comparar con el número de violaciones en cada estado por \(100,000\) individuos. Es importante escalar las variables antes de aplicar el método de PCA, ya que se puede ver que la mayoría de los componentes principales estarían impulsadas por la variable Assault ya que tiene una media muy alta y una varianza muy grande. Por lo que primero se estandarizarán las variables para que tengan media igual a \(0\) y desviación estandar igual a \(1\) antes de aplicar PCA.

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

Por defecto, la función prcomp() centra las variables para que tengan una media igual a \(0\). Usando el parámetro scale = TRUE se escalan las variables para que tengan desviación estandar igual a \(1\). El resultado de la función prcomp() contiene unas cantidades utiles.

names(pr.out)

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

Los componentes center y scale corresponden a las medias y a las desviaciones estandar de las variables que fueron usadas para hacer la escala para poder implementar PCA.

pr.out$center

  Murder  Assault UrbanPop     Rape 
   7.788  170.760   65.540   21.232 

pr.out$center

  Murder  Assault UrbanPop     Rape 
   7.788  170.760   65.540   21.232 

rotation es una matriz que provee los componentes principales. Cada columna de pr.out$rotation contiene el componente principal del vector.

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

Se puede ver que hay cuatro componentes principales. Usando la función prcomp() no se necesita multiplicar los datos por los vectores de los componentes principales para obtener los puntajes de los componentes.

dim(pr.out$x)

[1] 50  4

Se pueden graficar los dos primeros componentes principales:

biplot(pr.out, scale = 0)

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

La función prcomp() también devuelve la desviación estandar de cada componente principal. Por lo que podemos ver las desviaciones estadar:

pr.out$sdev

[1] 1.5748783 0.9948694 0.5971291 0.4164494

La varianza por cada componente principal se obtiene elevando al cuadrado esas desviaciones estandar.

pr.var <- pr.out$sdev^2
pr.var

[1] 2.4802416 0.9897652 0.3565632 0.1734301

Ahora se dividirá la varianza entre cada componente principal por el total de varianza para los cuatro componentes principales.

pve <- pr.var/sum(pr.var)
pve

[1] 0.62006039 0.24744129 0.08914080 0.04335752

Aqui se puede ver que el primer componente principal explica \(62\%\) de la varianza de la data. El siguiente componente principal explica el \(24.7\%\) de la varianza. Se puede graficar el PVE por cada componente.

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 computa la suma acumulativa de los elementos de un vector numérico.

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

[1]  1  3 11  8