library(faraway)
data(meatspec)
dim(meatspec)[1] 215 1011 Introducción
Este ejemplo está basado en una excelente publicación de Joaquín Amat Rodrigo y que puede ser consultada en este enlace.
El objetivo de este ejemplo es mostrar la utilidad que tiene en análisis de componentes principales (PCA) al crear un modelo de regresión.
2 Enunciado del problema
El contenido en grasa de la carne pude medirse 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.
El objetivo es crear los siguientes modelos:
- Modelo de regresión usando las variables originales \(X_1, X_2, \ldots, X_{100}\).
- Modelo de regresión usando las componentes principales \(CP_1, CP_2, \ldots, CP_{k}\).
La comparación del desempeño de los modelos se hará por medio del
\[ MSE = \frac{1}{n}\displaystyle\sum_{i=1}^n ( \hat{y_i} - y_i)^2 \]
Los datos para el ejemplo se pueden obtener de la siguiente manera.
Para entrenar los modelos y comparar los desempeños vamos a particionar los datos así:
training <- meatspec[1:172, ] # Datos de entrenamiento
test <- meatspec[173:215, ] # Datos de prueba3 Modelo de regresión usando las variables originales
Vamos a ajustar el modelo 1 así:
modelo1 <- lm(fat ~ ., data = training)
summary(modelo1)
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-16MSE empleando los datos de entrenamiento.
predicciones <- predict(modelo1, newdata = training)
training_mse_1 <- mean((predicciones - training$fat)^2)
training_mse_1[1] 0.4765372MSE empleando los datos de prueba.
predicciones <- predict(modelo1, newdata = test)
test_mse_1 <- mean((predicciones - test$fat)^2)
test_mse_1[1] 14.546594 Modelo de regresión usando las componentes principales
Vamos a crear un modelo de regresión usando \(CP_1, CP_2, \ldots, CP_{k}\). Para averiguar el valor apropiado de \(k\) usaremos cross-validation de la siguiente forma.
library(pls)
Attaching package: 'pls'The following object is masked from 'package:stats':
loadingsset.seed(123)
modelo2 <- pcr(formula = fat ~ ., data = training,
scale = TRUE, center = TRUE,
validation = "CV")
modelo_2_CV <- MSEP(modelo2, estimate = "CV")
which.min(modelo_2_CV$val)[1] 22De la salida anterior se tiene que el número \(k=22\), es decir, conviene utilizar 22 componentes principales.
En la siguiente figura se muestra el valor de MSE para todos los valores posibles de \(k\).
par(mfrow = c(1,2))
plot(modelo_2_CV$val, main = "MSE vs nº componentes",
type = "l", ylab = "MSE", las=1,
col = "tomato", xlab = "Componentes")
plot(modelo_2_CV$val, main = "zoom",
type = "l", ylab = "MSE", las=1,
xlab = "Componentes", col = "tomato", ylim = c(0,20))
Vamos ahora a calcular el MSE con el modelo2.
MSE empleando los datos de entrenamiento.
predicciones <- predict(modelo2, newdata = training, ncomp = 22)
training_mse_2 <- mean((predicciones - training$fat)^2)
training_mse_2[1] 3.748703MSE empleando los datos de prueba.
predicciones <- predict(modelo2, newdata = test, ncomp = 22)
test_mse_2 <- mean((predicciones - test$fat)^2)
test_mse_2[1] 4.6540595 Comparación de los resultados
En la siguiente tabla se muestran los resultados de ambos modelos.
| Modelo | Conjunto de datos | MSE |
|---|---|---|
| 1 | Train | 0.4765372 |
| 1 | Test | 14.5465923 |
| 2 | Train | 3.7487028 |
| 2 | Test | 4.6540594 |