En primer lugar, cargamos las librerías necesarias

if (!require(tidyverse)) install.packages('tidyverse', dependencies = T)
## Loading required package: tidyverse
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.1     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.2     ✔ tibble    3.2.1
## ✔ lubridate 1.9.2     ✔ tidyr     1.3.0
## ✔ purrr     1.0.1     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
if (!require(caret)) install.packages('caret', dependencies = T)
## Loading required package: caret
## Loading required package: lattice
## 
## Attaching package: 'caret'
## 
## The following object is masked from 'package:purrr':
## 
##     lift
if (!require(MASS)) install.packages('MASS', dependencies = T)
## Loading required package: MASS
## 
## Attaching package: 'MASS'
## 
## The following object is masked from 'package:dplyr':
## 
##     select
if (!require(boot)) install.packages('boot', dependencies = T)
## Loading required package: boot
## 
## Attaching package: 'boot'
## 
## The following object is masked from 'package:lattice':
## 
##     melanoma
library(tidyverse) 
library(caret)
library(MASS)
library(boot)

A continuación, examinamos la base de datos

data("Boston")
sample_n(Boston, 3)

Creación del modelo

En primer lugar, creamos un modelo de regresión con todas las variables de la base de datos

lm.fit1 <- lm(medv ~., data = Boston)
summary (lm.fit1)
## 
## Call:
## lm(formula = medv ~ ., data = Boston)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.595  -2.730  -0.518   1.777  26.199 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.646e+01  5.103e+00   7.144 3.28e-12 ***
## crim        -1.080e-01  3.286e-02  -3.287 0.001087 ** 
## zn           4.642e-02  1.373e-02   3.382 0.000778 ***
## indus        2.056e-02  6.150e-02   0.334 0.738288    
## chas         2.687e+00  8.616e-01   3.118 0.001925 ** 
## nox         -1.777e+01  3.820e+00  -4.651 4.25e-06 ***
## rm           3.810e+00  4.179e-01   9.116  < 2e-16 ***
## age          6.922e-04  1.321e-02   0.052 0.958229    
## dis         -1.476e+00  1.995e-01  -7.398 6.01e-13 ***
## rad          3.060e-01  6.635e-02   4.613 5.07e-06 ***
## tax         -1.233e-02  3.760e-03  -3.280 0.001112 ** 
## ptratio     -9.527e-01  1.308e-01  -7.283 1.31e-12 ***
## black        9.312e-03  2.686e-03   3.467 0.000573 ***
## lstat       -5.248e-01  5.072e-02 -10.347  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.745 on 492 degrees of freedom
## Multiple R-squared:  0.7406, Adjusted R-squared:  0.7338 
## F-statistic: 108.1 on 13 and 492 DF,  p-value: < 2.2e-16

Con la técnica forward selection ajustamos el modelo

Boston2 <- Boston [,-3]
lm.fit2 <- lm(medv ~., data = Boston2)
summary (lm.fit2)
## 
## Call:
## lm(formula = medv ~ ., data = Boston2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.587  -2.737  -0.506   1.742  26.212 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.636e+01  5.091e+00   7.143 3.30e-12 ***
## crim        -1.084e-01  3.281e-02  -3.304 0.001022 ** 
## zn           4.593e-02  1.364e-02   3.368 0.000816 ***
## chas         2.716e+00  8.562e-01   3.173 0.001605 ** 
## nox         -1.743e+01  3.681e+00  -4.735 2.87e-06 ***
## rm           3.797e+00  4.158e-01   9.132  < 2e-16 ***
## age          6.971e-04  1.320e-02   0.053 0.957898    
## dis         -1.490e+00  1.948e-01  -7.648 1.08e-13 ***
## rad          2.999e-01  6.367e-02   4.710 3.22e-06 ***
## tax         -1.178e-02  3.378e-03  -3.489 0.000529 ***
## ptratio     -9.471e-01  1.296e-01  -7.308 1.10e-12 ***
## black        9.282e-03  2.682e-03   3.461 0.000586 ***
## lstat       -5.235e-01  5.052e-02 -10.361  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.741 on 493 degrees of freedom
## Multiple R-squared:  0.7406, Adjusted R-squared:  0.7343 
## F-statistic: 117.3 on 12 and 493 DF,  p-value: < 2.2e-16
Boston3 <- Boston2 [,-6]
lm.fit3 <- lm(medv ~., data = Boston3)
summary (lm.fit3)
## 
## Call:
## lm(formula = medv ~ ., data = Boston3)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.5984  -2.7386  -0.5046   1.7273  26.2373 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  36.341145   5.067492   7.171 2.73e-12 ***
## crim         -0.108413   0.032779  -3.307 0.001010 ** 
## zn            0.045845   0.013523   3.390 0.000754 ***
## chas          2.718716   0.854240   3.183 0.001551 ** 
## nox         -17.376023   3.535243  -4.915 1.21e-06 ***
## rm            3.801579   0.406316   9.356  < 2e-16 ***
## dis          -1.492711   0.185731  -8.037 6.84e-15 ***
## rad           0.299608   0.063402   4.726 3.00e-06 ***
## tax          -0.011778   0.003372  -3.493 0.000521 ***
## ptratio      -0.946525   0.129066  -7.334 9.24e-13 ***
## black         0.009291   0.002674   3.475 0.000557 ***
## lstat        -0.522553   0.047424 -11.019  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.736 on 494 degrees of freedom
## Multiple R-squared:  0.7406, Adjusted R-squared:  0.7348 
## F-statistic: 128.2 on 11 and 494 DF,  p-value: < 2.2e-16

El modelo propuesto está compuesto por el intercepto y las variables crim, zn, chas, nox, rm, dis, rad, tax, ptratio, black y lstat explica un 73.48% de la varianza. En este sentido, observamos variables que afectan de manera positiva ( zn, chas, rm, rad y black) y otras de manera negativa ( crim, nox, dis, tax, ptratio y lstat)

Así pues, el modelo de regresión lineal propuesto es:

medv = 36.341145 - (crim x 0.108413) + (zn x 0.045845) + (chas x 2.718716) - (nox x 17.376023) + (rm x 3.801579) - (dis x 1.492711) + (rad x 0.299608) - (tax x 0.011778) - (ptratio x 0.946525) + (black x 0.009291) - (lstat x 0.522553)

En concreto, y tal como podemos observar, la variable crim tiene un peso negativo de -0.108413 y un p-valor de 0.001010.

A continuación renombramos la base de datos y la regresión del modelo con el que trabajaremos

DF.Modelo <- Boston3
lm.Modelo <- lm(medv ~ ., data = DF.Modelo)

Posteriormente realizamos las predicciones y calculamos R2, RMSE y MAE

predictions <- lm.Modelo %>% predict(DF.Modelo)
data.frame(R2 = R2(predictions, DF.Modelo$medv), 
            RMSE = RMSE(predictions,DF.Modelo$medv), 
            MAE = MAE(predictions, DF.Modelo$medv))

Estimamos la tasa de error de predicción. Esta ha de ser lo más baja posible

RMSE(predictions, DF.Modelo$medv)/mean(DF.Modelo$medv)
## [1] 0.2076854

Método de validación Leave One Out Cross-Validation (LOOCV)

En primer lugar, especificamos el modelo de validación Leave One Out Cross-Validation (LOOCV) y, posteriormente, realizaremos la validación del modelo.

train.control <- trainControl(method = "LOOCV")
model_LOOCV <- train(medv ~ . , data = DF.Modelo, method = "lm",
               trControl = train.control)
print(model_LOOCV)
## Linear Regression 
## 
## 506 samples
##  11 predictor
## 
## No pre-processing
## Resampling: Leave-One-Out Cross-Validation 
## Summary of sample sizes: 505, 505, 505, 505, 505, 505, ... 
## Resampling results:
## 
##   RMSE      Rsquared   MAE     
##   4.849046  0.7216232  3.370289
## 
## Tuning parameter 'intercept' was held constant at a value of TRUE

Tal como podemos observar, el modelo con el método de validación Leave One Out Cross-Validation (LOOCV) muestra una R2 de 0.72162323. Esto es ligeramente inferior al modelo orginal. Tambiés comprombamos que es superior en el RMSE.

A continuación, observamos como el peso de la variable crim es exctamente el mismo.

summary(model_LOOCV)
## 
## Call:
## lm(formula = .outcome ~ ., data = dat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.5984  -2.7386  -0.5046   1.7273  26.2373 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  36.341145   5.067492   7.171 2.73e-12 ***
## crim         -0.108413   0.032779  -3.307 0.001010 ** 
## zn            0.045845   0.013523   3.390 0.000754 ***
## chas          2.718716   0.854240   3.183 0.001551 ** 
## nox         -17.376023   3.535243  -4.915 1.21e-06 ***
## rm            3.801579   0.406316   9.356  < 2e-16 ***
## dis          -1.492711   0.185731  -8.037 6.84e-15 ***
## rad           0.299608   0.063402   4.726 3.00e-06 ***
## tax          -0.011778   0.003372  -3.493 0.000521 ***
## ptratio      -0.946525   0.129066  -7.334 9.24e-13 ***
## black         0.009291   0.002674   3.475 0.000557 ***
## lstat        -0.522553   0.047424 -11.019  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.736 on 494 degrees of freedom
## Multiple R-squared:  0.7406, Adjusted R-squared:  0.7348 
## F-statistic: 128.2 on 11 and 494 DF,  p-value: < 2.2e-16

Método de validación K-Fold

En primer lugar, especificamos el modelo de validación K-Fold y, posteriormente, realizaremos la validación del modelo.

set.seed(123)
train.control <- trainControl(method = "cv", number = 10)
model_K_Fold <- train(medv ~ . , data = DF.Modelo, method = "lm", 
               trControl = train.control)
print(model_K_Fold)
## Linear Regression 
## 
## 506 samples
##  11 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold) 
## Summary of sample sizes: 455, 456, 456, 456, 456, 456, ... 
## Resampling results:
## 
##   RMSE      Rsquared   MAE     
##   4.809318  0.7249622  3.359429
## 
## Tuning parameter 'intercept' was held constant at a value of TRUE

A continuación, observamos como el peso de la variable crim es exctamente el mismo.

summary(model_K_Fold)
## 
## Call:
## lm(formula = .outcome ~ ., data = dat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.5984  -2.7386  -0.5046   1.7273  26.2373 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  36.341145   5.067492   7.171 2.73e-12 ***
## crim         -0.108413   0.032779  -3.307 0.001010 ** 
## zn            0.045845   0.013523   3.390 0.000754 ***
## chas          2.718716   0.854240   3.183 0.001551 ** 
## nox         -17.376023   3.535243  -4.915 1.21e-06 ***
## rm            3.801579   0.406316   9.356  < 2e-16 ***
## dis          -1.492711   0.185731  -8.037 6.84e-15 ***
## rad           0.299608   0.063402   4.726 3.00e-06 ***
## tax          -0.011778   0.003372  -3.493 0.000521 ***
## ptratio      -0.946525   0.129066  -7.334 9.24e-13 ***
## black         0.009291   0.002674   3.475 0.000557 ***
## lstat        -0.522553   0.047424 -11.019  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.736 on 494 degrees of freedom
## Multiple R-squared:  0.7406, Adjusted R-squared:  0.7348 
## F-statistic: 128.2 on 11 and 494 DF,  p-value: < 2.2e-16

Podemos observar que se mejora ligeramewnte el resultado obtenido respecto a LOOCV ya que el RMSE es ligeramente inferior y a la vez se mejora el grado de calidad de ajuste al aincrementarse el R2.

Método de validación K-Fold Repetido

En primer lugar, especificamos el modelo de validación K-Fold y, posteriormente, realizaremos la validación del modelo.

set.seed(123)
train.control <- trainControl(method = "repeatedcv", number = 10, 
                              repeats = 3)
model_K_Fold_Repetido <- train(medv ~ . , data = DF.Modelo, method = "lm", 
               trControl = train.control)
print(model_K_Fold_Repetido)
## Linear Regression 
## 
## 506 samples
##  11 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold, repeated 3 times) 
## Summary of sample sizes: 455, 456, 456, 456, 456, 456, ... 
## Resampling results:
## 
##   RMSE      Rsquared   MAE     
##   4.777001  0.7315143  3.373325
## 
## Tuning parameter 'intercept' was held constant at a value of TRUE

Nuevamente se produce una ligera mejoría del RMSE ya se gana calidad de ajuste al incrementarse el R2

summary(model_K_Fold_Repetido)
## 
## Call:
## lm(formula = .outcome ~ ., data = dat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.5984  -2.7386  -0.5046   1.7273  26.2373 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  36.341145   5.067492   7.171 2.73e-12 ***
## crim         -0.108413   0.032779  -3.307 0.001010 ** 
## zn            0.045845   0.013523   3.390 0.000754 ***
## chas          2.718716   0.854240   3.183 0.001551 ** 
## nox         -17.376023   3.535243  -4.915 1.21e-06 ***
## rm            3.801579   0.406316   9.356  < 2e-16 ***
## dis          -1.492711   0.185731  -8.037 6.84e-15 ***
## rad           0.299608   0.063402   4.726 3.00e-06 ***
## tax          -0.011778   0.003372  -3.493 0.000521 ***
## ptratio      -0.946525   0.129066  -7.334 9.24e-13 ***
## black         0.009291   0.002674   3.475 0.000557 ***
## lstat        -0.522553   0.047424 -11.019  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.736 on 494 degrees of freedom
## Multiple R-squared:  0.7406, Adjusted R-squared:  0.7348 
## F-statistic: 128.2 on 11 and 494 DF,  p-value: < 2.2e-16

Método Bootstrap

Utilizarmeos ahora el método de bootstrap

train.control <- trainControl(method = "boot", number = 506)
model_Bootstrap <- train(medv ~ . , data = DF.Modelo, method = "lm", 
               trControl = train.control)
print(model_Bootstrap)
## Linear Regression 
## 
## 506 samples
##  11 predictor
## 
## No pre-processing
## Resampling: Bootstrapped (506 reps) 
## Summary of sample sizes: 506, 506, 506, 506, 506, 506, ... 
## Resampling results:
## 
##   RMSE      Rsquared   MAE     
##   4.921285  0.7183245  3.441539
## 
## Tuning parameter 'intercept' was held constant at a value of TRUE

En este caso el resultado es peor que en todos los métodos anteiores, tanto en _RMSE como en R2

A continuación observamos el peso de cada variable para analizar crim. Tal como podemos observar el peso es de -2.013307245, que es bastante superior al modelo original.

model_coef <- function(data, index)
  {coef(lm(medv ~., data = DF.Modelo, subset = index))}
model_coef(DF.Modelo, 1:506)
##   (Intercept)          crim            zn          chas           nox 
##  36.341145004  -0.108413345   0.045844929   2.718716303 -17.376023429 
##            rm           dis           rad           tax       ptratio 
##   3.801578840  -1.492711460   0.299608454  -0.011777973  -0.946524570 
##         black         lstat 
##   0.009290845  -0.522553457

A contiunación, calcularemos los errores estándar usando 5000 estimaciones de los coeficientes mediante bootstrap:

boot(DF.Modelo, model_coef, 5000)
## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = DF.Modelo, statistic = model_coef, R = 5000)
## 
## 
## Bootstrap Statistics :
##           original        bias    std. error
## t1*   36.341145004 -6.070494e-01 7.852530488
## t2*   -0.108413345  4.381154e-03 0.033974838
## t3*    0.045844929 -1.024345e-03 0.013475915
## t4*    2.718716303 -3.033303e-02 1.308624982
## t5*  -17.376023429  2.033699e-01 3.322185385
## t6*    3.801578840  5.308337e-02 0.796774251
## t7*   -1.492711460  1.347463e-02 0.220358365
## t8*    0.299608454 -3.201716e-03 0.062601655
## t9*   -0.011777973 -5.010388e-05 0.002910928
## t10*  -0.946524570  5.702314e-03 0.114809286
## t11*   0.009290845  5.557791e-05 0.002800442
## t12*  -0.522553457  1.022857e-03 0.085371309

Con estos valores, podemos construir intervalos de confianza para los coeficientes. Por ejemplo, el intervalo de confianza del 95% para un crim se define como Crim±2*SE(crim):

• Límite inferior = -0.108413345 – (2*0.034977083) = -0.178367511

• Límite superior = -0.108413345 + (2*0.034977083) = -0.038459179

Es decir, hay aproximadamente un 95% de probabilidad de que el intervalo [-0.178367511, -0.038459179] contenga el verdadero valor de la variables crim.

Conclusiones

A continuación mostramos una tabla en la que exponemos los datos obtenidos

Modelo RMSE R Squared Mae
LOOCV 4.849046 0.7216232 3.370289
K Fold 4.809318 0.7249622 3.359429
K Fold repetido 4.777001 0.7315143 3.373325
Bootstrap 4.921285 0.7183245 3.441539

Tal como hemos venido desarrollando, cada método de validación nos ha reportado reusltados diferentes. En este sentido, el método de validación K-Fold Repetido ofrecía el RMSE más bajo y,a su vez, el R2 más alto. Por tanto, podemos concluir que, en este caso, es el método más adecuado.