Analisis De Regresion Multiple
Samuel Coronado
2024-08-29
#install.packages("devtools") solo la primera vez
#devtools::install_github("centromagis/paqueteMODELOS", force =TRUE)
library(paqueteMODELOS)
library(dplyr)
library(plotly)
library(naniar)
library(mice)
library(plotly)
library(corrplot)
library(car)
library(caret)
library(Metrics)
library(tseries)
library(lmtest)
1 Punto 1
data("vivienda")
str(vivienda)## spc_tbl_ [8,322 × 13] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
## $ id : num [1:8322] 1147 1169 1350 5992 1212 ...
## $ zona : chr [1:8322] "Zona Oriente" "Zona Oriente" "Zona Oriente" "Zona Sur" ...
## $ piso : chr [1:8322] NA NA NA "02" ...
## $ estrato : num [1:8322] 3 3 3 4 5 5 4 5 5 5 ...
## $ preciom : num [1:8322] 250 320 350 400 260 240 220 310 320 780 ...
## $ areaconst : num [1:8322] 70 120 220 280 90 87 52 137 150 380 ...
## $ parqueaderos: num [1:8322] 1 1 2 3 1 1 2 2 2 2 ...
## $ banios : num [1:8322] 3 2 2 5 2 3 2 3 4 3 ...
## $ habitaciones: num [1:8322] 6 3 4 3 3 3 3 4 6 3 ...
## $ tipo : chr [1:8322] "Casa" "Casa" "Casa" "Casa" ...
## $ barrio : chr [1:8322] "20 de julio" "20 de julio" "20 de julio" "3 de julio" ...
## $ longitud : num [1:8322] -76.5 -76.5 -76.5 -76.5 -76.5 ...
## $ latitud : num [1:8322] 3.43 3.43 3.44 3.44 3.46 ...
## - attr(*, "spec")=List of 3
## ..$ cols :List of 13
## .. ..$ id : list()
## .. .. ..- attr(*, "class")= chr [1:2] "collector_double" "collector"
## .. ..$ zona : list()
## .. .. ..- attr(*, "class")= chr [1:2] "collector_character" "collector"
## .. ..$ piso : list()
## .. .. ..- attr(*, "class")= chr [1:2] "collector_character" "collector"
## .. ..$ estrato : list()
## .. .. ..- attr(*, "class")= chr [1:2] "collector_double" "collector"
## .. ..$ preciom : list()
## .. .. ..- attr(*, "class")= chr [1:2] "collector_double" "collector"
## .. ..$ areaconst : list()
## .. .. ..- attr(*, "class")= chr [1:2] "collector_double" "collector"
## .. ..$ parqueaderos: list()
## .. .. ..- attr(*, "class")= chr [1:2] "collector_double" "collector"
## .. ..$ banios : list()
## .. .. ..- attr(*, "class")= chr [1:2] "collector_double" "collector"
## .. ..$ habitaciones: list()
## .. .. ..- attr(*, "class")= chr [1:2] "collector_double" "collector"
## .. ..$ tipo : list()
## .. .. ..- attr(*, "class")= chr [1:2] "collector_character" "collector"
## .. ..$ barrio : list()
## .. .. ..- attr(*, "class")= chr [1:2] "collector_character" "collector"
## .. ..$ longitud : list()
## .. .. ..- attr(*, "class")= chr [1:2] "collector_double" "collector"
## .. ..$ latitud : list()
## .. .. ..- attr(*, "class")= chr [1:2] "collector_double" "collector"
## ..$ default: list()
## .. ..- attr(*, "class")= chr [1:2] "collector_guess" "collector"
## ..$ delim : chr ";"
## ..- attr(*, "class")= chr "col_spec"
## - attr(*, "problems")=<externalptr>summarytools::freq(vivienda$tipo)## Frequencies
## vivienda$tipo
## Type: Character
##
## Freq % Valid % Valid Cum. % Total % Total Cum.
## ----------------- ------ --------- -------------- --------- --------------
## Apartamento 5100 61.305 61.305 61.283 61.283
## Casa 3219 38.695 100.000 38.681 99.964
## <NA> 3 0.036 100.000
## Total 8322 100.000 100.000 100.000 100.000Antes de aplicar el filtro y quedarnos solo con aquellas viviendas que sean solo apartamentos, primero verificamos que no haya problemas de digitalizacion en la base de datos, para ello se hizo una tabulacion de la variable tipo y se encontraron que hay 5100 apartamentos y el nombre se encuentra escrito uniformemente.
vivienda1 <- vivienda %>% filter(tipo == "Apartamento")
head(vivienda1,3)## # A tibble: 3 × 13
## id zona piso estrato preciom areaconst parqueaderos banios habitaciones
## <dbl> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1212 Zona N… 01 5 260 90 1 2 3
## 2 1724 Zona N… 01 5 240 87 1 3 3
## 3 2326 Zona N… 01 4 220 52 2 2 3
## # ℹ 4 more variables: tipo <chr>, barrio <chr>, longitud <dbl>, latitud <dbl>A continuación, se revisara si la base de datos tiene valores vacios, los cuales pueden afectar posteriormente el analisis.
gg_miss_var(vivienda1) # grafico de datos faltantes
Debido a la gran cantidad de datos faltantes que contienen las variables Piso y Parqueadero, no es conveniente eliminar todas la filas que contengas valores vacios, ya que se perderia mucha informacion importante para nuestro analisis, por lo tanto se procedera a reemplazar esa informacion por la mediana correspondiente a cada variable.
# Rellenar los valores faltantes con la mediana de cada columna
vivienda1 <- vivienda1 %>%
mutate(across(everything(), ~ ifelse(is.na(.), median(., na.rm = TRUE), .)))
gg_miss_var(vivienda1) # grafico de datos faltantes
Al haber realizado dicho procedimiento, como se puede verificar en el grafico anterior, ninguna variable cuenta con valores faltantes.
2 Punto 2
plot_ly(data = vivienda1, x = ~areaconst, y = ~preciom, type = 'scatter', mode = 'markers') %>%
layout(title = "Scatter Plot: Precio vs Area Construida",
xaxis = list(title = "Area Construida"),
yaxis = list(title = "Precio de los apartamentos"))Como se puede apreciar en el grafico anterior, se evidencia que hay una tendencia positiva entre el tamaño que tenga el apartamento, representado por los metros construidos y el precio que estos logran alcanzar, ya que entre mayor sea el area construida, mayor sera el precio al que se logran vender, mas adelante a traves de la matriz de correlacion obtendremos un valor numerico del grado de correlacion que tienen estas 2 variables.
plot_ly(data = vivienda1, x = ~as.factor(estrato), y = ~preciom, type = 'box', color = ~as.factor(estrato)) %>%
layout(title = "Boxplot: Precio de la Vivienda por Estrato",
xaxis = list(title = "Estrato"),
yaxis = list(title = "Precio de la Vivienda"))Al hacer el analisis bivariado entre el precio de los apartamentos segun el estrato, se logra observar que entre mas alto es el estrato en donde se encuentre el aparatmento, mayor precio de venta estos alcanzan.
plot_ly(data = vivienda1, x = ~as.factor(banios), y = ~preciom, type = 'violin', color = ~as.factor(banios),
box = list(visible = TRUE), # Incluye el boxplot dentro del violin plot
meanline = list(visible = TRUE)) %>%
layout(title = "Violin Plot: Precio de los Aparatmentos por Numero de Baños",
xaxis = list(title = "Numero de Baños"),
yaxis = list(title = "Precio de los apartamentos"))Del anterior grafico de violinplot se puede observar que hay una relacion directamente proporcional, entre el precio de los apartamentos y el numero de baños, ya que entre mas numeros de baños estos tengan, mas alto sera su precio de venta.
plot_ly(data = vivienda1, x = ~as.factor(parqueaderos), y = ~preciom, type = 'violin', color = ~as.factor(parqueaderos),
box = list(visible = TRUE), # Incluye el boxplot dentro del violin plot
meanline = list(visible = TRUE)) %>%
layout(title = "Violin Plot: Precio de los apartamentos vs Numero de Parqueaderos",
xaxis = list(title = "Numero de Parqueaderos"),
yaxis = list(title = "Precio de los apartamentos"))Al hacer el analisis entre el precio de los apartamentos vs el numero de parqueaderos, al iagual que ocurrio con el numero de baños, se logra ver una tendencia positiva entre estas 2 variables, sin embargo a partir de 10 parqueaderos, el precio de los apartamentos cae bruscamente, lo que amerita encontrar posibles causas de este cambio en los precios.
plot_ly(data = vivienda1, x = ~as.factor(habitaciones), y = ~preciom, type = 'violin', color = ~as.factor(habitaciones),
box = list(visible = TRUE), # Incluye el boxplot dentro del violin plot
meanline = list(visible = TRUE)) %>%
layout(title = "Violin Plot: Precio de los apartamentos vs Numero de Habitaciones",
xaxis = list(title = "Numero de Habitaciones"),
yaxis = list(title = "Precio de los apartamentos"))Al ver el siguiente grafico, de la relacion de precios con el numero de habitaciones, no se logra observar una clara tendencia, puesto que el precio de los apartamentos va fluctuando en diferentes tramos, a medida que aumenta el numero de habitaciones.
plot_ly(data = vivienda1, x = ~as.factor(zona), y = ~preciom, type = 'box', color = ~as.factor(zona)) %>%
layout(title = "Boxplot: Precio de los apartamentos vs la zona",
xaxis = list(title = "Zona de la ciudad"),
yaxis = list(title = "Precio de los apartamentos"))Al analizar los precios vs la zona en donde se encuentra ubicado el apartamentos, se logra ver claramente que la zona oeste de la ciudad de cali es el lugar donde a mejor precio de venta alcanzan los apartamentos seguido de la zona sur, en su contra parte la zona oriental y la zona centro son las zonas donde los apartamentos se venden mas barato.
3 Punto 3
Para realizar la estimacion del modelo, primero realizaremos una matriz de correlacion, para ver el grado de correlacion que tienen las variables independientes, con nuestra variable dependiente precio.
# Calcular la matriz de correlación usando Spearman
cor_matrix <- cor(vivienda1[, c(5, 4, 6, 7, 8, 9)], method = "spearman", use = "complete.obs")
# Visualizar la matriz de correlación con colores
corrplot(cor_matrix, method = "color", col = colorRampPalette(c("blue", "white", "red"))(200),
type = "lower", tl.col = "black", tl.srt = 45, addCoef.col = "black", number.cex = 0.7)
Como se puede observar en la matriz de correlacion, realizada por el metodo de spearman, la mayoria de las variables independientes estan altamente correlacionadas con la variable dependiente precio, a excepcion de la variable habitaciones, cuya correlacion es baja sin embargo, se incluira dentro del modelo para evaluar el efecto que esta tiene sobre el precio de los apartamentos y ver si es significativa dentro del modelo.
modelo<-lm(preciom ~ estrato + areaconst + parqueaderos + banios + habitaciones, data=vivienda1)
summary(modelo)##
## Call:
## lm(formula = preciom ~ estrato + areaconst + parqueaderos + banios +
## habitaciones, data = vivienda1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1706.19 -50.41 0.97 43.66 1030.31
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -269.2153 12.6949 -21.21 <2e-16 ***
## estrato 53.6219 2.5317 21.18 <2e-16 ***
## areaconst 2.0141 0.0423 47.61 <2e-16 ***
## parqueaderos 92.8962 3.6902 25.17 <2e-16 ***
## banios 48.1824 2.9578 16.29 <2e-16 ***
## habitaciones -36.0841 3.2135 -11.23 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 130.8 on 5094 degrees of freedom
## Multiple R-squared: 0.7958, Adjusted R-squared: 0.7956
## F-statistic: 3971 on 5 and 5094 DF, p-value: < 2.2e-16Analizando la siguiente salida podemos decir que todas las variables son significativas a nivel individual ya que si nos fijamos en el p-valor que es 2e-16 esta por debajo del p critico del 5% o 0.05 lo que me quiere decir que se rechaza la hipotesis nula de no significancia, por tanto todas son significativas, si analizamos el R2 este arroja un valor de 0.78 lo que me indica que las variables independientes explican en un 78% las variaciones que se presentan en el precio de los apartamentos, este es un muy buen indicador de ajuste del modelo.
3.0.1 Interpretaciones de los coeficientes
Acontinuacion se presenta la ecuacion de regresion
\[ Precio = B_1 + B_2 Estrato + B_3Areaconst + B_4 Parqueaderos + B_5 Baños + B_6 Habitacion \]
La variable estrato el cual es el \(B_2\) arrojo un valor de 56.24 lo que me indica que a medica que el estrato aumenta en una unidad, el precio de los apartamentos aumenta en promedio 56 millones de pesos.
Area Construida tuvo un valor en su coeficiente \(B_3\) de 2.004 lo que me indica que por cada metro construido adicional que tenga el apartamento su valor aumentara en promedio 2 millones de pesos.
La variable parqueadero que es el \(B_4\) obtuvo un valor de 90.42 lo que indica que a medida que aumenten el numero de parqueaderos en una unidad, el precio de los apartamentos aumentara en promedio 90 millones pesos.
La Variable baños que es el \(B_5\) obtuvo un valor de 54.84 lo que indica que por cada baño adicional que tenga el apartamento su precio aumentara en promedio 54 millones de pesos.
Por ultimo la variable habitacion que es el \(B_6\) obtuvo un valor de -42.66 lo que me indica que por cada habitacion adicional que tenga el apartamento su precio disminuye en promedio 42 millones de pesos.
Despues del analisis anterior, la unica variable que no se ajusta a los resultados esperados de manera aprioris, es la variable habitacion, ya que entre mas habitaciones, quiere decir que el apartamento es mas grande por tanto deberia tomar un mayor valor, sin embargo segun los resultados la relacion es inversa, por lo que se propone que esta variable no se incluya en el modelo ya que esta generando inconsistencias.
4 Punto 4
4.0.1 Supuesto de Normalidad
# Histograma de los Residuos
hist(modelo$residuals,col="cadetblue2",main = "Histograma de los Residuos")
qqnorm(modelo$residuals,col="cadetblue2")
qqline(modelo$residuals,col="blue")
segun el Histograma y el Grafico Quantil Quantil parece que los residuos siguen una distribucion normal, sin embargo se aplicara la prueba Shapiro Wilk para salir de dudas.
jarque.bera.test(modelo$residuals)##
## Jarque Bera Test
##
## data: modelo$residuals
## X-squared = 57604, df = 2, p-value < 2.2e-16Segun el test de Jarque Bera este rechaza la hipotesis nula de normalidad en los residuos, por lo tanto se puede decir que el modelo no cumple con este supuesto.
4.0.2 Supuesto de que los errores tienen media cero
# se aplica el test de Reset
options(scipen = 999)
mean(modelo$residuals)## [1] 0.00000000000001026864Como se puede observar en los resultados, en promedio los errores son cero, por tanto se cumple con este supuesto.
4.0.3 Supuesto de Varianza Constante en los errores
Para ello se aplicara la prueba de Breusch Pagan para evaluar si se cumple con el supuesto de Homocedasticidad en los errores.
bptest(modelo)##
## studentized Breusch-Pagan test
##
## data: modelo
## BP = 1424.4, df = 5, p-value < 0.00000000000000022Como se puede observar en los resultados del p-valor 0,00000000000000022 esta por debajo del p-critico de 0,05 lo que indica que rechaza la hipotesis nula de los que los residuos son homocedasticos, por tanto el modelo no cumple con este supuesto y tiene problemas de heterocedasticidad.
4.0.4 Supuesto de Independencia de los errores o no autocorrelacion
plot(modelo)
Tal como se puede evidenciar en la imagen de los residuos, se logra ver un patron lo que indica que los residuos pueden estar correlacionados entre si, para salir de dudas se aplicara la prueba de Durbin Watson.
dwtest(modelo)##
## Durbin-Watson test
##
## data: modelo
## DW = 1.6329, p-value < 0.00000000000000022
## alternative hypothesis: true autocorrelation is greater than 0Segun el resultado de la prueba de durbin watson este rechaza la hipotesis nula de no autocorrelacion serial de los errores, por tanto el modelo tiene problemas de autocorrelacion.
4.0.5 Validacion de Multicolinealidad
Para ello se aplicara la prueba conocida como VIF o factor de inflacion de varianza.
vif(modelo)## estrato areaconst parqueaderos banios habitaciones
## 1.826368 2.567706 2.042818 2.979718 1.407247Segun los resultados de la Prueba de Factor de inflacion de varianza, ninguna de las variables supera el valor de 10 o incluso el de 5, el cual es el valor critico como referencia para saber si hay problemas de multicolinealidad, por tanto se puede decir que el modelo no tiene este problema.
5 Punto 5
# Configurar una semilla para reproducibilidad
set.seed(123)
# Realizar la partición (70% entrenamiento, 30% prueba)
trainIndex <- createDataPartition(vivienda1$preciom, p = 0.7, list = FALSE)
# Dividir los datos en entrenamiento y prueba
train_data <- vivienda1[trainIndex, ]
test_data <- vivienda1[-trainIndex, ]
# Estimar el modelo usando solo el conjunto de entrenamiento
modelo_train <- lm(preciom ~ estrato + areaconst + parqueaderos + banios + habitaciones, data=train_data)
# Mostrar los resultados del modelo entrenado
summary(modelo_train)##
## Call:
## lm(formula = preciom ~ estrato + areaconst + parqueaderos + banios +
## habitaciones, data = train_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1670.51 -50.45 0.20 43.45 1043.15
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -281.00776 15.14240 -18.558 <0.0000000000000002 ***
## estrato 57.54070 3.04266 18.911 <0.0000000000000002 ***
## areaconst 1.97309 0.05047 39.091 <0.0000000000000002 ***
## parqueaderos 93.22171 4.50302 20.702 <0.0000000000000002 ***
## banios 43.55958 3.55251 12.262 <0.0000000000000002 ***
## habitaciones -33.33537 3.81823 -8.731 <0.0000000000000002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 131.6 on 3566 degrees of freedom
## Multiple R-squared: 0.7861, Adjusted R-squared: 0.7858
## F-statistic: 2621 on 5 and 3566 DF, p-value: < 0.00000000000000022Al haber hecho la estimacion con el 70% de los datos entrenados, se analiza el R2 que da 0.77 o 77% lo que me indica que el modelo tiene un buen ajuste y ademas se obtuvieron resultados muy similares, a los que arrojo el modelo hecho con todos los datos, sin hacer particion.
6 Punto 6
# Realizar predicciones con el modelo usando el conjunto de prueba
predicciones <- predict(modelo_train, newdata = test_data)
# Crear un dataframe con los valores reales y predichos
resultados <- data.frame(Real = test_data$preciom, Prediccion = predicciones)
# Gráfico interactivo usando plotly
plot_ly(resultados, x = ~Real, y = ~Prediccion, type = 'scatter', mode = 'markers',
marker = list(color = 'blue', size = 10)) %>%
layout(title = "Comparación de Valores Reales vs Predicciones",
xaxis = list(title = "Valores Reales (Precio)"),
yaxis = list(title = "Predicciones (Precio)"),
shapes = list(
list(type = "line",
x0 = min(resultados$Real), x1 = max(resultados$Real),
y0 = min(resultados$Real), y1 = max(resultados$Real),
line = list(color = "red", dash = "dash"))
))head(resultados,10)## Real Prediccion
## 1 240 302.24897
## 2 220 225.31222
## 3 385 383.48057
## 4 100 58.40634
## 5 170 187.29937
## 6 130 48.74805
## 7 92 113.44571
## 8 225 215.28139
## 9 370 314.12521
## 10 150 215.28139
7 Punto 7
# Eliminar las filas que tengan valores NA en alguna de las dos columnas
resultados <- resultados[complete.cases(resultados), ]
# Calcular de nuevo las métricas usando los datos limpios
mse <- mean((resultados$Real - resultados$Prediccion)^2)
mae <- mean(abs(resultados$Real - resultados$Prediccion))
r2 <- cor(resultados$Real, resultados$Prediccion)^2
# Mostrar los resultados
cat("Error Cuadrático Medio (MSE):", mse, "\n")## Error Cuadrático Medio (MSE): 16669.2cat("Error Absoluto Medio (MAE):", mae, "\n")## Error Absoluto Medio (MAE): 78.52596cat("Coeficiente de Determinación (R²):", r2, "\n")## Coeficiente de Determinación (R²): 0.816364Segun los resultados de las metricas, se obtuvo un error absoluto medio (MAE) de 85.18 este valor se considera bajo, lo que quiere decir que las predicciones estan muy cerca de los valores reales, tambien si se analiza el R2 el valor que se obtuvo fue bastante alto de 0.80 es decir que el modelo tiene un buen ajuste del 80% o en otras palabras las variables independientes influyen en forma determinante en la varianza de los precios de los apartamentos.