Regresion_Lineal_Ajustada
Irving Tolosa
2024
paso 1 - Cargar y explorar datos
Cargamos las biblioteca necesaria y la base de datos Boston, con la cual trabajamos la semana pasada y habiamos visto que no cumpliamos con todos los crterios de linealidad.
library(MASS) # Conjunto de datos Boston
library(carData) #paquete necesario para usar car
library(car) # Para la validación de supuestos
library(ggplot2) # Para visualización
library(lattice) #paquete necesario para usar caret
library(caret) # Validación cruzada y división de datos
datos <- Boston
head(datos)## crim zn indus chas nox rm age dis rad tax ptratio black lstat
## 1 0.00632 18 2.31 0 0.538 6.575 65.2 4.0900 1 296 15.3 396.90 4.98
## 2 0.02731 0 7.07 0 0.469 6.421 78.9 4.9671 2 242 17.8 396.90 9.14
## 3 0.02729 0 7.07 0 0.469 7.185 61.1 4.9671 2 242 17.8 392.83 4.03
## 4 0.03237 0 2.18 0 0.458 6.998 45.8 6.0622 3 222 18.7 394.63 2.94
## 5 0.06905 0 2.18 0 0.458 7.147 54.2 6.0622 3 222 18.7 396.90 5.33
## 6 0.02985 0 2.18 0 0.458 6.430 58.7 6.0622 3 222 18.7 394.12 5.21
## medv
## 1 24.0
## 2 21.6
## 3 34.7
## 4 33.4
## 5 36.2
## 6 28.7summary(datos)## crim zn indus chas
## Min. : 0.00632 Min. : 0.00 Min. : 0.46 Min. :0.00000
## 1st Qu.: 0.08205 1st Qu.: 0.00 1st Qu.: 5.19 1st Qu.:0.00000
## Median : 0.25651 Median : 0.00 Median : 9.69 Median :0.00000
## Mean : 3.61352 Mean : 11.36 Mean :11.14 Mean :0.06917
## 3rd Qu.: 3.67708 3rd Qu.: 12.50 3rd Qu.:18.10 3rd Qu.:0.00000
## Max. :88.97620 Max. :100.00 Max. :27.74 Max. :1.00000
## nox rm age dis
## Min. :0.3850 Min. :3.561 Min. : 2.90 Min. : 1.130
## 1st Qu.:0.4490 1st Qu.:5.886 1st Qu.: 45.02 1st Qu.: 2.100
## Median :0.5380 Median :6.208 Median : 77.50 Median : 3.207
## Mean :0.5547 Mean :6.285 Mean : 68.57 Mean : 3.795
## 3rd Qu.:0.6240 3rd Qu.:6.623 3rd Qu.: 94.08 3rd Qu.: 5.188
## Max. :0.8710 Max. :8.780 Max. :100.00 Max. :12.127
## rad tax ptratio black
## Min. : 1.000 Min. :187.0 Min. :12.60 Min. : 0.32
## 1st Qu.: 4.000 1st Qu.:279.0 1st Qu.:17.40 1st Qu.:375.38
## Median : 5.000 Median :330.0 Median :19.05 Median :391.44
## Mean : 9.549 Mean :408.2 Mean :18.46 Mean :356.67
## 3rd Qu.:24.000 3rd Qu.:666.0 3rd Qu.:20.20 3rd Qu.:396.23
## Max. :24.000 Max. :711.0 Max. :22.00 Max. :396.90
## lstat medv
## Min. : 1.73 Min. : 5.00
## 1st Qu.: 6.95 1st Qu.:17.02
## Median :11.36 Median :21.20
## Mean :12.65 Mean :22.53
## 3rd Qu.:16.95 3rd Qu.:25.00
## Max. :37.97 Max. :50.00Paso 2 - Aplicar las transformaciones de variables
En estra primera etapa aplicaremos una transformación logaritmica y una raíz cuadrática, intentando disminuir la curva de homocestacidad y la dispensión de los residuos observados en el gráfico QQ-plot de la semana pasada.
# Aplicamos la transformación logarítmica a 'lstat'
datos$log_lstat <- log(datos$lstat)
# Verificar que la nueva variable fue creada correctamente
head(datos$log_lstat)## [1] 1.605430 2.212660 1.393766 1.078410 1.673351 1.650580# Aplicamos la transformación logarítmica a 'lstat'
datos$sqrt_lstat <- sqrt(datos$lstat)
# Verificar que la nueva variable fue creada correctamente
head(datos$sqrt_lstat)## [1] 2.231591 3.023243 2.007486 1.714643 2.308679 2.282542Paso 3 - División del conjunto de datos
Procedemos a dividir el conjunto de datos en una proporsión que es muy usada en este tipo de estudios, 70 % para entrenamiento y un 30% para pruebas y validadción.
set.seed(123)
indice <- createDataPartition(datos$medv, p = 0.7, list = FALSE)
entrenamiento <- datos[indice, ]
prueba <- datos[-indice, ]Paso 4 - Generación de modelos con las variables transformadas
modelo_ajustado <- lm(medv ~ rm + lstat, data = entrenamiento)
summary(modelo_ajustado)##
## Call:
## lm(formula = medv ~ rm + lstat, data = entrenamiento)
##
## Residuals:
## Min 1Q Median 3Q Max
## -17.5175 -3.3191 -0.8718 1.9854 27.2216
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.6826 3.4577 -0.197 0.844
## rm 4.9476 0.4822 10.261 <2e-16 ***
## lstat -0.6313 0.0500 -12.626 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.355 on 353 degrees of freedom
## Multiple R-squared: 0.6494, Adjusted R-squared: 0.6474
## F-statistic: 326.9 on 2 and 353 DF, p-value: < 2.2e-16modelo_log <- lm(medv ~ rm + log_lstat, data = entrenamiento)
summary(modelo_log)##
## Call:
## lm(formula = medv ~ rm + log_lstat, data = entrenamiento)
##
## Residuals:
## Min 1Q Median 3Q Max
## -16.5637 -2.8589 -0.4519 2.4528 27.3552
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 20.6696 3.9341 5.254 2.58e-07 ***
## rm 3.7891 0.4619 8.204 4.39e-15 ***
## log_lstat -9.2892 0.5759 -16.130 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.896 on 353 degrees of freedom
## Multiple R-squared: 0.707, Adjusted R-squared: 0.7054
## F-statistic: 425.9 on 2 and 353 DF, p-value: < 2.2e-16modelo_sqrt <- lm(medv ~ rm + sqrt_lstat, data = entrenamiento)
summary(modelo_sqrt)##
## Call:
## lm(formula = medv ~ rm + sqrt_lstat, data = entrenamiento)
##
## Residuals:
## Min 1Q Median 3Q Max
## -17.0054 -3.1943 -0.5723 2.0991 27.0488
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 13.2628 3.8391 3.455 0.000618 ***
## rm 4.2786 0.4703 9.098 < 2e-16 ***
## sqrt_lstat -5.1842 0.3536 -14.660 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.087 on 353 degrees of freedom
## Multiple R-squared: 0.6837, Adjusted R-squared: 0.6819
## F-statistic: 381.5 on 2 and 353 DF, p-value: < 2.2e-16Paso 5 - Evaluación de los modelos
Evaluamos el rendimiento de ambos modelos con el conjunto de prueba utilizando el MAE(Error absoluto medio) y el MSE(Error cuadrático medio).
pred_ajustado <- predict(modelo_ajustado, prueba)
pred_log <- predict(modelo_log, prueba)
pred_sqrt <- predict(modelo_sqrt, prueba)
mae_ajustado <- mean(abs(pred_ajustado - prueba$medv))
mse_ajustado <- mean((pred_ajustado - prueba$medv)^2)
mae_log <- mean(abs(pred_log - prueba$medv))
mse_log <- mean((pred_log - prueba$medv)^2)
mae_sqrt <- mean(abs(pred_sqrt - prueba$medv))
mse_sqrt <- mean((pred_sqrt - prueba$medv)^2)
cat("MAE Log:", mae_ajustado, "MSE Log:", mse_ajustado, "\n")## MAE Log: 4.092162 MSE Log: 35.56467cat("MAE Log:", mae_log, "MSE Log:", mse_log, "\n")## MAE Log: 3.715681 MSE Log: 27.12551cat("MAE Sqrt:", mae_sqrt, "MSE Sqrt:", mse_sqrt, "\n")## MAE Sqrt: 3.865623 MSE Sqrt: 31.00546Dado que el que tiene un menor MAE y MSE es el modelo con la transformación logaritmica, nos quearemos con este modelo, ahora como siguiente paso hay que validar los supuestos de linealidad.
Paso 6 - Validando la linealidad del modelo
ggplot(entrenamiento, aes(x = log_lstat, y = medv)) +
geom_point() + geom_smooth(method = "lm") +
labs(title = "Relación lineal entre lstat y medv")## `geom_smooth()` using formula = 'y ~ x'
#Prueba de independencia de errores
durbinWatsonTest(modelo_log)## lag Autocorrelation D-W Statistic p-value
## 1 0.510388 0.9556693 0
## Alternative hypothesis: rho != 0# Gráfico de homocedasticidad
plot(modelo_log, which = 3)
# Gráfico QQ-plot para normalidad de residuos
qqPlot(modelo_log)
## 365 373
## 254 259# Calcular VIF para detectar multicolinealidad
vif(modelo_log)## rm log_lstat
## 1.710332 1.710332Con las pruebas obtenidas podemos verificar que el modelo si mejoró en función de nuestro modelo de la semana pasada, sin embargo aún no cumple con locriterios de linealidad, en este punto tenemos varias opciones, podemos agregar más variables, ocupar métodos más robusto de modelado o intentar un modelado no lineal. En nuestro caso y porque así lo marca el temario pasaremos directo a usar un modelo polinomial, ya que indagar en técnicas más especificas de regresión lineal pudiese ser objeto de estudio de un curso completo.
Paso 7 - APlicando un modelo polinomial
Para este ejemplo aplicaremos una regresión polinomial de grado 2, si observamos mejora podemos aplicar una regresión de un grado mayor, pero hay que tener cuidado con esto para evitar el sobreajuste.
modelo_polinomial <- lm(medv ~ rm + poly(lstat, 2), data = entrenamiento)
summary(modelo_polinomial)##
## Call:
## lm(formula = medv ~ rm + poly(lstat, 2), data = entrenamiento)
##
## Residuals:
## Min 1Q Median 3Q Max
## -18.0139 -3.0003 -0.3518 2.4446 27.2384
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.8653 2.8165 -1.727 0.085 .
## rm 4.3445 0.4462 9.737 < 2e-16 ***
## poly(lstat, 2)1 -89.3947 6.1338 -14.574 < 2e-16 ***
## poly(lstat, 2)2 41.7906 4.9557 8.433 8.81e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.892 on 352 degrees of freedom
## Multiple R-squared: 0.7083, Adjusted R-squared: 0.7058
## F-statistic: 284.9 on 3 and 352 DF, p-value: < 2.2e-16Ya que tenemos el modelo tenemos que proceder a evaluarlo y a realizar predicciones para saber si mejoró o no con respecto a nuestros intentos previos.
# Realizar predicciones con el conjunto de prueba
predicciones_polinomiales <- predict(modelo_polinomial, newdata = prueba)
# Calcular el MAE y MSE para evaluar el modelo
mae_polinomiales <- mean(abs(predicciones_polinomiales - prueba$medv))
mse_polinomiales <- mean((predicciones_polinomiales - prueba$medv)^2)
cat("MAE del modelo polinomial:", mae_polinomiales, "\n")## MAE del modelo polinomial: 3.789103cat("MSE del modelo polinomial:", mse_polinomiales, "\n")## MSE del modelo polinomial: 28.59745Finalmente procedemos a evaluar de nueva cuenta la linealidad.
#Prueba de independencia de errores
durbinWatsonTest(modelo_polinomial)## lag Autocorrelation D-W Statistic p-value
## 1 0.4975474 0.9802912 0
## Alternative hypothesis: rho != 0# Gráfico de homocedasticidad
plot(modelo_polinomial, which = 3)
# Gráfico QQ-plot para normalidad de residuos
qqPlot(modelo_polinomial)
## 365 373
## 254 259# Calcular VIF para detectar multicolinealidad
vif(modelo_polinomial)## GVIF Df GVIF^(1/(2*Df))
## rm 1.59874 1 1.264413
## poly(lstat, 2) 1.59874 2 1.124461Conclusiones
Después de hacer varios intentos seguimos teniendo los mismos inconvenientes que con el modelo optimizado, lo que nos indica que una regresión lineal no es una buena técnica para caracterizar y predecir la variable mdev. La siguiente semana aplicaremos la técnica de redes neuronales para obtener un mejor modelo de predicción sobre esta base de datos.
Finalmente, para concluir, imaginemos que no tuvisemos esos problemas con la linealidad, ¿cómo hubiesemos obtenido la función?
# Coeficientes del modelo logarítmico
summary(modelo_log)$coefficients## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 20.669646 3.9341277 5.253934 2.580952e-07
## rm 3.789092 0.4618666 8.203867 4.392251e-15
## log_lstat -9.289206 0.5758907 -16.130155 3.054764e-44Para el caso del modelo de logaritmo la función de acuerdo a sus coeficiebtes es: \[ \text{medv} = 20.6696 + 3.7891 \times \text{rm} - 9.2892 \times \log(\text{lstat}) \]
# Coeficientes del modelo polinomial
summary(modelo_polinomial)$coefficients## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.865260 2.8165457 -1.727385 8.497564e-02
## rm 4.344542 0.4461758 9.737286 5.365542e-20
## poly(lstat, 2)1 -89.394709 6.1338427 -14.574014 5.629972e-38
## poly(lstat, 2)2 41.790609 4.9557077 8.432824 8.810803e-16Y para el caso del modelo polinomial la función de acuerdo a sus coeficientes es:
\[ \text{medv} = -4.8653 + 4.3445 \times \text{rm} - 89.3947 \times \text{lstat} + 41.7906 \times \text{lstat}^2 \]
Finalmente no olviden realizar la actividad de la semana, la cual consiste en usar la técnica de regresión lineal a partir de una base de datos de su elección o emplear la proporcionada en la actividad en formato .csv para intentar predecir una variable.