Recuerda que al refererinos a Aprendizaje No Supervisado nos interesa entender la estructura y las relaciones subyacentes en un conjunto de datos, las predicciones no forman parte de este tipo de análisis.

1. Análisis de Componentes Principales (PCA)

PCA es una técnica utilizada para describir un conjunto de datos en términos de nuevas componentes no correlacionadas. Los componentes se ordenan por la cantidad de varianza original que describen.

En este laboratorio, realizamos PCA en el conjunto de datos de USArrests, que es parte de R. Las filas del conjunto de datos contienen los 50 estados, en orden alfabético.

states <- row.names(USArrests )
states
##  [1] "Alabama"        "Alaska"         "Arizona"        "Arkansas"      
##  [5] "California"     "Colorado"       "Connecticut"    "Delaware"      
##  [9] "Florida"        "Georgia"        "Hawaii"         "Idaho"         
## [13] "Illinois"       "Indiana"        "Iowa"           "Kansas"        
## [17] "Kentucky"       "Louisiana"      "Maine"          "Maryland"      
## [21] "Massachusetts"  "Michigan"       "Minnesota"      "Mississippi"   
## [25] "Missouri"       "Montana"        "Nebraska"       "Nevada"        
## [29] "New Hampshire"  "New Jersey"     "New Mexico"     "New York"      
## [33] "North Carolina" "North Dakota"   "Ohio"           "Oklahoma"      
## [37] "Oregon"         "Pennsylvania"   "Rhode Island"   "South Carolina"
## [41] "South Dakota"   "Tennessee"      "Texas"          "Utah"          
## [45] "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. Para esto utilizaremos la función apply() que nos permite aplicar una función -en este caso, la función mean()- a cada fila o columna del conjunto de datos. El segundo argumento denota si queremos calcular la media de las filas, 1, o de las columnas, 2.

apply(USArrests, 2, mean)
##   Murder  Assault UrbanPop     Rape 
##    7.788  170.760   65.540   21.232

Vemos que hay en promedio tres veces más violaciones que asesinatos, y más de ocho veces más asaltos que violaciones. Notamos que las variables tienen medias muy diferentes. También podemos examinar las varianzas de las cuatro variables utilizando 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 también tengan varianzas muy diferentes: la variable UrbanPop mide el porcentaje de la población de cada estado que vive en una zona urbana, que no es un número comparable al número de violaciones en cada estado por cada 100,000 individuos. Si no pudiéramos escalar las variables antes de realizar el PCA, entonces la mayoría de los principales componentes que observamos estarían impulsados por la variable Assault, ya que tiene, por mucho, la mayor media y varianza. Por lo tanto, es importante estandarizar las variables para tener una media cero y una desviación estándar de uno antes de realizar el PCA.

Ahora realizamos el análisis de los componentes principales usando la función prcomp(), que es una de las varias funciones en R que realizan PCA.

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

Por defecto, la función prcomp() centra las variables para que tengan un promedio de cero. Usando la opción scale=TRUE, escalamos las variables para tener una desviación estándar de uno. La salida de prcomp() contiene una serie de cantidades útiles que vemos utilizando la función names() en el objeto.

names(pr.out)
## [1] "sdev"     "rotation" "center"   "scale"    "x"

Los componentes de center y scale corresponden a la media y a las desviaciones estándar de las variables que se utilizaron para la escala antes de aplicar la 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 rotation proporciona las cargas de los componentes principales; cada columna de pr.out$rotation contiene el correspondiente vector de carga de los componentes principales.

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

Vemos que hay cuatro componentes principales distintos.

Utilizando la función prcomp(), no es necesario multiplicar explícitamente los datos por los vectores de carga del componente principal para obtener los vectores de puntuación del componente principal. Más bien la matriz de 50 × 4 x tiene como columnas los vectores de puntuación del componente principal.

pr.out$x
##                        PC1         PC2         PC3          PC4
## Alabama        -0.97566045  1.12200121 -0.43980366  0.154696581
## Alaska         -1.93053788  1.06242692  2.01950027 -0.434175454
## Arizona        -1.74544285 -0.73845954  0.05423025 -0.826264240
## Arkansas        0.13999894  1.10854226  0.11342217 -0.180973554
## California     -2.49861285 -1.52742672  0.59254100 -0.338559240
## Colorado       -1.49934074 -0.97762966  1.08400162  0.001450164
## Connecticut     1.34499236 -1.07798362 -0.63679250 -0.117278736
## Delaware       -0.04722981 -0.32208890 -0.71141032 -0.873113315
## Florida        -2.98275967  0.03883425 -0.57103206 -0.095317042
## Georgia        -1.62280742  1.26608838 -0.33901818  1.065974459
## Hawaii          0.90348448 -1.55467609  0.05027151  0.893733198
## Idaho           1.62331903  0.20885253  0.25719021 -0.494087852
## Illinois       -1.36505197 -0.67498834 -0.67068647 -0.120794916
## Indiana         0.50038122 -0.15003926  0.22576277  0.420397595
## Iowa            2.23099579 -0.10300828  0.16291036  0.017379470
## Kansas          0.78887206 -0.26744941  0.02529648  0.204421034
## Kentucky        0.74331256  0.94880748 -0.02808429  0.663817237
## Louisiana      -1.54909076  0.86230011 -0.77560598  0.450157791
## Maine           2.37274014  0.37260865 -0.06502225 -0.327138529
## Maryland       -1.74564663  0.42335704 -0.15566968 -0.553450589
## Massachusetts   0.48128007 -1.45967706 -0.60337172 -0.177793902
## Michigan       -2.08725025 -0.15383500  0.38100046  0.101343128
## Minnesota       1.67566951 -0.62590670  0.15153200  0.066640316
## Mississippi    -0.98647919  2.36973712 -0.73336290  0.213342049
## Missouri       -0.68978426 -0.26070794  0.37365033  0.223554811
## Montana         1.17353751  0.53147851  0.24440796  0.122498555
## Nebraska        1.25291625 -0.19200440  0.17380930  0.015733156
## Nevada         -2.84550542 -0.76780502  1.15168793  0.311354436
## New Hampshire   2.35995585 -0.01790055  0.03648498 -0.032804291
## New Jersey     -0.17974128 -1.43493745 -0.75677041  0.240936580
## New Mexico     -1.96012351  0.14141308  0.18184598 -0.336121113
## New York       -1.66566662 -0.81491072 -0.63661186 -0.013348844
## North Carolina -1.11208808  2.20561081 -0.85489245 -0.944789648
## North Dakota    2.96215223  0.59309738  0.29824930 -0.251434626
## Ohio            0.22369436 -0.73477837 -0.03082616  0.469152817
## Oklahoma        0.30864928 -0.28496113 -0.01515592  0.010228476
## Oregon         -0.05852787 -0.53596999  0.93038718 -0.235390872
## Pennsylvania    0.87948680 -0.56536050 -0.39660218  0.355452378
## Rhode Island    0.85509072 -1.47698328 -1.35617705 -0.607402746
## South Carolina -1.30744986  1.91397297 -0.29751723 -0.130145378
## South Dakota    1.96779669  0.81506822  0.38538073 -0.108470512
## Tennessee      -0.98969377  0.85160534  0.18619262  0.646302674
## Texas          -1.34151838 -0.40833518 -0.48712332  0.636731051
## Utah            0.54503180 -1.45671524  0.29077592 -0.081486749
## Vermont         2.77325613  1.38819435  0.83280797 -0.143433697
## Virginia        0.09536670  0.19772785  0.01159482  0.209246429
## Washington      0.21472339 -0.96037394  0.61859067 -0.218628161
## West Virginia   2.08739306  1.41052627  0.10372163  0.130583080
## Wisconsin       2.05881199 -0.60512507 -0.13746933  0.182253407
## Wyoming         0.62310061  0.31778662 -0.23824049 -0.164976866
dim(pr.out$x)
## [1] 50  4

Podemos plotear los dos primeros componentes principales de la siguiente manera:

biplot(pr.out, scale=0)

El argumento scale=0 de biplot() asegura que las flechas se escalen para representar las cargas; otros valores de scale dan biplots ligeramente diferentes con interpretaciones diferentes.

En la figura que acabamos de elaborar tenemos la mayor parte de los componentes principales en sentido negativo, si cambiamos el signo de todos los valores podemos obtener un plot más interpretable con componentes principales mayormente positivos. Para esto, hacemos lo siguiente:

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

La función prcomp() también produce la desviación estándar de cada componente principal. Por ejemplo, en el conjunto de datos USArrests, podemos acceder a estas desviaciones estándar usando $sdev o simplemente utilizando la función print():

pr.out$sdev
## [1] 1.5748783 0.9948694 0.5971291 0.4164494
print(pr.out)
## Standard deviations (1, .., p=4):
## [1] 1.5748783 0.9948694 0.5971291 0.4164494
## 
## Rotation (n x k) = (4 x 4):
##                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

La varianza explicada por cada componente principal se obtiene mediante la cuadrática de éstos:

pr.var <- pr.out$sdev^2
pr.var
## [1] 2.4802416 0.9897652 0.3565632 0.1734301

Para calcular la proporción de la varianza explicada por cada componente principal, simplemente dividimos la varianza explicada por cada componente principal entre la varianza total explicada por los cuatro componentes principales o directamente utilizando la función summary() que además nos proporciona la varianza acumulada:

pve <- pr.var/sum(pr.var)
pve
## [1] 0.62006039 0.24744129 0.08914080 0.04335752
summary(pr.out)
## Importance of components:
##                           PC1    PC2     PC3     PC4
## Standard deviation     1.5749 0.9949 0.59713 0.41645
## Proportion of Variance 0.6201 0.2474 0.08914 0.04336
## Cumulative Proportion  0.6201 0.8675 0.95664 1.00000

Vemos 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 así sucesivamente. Podemos graficar la proporción de la varianza explicada (PVE), explicado por cada componente, así como el PVE acumulado, de la siguiente manera:

plot(pve, xlab="Componente Principal", ylab="Proporción de la Varianza Explicada", ylim=c(0,1), type='b')
text(pve,
     labels = round(pve,3),
     cex = 0.8, pos = 1, col = "black")

plot(cumsum(pve), xlab="Componente Principal", ylab="Proporción Acumulada de la Varianza Explicada", ylim=c(0,1), type='b')
text(cumsum(pve),
     labels = round(cumsum(pve),3),
     cex = 0.8, pos = 1, col = "black")

Obsérvese que la función cumsum() calcula la suma acumulativa de los elementos de un vector numérico. Por ejemplo:

a <- c(1,2,8,-3)
cumsum(a)
## [1]  1  3 11  8

2. Resuelve

Utilice la base de datos Boston de la libreria MASS. Cree un objeto pr.out.Boston que contenga el análisis de componentes principales de esta base de datos, no olvide incluir el argumento para escalar las variables.

library(MASS)

apply(Boston, 2, mean)
##         crim           zn        indus         chas          nox           rm 
##   3.61352356  11.36363636  11.13677866   0.06916996   0.55469506   6.28463439 
##          age          dis          rad          tax      ptratio        black 
##  68.57490119   3.79504269   9.54940711 408.23715415  18.45553360 356.67403162 
##        lstat         medv 
##  12.65306324  22.53280632
apply(Boston, 2, var)
##         crim           zn        indus         chas          nox           rm 
## 7.398658e+01 5.439368e+02 4.706444e+01 6.451297e-02 1.342764e-02 4.936709e-01 
##          age          dis          rad          tax      ptratio        black 
## 7.923584e+02 4.434015e+00 7.581637e+01 2.840476e+04 4.686989e+00 8.334752e+03 
##        lstat         medv 
## 5.099476e+01 8.458672e+01
pr.out.Boston <-  prcomp(Boston, scale=TRUE)
names(pr.out.Boston)
## [1] "sdev"     "rotation" "center"   "scale"    "x"
pr.out.Boston$center
##         crim           zn        indus         chas          nox           rm 
##   3.61352356  11.36363636  11.13677866   0.06916996   0.55469506   6.28463439 
##          age          dis          rad          tax      ptratio        black 
##  68.57490119   3.79504269   9.54940711 408.23715415  18.45553360 356.67403162 
##        lstat         medv 
##  12.65306324  22.53280632
pr.out.Boston$scale
##        crim          zn       indus        chas         nox          rm 
##   8.6015451  23.3224530   6.8603529   0.2539940   0.1158777   0.7026171 
##         age         dis         rad         tax     ptratio       black 
##  28.1488614   2.1057101   8.7072594 168.5371161   2.1649455  91.2948644 
##       lstat        medv 
##   7.1410615   9.1971041
pr.out.Boston$rotation
##                  PC1          PC2          PC3          PC4          PC5
## crim     0.242284451 -0.065873108  0.395077419 -0.100366211  0.004957659
## zn      -0.245435005 -0.148002653  0.394545713 -0.342958421  0.114495002
## indus    0.331859746  0.127075668 -0.066081913  0.009626936 -0.022583692
## chas    -0.005027133  0.410668763 -0.125305293 -0.700406497 -0.535197817
## nox      0.325193880  0.254276363 -0.046475549 -0.053707583  0.194605570
## rm      -0.202816554  0.434005810  0.353406095  0.293357309 -0.008320481
## age      0.296976574  0.260303205 -0.200823078  0.078426326  0.149750092
## dis     -0.298169809 -0.359149977  0.157068710 -0.184747787 -0.106219480
## rad      0.303412754  0.031149596  0.418510334  0.051374381 -0.230352185
## tax      0.324033052  0.008851406  0.343232194  0.026810695 -0.163425820
## ptratio  0.207679535 -0.314623061  0.000399092  0.342036328 -0.615707380
## black   -0.196638358  0.026481032 -0.361375914  0.201741185 -0.367460674
## lstat    0.311397955 -0.201245177 -0.161060336 -0.242621217  0.178358870
## medv    -0.266636396  0.444924411  0.163188735  0.180297553 -0.050659893
##                 PC6          PC7         PC8          PC9         PC10
## crim    -0.22462703  0.777083366 -0.15740140  0.254211798 -0.071384615
## zn      -0.33574694 -0.274178365  0.38031404  0.382899480  0.245579673
## indus   -0.08082495 -0.340273839 -0.17174578  0.627048264 -0.254827026
## chas     0.16264906  0.074075775  0.03292700 -0.018642967 -0.041706916
## nox     -0.14899191 -0.198092965 -0.04745838 -0.043024391 -0.211620349
## rm       0.13108056  0.074084938  0.43761566 -0.003666947 -0.526133916
## age     -0.06086960  0.118580363  0.58810569 -0.043265822  0.245647942
## dis      0.01162335 -0.104397844  0.12823060 -0.175802196 -0.299412026
## rad     -0.13493732 -0.137080107 -0.07464872 -0.463439397  0.115793486
## tax     -0.18847146 -0.313984433 -0.07099212 -0.179446555 -0.008366413
## ptratio  0.27901731  0.001485608  0.28346960  0.274525949  0.160474164
## black   -0.78590728  0.074842780  0.04444175 -0.060975651 -0.146292237
## lstat   -0.09197211  0.083213083  0.35748247 -0.171810921  0.066647267
## medv    -0.05402804 -0.009964973 -0.15230879  0.070751083  0.575547284
##                 PC11        PC12         PC13         PC14
## crim    -0.071068781  0.06327612  0.097032312  0.059114176
## zn      -0.127709065 -0.22112210 -0.132375830 -0.096296807
## indus    0.273797614  0.34840828  0.083716854 -0.235472877
## chas    -0.009968402 -0.01903975 -0.049917454  0.023488966
## nox     -0.437475550 -0.44909357  0.524974687  0.087649148
## rm       0.223951923 -0.12560554 -0.049893596  0.007190515
## age     -0.329630928  0.48633905 -0.051462562 -0.038227027
## dis     -0.114600078  0.49356822  0.552292172  0.047124029
## rad      0.042213348  0.01863641 -0.006278474 -0.634975332
## tax      0.042794054  0.17042179 -0.242987756  0.698822190
## ptratio -0.099991841 -0.23214842  0.188347079  0.055738160
## black    0.039194858 -0.04152885 -0.021078199 -0.016165280
## lstat    0.683032690 -0.18189209  0.249489863  0.083143795
## medv     0.242001064  0.09828580  0.469629324  0.134127182

Muestre la desviación estándar de cada componente principal.

pr.out.Boston$sdev
##  [1] 2.5585132 1.2843410 1.1614241 0.9415625 0.9224421 0.8124105 0.7317177
##  [8] 0.6348831 0.5265582 0.5022524 0.4612919 0.4277704 0.3660733 0.2456149
print(pr.out.Boston)
## Standard deviations (1, .., p=14):
##  [1] 2.5585132 1.2843410 1.1614241 0.9415625 0.9224421 0.8124105 0.7317177
##  [8] 0.6348831 0.5265582 0.5022524 0.4612919 0.4277704 0.3660733 0.2456149
## 
## Rotation (n x k) = (14 x 14):
##                  PC1          PC2          PC3          PC4          PC5
## crim     0.242284451 -0.065873108  0.395077419 -0.100366211  0.004957659
## zn      -0.245435005 -0.148002653  0.394545713 -0.342958421  0.114495002
## indus    0.331859746  0.127075668 -0.066081913  0.009626936 -0.022583692
## chas    -0.005027133  0.410668763 -0.125305293 -0.700406497 -0.535197817
## nox      0.325193880  0.254276363 -0.046475549 -0.053707583  0.194605570
## rm      -0.202816554  0.434005810  0.353406095  0.293357309 -0.008320481
## age      0.296976574  0.260303205 -0.200823078  0.078426326  0.149750092
## dis     -0.298169809 -0.359149977  0.157068710 -0.184747787 -0.106219480
## rad      0.303412754  0.031149596  0.418510334  0.051374381 -0.230352185
## tax      0.324033052  0.008851406  0.343232194  0.026810695 -0.163425820
## ptratio  0.207679535 -0.314623061  0.000399092  0.342036328 -0.615707380
## black   -0.196638358  0.026481032 -0.361375914  0.201741185 -0.367460674
## lstat    0.311397955 -0.201245177 -0.161060336 -0.242621217  0.178358870
## medv    -0.266636396  0.444924411  0.163188735  0.180297553 -0.050659893
##                 PC6          PC7         PC8          PC9         PC10
## crim    -0.22462703  0.777083366 -0.15740140  0.254211798 -0.071384615
## zn      -0.33574694 -0.274178365  0.38031404  0.382899480  0.245579673
## indus   -0.08082495 -0.340273839 -0.17174578  0.627048264 -0.254827026
## chas     0.16264906  0.074075775  0.03292700 -0.018642967 -0.041706916
## nox     -0.14899191 -0.198092965 -0.04745838 -0.043024391 -0.211620349
## rm       0.13108056  0.074084938  0.43761566 -0.003666947 -0.526133916
## age     -0.06086960  0.118580363  0.58810569 -0.043265822  0.245647942
## dis      0.01162335 -0.104397844  0.12823060 -0.175802196 -0.299412026
## rad     -0.13493732 -0.137080107 -0.07464872 -0.463439397  0.115793486
## tax     -0.18847146 -0.313984433 -0.07099212 -0.179446555 -0.008366413
## ptratio  0.27901731  0.001485608  0.28346960  0.274525949  0.160474164
## black   -0.78590728  0.074842780  0.04444175 -0.060975651 -0.146292237
## lstat   -0.09197211  0.083213083  0.35748247 -0.171810921  0.066647267
## medv    -0.05402804 -0.009964973 -0.15230879  0.070751083  0.575547284
##                 PC11        PC12         PC13         PC14
## crim    -0.071068781  0.06327612  0.097032312  0.059114176
## zn      -0.127709065 -0.22112210 -0.132375830 -0.096296807
## indus    0.273797614  0.34840828  0.083716854 -0.235472877
## chas    -0.009968402 -0.01903975 -0.049917454  0.023488966
## nox     -0.437475550 -0.44909357  0.524974687  0.087649148
## rm       0.223951923 -0.12560554 -0.049893596  0.007190515
## age     -0.329630928  0.48633905 -0.051462562 -0.038227027
## dis     -0.114600078  0.49356822  0.552292172  0.047124029
## rad      0.042213348  0.01863641 -0.006278474 -0.634975332
## tax      0.042794054  0.17042179 -0.242987756  0.698822190
## ptratio -0.099991841 -0.23214842  0.188347079  0.055738160
## black    0.039194858 -0.04152885 -0.021078199 -0.016165280
## lstat    0.683032690 -0.18189209  0.249489863  0.083143795
## medv     0.242001064  0.09828580  0.469629324  0.134127182

Obtenga la varianza explicada por cada componente principal, recuerde que se obtiene mediante la cuadrática de la desviación estándar.

pr.var.Boston <- pr.out.Boston$sdev^2
pr.var.Boston
##  [1] 6.54598958 1.64953191 1.34890592 0.88653987 0.85089944 0.66001077
##  [7] 0.53541080 0.40307658 0.27726358 0.25225744 0.21279025 0.18298750
## [13] 0.13400970 0.06032666

Calcule la proporción de la varianza explicada por cada componente principal.

pve.Boston <- pr.var.Boston/sum(pr.var.Boston)
pve.Boston
##  [1] 0.467570684 0.117823708 0.096350423 0.063324276 0.060778531 0.047143627
##  [7] 0.038243629 0.028791184 0.019804542 0.018018389 0.015199304 0.013070536
## [13] 0.009572121 0.004309047
summary(pr.out.Boston)
## Importance of components:
##                           PC1    PC2     PC3     PC4     PC5     PC6     PC7
## Standard deviation     2.5585 1.2843 1.16142 0.94156 0.92244 0.81241 0.73172
## Proportion of Variance 0.4676 0.1178 0.09635 0.06332 0.06078 0.04714 0.03824
## Cumulative Proportion  0.4676 0.5854 0.68174 0.74507 0.80585 0.85299 0.89123
##                            PC8    PC9    PC10   PC11    PC12    PC13    PC14
## Standard deviation     0.63488 0.5266 0.50225 0.4613 0.42777 0.36607 0.24561
## Proportion of Variance 0.02879 0.0198 0.01802 0.0152 0.01307 0.00957 0.00431
## Cumulative Proportion  0.92003 0.9398 0.95785 0.9730 0.98612 0.99569 1.00000

Grafique la Proporción Acumulada de la Varianza Explicada.

plot(pve.Boston, xlab="Componente Principal", ylab="Proporción de la Varianza Explicada", ylim=c(0,1), type='b')
text(pve.Boston,
     labels = round(pve,3),
     cex = 0.8, pos = 1, col = "black")

plot(cumsum(pve.Boston), xlab="Componente Principal", ylab="Proporción Acumulada de la Varianza Explicada", ylim=c(0,1), type='b')
text(cumsum(pve.Boston),
     labels = round(cumsum(pve.Boston),3),
     cex = 0.8, pos = 1, col = "black")

¿Cuántos componentes principales explican por lo menos el 80% de la varianza?

5