INFORME A&C

CONTEXTO

El siguiente informe es dirigido a la inmobiliaria A&C, la cual ha solicitado el apoyo de un científico de datos para construir un modelo que les oriente sobre los precios de los inmuebles. Este modelo se basará en los datos de ofertas de vivienda descargados del portal Fincaraiz para apartamentos de estrato \(4\) con área construida menor a \(200\) \(m^2\) (vivienda4.RDS).

El objetivo final es intentar ajustar un modelo robusto y preciso que pueda ayudar a la inmobiliaria A&C a determinar los precios de los inmuebles a negociar en el futuro.

ANÁLISIS EXPLORATORIO

Se inicia cargando las librerias que serán de utilidad para este análisis:

library(devtools)
library(paqueteMETODOS)
library(skimr)
library(car)
library(tidyverse)
library(ggpubr)
library(gridExtra)
library(GGally)
library(car)
library(bestNormalize)
library(stargazer)
library(lmtest)

Luego se cargan los datos que se encuentran en el archivo (vivienda4.RDS)

data(vivienda4)
vivienda <- as.data.frame(vivienda4)

Se hace un primer resumen de estos datos para entender si es necesario realizar algún tipo de preprocesamiento antes de avanzar con su posterior análisis y modelación.

summary(vivienda)

##            zona      estrato     preciom        areaconst     
##  Zona Centro :   8   3:   0   Min.   : 78.0   Min.   : 40.00  
##  Zona Norte  : 288   4:1706   1st Qu.:160.0   1st Qu.: 60.00  
##  Zona Oeste  :  60   5:   0   Median :210.0   Median : 75.00  
##  Zona Oriente:   6   6:   0   Mean   :225.4   Mean   : 87.63  
##  Zona Sur    :1344            3rd Qu.:265.0   3rd Qu.: 98.00  
##                               Max.   :760.0   Max.   :200.00  
##           tipo     
##  Apartamento:1363  
##  Casa       : 343  
##                    
##                    
##                    
##

Se observa que solo hay un estrato (el \(4\)) pero se mencionan más con valor de \(0\), también que hay \(5\) zonas donde la que posee la mayor cantidad de datos es la Zona Sur, seguida de la Zona Norte, pero antes de intentar calcular algún porcentaje se observa que se tienen \(343\) viviendas bajo la categoría de Casa, esto se debe quitar ya que no hace parte del interés de la inmobiliaria.

Se procede entonces con la depuración de los datos:

# Filtramos el dataframe
  vivienda_filtrada <- vivienda %>%
# Quitamos la variable estrato
  select(-estrato) %>%
# Filtramos para solo usar los Apartamentos
  filter(tipo == "Apartamento") %>%
# Corroboramos que la variable areaconst sean valores menores o iguales a 200
  filter(areaconst <= 200)
# Filtramos vivienda tipo
vivienda_filtrada$tipo <- droplevels(vivienda_filtrada$tipo)

Se realiza un nuevo resumen ahora sobre los datos de vivienda_filtrada

# Para todas las variables
summary(vivienda_filtrada)

##            zona         preciom        areaconst               tipo     
##  Zona Centro :   7   Min.   : 78.0   Min.   : 40.00   Apartamento:1363  
##  Zona Norte  : 237   1st Qu.:153.5   1st Qu.: 60.00                     
##  Zona Oeste  :  52   Median :185.0   Median : 70.00                     
##  Zona Oriente:   2   Mean   :202.4   Mean   : 75.48                     
##  Zona Sur    :1065   3rd Qu.:240.0   3rd Qu.: 84.00                     
##                      Max.   :645.0   Max.   :200.00

Respecto a la primera revisión de los datos, se observa que se retira la variable estrato, ya que solo aporta contexto al análisis, luego se observa como ahora, la variable tipo solo tiene información sobre apartamentos. Respecto al dataframe en general, se tienen \(1363\) registros, y \(4\) variables, donde:

Zona (zona), es una variable categórica nominal con 5 categorías de las cuales, la Zona Sur y Norte poseen la mayoría de los apartamentos de estrato \(4\) con el \(78.1\%\) y el \(17.4\%\) del total de registros de apartamentos respectivamente, como lo muestra la siguiente figura.

# Gráfico de barras para la variable zona
vivienda_filtrada %>%
  group_by(zona) %>%
  summarise(n = n()) %>%
  mutate(porcentaje = n / sum(n) * 100) %>%
  ggplot(aes(x=zona, y=porcentaje, fill=zona)) +
  geom_bar(stat="identity", color="black") +
  geom_text(aes(label=paste0(round(porcentaje,1),"%")), vjust=-0.3) +
  scale_fill_brewer(palette="BrBG") +
  theme_minimal() +
  labs(title="Gráfico de barras zona",
       x="zona",
       y="Porcentaje (%)") +
  theme(plot.title = element_text(hjust = 0.5))

Precio (preciom), es una variable continua, para este caso concreto, el precio de los apartamentos varía desde un mínimo de \(78\) millones COP hasta un máximo de \(645\) millones COP. El precio medio es de aproximadamente \(202.4\) millones COP, y el \(50\%\) de los apartamentos tienen un precio entre \(153.5\) y \(240\) millones COP (según el primer y tercer cuartil).
Área (areaconst), es también una variable continua. El área construida de los apartamentos varía desde un mínimo de \(40\) \(m^2\) hasta un máximo de \(200\) \(m^2\). En promedio, los apartamentos tienen un área construida de aproximadamente \(75.5\) \(m^2\), y el \(50\%\) de los apartamentos tienen un área construida entre \(60\) y \(84\) \(m^2\) (según el primer y tercer cuartil).

# Histograma para la variable preciom
p1 <- ggplot(vivienda_filtrada, aes(x=preciom)) +
  geom_histogram(binwidth=20, fill="#E0EEEE", color="black") +
  theme_minimal() +
  labs(title="Histograma de precio",
       x="Precio en millones COP",
       y="Frecuencia") +
  theme(plot.title = element_text(hjust = 0.5))

# Histograma para la variable área construida
p2 <- ggplot(vivienda_filtrada, aes(x=areaconst)) +
  geom_histogram(binwidth=10, fill="#EED5B7", color="black") +
  theme_minimal() +
  labs(title="Histograma de área construida",
       x="Área en metros cuadrados",
       y="") +
  theme(plot.title = element_text(hjust = 0.5))

grid.arrange(p1, p2, ncol=2)

Ahora, tanto la variable preciom (precio) como la variable areaconst (área construida) muestran una distribución con cola derecha pesada. Esto indica que, aunque la mayoría de los apartamentos tienen precios y áreas más bajos, hay algunos apartamentos con precios y áreas extremadamente altos. Estos valores extremos arrastran la media hacia la derecha, creando una cola derecha pesada en ambos histogramas.

En el caso del precio, la mayoría de los apartamentos tienen un precio entre \(78\) millones COP y \(240\) millones COP, pero algunos apartamentos tienen precios tan altos como \(645\) millones COP. Similarmente, para el área construida, la mayoría de los apartamentos tienen un área entre \(40 m^2\) y \(84\) \(m^2\), pero algunos apartamentos tienen áreas tan grandes como \(200\) \(m^2\).

Lo anterior sugiere que, aunque la mayoría de los apartamentos en el mercado son más asequibles y más pequeños, hay una minoría de apartamentos que son significativamente más caros y más grandes.

Para profundizar en la relación entre la variable respuesta preciom en función de la variable predictora areaconst se continua con un análisis bivariado

ggplot(vivienda_filtrada, aes(x=areaconst, y=preciom)) +
  geom_point() +
  theme_minimal() +
  labs(title="Gráfico de dispersión",
       x="Área construida en metro cuadrado",
       y="Precio en millones COP") +
  theme(plot.title = element_text(hjust = 0.5))

Resulta evidente que entre \(50\) y \(100\) metros cuadrados se concentran los apartamentos con precios entre los \(200\) y \(400\) millones COP, también parece haber una relación lineal positiva entre ambas variables por lo que surge la necesidad de revisar mediante el coeficiente de pearson, al tratarse de dos variables continuas, cual es su correlación lineal.

ggpairs(vivienda_filtrada[,2:3], title="Precio vs Área") +
  theme(plot.title = element_text(hjust = 0.5))

Se observa una correlación fuerte positiva de \(0.748\) entre la variable precio y área. Esto significa que a medida que el área de un apartamento aumenta, también lo hace su precio, en general. Vale aclarar que la correlación no expresa causalidad por lo que el aumento de precio puede estar ligado a otras variables que no se estarían considerando actualmente.

cor.test(vivienda_filtrada$areaconst,vivienda_filtrada$preciom, method = "pearson")

## 
##  Pearson's product-moment correlation
## 
## data:  vivienda_filtrada$areaconst and vivienda_filtrada$preciom
## t = 41.595, df = 1361, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.7237938 0.7706237
## sample estimates:
##       cor 
## 0.7481389

Finalmente, se realiza una prueba de hipotesis para determinar si la correlación es significativa para estas varibles.

El valor p es menor que \(2.2e^{-16}\), lo que es esencialmente cero. Esto indica que la correlación de \(0.748\) entre areaconst y preciom es estadísticamente significativa. En otras palabras, se puede rechazar la hipótesis nula de que no hay correlación entre estas dos variables.

Además, el intervalo de confianza del \(95\%\) para la correlación está entre \(0.7237938\) y \(0.7706237\), lo que indica que si se repitiera el estudio muchas veces, se esperaría que el \(95\%\) de los intervalos de confianza calculados contengan la verdadera correlación.

MODELACIÓN

Se estima el siguiente modelo

\[ precio = f(area) + ε \]

# Regresión
modelo <- lm(preciom ~ areaconst, data = vivienda_filtrada)
summary(modelo)

## 
## Call:
## lm(formula = preciom ~ areaconst, data = vivienda_filtrada)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -225.404  -23.902   -4.754   25.763  209.021 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 39.04679    4.09977   9.524   <2e-16 ***
## areaconst    2.16473    0.05204  41.595   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 43.34 on 1361 degrees of freedom
## Multiple R-squared:  0.5597, Adjusted R-squared:  0.5594 
## F-statistic:  1730 on 1 and 1361 DF,  p-value: < 2.2e-16

El modelo estimado se puede escribir de la siguiente manera:

\[ \hat{preciom_i} = 39.04679 + 2.16473 \cdot areaconst_i \]

Donde:

\(\hat{preciom}\) es la variable dependiente estimada, es decir, la predicción del modelo para el precio de un apartamento.
\(areaconst\) es la variable independiente, es decir, el área construida del apartamento.
\(39.04679\) es el intercepto estimado del modelo (\(\hat{β}_0\)), que representa la predicción del modelo para el precio de un apartamento cuando \(areaconst\) es cero. Se interpreta que si un apartamento no tuviera área este costaría \(39.1\) millones COP. Esto es ilogico pero es el resultado que se ha obtenido con la estimación.
\(2.16473\) es el coeficiente estimado para \(areaconst\) (\(\hat{β}_1\)), que representa el cambio esperado en \(\hat{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 \(2.2\) millones COP.

ggscatter(vivienda_filtrada, x = "areaconst", y = "preciom",
          add = "reg.line", conf.int = TRUE, 
          xlab = "Área construida en metro cuadrado", ylab = "Precio en millones COP") +
  stat_smooth(method = lm, se = FALSE, color = "red", linetype = "solid") +
  theme_minimal() +
  ggtitle("Gráfico de dispersión con línea de regresión") +
  theme(plot.title = element_text(hjust = 0.5))

La figura muestra el modelo ajustado en la nube de dispersión.

INFERENCIAS SOBRE EL MODELO

Se plantean las siguientes hipotesis sobre el coeficiente \(\beta_1\)

\(H_0:\) La variable ‘areaconst’ no tiene un efecto significativo en ‘preciom’, es decir, el coeficiente β1 es igual a cero. \(H_1:\)La variable ‘areaconst’ tiene un efecto significativo en ‘preciom’, es decir, el coeficiente β1 es distinto de cero.

Se construye un intervalo de confianza (\(95\%\)) para el coeficiente \(β_1\)

# Calcula el intervalo de confianza del 95% para los coeficientes del modelo
ic <- confint(modelo, level = 0.95)

# Muestra el intervalo de confianza para β1
ic["areaconst", ]

##    2.5 %   97.5 % 
## 2.062640 2.266826

Observando los límites inferior y superior del intervalo de confianza, se observa que el intervalo no incluye el valor cero, por tanto se puede rechazar la hipótesis nula de que \(β_1\) es igual a cero al nivel de confianza del \(95\%\). Por tanto se concluye que la variable areaconst tiene un efecto significativo sobre la variable preciom.

Adicional, este resultado se confirma realizando la correspondiente \(prueba\) \(t\), sobre el coeficiente \(β_1\)

# Asume que 'modelo' es tu modelo ajustado
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|) 
##  4.159515e+01 1.056327e-244

Se observa como el valor \(p\) es aproximadamente cero, por tanto es menor al valor de significancia \(\alpha=0.05\), lo que indica que se puede rechazar la hipótesis nula.

Para calcular e interpretar el indicador de bondad \(R^2\), se puede usar la función summary() en R:

# Obtenemos el resumen del modelo
resumen <- summary(modelo)

# Obtenemos el valor de R cuadrado del modelo
R2 <- resumen$r.squared

# Imprimimos el valor de R cuadrado
R2

## [1] 0.5597117

El valor de \(R^2\) representa la proporción de la varianza total de preciom que se explica por el modelo. Un valor de \(R^2\) cercano a \(1\) indica que una gran proporción de la variabilidad en preciom puede ser explicada por areaconst, mientras que un valor de \(R^2\) cercano a 0 indica lo contrario. Para este caso particular, la variable areaconst explica aproximadamente el \(56\%\) de la variabilidad de la variable preciom.

CASO DE USO

Para este caso resulta interesante responder la siguiente pregunta. ¿Cuál sería el precio promedio estimado para un apartamento de 110 metros cuadrados?

predict(modelo, data.frame(areaconst = 110), interval = "prediction", level = 0.95)

##        fit      lwr      upr
## 1 277.1674 192.0449 362.2899

Para un apartamento en estrato \(4\) de \(110\) \(m^2\) el precio promedio estimado es aproximadamente de \(277\) millones COP. Además con una confianza del \(95\%\), el precio promedio de un apartamento que tiene un área de \(110\) \(m^2\) está entre los \(192\) y los \(362\) millones COP.

Ahora, se analiza si un apartamento con este mismo tamaño de área (\(110\) \(m^2\)) con un precio de \(200\) millones COP se puede considerar una oferta atractiva de tomar.

predict(modelo, data.frame(areaconst = 110), interval = "confidence", level = 0.95)

##        fit      lwr      upr
## 1 277.1674 272.9573 281.3775

El intervalo evidencia significancia al no incluir el valor de cero. Con una confianza del \(95\%\), el precio promedio de los apartamentos que tienen en promedio \(100\) \(m^2\) de área construida en estrato \(4\), está entre los \(273\) y los \(281\) millones COP. Por ende, se puede concluir que un apartamento por \(200\) millones COP, sí resulta en una oferta atractiva.

Sin embargo, solo el precio no debería evidenciar si la oferta es atractiva o no, otros aspectos que se deberían tener en consideración son la ubicación espacial, la condición actual del apartamento, el piso en el que se encuentra el apartamento, valor de administración, ascensor, parqueadero, cercanía al trabajo, número de habitaciones y baños, entre otros factores conformes a la situación actual de vida de la persona o personas que desean habitarlo.

VALIDACIÓN DE SUPUESTOS

Para validar los diferentes supuestos del modelo se realizan las siguiente pruebas sobre el modelo actual:

# Residuales
residuos <- residuals(modelo)

# Creamos un dataframe con los residuos
dfr <- data.frame(residuos = residuals(modelo))

# Creamos el gráfico de caja
ggplot(dfr, aes(y = residuos)) +
  geom_boxplot() +
  theme_minimal() +
  labs(title = "Gráfico de caja de los residuos", y = "Residuos") +
  theme(plot.title = element_text(hjust = 0.5))

# Creamos el gráfico de dispersión
ggplot(dfr, aes(x = 1:length(residuos), y = residuos)) +
  geom_point() +
  theme_minimal() +
  labs(title = "Gráfico de dispersión de los residuos", x = "Observación", y = "Residuos") +
  theme(plot.title = element_text(hjust = 0.5))

Tanto el scatterplot, como el boxplot, evidencian outliers, ya que el boxplot muestra como minimo y máximo, valores entre \(-100\) y \(100\) por tanto, todo valor que sea superior o inferior a estos límites se van a entender como outliers.

El gráfico de residuales también es muy amplio y la diferencia entre los datos y los valores estimados por el modelo, no deberían tener diferencias tan grandes, por lo que se garantiza que si existen outliers.

Antes de continuar con la validación de supuesto se identificarán los registros que corresponden a outliers.

# Cálculo del IQR
Q1 <- quantile(residuos, 0.25)
Q3 <- quantile(residuos, 0.75)
IQR <- Q3 - Q1

# Identificación de outliers
indices_outliers <- which((residuos < (Q1 - 1.5 * IQR)) | (residuos > (Q3 + 1.5 * IQR)))

# Extraemos las observaciones que generan outliers
outliers <- vivienda_filtrada[indices_outliers, ]

# Observaciones que generan outliers
outliers

##             zona preciom areaconst        tipo
## 118   Zona Norte     340       200 Apartamento
## 126     Zona Sur     420       124 Apartamento
## 129     Zona Sur     295        50 Apartamento
## 183     Zona Sur     230       175 Apartamento
## 252     Zona Sur     250       143 Apartamento
## 265   Zona Oeste     430       124 Apartamento
## 273     Zona Sur     170       120 Apartamento
## 276     Zona Sur     165       122 Apartamento
## 280  Zona Centro     150       129 Apartamento
## 376     Zona Sur     160       160 Apartamento
## 383     Zona Sur     340        90 Apartamento
## 384     Zona Sur     222       145 Apartamento
## 395     Zona Sur     175       110 Apartamento
## 479     Zona Sur     350        68 Apartamento
## 526   Zona Norte     420       100 Apartamento
## 528   Zona Norte     360        83 Apartamento
## 568     Zona Sur     320        79 Apartamento
## 572     Zona Sur     186       133 Apartamento
## 580   Zona Oeste     210       144 Apartamento
## 582   Zona Oeste     300        74 Apartamento
## 619   Zona Oeste     200       143 Apartamento
## 664     Zona Sur     645       184 Apartamento
## 669   Zona Norte     250       160 Apartamento
## 724   Zona Norte     314        68 Apartamento
## 739     Zona Sur     140       114 Apartamento
## 769     Zona Sur     500       123 Apartamento
## 771     Zona Sur     590       159 Apartamento
## 779   Zona Norte     107       102 Apartamento
## 788     Zona Sur     250       160 Apartamento
## 796   Zona Oeste     440       126 Apartamento
## 803   Zona Oeste     350        70 Apartamento
## 807   Zona Norte     280       173 Apartamento
## 831   Zona Norte     510       121 Apartamento
## 967     Zona Sur     300        62 Apartamento
## 1017    Zona Sur     340        92 Apartamento
## 1080    Zona Sur     375       107 Apartamento
## 1158    Zona Sur     355        94 Apartamento
## 1182    Zona Sur     355        94 Apartamento
## 1224    Zona Sur     380       107 Apartamento
## 1296    Zona Sur     342        78 Apartamento

Se observa que los outliers corresponden a \(40\) registros de la base de vivienda_filtrada

Posibles tratamientos a los datos para tratar con los outliers pueden ser, eliminarlos, transformar los datos y ver si mejora el ajuste del modelo, por último analizar alguno de ellos se puede considerar un punto influyente. Más adelante se desarrollarán posibles tratamientos para estos casos.

Ahora se prueba normalidad:

# Test de Shapiro - Wilk
shapiro.test(residuos)

## 
##  Shapiro-Wilk normality test
## 
## data:  residuos
## W = 0.96486, p-value < 2.2e-16

Dado que el valor p es muy pequeño (menor a \(2.2e-16\)), se rechaza la hipótesis nula y se concluye que los datos no presentan evidencia suficiente sobre seguir una distribución normal.

# Prueba Breusch-Pagan para la homocedasticidad
bptest <- bptest(modelo)
bptest

## 
##  studentized Breusch-Pagan test
## 
## data:  modelo
## BP = 292.99, df = 1, p-value < 2.2e-16

Dado que el valor p es muy pequeño (menor a \(2.2e-16\)), se puede rechazar la hipótesis nula y concluir que la varianza de los errores no es constante, lo que indica heterocedasticidad.

La heterocedasticidad puede ser problemática porque puede hacer que los estimadores de mínimos cuadrados ordinarios (MCO) sean ineficientes y puede sesgar las pruebas de hipótesis.

# Prueba Durbin-Watson para la no autocorrelación
dwtest <- dwtest(modelo)
dwtest

## 
##  Durbin-Watson test
## 
## data:  modelo
## DW = 1.443, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0

El valor p es menor a \(2.2e-16\), lo que indica que hay evidencia significativa de autocorrelación positiva en los residuos del modelo.

Por lo anterior se concluye que el modelo, incumple todos los supuestos del modelo. Para resolver esto, se puede probar realizar transformaciones a la variable dependiente, o incluso a la independiente, otra opción puede ser revisar puntos influyentes y ver si estos afectan al modelo quitandolos y revisando si los supuestos y ajuste mejoran, revisar el tratamiento de outliers es otra opción para intentar mejorar el modelo y sus supuestos, revisar la inclusión de otras variables que se asocien al precio de un apartamento estrato 4, revisar otro tipo de modelos puesto que los modelos de regresión no simpre son la mejor opción para modelar este tipo de datos que al parecer no siguen una distribución normal.

TRANSFORMACIÓN BOX COX

Revisando la parte de transformaciones, se estima y construye un nuevo modelo, a partir de la metodología de boxcox

MASS::boxcox(modelo)

# Calcula la transformación Box-Cox
resultado_boxcox <- MASS::boxcox(modelo, plotit = FALSE)

# Encuentra el valor de lambda que maximiza la log-verosimilitud
lambda_optimo <- resultado_boxcox$x[which.max(resultado_boxcox$y)]
lambda_optimo

## [1] -0.1

Por lo anterior se entiende que el valor de la trasnformación a realizar, corresponde a \(-0.1\) por lo que usando la función bcPower se procede a realizar dicha trasnformación

# Asume que 'df' es tu dataframe y 'lambda_optimo' es el valor de lambda que obtuviste
vivienda_filtrada$preciom_transformado <- bcPower(vivienda_filtrada$preciom, lambda_optimo)

El nuevo modelo estimado con la transformación aplicada es

modelo2 <- lm(preciom_transformado ~ areaconst, data = vivienda_filtrada)
summary(modelo2)

## 
## Call:
## lm(formula = preciom_transformado ~ areaconst, data = vivienda_filtrada)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.57687 -0.07698 -0.00657  0.09361  0.38737 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 3.6741110  0.0114583  320.65   <2e-16 ***
## areaconst   0.0055180  0.0001455   37.94   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1211 on 1361 degrees of freedom
## Multiple R-squared:  0.514,  Adjusted R-squared:  0.5136 
## F-statistic:  1439 on 1 and 1361 DF,  p-value: < 2.2e-16

A este modelo también se decide probar sus supuestos, pese a que su ajuste al usar la transformación en realidad empeoró frente al modelo inicial.

# Residuales
residuos2 <- residuals(modelo2)

# Creamos un dataframe con los residuos
dfr2 <- data.frame(residuos2 = residuals(modelo2))

# Creamos el gráfico de caja
ggplot(dfr2, aes(y = residuos2)) +
  geom_boxplot() +
  theme_minimal() +
  labs(title = "Gráfico de caja de los residuos modelo 2", y = "Residuos") +
  theme(plot.title = element_text(hjust = 0.5))

# Creamos el gráfico de dispersión
ggplot(dfr2, aes(x = 1:length(residuos2), y = residuos2)) +
  geom_point() +
  theme_minimal() +
  labs(title = "Gráfico de dispersión de los residuos modelo 2", x = "Observación", y = "Residuos") +
  theme(plot.title = element_text(hjust = 0.5))

La transformación evidencia una escala disminuida (entre \(-0.6\) y \(0.4\)), lo que cambia la interpretación de los coeficientes del modelos. También aparenta mostrar menos cantidad de outliers. Se procede a identificar el número de outliers para confirmar si esto es real.

# Cálculo del IQR
Q1_2 <- quantile(residuos2, 0.25)
Q3_2 <- quantile(residuos2, 0.75)
IQR_2 <- Q3_2 - Q1_2

# Identificación de outliers
indices_outliers_2 <- which((residuos2 < (Q1_2 - 1.5 * IQR_2)) | (residuos2 > (Q3_2 + 1.5 * IQR_2)))

# Extraemos las observaciones que generan outliers
outliers_2 <- vivienda_filtrada[indices_outliers_2, ]

# Observaciones que generan outliers
outliers_2

##            zona preciom areaconst        tipo preciom_transformado
## 118  Zona Norte     340       200 Apartamento             4.417200
## 129    Zona Sur     295        50 Apartamento             4.337375
## 183    Zona Sur     230       175 Apartamento             4.194666
## 276    Zona Sur     165       122 Apartamento             3.998613
## 280 Zona Centro     150       129 Apartamento             3.941141
## 376    Zona Sur     160       160 Apartamento             3.980118
## 479    Zona Sur     350        68 Apartamento             4.433359
## 572    Zona Sur     186       133 Apartamento             4.070082
## 592    Zona Sur      90        53 Apartamento             3.623597
## 593    Zona Sur      78        46 Apartamento             3.531694
## 619  Zona Oeste     200       143 Apartamento             4.112960
## 739    Zona Sur     140       114 Apartamento             3.899194
## 779  Zona Norte     107       102 Apartamento             3.732972
## 803  Zona Oeste     350        70 Apartamento             4.433359

Se evidencia gráficamente que los outliers con la transformación box-cox se reducen de \(40\) a \(14\) valores.

Se prueban nuevamente los supuestos sobre el modelo transformado para confirmar si los supuestos se cumplen:

# Prueba de Shapiro-Wilk para la normalidad
shapiro_test <- shapiro.test(residuos2)
shapiro_test

## 
##  Shapiro-Wilk normality test
## 
## data:  residuos2
## W = 0.99002, p-value = 5.194e-08

# Prueba Breusch-Pagan para la homocedasticidad
bptest <- bptest(modelo2)
bptest

## 
##  studentized Breusch-Pagan test
## 
## data:  modelo2
## BP = 135.15, df = 1, p-value < 2.2e-16

# Prueba Durbin-Watson para la no autocorrelación
dwtest <- dwtest(modelo2)
dwtest

## 
##  Durbin-Watson test
## 
## data:  modelo2
## DW = 1.312, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0

El supuesto de normalidad, homocedasticidad e independencia de los errores como evidencian las preubas anteriores siguen sin cumplirse, pero parecen haber mejorado un poco, respecto al modelo inicial.

COMPARACIÓN DE MODELOS

Ya que los modelos siguen sin cumplir los supuestos, compararlos no tiene mucho sentido, sin embargo el modelo inicial posee un \(R^2=0.5597\), y el modelo con la transformación de box-cox tiene un \(R^2=0.514\). Se sugiere seguir revisando transformaciones para intentar encontrar un modelo que mejore tanto el ajuste y por ende la explicabilidad del área sobre el precio, como que permita realizar inferencia validas al cumplir con los diferentes supuestos de los modelos de regresión.

PROBANDO MÁS MODELOS

modelo1_1=lm(preciom ~ areaconst, data=vivienda_filtrada)           # Lin - Lin
modelo2_1=lm(preciom ~ log(areaconst), data=vivienda_filtrada)      # Lin - Log
modelo3_1=lm(log(preciom) ~ areaconst, data=vivienda_filtrada)      # Log - Lin
modelo4_1=lm(log(preciom) ~ log(areaconst), data=vivienda_filtrada) # Log - Log
modelo5_1=lm(sqrt(preciom) ~ areaconst, data=vivienda_filtrada)
modelo6_1=lm((1/preciom) ~ areaconst, data=vivienda_filtrada)


stargazer(modelo1_1, modelo2_1, modelo3_1, modelo4_1, modelo5_1, modelo6_1, type="text", df=FALSE)

## 
## ==================================================================================================
##                                                  Dependent variable:                              
##                     ------------------------------------------------------------------------------
##                              preciom                log(preciom)        sqrt(preciom) (1/preciom) 
##                         (1)          (2)          (3)          (4)           (5)          (6)     
## --------------------------------------------------------------------------------------------------
## areaconst             2.165***                  0.009***                  0.071***    -0.00005*** 
##                       (0.052)                   (0.0002)                   (0.002)     (0.00000)  
##                                                                                                   
## log(areaconst)                    195.419***                 0.882***                             
##                                    (4.445)                   (0.020)                              
##                                                                                                   
## Constant             39.047***   -635.532***    4.551***     1.484***     8.728***      0.009***  
##                       (4.100)      (19.092)     (0.019)      (0.087)       (0.138)      (0.0001)  
##                                                                                                   
## --------------------------------------------------------------------------------------------------
## Observations           1,363        1,363        1,363        1,363         1,363        1,363    
## R2                     0.560        0.587        0.520        0.582         0.545        0.453    
## Adjusted R2            0.559        0.587        0.519        0.582         0.545        0.453    
## Residual Std. Error    43.339       41.982       0.205        0.191         1.458        0.001    
## F Statistic         1,730.157*** 1,933.199*** 1,473.424*** 1,894.288*** 1,629.340***  1,127.741***
## ==================================================================================================
## Note:                                                                  *p<0.1; **p<0.05; ***p<0.01

Las transformaciones apenas si mejoran el ajuste del modelo inicial. Se decide probar un modelo sin intercepto.

# Ajusta el modelo sin intercepto
modelo_sin_intercepto <- lm(preciom ~ 0 + areaconst, data = vivienda_filtrada)
summary(modelo_sin_intercepto)

## 
## Call:
## lm(formula = preciom ~ 0 + areaconst, data = vivienda_filtrada)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -262.343  -16.933    3.018   28.829  190.603 
## 
## Coefficients:
##           Estimate Std. Error t value Pr(>|t|)    
## areaconst  2.63964    0.01538   171.6   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 44.74 on 1362 degrees of freedom
## Multiple R-squared:  0.9558, Adjusted R-squared:  0.9557 
## F-statistic: 2.944e+04 on 1 and 1362 DF,  p-value: < 2.2e-16

Al quitar el intercepto, el ajuste del modelo incrementa de forma abismal \(R^2=0.9558\) comparado con todos los demás modelos que se habían revisado (hasta el momento se han revisado \(7\) modelos). Por lo que se decide revisar los supuestos del modelo sin intercepto.

# Residuales
residuos111 <- residuals(modelo_sin_intercepto)

# Prueba de Shapiro-Wilk para la normalidad
shapiro_test <- shapiro.test(residuos111)
print(shapiro_test)

## 
##  Shapiro-Wilk normality test
## 
## data:  residuos111
## W = 0.94572, p-value < 2.2e-16

# Calcula los valores ajustados
valores_ajustados <- fitted(modelo_sin_intercepto)

# Crea un gráfico de residuos vs valores ajustados
plot(valores_ajustados, residuos111,
     xlab = "Valores ajustados",
     ylab = "Residuos",
     main = "Gráfico de residuos vs valores ajustados")

# Añade una línea horizontal en y = 0
abline(h = 0, lty = 2)

# Prueba Durbin-Watson para la no autocorrelación
dwtest <- dwtest(modelo_sin_intercepto)
print(dwtest)

## 
##  Durbin-Watson test
## 
## data:  modelo_sin_intercepto
## DW = 1.433, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0

# Visualización gráfica para detectar outliers
boxplot(residuos111)

plot(residuos111)

La prueba de homocedasticidad Breush-Pagan mediante su función bptest() no permite evaluar modelos sin intercepto por lo que se construye una gráfica en su lugar para hacer la respectiva revisión de este supuesto. Sin embargo se observa que al igual que los modelos anteriores los supuestos del modelo siguen mostrar evidencia suficiente de que se logren cumplir.

Revisando un poco más de literatura, Minitab en su web menciona que si se incumple el supuesto de homocedasticidad en un modelo, es posible usar modelos de regresión ponderados para corregir este supuesto. A continuación se construye un código, con el que se espera probar si es cierto lo que menciona esta web. https://support.minitab.com/es-mx/minitab/21/help-and-how-to/statistical-modeling/regression/supporting-topics/basics/weighted-regression/.

# Cargar la librería necesaria
library(lmtest)

# Ajustar el modelo de regresión lineal
modeloP <- lm(preciom ~ areaconst, data = vivienda_filtrada)

# Realizar la prueba de Breusch-Pagan para detectar heterocedasticidad
bp_test <- bptest(modeloP)

# Si el p-valor es menor que 0.05, entonces hay evidencia de heterocedasticidad
if (bp_test$p.value < 0.05) {
  # En este caso, podríamos intentar ajustar un modelo de regresión ponderada
  pesos <- 1 / fitted(modeloP)^2
  modelo_ponderado <- lm(preciom ~ areaconst, data = vivienda_filtrada, weights = pesos)
}

summary(modelo_ponderado)

## 
## Call:
## lm(formula = preciom ~ areaconst, data = vivienda_filtrada, weights = pesos)
## 
## Weighted Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.65562 -0.11959 -0.01259  0.12355  1.05325 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  15.8714     4.2244   3.757 0.000179 ***
## areaconst     2.4800     0.0601  41.263  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.19 on 1361 degrees of freedom
## Multiple R-squared:  0.5558, Adjusted R-squared:  0.5554 
## F-statistic:  1703 on 1 and 1361 DF,  p-value: < 2.2e-16

Respecto al modelo inicial lm(preciom ~ areaconst, data = vivienda_filtrada) su \(R^2\), se podría decir que ambos modelos explican el mismo porcentaje del precio de los apartamentos, o ajusta de forma similar, ahora se hace la prueba de homocedasticidad la cual se espera, sea reparada por el método.

# Realizar la prueba de Breusch-Pagan para detectar heterocedasticidad
bp_test <- bptest(modelo_ponderado)
bp_test

## 
##  studentized Breusch-Pagan test
## 
## data:  modelo_ponderado
## BP = 0.0053968, df = 1, p-value = 0.9414

Se observa que no se rechaza la hipótesis nula, por tanto hay evidencia suficiente para afirmar que el modelo es homocedastico. Se revisa los otros supuestos.

# Prueba de Shapiro-Wilk
shapiro_test <- shapiro.test(modelo_ponderado$residuals)
print(shapiro_test)

## 
##  Shapiro-Wilk normality test
## 
## data:  modelo_ponderado$residuals
## W = 0.9558, p-value < 2.2e-16

library(nortest)
# Prueba de Anderson-Darling
anderson_test <- ad.test(modelo_ponderado$residuals)
print(anderson_test)

## 
##  Anderson-Darling normality test
## 
## data:  modelo_ponderado$residuals
## A = 11.336, p-value < 2.2e-16

# QQPlot
qqnorm(modelo_ponderado$residuals)
qqline(modelo_ponderado$residuals, col = "steelblue", lwd = 2)

No se cumple el supuesto de normalidad. El algoritmo de autocorrelación que se usa para la prueba Durbin Watson no funciona por lo que se revisa de forma gráfica los residuales.

# Calcula los valores ajustados
valores_ajustados2 <- fitted(modelo_ponderado)

# Crea un gráfico de residuos vs valores ajustados
plot(valores_ajustados2, modelo_ponderado$residuals,
     xlab = "Valores ajustados",
     ylab = "Residuos",
     main = "Gráfico de residuos vs valores ajustados")

# Añade una línea horizontal en y = 0
abline(h = 0, lty = 2)

El qqplot muestra que el modelo posee colas pesadas, tanto inferior como superior, por lo que debe serguirse revisando que técnicas pueden servir para corregir este tipo de comportamiento natural de los datos, así que aunque se cumple el supuesto de homocedasticidad con esta forma de modelación, no se corrigen los demás supuestos.

CONCLUSIONES

Es evidente que hay relación entre las variables precio y área de los apartamentos de estrato \(4\) con menos de \(200\) \(m^2\), como lo evidencia la prueba de correlación, sin embargo se hace notable la presencia de outliers en los datos, por lo que resulta de interés realizar algún tipo de tratamiento que permita trabajar con ellos, o hasta excluírlos si esto permite mejorar el modelo, sin embargo, algo que también se revisó fue la extracción de puntos influyentes estos arrojaron una mejora en el ajuste, sin embargo no se dejan consignados en el informe ya que existia una buena cantidad de modelos ya probados.

La modelación quitando el intercepto sorprende mostrando un \(R^2=0.95\), y en realidad resulta bastante lógico quitarlo, ya que un apartamento sin área (no existe), no debería de poseer un precio, cosa que pasa con todos los demás modelos, sin embargo al no cumplir con los supuestos no evidencia ser tampoco la mejor opción a la hora de hacer la elección de un modelo, que sirva para extrapolar los precios de los apartamentos.

Las transformaciones no fueron realmente útiles ya que además de que afectan directamente la interpretación de los coeficientes del modelo, no lograron solucionar el problema de incumplimiento de los supuestos que es para lo que se supone sirve realizarlas, la recomendación respecto a esto sería pensar en realizar otro tipo de modelación o conseguir más covariables si es posible para mejorar el modelo en cuestión antes de pensar en transformar variables que no tienen interpretación al transformarse.

Por último, el modelo ponderado, tendría que revisarse su interpretación en los coeficientes como cambia frente a un modelo sin ponderar, pero aparentemente logra arreglar el supuesto de homocedasticidad, algo que resulta interesante de seguir estudiando como revisar si existen otras ponderaciones que pudiesen reparar los demás supuestos, combinar esta metodología con transformaciones y hasta revisar si es posible ponderar el modelo sin intercepto.

Por ahora, el modelo base, el modelo sin intercepto y el modelo ponderado, son los modelos más interesantes de seguir revisando ya que las transformaciones alteran la interpretación de los coeficientes, y obtener una predicción o extrapolación de la variable precio necesitará de un procedimiento adicional para recuperar el precio real que se estima con este tipo de modelos. Para continuar con el estudio hace falta realizar pruebas que involucren validación cruzadas y ayuden a complementar el \(R^2\), esto en caso que no se pudieran agregar más covariables que mejoren el modelo actual o se agreguen indexaciones como el espacio mediante coordenadas, para hacer otro tipo de modelos que involucren la espacialidad.