Crear, evaluar y comparar un modelo de regresión lineal simple y un moldeo polinómico para predecir ???
La regresión lineal simple apoya a la correlación estimando coeficientes para hacer predicciones.
\[ Y = a + bx \]
\[ Y = \beta_0 + \beta_1\cdot x_i + \epsilon \]
\[ Y = \beta0 + \beta_1\cdot x^2+\beta_2\cdot x^3 + ... \beta_n \cdot x^n \]
library(ggplot2)
# library(plotly) # no se está usando
library(knitr)
library(PerformanceAnalytics) # Para coorelaciones gráficas
library(caret) # Para particionar
library(Metrics) # Para determinar rmse
datos <- read.csv("https://raw.githubusercontent.com/rpizarrog/Machine-Learning-con-R/main/datos/Advertising.csv")
str(datos)
## 'data.frame': 200 obs. of 5 variables:
## $ X : int 1 2 3 4 5 6 7 8 9 10 ...
## $ TV : num 230.1 44.5 17.2 151.5 180.8 ...
## $ Radio : num 37.8 39.3 45.9 41.3 10.8 48.9 32.8 19.6 2.1 2.6 ...
## $ Newspaper: num 69.2 45.1 69.3 58.5 58.4 75 23.5 11.6 1 21.2 ...
## $ Sales : num 22.1 10.4 9.3 18.5 12.9 7.2 11.8 13.2 4.8 10.6 ...
summary(datos)
## X TV Radio Newspaper
## Min. : 1.00 Min. : 0.70 Min. : 0.000 Min. : 0.30
## 1st Qu.: 50.75 1st Qu.: 74.38 1st Qu.: 9.975 1st Qu.: 12.75
## Median :100.50 Median :149.75 Median :22.900 Median : 25.75
## Mean :100.50 Mean :147.04 Mean :23.264 Mean : 30.55
## 3rd Qu.:150.25 3rd Qu.:218.82 3rd Qu.:36.525 3rd Qu.: 45.10
## Max. :200.00 Max. :296.40 Max. :49.600 Max. :114.00
## Sales
## Min. : 1.60
## 1st Qu.:10.38
## Median :12.90
## Mean :14.02
## 3rd Qu.:17.40
## Max. :27.00
Diagrama de dispersión
source("https://raw.githubusercontent.com/rpizarrog/Machine-Learning-con-R/main/funciones/simular_women")
f_dispersion(data.frame(datos$TV, datos$Sales))
Seleccionar las columnas de interes
datos <- select(datos, TV, Sales)
En caso necesario
Aleatoriamente se repate las observaciones con el 70% para datos de entrenamiento y el 30% para datos de validación.
Sembrar una semilla con set.seed()
set.seed(2022)
n <- nrow(datos) # cantidad de observaciones
entrena <- createDataPartition(y = datos$Sales, p = 0.70, list = FALSE, times = 1)
# Datos entrenamiento
datos.entrenamiento <- datos[entrena, ] # [renglones, columna]
# Datos validación
datos.validacion <- datos[-entrena, ]
datos.entrenamiento
## TV Sales
## 1 230.1 22.1
## 2 44.5 10.4
## 3 17.2 9.3
## 4 151.5 18.5
## 5 180.8 12.9
## 6 8.7 7.2
## 7 57.5 11.8
## 8 120.2 13.2
## 9 8.6 4.8
## 11 66.1 8.6
## 12 214.7 17.4
## 13 23.8 9.2
## 14 97.5 9.7
## 15 204.1 19.0
## 16 195.4 22.4
## 17 67.8 12.5
## 18 281.4 24.4
## 19 69.2 11.3
## 23 13.2 5.6
## 25 62.3 9.7
## 28 240.1 15.9
## 29 248.8 18.9
## 32 112.9 11.9
## 37 266.9 25.4
## 38 74.7 14.7
## 39 43.1 10.1
## 40 228.0 21.5
## 41 202.5 16.6
## 43 293.6 20.7
## 44 206.9 12.9
## 45 25.1 8.5
## 46 175.1 14.9
## 47 89.7 10.6
## 49 227.2 14.8
## 51 199.8 11.4
## 52 100.4 10.7
## 53 216.4 22.6
## 56 198.9 23.7
## 57 7.3 5.5
## 58 136.2 13.2
## 59 210.8 23.8
## 61 53.5 8.1
## 62 261.3 24.2
## 64 102.7 14.0
## 65 131.1 18.0
## 68 139.3 13.4
## 70 216.8 22.3
## 71 199.1 18.3
## 72 109.8 12.4
## 73 26.8 8.8
## 74 129.4 11.0
## 75 213.4 17.0
## 77 27.5 6.9
## 78 120.5 14.2
## 79 5.4 5.3
## 80 116.0 11.0
## 81 76.4 11.8
## 82 239.8 12.3
## 85 213.5 21.7
## 86 193.2 15.2
## 88 110.7 16.0
## 89 88.3 12.9
## 91 134.3 11.2
## 92 28.6 7.3
## 94 250.9 22.2
## 96 163.3 16.9
## 97 197.6 11.7
## 98 184.9 15.5
## 99 289.7 25.4
## 100 135.2 17.2
## 102 296.4 23.8
## 103 280.2 14.8
## 105 238.2 20.7
## 106 137.9 19.2
## 108 90.4 8.7
## 110 255.4 19.8
## 113 175.7 14.1
## 115 78.2 14.6
## 117 139.2 12.2
## 119 125.7 15.9
## 121 141.3 15.5
## 122 18.8 7.0
## 123 224.0 11.6
## 124 123.1 15.2
## 126 87.2 10.6
## 127 7.8 6.6
## 129 220.3 24.7
## 131 0.7 1.6
## 132 265.2 12.7
## 134 219.8 19.6
## 135 36.9 10.8
## 136 48.3 11.6
## 137 25.6 9.5
## 138 273.7 20.8
## 141 73.4 10.9
## 142 193.7 19.2
## 143 220.5 20.1
## 145 96.2 11.4
## 146 140.3 10.3
## 147 240.1 13.2
## 148 243.2 25.4
## 150 44.7 10.1
## 152 121.0 11.6
## 154 171.3 19.0
## 155 187.8 15.6
## 156 4.1 3.2
## 157 93.9 15.3
## 158 149.8 10.1
## 159 11.7 7.3
## 160 131.7 12.9
## 161 172.5 14.4
## 162 85.7 13.3
## 163 188.4 14.9
## 164 163.5 18.0
## 165 117.2 11.9
## 166 234.5 11.9
## 167 17.9 8.0
## 168 206.8 12.2
## 169 215.4 17.1
## 170 284.3 15.0
## 173 19.6 7.6
## 174 168.4 11.7
## 175 222.4 11.5
## 176 276.9 27.0
## 179 276.7 11.8
## 180 165.6 12.6
## 181 156.6 10.5
## 183 56.2 8.7
## 184 287.6 26.2
## 185 253.8 17.6
## 186 205.0 22.6
## 187 139.5 10.3
## 188 191.1 17.3
## 189 286.0 15.9
## 190 18.7 6.7
## 192 75.5 9.9
## 193 17.2 5.9
## 195 149.7 17.3
## 196 38.2 7.6
## 197 94.2 9.7
## 198 177.0 12.8
## 200 232.1 13.4
datos.validacion
## TV Sales
## 10 199.8 10.6
## 20 147.3 14.6
## 21 218.4 18.0
## 22 237.4 12.5
## 24 228.3 15.5
## 26 262.9 12.0
## 27 142.9 15.0
## 30 70.6 10.5
## 31 292.9 21.4
## 33 97.2 9.6
## 34 265.6 17.4
## 35 95.7 9.5
## 36 290.7 12.8
## 42 177.0 17.1
## 48 239.9 23.2
## 50 66.9 9.7
## 54 182.6 21.2
## 55 262.7 20.2
## 60 210.7 18.4
## 63 239.3 15.7
## 66 69.0 9.3
## 67 31.5 9.5
## 69 237.4 18.9
## 76 16.9 8.7
## 83 75.3 11.3
## 84 68.4 13.6
## 87 76.3 12.0
## 90 109.8 16.7
## 93 217.7 19.4
## 95 107.4 11.5
## 101 222.4 11.7
## 104 187.9 14.7
## 107 25.0 7.2
## 109 13.1 5.3
## 111 225.8 13.4
## 112 241.7 21.8
## 114 209.6 15.9
## 116 75.1 12.6
## 118 76.4 9.4
## 120 19.4 6.6
## 125 229.5 19.7
## 128 80.2 8.8
## 130 59.6 9.7
## 133 8.4 5.7
## 139 43.0 9.6
## 140 184.9 20.7
## 144 104.6 10.4
## 149 38.0 10.9
## 151 280.7 16.1
## 153 197.6 16.6
## 171 50.0 8.4
## 172 164.5 14.5
## 177 248.4 20.2
## 178 170.2 11.7
## 182 218.5 12.2
## 191 39.5 10.8
## 194 166.8 19.6
## 199 283.6 25.5
El modelo se construye con datos de entrenamiento ### Modelo de Regresión lineal simple
modelo.rs <- lm(data = datos.entrenamiento,
formula = Sales ~ TV)
summary(modelo.rs)
##
## Call:
## lm(formula = Sales ~ TV, data = datos.entrenamiento)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.8752 -1.7706 -0.2902 2.1879 6.9256
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.801923 0.551158 12.34 <2e-16 ***
## TV 0.050138 0.003278 15.30 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.318 on 140 degrees of freedom
## Multiple R-squared: 0.6256, Adjusted R-squared: 0.623
## F-statistic: 234 on 1 and 140 DF, p-value: < 2.2e-16
modelo.rpoly2 <- lm(data = datos.entrenamiento,
formula = Sales ~ poly(TV, degree = 2))
summary(modelo.rpoly2)
##
## Call:
## lm(formula = Sales ~ poly(TV, degree = 2), data = datos.entrenamiento)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.3015 -1.7900 -0.2512 2.1656 6.8913
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 14.0775 0.2779 50.664 <2e-16 ***
## poly(TV, degree = 2)1 50.7569 3.3111 15.330 <2e-16 ***
## poly(TV, degree = 2)2 -4.1996 3.3111 -1.268 0.207
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.311 on 139 degrees of freedom
## Multiple R-squared: 0.6299, Adjusted R-squared: 0.6246
## F-statistic: 118.3 on 2 and 139 DF, p-value: < 2.2e-16
modelo.rpoly5 <- lm(data = datos.entrenamiento,
formula = Sales ~ poly(TV, degree = 5))
summary(modelo.rpoly5)
##
## Call:
## lm(formula = Sales ~ poly(TV, degree = 5), data = datos.entrenamiento)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.2657 -1.7632 -0.1326 1.8232 6.9183
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 14.078 0.276 51.001 <2e-16 ***
## poly(TV, degree = 5)1 50.757 3.289 15.431 <2e-16 ***
## poly(TV, degree = 5)2 -4.200 3.289 -1.277 0.204
## poly(TV, degree = 5)3 4.261 3.289 1.295 0.197
## poly(TV, degree = 5)4 -3.294 3.289 -1.002 0.318
## poly(TV, degree = 5)5 4.845 3.289 1.473 0.143
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.289 on 136 degrees of freedom
## Multiple R-squared: 0.6427, Adjusted R-squared: 0.6295
## F-statistic: 48.92 on 5 and 136 DF, p-value: < 2.2e-16
En cuanto a a métrica R Square que significa la represntatividad estadística de la variable independiete sobre la variable dependiente (TV sobre Sales), el mejor modelo sería el polinómioc de segndo nivel dado que :
Se usa la función predcit() para hacer predicciones con los datos de validación luego se compara con los datos reales.
predicciones <- predict(modelo.rs, newdata = datos.validacion)
predicciones
## 10 20 21 22 24 26 27 30
## 16.819528 14.187274 17.752098 18.704723 18.248466 19.983246 13.966667 10.341677
## 31 33 34 35 36 42 48 50
## 21.487391 11.675352 20.118619 11.600145 21.377087 15.676378 18.830069 10.156166
## 54 55 60 63 66 67 69 76
## 15.957152 19.973219 17.366034 18.799986 10.261456 8.381275 18.704723 7.649258
## 83 84 87 90 93 95 101 104
## 10.577326 10.231373 10.627465 12.307093 17.717001 12.186762 17.952651 16.222884
## 107 109 111 112 114 116 118 120
## 8.055377 7.458733 18.123120 18.920317 17.310882 10.567299 10.632478 7.774603
## 125 128 130 133 139 140 144 149
## 18.308632 10.823003 9.790157 7.223083 8.957864 16.072469 12.046375 8.707173
## 151 153 171 172 177 178 182 191
## 20.875706 16.709224 9.308831 15.049651 19.256243 15.335438 17.757112 8.782380
## 194 199
## 15.164969 21.021106
Construir un data frame para comparar
comparaciones <- data.frame(datos.validacion, predicciones)
comparaciones
## TV Sales predicciones
## 10 199.8 10.6 16.819528
## 20 147.3 14.6 14.187274
## 21 218.4 18.0 17.752098
## 22 237.4 12.5 18.704723
## 24 228.3 15.5 18.248466
## 26 262.9 12.0 19.983246
## 27 142.9 15.0 13.966667
## 30 70.6 10.5 10.341677
## 31 292.9 21.4 21.487391
## 33 97.2 9.6 11.675352
## 34 265.6 17.4 20.118619
## 35 95.7 9.5 11.600145
## 36 290.7 12.8 21.377087
## 42 177.0 17.1 15.676378
## 48 239.9 23.2 18.830069
## 50 66.9 9.7 10.156166
## 54 182.6 21.2 15.957152
## 55 262.7 20.2 19.973219
## 60 210.7 18.4 17.366034
## 63 239.3 15.7 18.799986
## 66 69.0 9.3 10.261456
## 67 31.5 9.5 8.381275
## 69 237.4 18.9 18.704723
## 76 16.9 8.7 7.649258
## 83 75.3 11.3 10.577326
## 84 68.4 13.6 10.231373
## 87 76.3 12.0 10.627465
## 90 109.8 16.7 12.307093
## 93 217.7 19.4 17.717001
## 95 107.4 11.5 12.186762
## 101 222.4 11.7 17.952651
## 104 187.9 14.7 16.222884
## 107 25.0 7.2 8.055377
## 109 13.1 5.3 7.458733
## 111 225.8 13.4 18.123120
## 112 241.7 21.8 18.920317
## 114 209.6 15.9 17.310882
## 116 75.1 12.6 10.567299
## 118 76.4 9.4 10.632478
## 120 19.4 6.6 7.774603
## 125 229.5 19.7 18.308632
## 128 80.2 8.8 10.823003
## 130 59.6 9.7 9.790157
## 133 8.4 5.7 7.223083
## 139 43.0 9.6 8.957864
## 140 184.9 20.7 16.072469
## 144 104.6 10.4 12.046375
## 149 38.0 10.9 8.707173
## 151 280.7 16.1 20.875706
## 153 197.6 16.6 16.709224
## 171 50.0 8.4 9.308831
## 172 164.5 14.5 15.049651
## 177 248.4 20.2 19.256243
## 178 170.2 11.7 15.335438
## 182 218.5 12.2 17.757112
## 191 39.5 10.8 8.782380
## 194 166.8 19.6 15.164969
## 199 283.6 25.5 21.021106
rmse Root Mean Stándar Error, este valor normalmente se compara contra otro modelo y el que esté mas cerca de cero es mejor.
RMSE es una forma útil de ver qué tan bien un modelo de regresión puede ajustarse a un conjunto de datos.
Cuanto mayor sea el rmse, mayor será la diferencia entre los valores predichos y reales, lo que significa que peor se ajusta un modelo de regresión a los datos. Por el contrario, cuanto más pequeño sea el rmse, mejor podrá un modelo ajustar los datos.
\[ rmse = \sqrt{\frac{\sum(predicho_i - real_i)^{2}}{n}} \]
rmse <- rmse(actual = comparaciones$TV, predicted = comparaciones$predicciones)
rmse
## [1] 160.3604
Se usa la función predcit() para hacer predicciones con los datos de validación luego se compara con los datos reales.
predicciones <- predict(modelo.rpoly2, newdata = datos.validacion)
predicciones
## 10 20 21 22 24 26 27 30
## 17.028317 14.569565 17.829207 18.609391 18.240507 19.596234 14.350206 10.451326
## 31 33 34 35 36 42 48 50
## 20.668829 11.950308 19.696681 11.867778 20.593419 15.996468 18.709193 10.236869
## 54 55 60 63 66 67 69 76
## 16.255018 19.588763 17.502112 18.685301 10.358766 8.111550 18.609391 7.196249
## 83 84 87 90 93 95 101 104
## 10.721648 10.323986 10.778861 12.634129 17.799731 12.505177 17.996642 16.496651
## 107 109 111 112 114 116 118 120
## 7.706848 6.954308 18.137625 18.780640 17.454870 10.710193 10.784577 7.354584
## 125 128 130 133 139 140 144 149
## 18.289654 11.000977 9.809488 6.652944 8.816571 16.360244 12.353960 8.511765
## 151 153 171 172 177 178 182 191
## 20.244170 16.931159 9.238838 15.407334 19.043556 15.678037 17.833413 8.603486
## 194 199
## 15.516980 20.346545
Construir un data frame para comparar
comparaciones <- data.frame(datos.validacion, predicciones)
comparaciones
## TV Sales predicciones
## 10 199.8 10.6 17.028317
## 20 147.3 14.6 14.569565
## 21 218.4 18.0 17.829207
## 22 237.4 12.5 18.609391
## 24 228.3 15.5 18.240507
## 26 262.9 12.0 19.596234
## 27 142.9 15.0 14.350206
## 30 70.6 10.5 10.451326
## 31 292.9 21.4 20.668829
## 33 97.2 9.6 11.950308
## 34 265.6 17.4 19.696681
## 35 95.7 9.5 11.867778
## 36 290.7 12.8 20.593419
## 42 177.0 17.1 15.996468
## 48 239.9 23.2 18.709193
## 50 66.9 9.7 10.236869
## 54 182.6 21.2 16.255018
## 55 262.7 20.2 19.588763
## 60 210.7 18.4 17.502112
## 63 239.3 15.7 18.685301
## 66 69.0 9.3 10.358766
## 67 31.5 9.5 8.111550
## 69 237.4 18.9 18.609391
## 76 16.9 8.7 7.196249
## 83 75.3 11.3 10.721648
## 84 68.4 13.6 10.323986
## 87 76.3 12.0 10.778861
## 90 109.8 16.7 12.634129
## 93 217.7 19.4 17.799731
## 95 107.4 11.5 12.505177
## 101 222.4 11.7 17.996642
## 104 187.9 14.7 16.496651
## 107 25.0 7.2 7.706848
## 109 13.1 5.3 6.954308
## 111 225.8 13.4 18.137625
## 112 241.7 21.8 18.780640
## 114 209.6 15.9 17.454870
## 116 75.1 12.6 10.710193
## 118 76.4 9.4 10.784577
## 120 19.4 6.6 7.354584
## 125 229.5 19.7 18.289654
## 128 80.2 8.8 11.000977
## 130 59.6 9.7 9.809488
## 133 8.4 5.7 6.652944
## 139 43.0 9.6 8.816571
## 140 184.9 20.7 16.360244
## 144 104.6 10.4 12.353960
## 149 38.0 10.9 8.511765
## 151 280.7 16.1 20.244170
## 153 197.6 16.6 16.931159
## 171 50.0 8.4 9.238838
## 172 164.5 14.5 15.407334
## 177 248.4 20.2 19.043556
## 178 170.2 11.7 15.678037
## 182 218.5 12.2 17.833413
## 191 39.5 10.8 8.603486
## 194 166.8 19.6 15.516980
## 199 283.6 25.5 20.346545
rmse Root Mean Stándar Error, este valor normalmente se compara contra otro modelo y el que esté mas cerca de cero es mejor.
RMSE es una forma útil de ver qué tan bien un modelo de regresión puede ajustarse a un conjunto de datos.
Cuanto mayor sea el rmse, mayor será la diferencia entre los valores predichos y reales, lo que significa que peor se ajusta un modelo de regresión a los datos. Por el contrario, cuanto más pequeño sea el rmse, mejor podrá un modelo ajustar los datos.
rmse <- rmse(actual = comparaciones$TV, predicted = comparaciones$predicciones)
rmse
## [1] 160.4068
Unicamente se reflejan las tendencias lineal simple y polinómico de segundo nivel.
ggplot(data = datos.entrenamiento) +
geom_point(aes(x = TV, y = Sales), col='blue') +
geom_line(aes(x = TV, y = modelo.rs$fitted.values), col='red') +
geom_line(aes(x = TV, y = modelo.rpoly2$fitted.values), col='green')