Paso 1 - Carga de bibliotecas y conjunto de datos

En caso de que alguna biblioteca te marque error al intentar cargar, deberás de instalarla y posteriormente volver a ejecutar esta sección.

#cargar bibliotecas
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

data("Boston") # Carga el conjunto de datos Boston

Paso 2 - Exploración de datos

Imprimos la información del conjunto de datos de Boston, el cual contiene datos del valor medio de viviendas y diferentes características del área.

# Visualizamos la estructura y el resumen estadístico de los datos
str(Boston)
## 'data.frame':    506 obs. of  14 variables:
##  $ crim   : num  0.00632 0.02731 0.02729 0.03237 0.06905 ...
##  $ zn     : num  18 0 0 0 0 0 12.5 12.5 12.5 12.5 ...
##  $ indus  : num  2.31 7.07 7.07 2.18 2.18 2.18 7.87 7.87 7.87 7.87 ...
##  $ chas   : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ nox    : num  0.538 0.469 0.469 0.458 0.458 0.458 0.524 0.524 0.524 0.524 ...
##  $ rm     : num  6.58 6.42 7.18 7 7.15 ...
##  $ age    : num  65.2 78.9 61.1 45.8 54.2 58.7 66.6 96.1 100 85.9 ...
##  $ dis    : num  4.09 4.97 4.97 6.06 6.06 ...
##  $ rad    : int  1 2 2 3 3 3 5 5 5 5 ...
##  $ tax    : num  296 242 242 222 222 222 311 311 311 311 ...
##  $ ptratio: num  15.3 17.8 17.8 18.7 18.7 18.7 15.2 15.2 15.2 15.2 ...
##  $ black  : num  397 397 393 395 397 ...
##  $ lstat  : num  4.98 9.14 4.03 2.94 5.33 ...
##  $ medv   : num  24 21.6 34.7 33.4 36.2 28.7 22.9 27.1 16.5 18.9 ...
summary(Boston)
##       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.00

Después de observar la tabla, ¿tu qué variable seleccionarías como dependiente y cuáles como independientes?

Para nuestro caso, elegiremos como variable dependiente(a predecir), el valor medio de la vivienda, el cual esta etiquetado como medv. Mientras que las variables independientes serán: crim,rm, age, tax, lstat. Si no recuerdas el significado de las etiquetas puedes consultarla aa descripcíón de la base de datos en: http://www.cs.toronto.edu/~delve/data/boston/bostonDetail.html

Paso 3 - Selección de Variables Independientes

Ahora procederemos a construir el modelo con las variables seleccionadas.

# Modelo de regresión múltiple inicial, dado que estas variables son las que decidí contemplar.
modelo_completo <- lm(medv ~ crim + rm + age + tax + lstat, data = Boston)
summary(modelo_completo)
## 
## Call:
## lm(formula = medv ~ crim + rm + age + tax + lstat, data = Boston)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -16.411  -3.415  -1.037   1.915  29.613 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.858471   3.185415  -0.270  0.78766    
## crim        -0.060005   0.035489  -1.691  0.09150 .  
## rm           5.087041   0.447418  11.370  < 2e-16 ***
## age          0.020918   0.011415   1.833  0.06746 .  
## tax         -0.005904   0.001978  -2.985  0.00297 ** 
## lstat       -0.583764   0.056806 -10.276  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.445 on 500 degrees of freedom
## Multiple R-squared:  0.6529, Adjusted R-squared:  0.6495 
## F-statistic: 188.1 on 5 and 500 DF,  p-value: < 2.2e-16

Al hacer esto el software nos arrija un resultado de acuerdo a las variables que seleccionamos, ahora emplearemos el método paso a paso, para determinar cuales son las variables más significativas y poder descartar las que no lo son.

# Selección de variables usando Stepwise con AIC
modelo_stepwise <- step(modelo_completo, direction = "both")
## Start:  AIC=1721.04
## medv ~ crim + rm + age + tax + lstat
## 
##         Df Sum of Sq   RSS    AIC
## <none>               14825 1721.0
## - crim   1      84.8 14910 1721.9
## - age    1      99.6 14925 1722.4
## - tax    1     264.3 15090 1728.0
## - lstat  1    3131.3 17957 1816.0
## - rm     1    3833.0 18658 1835.4
summary(modelo_stepwise)
## 
## Call:
## lm(formula = medv ~ crim + rm + age + tax + lstat, data = Boston)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -16.411  -3.415  -1.037   1.915  29.613 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.858471   3.185415  -0.270  0.78766    
## crim        -0.060005   0.035489  -1.691  0.09150 .  
## rm           5.087041   0.447418  11.370  < 2e-16 ***
## age          0.020918   0.011415   1.833  0.06746 .  
## tax         -0.005904   0.001978  -2.985  0.00297 ** 
## lstat       -0.583764   0.056806 -10.276  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.445 on 500 degrees of freedom
## Multiple R-squared:  0.6529, Adjusted R-squared:  0.6495 
## F-statistic: 188.1 on 5 and 500 DF,  p-value: < 2.2e-16

Analiza los resultados, especificamente observa la columan de P, el cual indica la significancia de cada variable. Si el valor p es menor a 0.05(umbral típico), la variable se considera significativa, de lo contrario, podría ser eliminada porque no contibuye de manera importante. ¿Con qué variables debemos quedarnos? Si, es correcto, solo con rm y lstat.

Paso 4 - Optimización del modelo

Como primer paso dividamos el conjunto de datos en 70% para entrenamiento y un 30% para pruebas y validar su efectividad.

# Dividimos los datos en conjuntos de entrenamiento y prueba
set.seed(123)  # Fijar semilla para que tengamos los mismos resultados
indice <- createDataPartition(Boston$medv, p = 0.7, list = FALSE)
entrenamiento <- Boston[indice, ]
prueba <- Boston[-indice, ]

Posteriormente procedemos a crear el modelo optimizado, ahora solo utilizando las variables rm y lstat.

# Ajustamos el modelo con las variables seleccionadas
modelo_optimizacion <- lm(medv ~ rm + lstat, data = entrenamiento)
summary(modelo_optimizacion)
## 
## 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-16

Ahora es momento de medir el error cuadrático medio del conjunto de prueba, sección de los datos que estamos empleando para hacer las predicciones, dado que ya se tienen los valores, nos permite comparar la diferencia entre la predicción y el dato real.

# Realizar predicciones y calcular el MSE
predicciones <- predict(modelo_optimizacion, prueba)
mse <- mean((predicciones - prueba$medv)^2)
cat("Error Cuadrático Medio (MSE):", mse)
## Error Cuadrático Medio (MSE): 35.56467

Paso 5 - Validación de supuestos

1- Linealidad: Observa la relación lineal entre las variables predictoras seleccionadas y la variable dependiente.Al dibujar la gráfica se debe de observar una línea recta o una curva suave que se asemeje a una recta, de lo contrario esto no se cumpliría y tendríamos que revisar los parámetros del modelo empleado. En esta práctica se puede observar perfectamente la recta.

# Visualización de la relación lineal
ggplot(entrenamiento, aes(x = lstat, y = medv)) +
  geom_point() + geom_smooth(method = "lm") +
  labs(title = "Relación lineal entre lstat y medv")
## `geom_smooth()` using formula = 'y ~ x'

2- Independencia de los errores: Para esto aplicaremos la prueba Durbin-Watson para evaluar la autocorrelación de los errores. En esta salida esperamos números cercanos a 2, lo que nos indica independencia, valores muy bajos(cercanos a cero) o muy altos(cercanos a 4) nos indican correlación y tendríamos que revisar nuestor modelo. En nuestro caso nos arroja valores de 0.5 y 0.9, valores por debajo de 2 y más cercanos 0, con lo que concluimos la dependencia, lo que nos sugiere tener que ajustar nuestro modelo.

# Prueba de independencia de errores
durbinWatsonTest(modelo_optimizacion)
##  lag Autocorrelation D-W Statistic p-value
##    1       0.5162777     0.9503729       0
##  Alternative hypothesis: rho != 0

3- Homocedasticidad: Ahora verificaremos si los residuos tienen varianza constante a lo largo de los valores que se predijeron. Esta gráfica nos dibuja una líena roja indicando la tendencia, esta deber ser cercana a una horizontal, en nuestro caso se forma una U, lo que nos muestra que la varianza de los residuos no es constante a lo largo de los valores ajustados, concluyendo que no hay Homocedasticidad, con lo que se recomienda hacer una transformación a la variable dependiente mdev, esto lo estaremos viendo la próxima semana.

# Gráfico de homocedasticidad
plot(modelo_optimizacion, which = 3)

4- Normalidad de los residuos: Para esto realizaremos un gráfico QQ-plot para observar la distribución de los errores. En esta gráfica se dibuja una línea diagonal, todos los puntos deben estar cercanos a esta línea, especialmente hay que cuidar los extremos, ya que si hay separación de los puntos esto puede afectar la confiabilidad de los intervalos de confianza y pruebas de hipótesis. ¿Qué sucede en nuestro caso?

Si los valores extremos se nos separan, con lo que debemos ajustar nuestro modelo, una transformación de variables nos ayudaría, pero este tema lo veremos la próxima semana, donde retomaremos este ejercicio.

# Gráfico QQ-plot para normalidad de residuos
qqPlot(modelo_optimizacion)

## 370 373 
## 258 259

5- Multicolinealidad: Calculamos el factor de inflación de la varianza.

# Calcular VIF para detectar multicolinealidad
vif(modelo_optimizacion)
##       rm    lstat 
## 1.557671 1.557671

Observa los valores obtenidos, números superiores a 10 nos indican un problema de multicolinealidad y sugieren eliminar o transforma variables. Para este caso, nuestro valores son de 1.5, al estar por debajo de 10.

Paso 6 Conclusión

Aunque ya tenemos un modelo que puede realizar predicciones es muy arriesgado hacer uso del mismo así como está, ya que hay supuestos que no se han cumplido, es necesario hacer trasnformación de variables para ajustar el modelo y cumplir con estos supuestos. Tema que estaremos retomando la próxima semana para no hacer tan larga esta práctica, la cual retomaremos la próxima semana.

¿Cuál es nuestra función obtenida?

Para esto necesitamos consultar los coeficientes con la siguiente instrucción:

coef(modelo_optimizacion)
## (Intercept)          rm       lstat 
##  -0.6826325   4.9475826  -0.6313110

Con lo anterior podemos observar los valores, con lo que ya podemos escribir dicha ecuación, de acuerdo a nuestros resultados es: \[ \text{medv} = -0.6826 + 4.9476(\text{rm}) - 0.6313(\text{lstat}) \]

medv = -0.6826 + 4.9476(rm) - 0.6313(lstat)

Paso 7 - Uso del modelo esto para propósito de comparación

Realizamos predicciones con el conjunto de datos para este propósito y calculamos el error absoluto medio y el error cuadrático medio del modelo sin ajustes y trasnformación de varibles.

# Predicciones con el modelo ajustado
predicciones <- predict(modelo_optimizacion, Boston)
mae <- mean(abs(predicciones - Boston$medv))
mse <- mean((predicciones - Boston$medv)^2)
cat("Error absoluto medio (MAE):", mae, "\nError cuadrático medio (MSE):", mse, "\n")
## Error absoluto medio (MAE): 3.942629 
## Error cuadrático medio (MSE): 30.55139

El error absoluto medio mide el error promedio de las predicciones con respecto a los valores observados, sin considerar la dirección del error. Un valor bajo indica que en promedio las predicciones están cerca de los valores reales.

Por su parte, el error cuadrático medio penaliza los errores grandes al elevar al cuadrado la diferencia entre los valores predichos y reales, con lo que en este punto buscamos el menor número posible, ya que nos indicará un buen ajuste del modelo, si es grande puede indicar sobreajuste o errores significativos en las predicciones.

Para nuestro ejercicio, dado el número de error absoluto nos indica que tenemos un buen desempeño en general, sin embargo, es posible seguir optimizando para reducir el error cuadrático medio.