Modelación Estadística: Actividad 3
Daniel Arturo Sanchez Ortiz
2024-04-22
Contexto
La inmobiliara A&C requiere la construcción de un modelo estadistico que ayude a estimar los precios de las viviendas ofrecidas. Para este propósito, dispone de datos recopilados de apartamentos ubicados en estrato 4, obtenidos de la plataforma Fincaraíz. Estos datos incluyen variables relevantes como el área construida en metros cuadrados, la ubicación de la zona de la vivienda y su respectivo precio.
De acuerdo a este objetivo en este informe se pretende ajustar un modelo robusto y preciso que pueda ayudar a la inmobiliaria A&C a determinar los precios de los inmuebles.
Punto 1: Análisis exporatorio
Inicialmente cargamos las librerias utilizadas en la estimación de los modelos a probar y procedemos a cargar la base de datos, la cual se encuentra en “paqueteMETODOS” con el nombre de vivienda4.
'data.frame': 1706 obs. of 5 variables:
$ zona : Factor w/ 5 levels "Zona Centro",..: 2 2 2 5 2 2 2 2 2 2 ...
$ estrato : Factor w/ 4 levels "3","4","5","6": 2 2 2 2 2 2 2 2 2 2 ...
$ preciom : num 232 272 255 258 250 ...
$ areaconst: num 52 160 108 96 82 117 75 60 84 117 ...
$ tipo : Factor w/ 2 levels "Apartamento",..: 1 2 1 1 1 2 1 1 1 1 ...Realizamos un resumen de toda la base de datos:
zona estrato preciom areaconst
Zona Centro : 8 3: 0 Min. :207.4 Min. : 40.00
Zona Norte : 288 4:1706 1st Qu.:230.7 1st Qu.: 60.00
Zona Oeste : 60 5: 0 Median :238.8 Median : 75.00
Zona Oriente: 6 6: 0 Mean :243.7 Mean : 87.63
Zona Sur :1344 3rd Qu.:251.5 3rd Qu.: 98.00
Max. :309.7 Max. :200.00
tipo
Apartamento:1363
Casa : 343
Identificamos registros que corresponden a casas por lo que nos quedamos solo con los apartamentos:
bd_aptos <- subset(bd_viviendas, tipo == "Apartamento")Con los datos filtrados procedemos a realizar el análisis exploratorio de las variables precio de vivienda y área de la vivienda.
zona estrato preciom areaconst
Zona Centro : 7 3: 0 Min. :207.4 Min. : 40.00
Zona Norte : 237 4:1363 1st Qu.:228.8 1st Qu.: 60.00
Zona Oeste : 52 5: 0 Median :236.1 Median : 70.00
Zona Oriente: 2 6: 0 Mean :237.7 Mean : 75.48
Zona Sur :1065 3rd Qu.:243.6 3rd Qu.: 84.00
Max. :305.2 Max. :200.00
tipo
Apartamento:1363
Casa : 0
En la variable “zona”, se observa que la mayoría de los apartamentos se encuentran ubicados en la zona sur de la ciudad, con un total de 1065 viviendas ofertadas, mientras que la zona oriente muestra la menor oferta, con tan solo 2 viviendas disponibles. En cuanto al precio, se destaca que el valor promedio de la oferta es de 237.7 millones de pesos, con un valor máximo de 305.2 millones de pesos y un mínimo de 207.4 millones de pesos. En relación al área construida en metros cuadrados, se evidencia una variación desde un mínimo de \(40 m^{2}\) hasta un máximo de \(200 m^{2}\). El tamaño promedio de las viviendas es de \(70 m^{2}\), y el \(50\%\) de los apartamentos tienen un área construida comprendida entre \(60\) y \(84 m^{2}\), según el primer y tercer cuartil.
Ahora analizamos gráficamente la distribución de las variables:
par(mfrow=c(2, 2))
# Histograma del precio de la vivienda
hist(bd_aptos$preciom, main = "Distribución del Precio de la Vivienda",
xlab = "Precio de la Vivienda (millones COP)", ylab = "Frecuencia",xlim = c(200, 300))
# Boxplot del precio de la vivienda
boxplot(bd_aptos$preciom, main = "Boxplot del Precio de la Vivienda",
ylab = "Precio de la Vivienda (millones COP)", horizontal = TRUE, cex.lab = 0.7)
# Histograma del área de la vivienda
hist(bd_aptos$areaconst, main = "Distribución del Área de la Vivienda",
xlab = "Área de la Vivienda (m²)", ylab = "Frecuencia",xlim = c(40, 200))
# Boxplot del área de la vivienda
boxplot(bd_aptos$areaconst, main = "Boxplot del Área de la Vivienda",
ylab = "Área de la Vivienda (m²)" , horizontal = TRUE, cex.lab = 0.7)
Se observa que la variable precio de la vivienda presenta una distribución simétrica, aunque se destacan varios posibles datos atípicos en el extremo superior. En contraste, la variable área de la vivienda muestra una distribución asimétrica sesgada hacia la derecha, con una concentración significativa de registros entre 60 y 84 metros cuadrados. Además, se identifican posibles observaciones atípicas en el extremo superior de esta variable también.
Punto 2: Análisis exploratorio bivariado de datos
Para profundizar en la relación entre la variable respuesta “Precio de la vivienda” en función de la variable predictora “Área construida” se realiza un análisis bivariado:
Inicialmente vemos un gráfico de dispersión entre el precio y el área de la vivienda
plot(bd_aptos$areaconst, bd_aptos$preciom, main = "Dispersión entre Precio y Área de la Vivienda",
xlab = "Área de la Vivienda (m²)", ylab = "Precio de la Vivienda (millones COP)")
abline(lm(preciom ~ areaconst, data = bd_aptos), col = "red")
Gráficamente, se evidencia una relación positiva entre estas dos variables, dado que cuando el área de la vivienda es mayor, también aumenta el precio de las viviendas. Ahora validaremos esta observación con la estimación de la correlación de pearson.
correlation <- cor(bd_aptos$areaconst, bd_aptos$preciom)
cat("La correlación entre el área y el precio de la vivienda es:", correlation)## La correlación entre el área y el precio de la vivienda es: 0.8463271Se evidencia una correlación positiva fuerte de \(0.8463271\) entre la variable precio y área. Esto indica que conforme aumenta el área de un apartamento, también aumenta su precio.Sin embargo, a pesar de la fuerza de esta relación, es importante destacar que podemos realizar una prueba de hipótesis para determinar si la correlación entre estas variables es estadísticamente significativa.
cor.test(bd_aptos$areaconst,bd_aptos$preciom, method = "pearson")##
## Pearson's product-moment correlation
##
## data: bd_aptos$areaconst and bd_aptos$preciom
## t = 58.616, df = 1361, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.8305531 0.8607445
## sample estimates:
## cor
## 0.8463271El valor p obtenido es menor que \(2.2e -16\), prácticamente cero. Esto sugiere que la correlación de 0.8463271 entre el área construida y el precio es estadísticamente significativa. En otras palabras, podemos rechazar la hipótesis nula de que no hay correlación entre estas dos variables.
Ahora analizamos como es el comportamiento del precio de las viviendas segun la zona en la que se encuentran:
ggplot(bd_aptos, aes(x = zona, y = preciom)) +
geom_boxplot() +
labs(x = "Zona", y = "Precio de la vivienda", title = "Boxplot del precio de la vivienda por zona")
A pesar de que el numero de registros de cada zona es muy variable, teneniendo el mayor numero de registros en la zona sur y en la zona norte, se evidencia que los apartamentos de la zona centro son los que tienen los precios mas altos, seguidos por los apartamentos en la zona oeste.
Punto 3: Estimación modelo
Se estima el modelo de regresión lineal simple entre \(precio=f(area)+ \epsilon\):
Call:
lm(formula = preciom ~ areaconst, data = bd_aptos)
Residuals:
Min 1Q Median 3Q Max
-26.5139 -5.0886 -0.0031 4.6406 24.3309
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.001e+02 6.698e-01 298.67 <2e-16 ***
areaconst 4.984e-01 8.503e-03 58.62 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 7.081 on 1361 degrees of freedom
Multiple R-squared: 0.7163, Adjusted R-squared: 0.7161
F-statistic: 3436 on 1 and 1361 DF, p-value: < 2.2e-16(Intercept) areaconst
200.0634552 0.4984164 Obtenemos la siguiente ecuación:
\(\widehat{preciom}_i = 200.0634552 + 0.4984164 \times areaconst_i\)
Donde:
- \(\widehat{\textbf{preciom}}:\) es la variable dependiente estimada, es decir, la predicción del modelo para el precio de un apartamento.
- \(\textbf{areaconst:}\) es la variable independiente, es decir, el área construida del apartamento.
- \(\textbf{200.0634552:}\) es la constante o intercepto estimado del modelo \((\beta _0)\), que representa la predicción del modelo para el precio de un apartamento cuando areaconst es cero, lo cual puede resultar poco realista. Se interpreta que si un apartamento no tuviera área este costaría 200.06 millones de pesos. Al carecer de interpretación lógica este resultado nos podria indicar que posiblemente hay otras variables que no se estan teniendo en cuenta y que pudieran estar influyendo en el precio base de la vivienda.
- \(\textbf{0.4984164:}\) es la pendiente o coeficiente estimado para areaconst \((\beta _1)\), que representa el cambio esperado en \(preciom\) por cada unidad de aumento en \(areaconst\). Se interpreta que por cada \(m^{2}\) que aumenta el apartamento en área, su precio se incrementa en aproximadamente \(\$ 498.416\) Pesos COP.
Punto 4: Significancia Coeficiente \((\beta _1)\)
Se contruye un intervalo de confianza \((95\%)\) para el coeficiente \((\beta _1)\):
# Intervalo de confianza del 95% para los coeficientes del modelo
ic <- confint(modelo, level = 0.95)
# Intervalo de confianza para β1
ic["areaconst", ]## 2.5 % 97.5 %
## 0.4817357 0.5150970Al calcular el intervalo de confianza para la pendiente del modelo de regresión lineal se identifica que este no incluye el valor de cero, significa que la pendiente o \((\beta _1)\) es significativamente diferente de cero. Lo cual nos indica que existe evidencia estadística para afirmar que hay una relación significativa entre el área construida y el precio de la vivienda.
Ahora comparamos este resultado frente a una prueba de hipótesis \(t\).
Se plantean las siguientes hipótesis:
\(H_0:\) La variable “areaconst” no tiene un efecto significativo en “preciom”, es decir, el coeficiente \((\beta _1)\) es igual a cero.
\(H_1:\) La variable “areaconst” tiene un efecto significativo en “preciom”, es decir, el coeficiente \((\beta _1)\) es distinto de cero.
resumen <- summary(modelo)
# Muestra el valor t y el valor p para β1
print(resumen$coefficients["areaconst", c("t value", "Pr(>|t|)")])## t value Pr(>|t|)
## 58.61576 0.00000El valor p es aproximadamente cero, por lo tanto es menor al valor de significancia \(\alpha =0.05\), es decir que la evidencia sugiere que se rechaza la hipótesis nula y se puede indicar que la variable “areaconst” si tiene un efecto significativo en “preciom”, es decir, el coeficiente \((\beta _1)\) es distinto de cero.
Punto 5: Indicador de bondad \(R^{2}\)
Traemos el indicador \(R^{2}\) del summary del modelo:
# Obtenemos el valor de R cuadrado del modelo
R2 <- resumen$r.squared
# Imprimimos el valor de R cuadrado
cat("El indicador de bondad R² es:", R2)## El indicador de bondad R² es: 0.7162696El resultado obtenido significa que aproximadamente el 71.62% de la variabilidad en el precio de la vivienda puede ser explicada por el área construida según el modelo de regresión lineal utilizado.En otras palabras el modelo captura una cantidad significativa de la relación entre estas dos variables.
Punto 6: Caso de uso
A continuación vamos a estimar el precio promedio para un apartamento de 110 metros cuadrados usando el modelo de regresión anterior:
# Nuevo registro para probar el modelo
nuevo_registro <- data.frame(areaconst = 110) # Nuevo registro con el área de la vivienda
# Realizar predicciones sobre el nuevo registro utilizando el modelo
prediccion <- predict(modelo, nuevo_registro)
# Imprimir la predicción
cat("El precio estimado de la vivienda con 110 m² de area construida es: $", prediccion*1000000, "COP" )## El precio estimado de la vivienda con 110 m² de area construida es: $ 254889255 COPEl precio promedio estimado para un apartamento de 110 metros cuadrados es de $ 254’889.255 COP. Ahora validamos cual es el precio promedio de los apartamentos de 110 metros cuadrados en la base de datos en cada zona.
Inicialmente identificamos los registros
resultados <- bd_aptos %>%
filter(areaconst == 110) %>%
group_by(zona)
resultadosSe identifican 15 apartamentos en la zona sur los cuales tienen un precio promedio de 255 millones COP y 1 apartamento en la zona oeste el cual tiene un precio de 263 millones COP. Identificamos que nuestra estimación esta muy cercana al valor promedio de los apartamentos de la zona sur pero puede sugerir un valor más económico para si el apartamento esta ubicado en la zona oeste de la ciudad.
Respecto a el valor de referencia de 200 millones se considera que esta es una oferta muy atractiva ya que esta por debajo del precio de la estimación y de los mismos apartamentos con las caracteristicas similares.
Como consideraciones adicionales se propone incluir la variable zona ya que esta puede tener influencia en el valor de las viviendas.
Punto 7: Validación de supuestos
A continuación, se realiza una inspección gráfica de los supuestos del modelo:
par(mfrow=c(2, 2))
plot(modelo)
No autocorrelación de errores
En el grafico de Residual vs Fitted se identifica que los residuos pareces distribuirse aleatoriamente alrededor de la línea horizontal en 0. Esto indicaría que los residuos no muestran un patrón sistemático en función de los valores ajustados, lo cual es un indicador positivo para la validez del modelo.
Probamos mediante una prueba de hipotesis de Durbin-Watson para evaluar la autocorrelación de los residuos en el modelo de regresión.
\(H_0:\) no hay autocorrelación entre los residuos
\(H_1:\) hay autocorrelación en los residuos.
# No autocorrelación de errores
lmtest::dwtest(modelo)##
## Durbin-Watson test
##
## data: modelo
## DW = 2.0204, p-value = 0.6435
## alternative hypothesis: true autocorrelation is greater than 0El estadístico DW es cercano a 2 lo que sugiere ausencia de autocorrelación de los errores. Por otro lado, el valor p de la prueba es mayor a un nivel de significancia \(\alpha = 0,05\) por lo que no se rechaza la hipotesis nula indicando que no hay evidencia estadística de que los términos de error están autocorrelacionados positivamente.
Normalidad
Respecto a la normalidad de los residuales observamos en la gráfica Q-Q Residuals que los puntos estan distribuidos a lo largo de una línea diagonal. Esto indicaría que los cuantiles observados y los cuantiles esperados se alinean de buena manera, lo que sugiere que los residuos se ajustan bien a una distribución normal.
Probamos mediante una prueba de hipotesis de normalidad shapiro para evaluar la normalidad de los residuales.
\(H_0:\) Los residuales provienen de una población con una distribución normal.
\(H_1:\) Los residuales no provienen de una población con una distribución normal.
# Normalidad
shapiro.test(modelo$residuals)##
## Shapiro-Wilk normality test
##
## data: modelo$residuals
## W = 0.99885, p-value = 0.5419Dado que el valor p es mayor que el nivel de significancia \(\alpha = 0,05\), no hay suficiente evidencia para rechazar la hipótesis nula. Es decir, no hay suficiente evidencia para concluir que los residuales no siguen una distribución normal.
Varianza contante
En el grafico de Scale-Location se identifica que los puntos parecen distribuirse aleatoriamente alrededor de la línea horizontal. Esto sugiere que los residuales tienen una varianza constante a lo largo de los valores ajustados, lo que se conoce como homocedasticidad.
Ahora probamos mediante una prueba de hipotesis de Breusch-Pagan si la varianza es constante:
\(H_0:\) La varianza de los residuos es constante (homocedasticidad).
\(H_1:\) La varianza de los residuos no es constante (heterocedasticidad).
# Varianza contante
lmtest::bptest(modelo)##
## studentized Breusch-Pagan test
##
## data: modelo
## BP = 0.83288, df = 1, p-value = 0.3614Dado que el valor p es mayor que el nivel de significancia \(\alpha = 0,05\), no hay suficiente evidencia para rechazar la hipótesis nula. Es decir, la varianza de los residuos es constante (homocedasticidad).
Punto 8: Transformaciones
A pesar de que el $R^{2} de 71,62% sugiere que el modelo de regresión lineal explica una cantidad significativa de la variabilidad en el precio de la vivienda utilizando el área construida como predictor principal, aún queda una proporción de variabilidad sin explicar que puede ser explorada mediante la inclusión de otras variables en el modelo o algunas transformaciones. A continuación, se presentan algunas transformaciones realizadas para mejorar el ajuste del modelo:
modelo1=lm(preciom ~ areaconst, data = bd_aptos) # Lin - Lin
modelo2=lm(preciom ~ log(areaconst), data=bd_aptos) # Lin - Log
modelo3=lm(log(preciom) ~ areaconst, data=bd_aptos) # Log - Lin
modelo4=lm(log(preciom) ~ log(areaconst), data=bd_aptos) # Log - Log
W=1/bd_aptos$areaconst
modelo5=lm(preciom ~ W, data=bd_aptos) # Hiperbolica
stargazer(modelo1, modelo2, modelo3, modelo4, modelo5, type="text", df=FALSE)##
## =====================================================================================
## Dependent variable:
## -----------------------------------------------------------------
## preciom log(preciom) preciom
## (1) (2) (3) (4) (5)
## -------------------------------------------------------------------------------------
## areaconst 0.498*** 0.002***
## (0.009) (0.00004)
##
## log(areaconst) 42.878*** 0.174***
## (0.794) (0.003)
##
## W -3,182.604***
## (69.187)
##
## Constant 200.063*** 53.820*** 5.318*** 4.723*** 282.722***
## (0.670) (3.409) (0.003) (0.014) (1.005)
##
## -------------------------------------------------------------------------------------
## Observations 1,363 1,363 1,363 1,363 1,363
## R2 0.716 0.682 0.696 0.674 0.609
## Adjusted R2 0.716 0.682 0.695 0.674 0.608
## Residual Std. Error 7.081 7.496 0.030 0.031 8.317
## F Statistic 3,435.808*** 2,919.088*** 3,110.029*** 2,814.311*** 2,116.009***
## =====================================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01Punto 9: Comparación de ajustes
Se evidencia que ninguna transformación realizada supera la bondad del ajuste del primer modelo. Por lo tanto, no se realiza validación de supuestos adicionales ya que anteriormente probamos los supuestos del primer modelo.
Punto 10: Estimación varios modelos
Al identificar que las transformaciones no mejoraron el ajuste del modelo inicial, se realizan dos modelos adicionales, el primero incluyendo la variable zona y el segundo mediante minimos cuadrados ponderados, asignandole menos peso a los residuales mas grandes.
Modelo lineal multiple
Para este modelo agregamos la variable zona como variables dummy.
# Convertir la variable categórica 'zona' en variables ficticias (dummy variables)
bd_aptos2 <- cbind(bd_aptos, model.matrix(~ zona - 1, data = bd_aptos))
# Renombrar las variables ficticias para que sean más descriptivas
colnames(bd_aptos2)[6:10] <- c("ZonaCentro", "ZonaNorte", "ZonaOeste", "ZonaOriente", "ZonaSur")
# Realizar la regresión lineal
modeloMultiple <- lm(preciom ~ areaconst + ZonaSur + ZonaCentro + ZonaNorte + ZonaOeste + ZonaOriente, data = bd_aptos2)
# Ver los resultados del modelo
summary(modeloMultiple)##
## Call:
## lm(formula = preciom ~ areaconst + ZonaSur + ZonaCentro + ZonaNorte +
## ZonaOeste + ZonaOriente, data = bd_aptos2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -26.4354 -4.9696 -0.0298 4.6553 22.7673
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 197.07951 5.06208 38.933 <2e-16 ***
## areaconst 0.49788 0.00857 58.094 <2e-16 ***
## ZonaSur 2.97955 5.01269 0.594 0.552
## ZonaCentro 5.66451 5.68002 0.997 0.319
## ZonaNorte 2.83253 5.02845 0.563 0.573
## ZonaOeste 4.58753 5.10330 0.899 0.369
## ZonaOriente NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.081 on 1357 degrees of freedom
## Multiple R-squared: 0.7171, Adjusted R-squared: 0.7161
## F-statistic: 688 on 5 and 1357 DF, p-value: < 2.2e-16Se evidencia que el R cuadrado ajustado: no presenta mejora frente al primer modelo, por ende no vale la pena seguir enfatizando en este modelo, esto se puede deber al desbalance de las diferentes zonas.
Modelo minimos cuadrados ponderados
Para este modelo inicialmente calculamos los pesos que asignaremos:
# Calcular los pesos
residuos <- residuals(modelo)
Var_residuos <- var(residuos)
pesos <- 1 / residuals(modelo)^2
# Realizar la regresión ponderada por mínimos cuadrados
modeloPonderado <- lm(preciom ~ areaconst, data = bd_aptos, weights = pesos)
# Ver los resultados del modelo
summary(modeloPonderado)##
## Call:
## lm(formula = preciom ~ areaconst, data = bd_aptos, weights = pesos)
##
## Weighted Residuals:
## Min 1Q Median 3Q Max
## -1.0149 -1.0000 0.1317 1.0001 1.0286
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.001e+02 2.148e-02 9314 <2e-16 ***
## areaconst 4.984e-01 1.844e-04 2703 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1 on 1361 degrees of freedom
## Multiple R-squared: 0.9998, Adjusted R-squared: 0.9998
## F-statistic: 7.306e+06 on 1 and 1361 DF, p-value: < 2.2e-16# Calcular R²
summary(modelo)$r.squared## [1] 0.7162696summary(modeloPonderado)$r.squared## [1] 0.9998137# Coeficientes
modeloPonderado$coefficients## (Intercept) areaconst
## 200.0677614 0.4983516A pesar de que el R cuadrado mejoró notoriamente, se evidencia un posible sobreajuste de los datos al ser tan alto este indicador. Por otro lado, al revisar los coeficientes de este modelo son muy similares a los obtenidos con el primer modelo. En conclusión se sigue trabajando con el primer modelo obtenido ya que con este se cumplen los supuestos y se logra explicar gran parte de la variabilidad en el precio de la vivienda utilizando el área construida como predictor principal.