Análisis de Componentes Principales
2/8/2022
Ejemplo cálculo directo de PCA con R
El set de datos USArrests del paquete básico de
R contiene el porcentaje de asaltos (\(Assault\)), asesinatos (\(Murder\)) y secuestros (\(Rape\)) por cada 100,000 habitantes para
cada uno de los 50 estados de USA (1973). Además, también incluye el
porcentaje de la población de cada estado que vive en zonas rurales
(\(UrbanPoP\)).
data("USArrests")
head(USArrests)## Murder Assault UrbanPop Rape
## Alabama 13.2 236 58 21.2
## Alaska 10.0 263 48 44.5
## Arizona 8.1 294 80 31.0
## Arkansas 8.8 190 50 19.5
## California 9.0 276 91 40.6
## Colorado 7.9 204 78 38.7El promedio de los datos muestra que hay tres veces más secuestros que asesinatos y 8 veces más asaltos que secuestros.
apply(X = USArrests, MARGIN = 2, FUN = mean)## Murder Assault UrbanPop Rape
## 7.788 170.760 65.540 21.232La varianza es muy distinta entre las variables, en el caso de Assault, la varianza es varios órdenes de magnitud superior al resto.
apply(X = USArrests, MARGIN = 2, FUN = var)## Murder Assault UrbanPop Rape
## 18.97047 6945.16571 209.51878 87.72916Si no se estandarizan las variables para que tengan media cero y desviación estándar 1 antes de realizar el estudio \(\textbf{PCA}\), la variable Assault dominará la mayoría de las componentes principales.
La función prcomp() es una de las múltiples funciones en
R que realizan \(\textbf{PCA}\). Por defecto,
prcomp() centra las variables para que tengan media cero,
pero si se quiere además que su desviación estándar sea de uno, hay que
indicar scale = TRUE.
pca <- prcomp(USArrests, scale = TRUE)
names(pca)## [1] "sdev" "rotation" "center" "scale" "x"Los elementos center y scale almacenados en
el objeto pca contienen la media y desviación típica de las
variables previa estandarización (en la escala original).
pca$center## Murder Assault UrbanPop Rape
## 7.788 170.760 65.540 21.232pca$scale## Murder Assault UrbanPop Rape
## 4.355510 83.337661 14.474763 9.366385rotation contiene el valor de los loadings \(\phi\) para cada componente (\(eigenvector\)). El número máximo de
componentes principales se corresponde con el \(mínimo(n-1,p)\), que en este caso es \(min(49,4)=4\).
pca$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.08902432Analizar con detalle el vector de \(loadings\) que forma cada componente puede ayudar a interpretar que tipo de información recoge cada una de ellas. Por ejemplo, la primera componente es el resultado de la siguiente combinación lineal de las variables originales:
\[PC1=−0.5358995 Murder−0.5831836 Assault−0.2781909 UrbanPop−0.5434321 Rape\]
Los pesos asignados en la primera componente a las variables \(Assault\), \(Murder\) y \(Rape\) son aproximadamente iguales entre ellos y bastante superiores al asignado a \(UrbanPoP\), esto significa que la primera componente recoge mayoritariamente la información correspondiente a los delitos. En la segunda componente, es la variable \(UrbanPoP\) la que tiene con diferencia mayor peso, por lo que se corresponde principalmente con el nivel de urbanización del estado. Si bien en este ejemplo la interpretación de las componentes es bastante clara, no en todos los casos ocurre lo mismo.
La función prcomp() calcula automáticamente el valor de
las componentes principales para cada observación (\(principal\,\, component\,\, scores\))
multiplicando los datos por los vectores de \(loadings\). El resultado se almacena en la
matriz \(x\).
head(pca$x)## PC1 PC2 PC3 PC4
## Alabama -0.9756604 1.1220012 -0.43980366 0.154696581
## Alaska -1.9305379 1.0624269 2.01950027 -0.434175454
## Arizona -1.7454429 -0.7384595 0.05423025 -0.826264240
## Arkansas 0.1399989 1.1085423 0.11342217 -0.180973554
## California -2.4986128 -1.5274267 0.59254100 -0.338559240
## Colorado -1.4993407 -0.9776297 1.08400162 0.001450164dim(pca$x)## [1] 50 4Mediante la función biplot() se puede obtener una
representación bidimensional de las dos primeras componentes. Es
recomendable indicar el argumento scale = 0 para que las
flechas estén en la misma escala que las componentes.
biplot(x = pca, scale = 0, cex = 0.6, col = c("blue4", "brown3"))
Una vez calculadas las componentes principales, se puede conocer la varianza explicada por cada una de ellas, la proporción respecto al total y la proporción de varianza acumulada.
library(ggplot2)
pca$sdev^2## [1] 2.4802416 0.9897652 0.3565632 0.1734301prop_varianza <- pca$sdev^2 / sum(pca$sdev^2)
prop_varianza## [1] 0.62006039 0.24744129 0.08914080 0.04335752ggplot(data = data.frame(prop_varianza, pc = 1:4),
aes(x = pc, y = prop_varianza)) +
geom_col(width = 0.3) +
scale_y_continuous(limits = c(0,1)) +
theme_bw() +
labs(x = "Componente principal",
y = "Prop. de varianza explicada")
prop_varianza_acum <- cumsum(prop_varianza)
prop_varianza_acum## [1] 0.6200604 0.8675017 0.9566425 1.0000000ggplot(data = data.frame(prop_varianza_acum, pc = 1:4),
aes(x = pc, y = prop_varianza_acum, group = 1)) +
geom_point() +
geom_line() +
theme_bw() +
labs(x = "Componente principal",
y = "Prop. varianza explicada acumulada")
En este caso, la primera componente explica el 62\(\%\) de la varianza observada en los datos y la segunda el 24.7\(\%\). Las dos últimas componentes no superan por separado el 1\(\%\) de varianza explicada. Si se empleasen únicamente las dos primeras componentes se conseguiría explicar el 86.75\(\%\) de la varianza observada.
PCR: PCA aplicado a regresión lineal
La cuantificación del contenido en grasa de la carne pude hacerse
mediante técnicas de analítica química, sin embargo, este proceso es
costoso en tiempo y recursos. Una posible alternativa para reducir
costes y optimizar tiempo es emplear un espectrofotómetro (instrumento
capaz de detectar la absorbancia que tiene un material a diferentes
tipos de luz en función de sus características). Para comprobar su
efectividad se mide el espectro de absorbancia de 100 longitudes de onda
en 215 muestras de carne, cuyo contenido en grasa se obtiene también por
análisis químico para poder comparar los resultados. El set de datos
meatspec del paquete faraway contiene toda la
información.
library(faraway)
data(meatspec)
head(meatspec)## V1 V2 V3 V4 V5 V6 V7 V8 V9
## 1 2.61776 2.61814 2.61859 2.61912 2.61981 2.62071 2.62186 2.62334 2.62511
## 2 2.83454 2.83871 2.84283 2.84705 2.85138 2.85587 2.86060 2.86566 2.87093
## 3 2.58284 2.58458 2.58629 2.58808 2.58996 2.59192 2.59401 2.59627 2.59873
## 4 2.82286 2.82460 2.82630 2.82814 2.83001 2.83192 2.83392 2.83606 2.83842
## 5 2.78813 2.78989 2.79167 2.79350 2.79538 2.79746 2.79984 2.80254 2.80553
## 6 3.00993 3.01540 3.02086 3.02634 3.03190 3.03756 3.04341 3.04955 3.05599
## V10 V11 V12 V13 V14 V15 V16 V17 V18
## 1 2.62722 2.62964 2.63245 2.63565 2.63933 2.64353 2.64825 2.65350 2.65937
## 2 2.87661 2.88264 2.88898 2.89577 2.90308 2.91097 2.91953 2.92873 2.93863
## 3 2.60131 2.60414 2.60714 2.61029 2.61361 2.61714 2.62089 2.62486 2.62909
## 4 2.84097 2.84374 2.84664 2.84975 2.85307 2.85661 2.86038 2.86437 2.86860
## 5 2.80890 2.81272 2.81704 2.82184 2.82710 2.83294 2.83945 2.84664 2.85458
## 6 3.06274 3.06982 3.07724 3.08511 3.09343 3.10231 3.11185 3.12205 3.13294
## V19 V20 V21 V22 V23 V24 V25 V26 V27
## 1 2.66585 2.67281 2.68008 2.68733 2.69427 2.70073 2.70684 2.71281 2.71914
## 2 2.94929 2.96072 2.97272 2.98493 2.99690 3.00833 3.01920 3.02990 3.04101
## 3 2.63361 2.63835 2.64330 2.64838 2.65354 2.65870 2.66375 2.66880 2.67383
## 4 2.87308 2.87789 2.88301 2.88832 2.89374 2.89917 2.90457 2.90991 2.91521
## 5 2.86331 2.87280 2.88291 2.89335 2.90374 2.91371 2.92305 2.93187 2.94060
## 6 3.14457 3.15703 3.17038 3.18429 3.19840 3.21225 3.22552 3.23827 3.25084
## V28 V29 V30 V31 V32 V33 V34 V35 V36
## 1 2.72628 2.73462 2.74416 2.75466 2.76568 2.77679 2.78790 2.79949 2.81225
## 2 3.05345 3.06777 3.08416 3.10221 3.12106 3.13983 3.15810 3.17623 3.19519
## 3 2.67892 2.68411 2.68937 2.69470 2.70012 2.70563 2.71141 2.71775 2.72490
## 4 2.92043 2.92565 2.93082 2.93604 2.94128 2.94658 2.95202 2.95777 2.96419
## 5 2.94986 2.96035 2.97241 2.98606 3.00097 3.01652 3.03220 3.04793 3.06413
## 6 3.26393 3.27851 3.29514 3.31401 3.33458 3.35591 3.37709 3.39772 3.41828
## V37 V38 V39 V40 V41 V42 V43 V44 V45
## 1 2.82706 2.84356 2.86106 2.87857 2.89497 2.90924 2.92085 2.93015 2.93846
## 2 3.21584 3.23747 3.25889 3.27835 3.29384 3.30362 3.30681 3.30393 3.29700
## 3 2.73344 2.74327 2.75433 2.76642 2.77931 2.79272 2.80649 2.82064 2.83541
## 4 2.97159 2.98045 2.99090 3.00284 3.01611 3.03048 3.04579 3.06194 3.07889
## 5 3.08153 3.10078 3.12185 3.14371 3.16510 3.18470 3.20140 3.21477 3.22544
## 6 3.43974 3.46266 3.48663 3.51002 3.53087 3.54711 3.55699 3.55986 3.55656
## V46 V47 V48 V49 V50 V51 V52 V53 V54
## 1 2.94771 2.96019 2.97831 3.00306 3.03506 3.07428 3.11963 3.16868 3.21771
## 2 3.28925 3.28409 3.28505 3.29326 3.30923 3.33267 3.36251 3.39661 3.43188
## 3 2.85121 2.86872 2.88905 2.91289 2.94088 2.97325 3.00946 3.04780 3.08554
## 4 3.09686 3.11629 3.13775 3.16217 3.19068 3.22376 3.26172 3.30379 3.34793
## 5 3.23505 3.24586 3.26027 3.28063 3.30889 3.34543 3.39019 3.44198 3.49800
## 6 3.54937 3.54169 3.53692 3.53823 3.54760 3.56512 3.59043 3.62229 3.65830
## V55 V56 V57 V58 V59 V60 V61 V62 V63
## 1 3.26254 3.29988 3.32847 3.34899 3.36342 3.37379 3.38152 3.38741 3.39164
## 2 3.46492 3.49295 3.51458 3.53004 3.54067 3.54797 3.55306 3.55675 3.55921
## 3 3.11947 3.14696 3.16677 3.17938 3.18631 3.18924 3.18950 3.18801 3.18498
## 4 3.39093 3.42920 3.45998 3.48227 3.49687 3.50558 3.51026 3.51221 3.51215
## 5 3.55407 3.60534 3.64789 3.68011 3.70272 3.71815 3.72863 3.73574 3.74059
## 6 3.69515 3.72932 3.75803 3.78003 3.79560 3.80614 3.81313 3.81774 3.82079
## V64 V65 V66 V67 V68 V69 V70 V71 V72
## 1 3.39418 3.39490 3.39366 3.39045 3.38541 3.37869 3.37041 3.36073 3.34979
## 2 3.56045 3.56034 3.55876 3.55571 3.55132 3.54585 3.53950 3.53235 3.52442
## 3 3.18039 3.17411 3.16611 3.15641 3.14512 3.13241 3.11843 3.10329 3.08714
## 4 3.51036 3.50682 3.50140 3.49398 3.48457 3.47333 3.46041 3.44595 3.43005
## 5 3.74357 3.74453 3.74336 3.73991 3.73418 3.72638 3.71676 3.70553 3.69289
## 6 3.82258 3.82301 3.82206 3.81959 3.81557 3.81021 3.80375 3.79642 3.78835
## V73 V74 V75 V76 V77 V78 V79 V80 V81
## 1 3.33769 3.32443 3.31013 3.29487 3.27891 3.26232 3.24542 3.22828 3.21080
## 2 3.51583 3.50668 3.49700 3.48683 3.47626 3.46552 3.45501 3.44481 3.43477
## 3 3.07014 3.05237 3.03393 3.01504 2.99569 2.97612 2.95642 2.93660 2.91667
## 4 3.41285 3.39450 3.37511 3.35482 3.33376 3.31204 3.28986 3.26730 3.24442
## 5 3.67900 3.66396 3.64785 3.63085 3.61305 3.59463 3.57582 3.55695 3.53796
## 6 3.77958 3.77024 3.76040 3.75005 3.73929 3.72831 3.71738 3.70681 3.69664
## V82 V83 V84 V85 V86 V87 V88 V89 V90
## 1 3.19287 3.17433 3.15503 3.13475 3.11339 3.09116 3.06850 3.04596 3.02393
## 2 3.42465 3.41419 3.40303 3.39082 3.37731 3.36265 3.34745 3.33245 3.31818
## 3 2.89655 2.87622 2.85563 2.83474 2.81361 2.79235 2.77113 2.75015 2.72956
## 4 3.22117 3.19757 3.17357 3.14915 3.12429 3.09908 3.07366 3.04825 3.02308
## 5 3.51880 3.49936 3.47938 3.45869 3.43711 3.41458 3.39129 3.36772 3.34450
## 6 3.68659 3.67649 3.66611 3.65503 3.64283 3.62938 3.61483 3.59990 3.58535
## V91 V92 V93 V94 V95 V96 V97 V98 V99
## 1 3.00247 2.98145 2.96072 2.94013 2.91978 2.89966 2.87964 2.85960 2.83940
## 2 3.30473 3.29186 3.27921 3.26655 3.25369 3.24045 3.22659 3.21181 3.19600
## 3 2.70934 2.68951 2.67009 2.65112 2.63262 2.61461 2.59718 2.58034 2.56404
## 4 2.99820 2.97367 2.94951 2.92576 2.90251 2.87988 2.85794 2.83672 2.81617
## 5 3.32201 3.30025 3.27907 3.25831 3.23784 3.21765 3.19766 3.17770 3.15770
## 6 3.57163 3.55877 3.54651 3.53442 3.52221 3.50972 3.49682 3.48325 3.46870
## V100 fat
## 1 2.81920 22.5
## 2 3.17942 40.1
## 3 2.54816 8.4
## 4 2.79622 5.9
## 5 3.13753 25.5
## 6 3.45307 42.7dim(meatspec)## [1] 215 101El set de datos contiene 101 columnas. Las 100 primeras, nombradas como \(v_1, ..., v_{100}\) recogen el valor de absorbancia para cada una de las 100 longitudes de onda analizadas, y la columna fat el contenido en grasa medido por técnicas químicas.
Para poder evaluar la capacidad predictiva del modelo, se dividen las observaciones disponibles en dos grupos: uno de entrenamiento para ajustar el modelo (80\(\%\) de los datos) y uno de test (20\(\%\) de los datos).
training <- meatspec[1:172, ]
test <- meatspec[173:215, ]En primer lugar se ajusta un modelo incluyendo todas las longitudes de onda como predictores.
modelo <- lm(fat ~ ., data = training)
summary(modelo)##
## Call:
## lm(formula = fat ~ ., data = training)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.09837 -0.35779 0.04555 0.38080 2.33860
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.324 2.012 3.143 0.002439 **
## V1 12134.077 3659.798 3.316 0.001443 **
## V2 -12585.857 5971.891 -2.108 0.038605 *
## V3 -5107.556 9390.265 -0.544 0.588200
## V4 23880.493 17143.644 1.393 0.167977
## V5 -40509.555 22129.359 -1.831 0.071360 .
## V6 28469.416 19569.400 1.455 0.150134
## V7 -20901.082 12501.639 -1.672 0.098952 .
## V8 8369.465 7515.467 1.114 0.269193
## V9 -1539.328 5397.505 -0.285 0.776327
## V10 4706.267 7406.895 0.635 0.527217
## V11 7012.943 11720.620 0.598 0.551516
## V12 14891.444 20169.170 0.738 0.462749
## V13 -30963.902 26186.839 -1.182 0.240983
## V14 34338.612 22323.830 1.538 0.128444
## V15 -22235.237 13842.268 -1.606 0.112640
## V16 -7466.797 8558.172 -0.872 0.385890
## V17 6716.653 6561.805 1.024 0.309500
## V18 -2033.071 6741.330 -0.302 0.763851
## V19 8541.212 9419.998 0.907 0.367627
## V20 -1667.207 17300.433 -0.096 0.923500
## V21 -31972.494 24622.615 -1.299 0.198317
## V22 59526.389 27730.712 2.147 0.035244 *
## V23 -49241.388 23117.226 -2.130 0.036632 *
## V24 16184.597 16679.609 0.970 0.335180
## V25 12077.951 10751.912 1.123 0.265081
## V26 -12632.330 6774.573 -1.865 0.066361 .
## V27 -6298.837 7032.334 -0.896 0.373442
## V28 29625.988 9011.227 3.288 0.001573 **
## V29 -39374.835 13561.228 -2.903 0.004914 **
## V30 31251.427 18742.000 1.667 0.099829 .
## V31 -27238.189 21335.756 -1.277 0.205887
## V32 23009.543 19776.156 1.163 0.248522
## V33 -4584.373 14572.471 -0.315 0.753995
## V34 -5437.943 10344.728 -0.526 0.600754
## V35 -6128.931 8762.663 -0.699 0.486564
## V36 5599.605 6652.640 0.842 0.402776
## V37 -5569.160 6670.198 -0.835 0.406557
## V38 97.451 9291.480 0.010 0.991661
## V39 36021.407 12574.711 2.865 0.005488 **
## V40 -54273.400 17144.384 -3.166 0.002280 **
## V41 52084.876 21758.024 2.394 0.019318 *
## V42 -48458.089 23950.549 -2.023 0.046813 *
## V43 29334.488 20232.617 1.450 0.151500
## V44 -18282.834 13508.157 -1.353 0.180200
## V45 22110.934 9725.348 2.274 0.026020 *
## V46 -11735.692 6631.245 -1.770 0.081061 .
## V47 -514.521 3800.612 -0.135 0.892696
## V48 2551.480 6131.893 0.416 0.678592
## V49 3707.639 8970.401 0.413 0.680618
## V50 -25762.703 10934.783 -2.356 0.021236 *
## V51 46844.468 15367.852 3.048 0.003233 **
## V52 -47783.626 18069.344 -2.644 0.010065 *
## V53 26233.604 18822.491 1.394 0.167744
## V54 87.825 17403.836 0.005 0.995988
## V55 -8475.119 13232.005 -0.641 0.523908
## V56 3488.507 7228.428 0.483 0.630858
## V57 -1520.733 4988.093 -0.305 0.761355
## V58 2275.175 5495.630 0.414 0.680124
## V59 -5415.427 5721.475 -0.947 0.347099
## V60 7152.015 4754.317 1.504 0.136935
## V61 -4494.234 4512.937 -0.996 0.322702
## V62 3662.045 4811.634 0.761 0.449129
## V63 13993.987 7098.106 1.972 0.052563 .
## V64 -23252.133 8973.839 -2.591 0.011604 *
## V65 4373.731 10048.591 0.435 0.664695
## V66 4580.913 10146.146 0.451 0.653011
## V67 -837.676 10747.974 -0.078 0.938097
## V68 -7074.425 10852.430 -0.652 0.516587
## V69 9506.571 9739.256 0.976 0.332325
## V70 -2765.100 9519.031 -0.290 0.772295
## V71 -1125.135 8586.061 -0.131 0.896113
## V72 -7295.096 7489.488 -0.974 0.333341
## V73 17059.811 6522.093 2.616 0.010870 *
## V74 -9889.553 6543.945 -1.511 0.135162
## V75 -325.615 6125.973 -0.053 0.957759
## V76 782.219 5421.002 0.144 0.885677
## V77 8058.935 5793.416 1.391 0.168554
## V78 -15869.978 6448.208 -2.461 0.016282 *
## V79 21768.619 6435.678 3.382 0.001172 **
## V80 -28338.145 8180.874 -3.464 0.000906 ***
## V81 8523.317 10053.153 0.848 0.399384
## V82 22319.451 12098.046 1.845 0.069226 .
## V83 -17244.722 13991.685 -1.232 0.221829
## V84 -18325.836 14959.964 -1.225 0.224627
## V85 33345.457 13868.197 2.404 0.018808 *
## V86 -7955.157 14571.278 -0.546 0.586813
## V87 -7837.966 16141.553 -0.486 0.628762
## V88 -1815.552 17261.928 -0.105 0.916532
## V89 631.595 15684.751 0.040 0.967992
## V90 -2701.955 16187.612 -0.167 0.867911
## V91 4375.678 19400.005 0.226 0.822199
## V92 12925.188 16456.244 0.785 0.434816
## V93 -7441.235 12417.883 -0.599 0.550923
## V94 -2464.532 11815.234 -0.209 0.835366
## V95 -2090.635 9666.576 -0.216 0.829394
## V96 10912.352 9950.716 1.097 0.276505
## V97 -20331.405 11022.234 -1.845 0.069270 .
## V98 3948.443 8227.133 0.480 0.632753
## V99 6358.930 8652.372 0.735 0.464800
## V100 -263.365 4104.463 -0.064 0.949019
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.074 on 71 degrees of freedom
## Multiple R-squared: 0.997, Adjusted R-squared: 0.9928
## F-statistic: 237.5 on 100 and 71 DF, p-value: < 2.2e-16El valor \(R^2_{ajustado}\) obtenido es muy alto (0.9928) lo que indica que el modelo es capaz de predecir con gran exactitud el contenido en grasa de las observaciones con las que se ha entrenado. El hecho de que el modelo en conjunto sea significativo (p-value: < 2.2\(e^{-16}\)), pero que muy pocos de los predictores lo sean a nivel individual, es un indicativo de una posible redundancia entre los predictores (colinealidad).
¿Cómo de bueno es el modelo prediciendo nuevas observaciones que no han participado en ajuste? Al tratarse de un modelo de regresión, la estimación del error de predicción se obtiene mediante el \(Mean\,\, Square\,\, Error\,\, (MSE)\).
\[MSE = \frac{1}{2} \displaystyle\sum^n_{i=1} (\hat{y_i} - y_i)^2\]
# MSE empleando las observaciones de entrenamiento
training_mse <- mean((modelo$fitted.values - training$fat)^2)
training_mse## [1] 0.4765372# MSE empleando nuevas observaciones
predicciones <- predict(modelo, newdata = test)
test_mse <- mean((predicciones - test$fat)^2)
test_mse## [1] 14.54659Se observa que el modelo tiene un MSE muy bajo (0.48) cuando predice las mismas observaciones con las que se ha entrenado, pero 30 veces más alto (14.54) al predecir nuevas observaciones. Esto significa que el modelo no es útil, ya que el objetivo es aplicarlo para predecir el contenido en grasa de futuras muestras de carne. A este problema se le conoce como \(overfitting\). Una de las causas por las que un modelo puede sufrir \(overfitting\) es la incorporación de predictores innecesarios, que no aportan información o que la información que aportan es redundante.
Se recurre en primer lugar a la selección de predictores mediante \(stepwise\) selection empleando el AIC como criterio de evaluación:
modelo_step_selection <- step(object = modelo, trace = FALSE)
# Número de predictores del modelo resultante
length(modelo_step_selection$coefficients)## [1] 73# Training-MSE
training_mse <- mean((modelo_step_selection$fitted.values - training$fat)^2)
training_mse## [1] 0.5034001El proceso de \(stepwise\) selection devuelve como mejor modelo el formado por 73 de los 100 predictores disponibles. Al haber eliminado predictores del modelo, el \(training\,\, MSE\) siempre aumenta, en este caso de 0.48 a 0.05, pero el \(test-MSE\) se ha reducido a 12.88986.
Véase ahora el resultado si se ajusta el modelo empleando las componentes principales:
# Cálculo de componentes principales. Se excluye la columna con la variable
# respuesta *fat*
pca <- prcomp(training[, -101], scale. = TRUE)
# Se muestra la proporción de varianza explicada y acumulada de las 9 primeras
# componentes
summary(pca)$importance[, 1:9]## PC1 PC2 PC3 PC4 PC5
## Standard deviation 9.92492 1.043606 0.5357885 0.3312792 0.07898436
## Proportion of Variance 0.98504 0.010890 0.0028700 0.0011000 0.00006000
## Cumulative Proportion 0.98504 0.995930 0.9988000 0.9999000 0.99996000
## PC6 PC7 PC8 PC9
## Standard deviation 0.04974461 0.02700194 0.02059129 0.008603878
## Proportion of Variance 0.00002000 0.00001000 0.00000000 0.000000000
## Cumulative Proportion 0.99999000 0.99999000 1.00000000 1.000000000El estudio de la proporción de varianza explicada muestra que la primera componente recoge la mayor parte de la información (98.5\(\%\)), decayendo drásticamente la varianza en las sucesivas componentes.
Una vez obtenido el valor de las componentes para cada observación
(\(principal\,\, component\,\,
scores\)), puede ajustarse el modelo lineal empleando dichos
valores junto con la variable respuesta que le corresponde a cada
observación. Con la función pcr() del paquete
pls se evita tener que codificar cada uno de los pasos
intermedios.
Acorde a la proporción de varianza acumulada, emplear las 4 primeras componentes podría ser una buena elección, ya que en conjunto explican el 99.99100 \(\%\) de varianza.
library(pls)##
## Attaching package: 'pls'## The following object is masked from 'package:stats':
##
## loadingsmodelo_pcr <- pcr(formula = fat ~ ., data = training, scale. = TRUE, ncomp = 4)
# Test-MSE
predicciones <- predict(modelo_pcr, newdata = test, ncomp = 4)
test_mse <- mean((predicciones - test$fat)^2)
test_mse## [1] 20.55699El \(test-MSE\) obtenido (20.56) para el modelo que emplea como predictores las 4 primeras componentes es mucho mayor que el obtenido con el modelo generado por \(stepwise\) selection (12.89) e incluso que el obtenido incluyendo todos los predictores (14.54659). Esto significa que, o bien el hecho de emplear componentes principales como predictores no es útil para este caso, o que el número de componentes incluido no es el adecuado.
La función pcr() incluye la posibilidad de recurrir a
\(cross\,\, validation\) para
identificar el número óptimo de componentes con el que se minimiza el
MSE.
set.seed(123)
modelo_pcr <- pcr(formula = fat ~ ., data = training, scale. = TRUE,
validation = "CV")
modelo_pcr_CV <- MSEP(modelo_pcr, estimate = "CV")
which.min(modelo_pcr_CV$val)## [1] 20par(mfrow = c(1,2))
plot(modelo_pcr_CV$val, main = "MSE vs nº componentes", type = "l",
ylab = "MSE",
col = "blue", xlab = "Componentes")
plot(modelo_pcr_CV$val, main = "zoom", type = "l", ylab = "MSE",
xlab = "Componentes", col = "blue", ylim = c(0,20))
# Test-MSE
predicciones <- predict(modelo_pcr, newdata = test, ncomp = 18)
test_mse <- mean((predicciones - test$fat)^2)
test_mse## [1] 4.524698El número óptimo de componentes principales identificado por \(cross\,\, validation\) es de 18. Empleando este número en la \(PCR\) se consigue reducir el \(test-MSE\) a 4.52, un valor muy por debajo del conseguido con los otros modelos.