Enunciado

María comenzó como agente de bienes raíces en Cali hace 10 años. Después de laborar dos años para una empresa nacional, se trasladó a Bogotá y trabajó para otra agencia de bienes raíces. Sus amigos y familiares la convencieron de que con su experiencia y conocimientos del negocio debía abrir su propia agencia. Terminó por adquirir la licencia de intermediario y al poco tiempo fundó su propia compañía, C&A (Casas y Apartamentos) en Cali. Santiago y Lina, dos vendedores de la empresa anterior aceptaron trabajar en la nueva compaña. En la actualidad ocho agentes de bienes raíces colaboran con ella en C&A.

Actualmente las ventas de bienes raíces en Cali se han visto disminuidas de manera significativa en lo corrido del año. Durante este periodo muchas instituciones bancarias de ahorro y vivienda están prestando grandes sumas de dinero para la industria y la construcción comercial y residencial. Cuando el efecto producto de las tensiones políticas y sociales disminuya, se espera que la actividad económica de este sector se reactive.

Hace dos días, María recibió una carta solicitando asesoría para la compra de dos viviendas por parte de una compañía internacional que desea ubicar a dos de sus empleados con sus familias en la ciudad. Las solicitudes incluyen las siguientes condiciones:

Table 1: Casos para compra de viviendas.
Caracteristicas Vivienda_1 Vivienda_2
Tipo Casa Apartamento
área_construida 200 300
parqueaderos 1 3
baños 2 3
habitaciones 4 5
estrato 4 o 5 5 o 6
zona Norte Sur
crédito_preaprobado 350 850

Ayude a María a responder la solicitud, mediante técnicas modelación que usted conoce. Ella requiere le envíe un informe ejecutivo donde analice los dos casos y sus recomendaciones. Como soporte del informe debe anexar las estimaciones, validaciones y comparación de modelos requeridos.

Análisis del dataset

data("vivienda")
datos <- vivienda

Como ya conocemos el dataset, realizaremos el mismo proceso de análisis respecto a sus atributos:

# Diferenciamos las variables categóricas de las numéricas
datos_numericos <- datos %>% select(where(is.numeric), -id, -estrato)
datos_categoricos <- datos %>% select(where(is.factor), where(is.character), estrato) %>%
  mutate_if(is.character, as.factor) %>%
  mutate(estrato = as.ordered(estrato))

# imputamos con el método pmm
datos_numericosTemp <- mice(datos_numericos, m=5, maxit=50, meth='pmm', seed=500)
datos_numericos <- complete(datos_numericosTemp, 1)

# Evitamos por un momento la variable barrio, por sus muchos niveles
datos_categoricos_sinbarrio <- datos_categoricos %>% select(-barrio)

# Hacemos la imputación a las varaibles categóricas
datos_categoricosTemp <- mice(datos_categoricos_sinbarrio, m=3, maxit=10, meth='polyreg', seed=500)
datos_categoricos_imputados <- complete(datos_categoricosTemp, 1)

# Combinamos los datos nuevamente
datos_categoricos <- cbind(datos_categoricos_imputados, barrio =  datos_categoricos$barrio)

# Encontramos la moda de la variable barrio
moda_barrio <- names(sort(table(datos_categoricos$barrio), decreasing = TRUE)[1])

# imputamos la moda en los valores faltantes
datos_categoricos$barrio[is.na(datos_categoricos$barrio)] <- moda_barrio

# Volvemos a unir todos los datos
datos_manejados <- cbind(datos_categoricos, datos_numericos)

Caso 1

Paso 1:

Realice un filtro a la base de datos e incluya solo las ofertas de \(base1\): casas de la zona norte de la ciudad. Presente los primeros \(3\) registros de las bases y algunas tablas que comprueben la consulta. (Adicional un mapa con los puntos de las bases. Discutir si todos los puntos se ubican en la zona correspondiente o se presentan valores en otras zonas, ¿por qué?).

Para este paso, realizaremos un filtro respecto al primer caso.

datos_filtrados <- datos_manejados %>%
  filter(tipo == "Casa" & zona == "Zona Norte")

knitr::kable(head(datos_filtrados))
zona piso tipo estrato barrio preciom areaconst parqueaderos banios habitaciones longitud latitud
Zona Norte 02 Casa 5 acopi 320 150 2 4 6 -76.51341 3.47968
Zona Norte 02 Casa 5 acopi 780 380 2 3 3 -76.51674 3.48721
Zona Norte 02 Casa 6 acopi 750 445 4 7 6 -76.52950 3.38527
Zona Norte 02 Casa 4 acopi 625 355 3 5 5 -76.53179 3.40590
Zona Norte 02 Casa 5 acopi 750 237 2 6 6 -76.54044 3.36862
Zona Norte 02 Casa 4 acopi 600 160 1 4 5 -76.55210 3.42125

Verificamos que se haya hecho el filtro correctamente:

fil_tipo <- datos_filtrados['tipo'] %>% unique()
fil_zona <- datos_filtrados['zona'] %>% unique()

knitr::kable(cbind(fil_tipo, fil_zona))
Table 2: Tabla del filtrado.
tipo zona
Casa Zona Norte

Y a través de una gráfica de barras, podemos analizar mejor el tipo de casas existente.

Gráfica de barras respecto al tipo de propiedad y zona.

Figure 1: Gráfica de barras respecto al tipo de propiedad y zona.

Ahora, que ya tenemos hecho el filtro, graficaremos la ubicación de la vivienda respecto a la longitud y latitud de cada una, para asegurarnos que cada propiedad corresponda al filtro realizado.

Figure 2: Mapa de viviendas filtradas para Caso 1.

Como se nota en la figura 2, no todas las viviendas se encuentra en lo que podemos denominar la zona norte de la ciudad, así que haremos uso de otros datos para mejorar este problema fundamental. En este caso, usaremos datos de Divipola, o División Política Administrativa de la ciudad de Cali.

# Accedemos a los datos descargados
ruta_geodatabase <- "C:\\Users\\lvhue\\OneDrive\\Maestría\\Segundo semestre\\Modelos Estadísticos\\Unidad 2\\Actividad\\divipola.gdb"

# Verificamos qué exista la capa requerida, la de los barrios
st_layers(ruta_geodatabase)
Driver: OpenFileGDB 
Available layers:
                         layer_name     geometry_type features fields
1                           barrios     Multi Polygon      337      5
2                           comunas     Multi Polygon       22      4
3                    corregimientos     Multi Polygon       15      5
4                   suelo_expansion     Multi Polygon       20      5
5                           veredas     Multi Polygon       98      7
6                    T_1_DirtyAreas     Multi Polygon       12      3
7                   T_1_PointErrors             Point        0      7
8                    T_1_LineErrors Multi Line String        5      8
9                    T_1_PolyErrors     Multi Polygon       57      9
10 modelo_localidades_cali_distrito     Multi Polygon       10      6
11         bcs_lim_perimetro_urbano     Multi Polygon        1      4
12              perimetro_distrital     Multi Polygon        1      6
13              clasificacion_suelo     Multi Polygon        7      4
14 conflicto_perim_urbano_vs_barrio     Multi Polygon        7      3
                         crs_name
1  MAGNA-SIRGAS / Cali urban grid
2  MAGNA-SIRGAS / Cali urban grid
3  MAGNA-SIRGAS / Cali urban grid
4  MAGNA-SIRGAS / Cali urban grid
5  MAGNA-SIRGAS / Cali urban grid
6  MAGNA-SIRGAS / Cali urban grid
7  MAGNA-SIRGAS / Cali urban grid
8  MAGNA-SIRGAS / Cali urban grid
9  MAGNA-SIRGAS / Cali urban grid
10 MAGNA-SIRGAS / Cali urban grid
11 MAGNA-SIRGAS / Cali urban grid
12 MAGNA-SIRGAS / Cali urban grid
13 MAGNA-SIRGAS / Cali urban grid
14 MAGNA-SIRGAS / Cali urban grid
# Leemos la capa
barrios_cali <- st_read(ruta_geodatabase, layer = "barrios")
Reading layer `barrios' from data source 
  `C:\Users\lvhue\OneDrive\Maestría\Segundo semestre\Modelos Estadísticos\Unidad 2\Actividad\divipola.gdb' 
  using driver `OpenFileGDB'
Simple feature collection with 337 features and 5 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 1054098 ymin: 860192.1 xmax: 1068492 ymax: 879000.7
Projected CRS: MAGNA-SIRGAS / Cali urban grid
# La transformamos a WGS84 para que no haya errores
barrios_cali_wgs84 <- st_transform(barrios_cali, "+proj=longlat +datum=WGS84")

# Queremos el centroide de cada barrio para colocar allí el barrio, si es necesario
centroides_barrios <- st_centroid(barrios_cali_wgs84)
Warning: st_centroid assumes attributes are constant over geometries
centroides_coords <- as.data.frame(st_coordinates(centroides_barrios))
names(centroides_coords) <- c("long_centro", "lat_centro")

barrios_con_centro <- cbind(barrios_cali_wgs84, centroides_coords)

# Creamos un mapa base solo con los polígonos de cada barrio
mapabase <- leaflet(barrios_con_centro) %>%
  addProviderTiles(providers$OpenStreetMap, group = "Mapa base") %>%
  addPolygons(
    fillColor = "white",
    fillOpacity = 0.5,
    weight = 1,
    color = "black",
    group = "Barrios",
    popup = ~BARRIO
  ) 

# Agregamos los nombres
mapabarrios <- mapabase %>%
  addLabelOnlyMarkers(
    lng = ~long_centro, 
    lat = ~lat_centro, 
    label = ~BARRIO,
    labelOptions = labelOptions(
      noHide = TRUE, 
      direction = "center",
      textOnly = TRUE,
      style = list(
        "font-weight" = "bold",
        "font-size" = "10px",
        "color" = "black",
        "opacity" = 0.2
      )
    ),
    group = "Nombres de Barrios"
  ) 

# Colocamos las viviendas 
mapabase %>%
  addCircleMarkers(
    data = datos_manejados,
    lng = ~longitud,
    lat = ~latitud,
    radius = 8,
    stroke = FALSE,
    fillOpacity = 0.5,
    color = "red",
    group = "Viviendas",
    popup = ~paste("Precio:", preciom, "millones", "<br>", "Área:", areaconst, "m²")
  ) %>% # Agregamos un control de cada capa para poder interactuar con ellas
  addLayersControl(
    baseGroups = c("Mapa base"),
    overlayGroups = c("Barrios", "Viviendas"),
    options = layersControlOptions(collapsed = FALSE)
  )

Figure 3: Mapa de las viviendas con Divipola.

Sin embargo, este paso no nos ayudan si los datos de barrios o zonas están equivocados (dando un salto de fe a que los datos de longitud y latitud de las propiedades sean correctos). El manejo de estos errores se dará a través de la data Divipola. Por eso, transformaremos los datos que tenemos en datos espaciales y haremos un \(join\) de los datos respecto a las coordenadas, para obtener el barrio correcto donde corresponde cada vivienda.

# Convertimos los datos como espaciales
viviendas_sf <- st_as_sf(datos_manejados, coords = c("longitud", "latitud"), crs = st_crs(barrios_cali_wgs84))

# Y unimos los datos
viviendas_con_barrio_correcto <- st_join(viviendas_sf, barrios_cali_wgs84 %>% select(BARRIO, COMUNA))

knitr::kable(head(viviendas_con_barrio_correcto))
Table 3: Tabla del join entre datos Divipola y data de viviendas.
zona piso tipo estrato barrio preciom areaconst parqueaderos banios habitaciones BARRIO COMUNA geometry
Zona Oriente 03 Casa 3 20 de julio 250 70 1 3 6 20 de Julio 11 POINT (-76.51168 3.43382)
Zona Oriente 01 Casa 3 20 de julio 320 120 1 2 3 20 de Julio 11 POINT (-76.51237 3.43369)
Zona Oriente 02 Casa 3 20 de julio 350 220 2 2 4 El Prado 11 POINT (-76.51537 3.43566)
Zona Sur 02 Casa 4 3 de julio 400 280 3 5 3 3 de Julio 19 POINT (-76.54 3.435)
Zona Norte 01 Apartamento 5 acopi 260 90 1 2 3 Porvenir 04 POINT (-76.5135 3.45891)
Zona Norte 01 Apartamento 5 acopi 240 87 1 3 3 Lili 17 POINT (-76.517 3.36971)

Como podemos darnos cuenta en la tabla 3, hemos decidido también quedarnos con los datos de la columna COMUNA, pues será la forma en que manejaremos los datos de la columna Zona.

datos_robustos <- datos_manejados %>%
  st_drop_geometry() %>%
  select(-barrio)
  
datos_robustos$barrio <-viviendas_con_barrio_correcto$BARRIO
datos_robustos$comuna <-viviendas_con_barrio_correcto$COMUNA

Antes de seguir, verificamos si existe algún dato que no puedo ser etiquetados respecto al barrio y la comuna.

datona <- datos_robustos %>% filter(is.na(viviendas_con_barrio_correcto$BARRIO))

mapabase %>%
  addCircleMarkers(
    data = datona,
    lng = ~longitud,
    lat = ~latitud,
    radius = 8,
    color = "red",
    fillOpacity = 0.8
  )

Figure 4: Mapa para la vivienda no clasificada.

Respecto a los datos de Divipola y la figura 4, la vivienda no se encuentra dentro de las periferias de la ciudad de Cali por lo que decidiremos eliminarlo de los datos.

datos_robustos <- datos_robustos %>% filter(!is.na(barrio))

Nos hemos quedado con la columna COMUNA para utilizarla sobre la columna Zona. Según datos de la División Administrativa de Cali y agendas regionales como la \(UNIMINUTO\), podemos catalogar la distribución de zonas geográficas respecto a las localidades, que sería una forma más exacta.

Table 4: TAbla de localidades y comunas respectivas.
Zonas Comunas
Norte 02, 04, 05, 06
Centro Histórico 03
Centro 09, 08, 12, 07, 10, 19, 11
Oriente 13, 14, 21, 16, 15
Sur 17, 22
Oeste 18, 20, 01

Por lo tanto, hacemos los cambios tomados de la página \(Mapa \:de\:Zonas\: Geográficas\) de la gobernación de Cali como referencia, pero debido a que no hay datos para importar, se debió analizar el mapa de manera manual, lo que podría generar algunos errores.

comuna_localidad <- data.frame(
  comuna = c(
    "02", "04", "05", "06",
    "03",
    "09", "08", "12", "07", "10", "19", "11",
    "13", "14", "21", "16", "15",
    "17", "22",
    "18", "20", "01"
  ),
  localidad = c(
    rep("Norte", 4),
    "Centro Historico",
    rep("Centro", 7),
    rep("Oriente", 5),
    rep("Sur", 2),
    rep("Oeste", 3)
  )
)

datos_finales <- datos_robustos %>%
  left_join(comuna_localidad, by = "comuna") %>%
  mutate(zona = localidad) %>%
  select(-localidad)

Veamos el comportamiento de las zonas más exactas:

num_zonas <- nlevels(datos_finales$zona)
paleta_zonas <- colorFactor(
  palette = brewer.pal(min(num_zonas, 9), "Set1"),
  domain = datos_finales$zona
)

mapabase %>%
  addCircleMarkers(
    data = datos_finales, 
    lng = ~longitud,
    lat = ~latitud,
    radius = 4,
    color = ~paleta_zonas(zona),
    stroke = FALSE,
    fillOpacity = 0.8,
    popup = ~paste0(
      "<b>Tipo:</b> ", tipo, "<br/>",
      "<b>Zona:</b> ", zona, "<br/>",
      "<b>Barrio:</b> ", barrio, "<br/>",
      "<b>Precio:</b> ", preciom, " millones"
    ),
    group = "Viviendas por Zona"
  ) %>%
  addLegend(
    pal = paleta_zonas,
    values = unique(datos_finales$zona),  
    title = "Zona de la Propiedad",
    opacity = 1
  )

Figure 5: Mapa con las zonas de Cali.

Sin embargo, respecto a los datos de zonas geográficos, figura 5, todavía hay algunas zonas, significativamente en la zona centro, que pertenecen a la zona oeste. Ya que nuestro enfoque es en las zonas Norte y Sur, podría dejar como está, puesto que es debido a que la distribución por zonas no es exacta respecto a las comunas, o podríamos cambiar la etiqueta de muchos barrios, pero es una tarea engorrosa que podría ser útil en otros casos, no en el actual.

Finalmente, realizados los cambios correspondientes, podemos hacer el análisis de los diversos casos. Entonces, hacemos nuevamente el filtro del \(caso 1\):

datos_caso1 <- datos_finales %>%
  filter(tipo == "Casa" & zona == "Norte")

knitr::kable(head(datos_caso1))
Table 5: Primeros datos del filtrado con datos manejados.
zona piso tipo estrato preciom areaconst parqueaderos banios habitaciones longitud latitud barrio comuna
Norte 02 Casa 5 320 150 2 4 6 -76.51341 3.47968 Urbanización La Merced 02
Norte 02 Casa 5 780 380 2 3 3 -76.51674 3.48721 Urbanización La Flora 02
Norte 01 Casa 3 180 120 1 3 3 -76.49768 3.47060 Los Andes 05
Norte 03 Casa 5 520 455 3 5 4 -76.49966 3.46284 El Sena 05
Norte 03 Casa 3 380 300 2 5 8 -76.50743 3.46566 Manzanares 04
Norte 03 Casa 5 395 165 2 4 4 -76.51797 3.47651 Vipasa 02 
mapabase %>% 
  addCircleMarkers(
    data = datos_caso1,
    lng = ~longitud,
    lat = ~latitud,
    radius = 8,
    color = "red", 
    fillOpacity = 0.8
  )

Figure 6: Mapa de la zona norte de Cali.

Ahora estamos seguros que las viviendas escogidas pertenecen a la zona norte de Cali, como se visualiza en la figura 6.

Paso 2:

Realice un análisis exploratorio de datos enfocado en la correlación entre la variable respuesta (precio de la casa) en función del área construida, estrato, numero de baños, número de habitaciones y zona donde se ubica la vivienda. Use gráficos interactivos con el paquete plotly e interprete los resultados.

datos_analisis_estrato <- datos_caso1 %>%
  select(preciom, areaconst, estrato, banios, habitaciones)

La función ggpairs nos brindará mucha información de la correlación entre variables.

corr1 <- ggpairs(datos_analisis_estrato,
        aes(color = estrato),
        title = "Análisis de correlación de estrato")

ggplotly(corr1)
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Figure 7: Gráfica de correlación de viviendas para el caso 1.

Como se puede observar en la gráfica @ref(fig:correlaciónCaso1), respecto a los diferentes estratos (\(3\) al \(6\)), se tiene diferentes tipo de relaciones con las otras variables. Por ejemplo, en general (columna 1, fila 2), el precio tiene un dispersión positiva en aumento con el área, y dependiendo del estrato (ordenado respecto a sus niveles), se ve que hay mayor dispersión, con una correlación importante de \(0,731\), obteniendo una correlación significativamente mayor en el estrato \(4\) y la más pequeña en el estrato \(6\). El precio también tiene relación con la cantidad de baños, con una correlación positiva de \(0,519\), muy superior para el caso del estrato \(3\) y la más pequeña para el estrato \(6\). Por otro lado, las habitaciones tienen una correlación débil, pero al segmentar por estratos, se tiene que la cantidad para el estrato \(6\) es muy baja. Asimismo, la cantidad de baños y habitaciones tiene correlación regular contra el área de la propiedad y los diversos estratos tiene una contribución bastante semejante para cada uno.

En la gráfica central se visualiza que en la zona Norte son pocas las viviendas del estrato \(6\), y en su mayoría pertenecen a los estrato \(3\) y \(5\), que corresponde con la contribución de cada uno en las correlaciones.

Respecto a los precios, se nota que estos aumentan con el estrato (columna 1, fila 1), pese a la distribución de densidad. Semejante sucede con el área. No obstante, con los baños y habitaciones sucede diferente, puesto que la cantidad fluctúa entre los estratos, y no está tan definida, evidente también en las gráficas de dispersión del precio y área respecto a los baños y habitaciones en las propiedades.

La relación precio vs estrato (columna 3, fila 1) sí es interesante, pero evidente, puesto que entre mayor es el estrato, mayor es el valor de las propiedades.

Paso 3:

Estime un modelo de regresión lineal múltiple con las variables del punto anterior (precio = f(área construida, estrato, número de cuartos, número de parqueaderos, número de baños )) e interprete los coeficientes si son estadísticamente significativos. Las interpretaciones deber están contextualizadas y discutir si los resultados son lógicos. Adicionalmente interprete el coeficiente R2 y discuta el ajuste del modelo e implicaciones (que podrían hacer para mejorarlo).

Para estimar el modelo de regresión lineal múltiple, debemos usar la función lm.

m1 <- lm(preciom ~ areaconst + estrato + habitaciones + parqueaderos + banios, data = datos_caso1)

summary(m1)

Call:
lm(formula = preciom ~ areaconst + estrato + habitaciones + parqueaderos + 
    banios, data = datos_caso1)

Residuals:
    Min      1Q  Median      3Q     Max 
-817.86  -69.15  -17.41   39.30 1031.59 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  186.40065   22.29308   8.361 4.42e-16 ***
areaconst      0.76051    0.04921  15.453  < 2e-16 ***
estrato.L    251.17817   21.21955  11.837  < 2e-16 ***
estrato.Q     88.35640   15.49061   5.704 1.85e-08 ***
estrato.C     51.61442   13.42554   3.844 0.000134 ***
habitaciones  -5.52483    5.02195  -1.100 0.271719    
parqueaderos  44.28853    5.80831   7.625 9.72e-14 ***
banios        13.28705    6.63276   2.003 0.045605 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 160.5 on 594 degrees of freedom
Multiple R-squared:  0.6924,    Adjusted R-squared:  0.6887 
F-statistic:   191 on 7 and 594 DF,  p-value: < 2.2e-16

El resumen del modelo nos ayuda a interpretar los coeficientes y la estimación respecto a cada variable.

Iniciaremos con el areaconst: Tiene un coeficiente de \(0.76051\), con una significación mucho menor al umbral (\(p < 0,05\)), que nos dice que hay una fuerte relación entre el área y el precio, entonces: Por cada metro cuadrado adicional de una vivienda, el precio se incrementa en \(0.76051\) millones de pesos en promedio. Lógicamente, una propiedad con mayor área tendrá un precio superior.

Para estrato, debido a que la estamos tratando como un factor con jerarquía, nos indica las variable con los sufijos L, Q y C, que representan polinomios ortogonales, que son los términos lineal (L), cuadrático (Q) y cúbico (C). Para estos valores podemos decir que el precio, al moverse entre estratos, puede ser más o menos significantes en ciertos niveles de la variable.

Para habitaciones, el coeficiente es de \(-5.52483\). Podríamos decir que por cada habitación que aumenta en una vivienda, el precio decrece en \(-5.52483\) millones de pesos, lo que no es algo lógico, pero podría tratarse por alguna vivienda (o un grupo de viviendas) que posee bastantes habitaciones, pero con valores de otras variables pequeñas, que hizo que el modelo tuviera ese comportamiento. Normalmente, la existencia de varias habitaciones debería ser directamente proporcional al precio, pero no es el caso.

Para parqueaderos, el precio aumenta \(44.28853\) por cada parqueadero adicional, lo que es una suposición bastante lógica. Los compradores, habitualmente, tienen en cuenta espacios de estacionamiento en zonas urbanas muy pobladas.

Para banios, por cada baño adicional, el precio aumenta \(13.28705\) millones de pesos, que se asemeja al comportamiento de los parqueaderos y los baños.

Para entender el comportamiento del modelo, veremos otras variables. El error residual estándar es de \(160.5\), que significa que la diferencia de precios de una vivienda estimada con una real es de \(160.5\) millones de pesos, que es un valor elevado y debe tenerse en cuenta sobre el modelo. El R2 indica que el \(69.24\:\%\) de la variabilidad del precio es explicado por estas variables, que podría sugerir la exploración de otros modelos. Además, el valor f-estadístico es grande, lo que nos dice que el modelo es globalmente significante.

Ahora haremos una estimación paso a paso.

m0 <- lm(preciom ~ 1, data = datos_caso1)

forward <- step(m0, direction='forward', scope=formula(m1), trace=0)

knitr::kable(forward$anova)
Table 6: Tabla ANOVA de estimación paso a paso.
Step Df Deviance Resid. Df Resid. Dev AIC
NA NA 601 49723449 6817.681
+ areaconst -1 26576878.92 600 23146570 6359.374
+ estrato -3 6016080.31 597 17130490 6184.180
+ parqueaderos -1 1730574.90 596 15399915 6122.068
+ banios -1 72193.17 595 15327722 6121.240

La tabla ANOVA nos muestra cómo fue el proceso. El modelo original \(m0\) tiene un AIC elevado. Al agregar la variable areaconst, decrece el AIC y una reducción significativa en el Deviance residual, que nos dice que área explica gran parte de la variabilidad de los precios. Al agregar estrato, parqueaderos y baños el modelo paró, diciéndonos que encontró que agregar otras adicionales no mejorarán el AIC.

Por lo tanto, podría explorar la opción de eliminar la variable habitaciones del modelo:

m2 <- lm(preciom ~ areaconst + estrato + parqueaderos + banios, data = datos_caso1)

summary(m2)

Call:
lm(formula = preciom ~ areaconst + estrato + parqueaderos + banios, 
    data = datos_caso1)

Residuals:
    Min      1Q  Median      3Q     Max 
-824.58  -71.05  -17.92   38.71 1019.87 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  181.01482   21.75273   8.321 5.96e-16 ***
areaconst      0.74542    0.04727  15.769  < 2e-16 ***
estrato.L    258.38100   20.18780  12.799  < 2e-16 ***
estrato.Q     88.10814   15.49171   5.687 2.02e-08 ***
estrato.C     53.06431   13.36306   3.971 8.03e-05 ***
parqueaderos  43.85715    5.79608   7.567 1.46e-13 ***
banios         9.36247    5.59272   1.674   0.0946 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 160.5 on 595 degrees of freedom
Multiple R-squared:  0.6917,    Adjusted R-squared:  0.6886 
F-statistic: 222.5 on 6 and 595 DF,  p-value: < 2.2e-16

Lamentablemente, no se aumentó el R2.

Paso 4:

La validación de los supuestos nos ayudará a interpretar de mejor manera el modelo. Para esto, usaremos una función con los supuesto más comúnes.

# Creamos una lista para almacenar los datos
lista_supuestos <- list()

# Un vector para representar la cantidad
j = 1:4

for (i in j){
  grafico <- graficos_modelo(m1, i)
  
  lista_supuestos[[i]] <- grafico
}

# Mostramos los gráficos en una cuadrícula de 2x2
grid.arrange(grobs = lista_supuestos[j], nrow = 2)
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
Supuestos del modelo 1.

Figure 8: Supuestos del modelo 1.

En la figura 8, la gráfica de \(Fitted values vs Residuals\) nos muestra que los valores de los errores no se distribuyen de manera clara (lineal, etc.) sino que en forma de embudo en el cual aumenta el radio de su cono entre más incrementan los valores ajustados, que nos indica una claro problema de heteroscedasticidad: La varianza de los residuos no es constante, confirmado por el test Breusch-Pagan, que redondeado a 3 cifras es 0. La consideración sería transformar la variable precio, tal vez con un logaritmo, para estabilizarla.

La Normalidad de los residuos nos indica que los puntos se desvían de la línea recta, mucho más en los extremos, lo que indica que los residuos no siguen una distribución normal. La parte inferior está por debajo de la línea y la superior por encima, sugiriendo una distribución sesgada en estos puntos. No podemos confiar plenamente la validez de los intervalos de confianza y las pruebas de hipótesis.

Respecto a la dispersión, parece que la línea no es un ajuste exacto para los valores debido a su altar dispersión, puesto que la nube de puntos se expande, reforzando el problema de heterocedasticidad.

Intentaremos con una transformación logarítmica:

m3 <- lm(log(preciom) ~ areaconst + estrato + habitaciones + parqueaderos + banios, data = datos_caso1)

summary(m3)

Call:
lm(formula = log(preciom) ~ areaconst + estrato + habitaciones + 
    parqueaderos + banios, data = datos_caso1)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.16480 -0.18176 -0.01318  0.17220  1.32215 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)   5.327e+00  4.049e-02 131.569  < 2e-16 ***
areaconst     1.260e-03  8.937e-05  14.096  < 2e-16 ***
estrato.L     5.512e-01  3.854e-02  14.304  < 2e-16 ***
estrato.Q    -2.721e-02  2.813e-02  -0.967   0.3339    
estrato.C     7.998e-02  2.438e-02   3.280   0.0011 ** 
habitaciones  5.814e-03  9.120e-03   0.637   0.5241    
parqueaderos  7.096e-02  1.055e-02   6.727 4.07e-11 ***
banios        5.239e-02  1.205e-02   4.349 1.61e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.2914 on 594 degrees of freedom
Multiple R-squared:  0.7579,    Adjusted R-squared:  0.755 
F-statistic: 265.6 on 7 and 594 DF,  p-value: < 2.2e-16

Con este modelo, hemos solucionado varios problemas, ahora el precio aumenta con la cantidad de habitaciones, y el error R2 es más grande, lo que nos indica que cerca del \(75.79\:\%\) de los valores son representados por este modelo. Igualmente, vemos las gráficas de supuestos:

# Creamos una lista para almacenar los datos
lista_supuestos <- list()

# Un vector para representar la cantidad
j = 1:4

for (i in j){
  grafico <- graficos_modelo(m3, i)
  
  lista_supuestos[[i]] <- grafico
}

# Mostramos los gráficos en una cuadrícula de 2x2
grid.arrange(grobs = lista_supuestos[j], nrow = 2)
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
Gráfica de supuestos del modelo logarítmico.

Figure 9: Gráfica de supuestos del modelo logarítmico.

Y, a partir de la figura 9, podemos afirmar que los residuos siguen una distribución normal. Los residuos ya no tienen esa forma de embudo, pero ahora los valores fitted versus los residuales no se acercan a la línea horizontal deseada. Sin embargo, podemos considerarlo un mejor modelo que el original, por lo que lo usaremos en adelante.

Paso 5:

Creamos un dataframe con los datos de la vivienda deseada y hacemos la predicción con el modelo logarítmico. Se debe transformar el precio predicho por la misma transformación hecha.

vivienda1_estrato4 <- data.frame(
  areaconst = 200,
  parqueaderos = 1,
  banios = 2,
  habitaciones = 4,
  estrato = factor(4, levels = levels(datos_finales$estrato))
)

precio_predicho1 <- predict(m3, newdata = vivienda1_estrato4)

# Transformamos el precio para obtener el valor real (debido a la transformación)
precio_predicho_millones1 <- exp(precio_predicho1)

print(sprintf("El precio predicho para la vivienda es: %.2f millones de pesos",  precio_predicho_millones1))
[1] "El precio predicho para la vivienda es: 305.42 millones de pesos"

Además, agregamos el supuesto de que sea una vivienda en estrato \(5\):

vivienda1_estrato5 <- data.frame(
  areaconst = 200,
  parqueaderos = 1,
  banios = 2,
  habitaciones = 4,
  estrato = factor(5, levels = levels(datos_finales$estrato))
)

precio_predicho2 <- predict(m3, newdata = vivienda1_estrato5)

precio_predicho_millones2 <- exp(precio_predicho2)


print(sprintf("El precio predicho para la vivienda es: %.2f millones de pesos",  precio_predicho_millones2))
[1] "El precio predicho para la vivienda es: 351.04 millones de pesos"

El valor de este último es ligeramente superior, pero razonable.

Paso 6:

Para ambos casos, se cumple el precio del crédito preaprobado. Entonces, hacemos una búsqueda de viviendas que cumplan con esta condición. Haremos un filtro donde consideremos todas las casas que cumplen con las condiciones y son inferiores al precio establecido.

precio_max = 350

# Eliminamos columnas no requeridas

datos_caso1 <- datos_caso1 %>% select(-zona, -tipo, -piso, -comuna)

caso1_filtrado <- datos_caso1 %>% 
  filter(
    areaconst == 200,
    parqueaderos == 1,
    banios == 2,
    habitaciones == 4,
    estrato %in% c("4", "5"),
    preciom <= precio_max
  )

knitr::kable(caso1_filtrado)
Table 7: Viviendas que corresponde a las características del caso1.
estrato preciom areaconst parqueaderos banios habitaciones longitud latitud barrio

Pero, infortunadamente, no existen viviendas que cumplan con estas condiciones. Esto es debido a que nuestro filtro es demasiado estricto, por eso haremos uno más laxo.

caso1_filtrado_rango <- datos_caso1 %>%
  filter(
    areaconst >= 200 & areaconst <= 240,
    parqueaderos >= 1 & parqueaderos <=2,
    banios >= 2 & banios <= 3, 
    habitaciones >= 4 & banios <= 5, 
    estrato %in% c("4", "5"),
    preciom <= precio_max
  )

knitr::kable(head(caso1_filtrado_rango))
Table 8: Viviendas que se acercan a los valores especificados.
estrato preciom areaconst parqueaderos banios habitaciones longitud latitud barrio
5 350 240 2 3 6 -76.53136 3.48635 El Bosque
5 350 203 2 2 5 -76.51448 3.48531 Urbanización La Flora
5 350 216 2 2 4 -76.51218 3.48181 Urbanización La Merced
5 320 210 2 3 5 -76.51200 3.47600 La Isla
5 340 203 2 3 4 -76.51803 3.48257 Vipasa
5 350 240 2 3 4 -76.51800 3.48100 Vipasa

Hacemos una mapa con estas propiedades:

mapabase %>% 
  addCircleMarkers(
    data = caso1_filtrado_rango,
    lng = ~longitud,
    lat = ~latitud,
    radius = 8,
    color = "steelblue", 
    fillOpacity = 0.8
  )

Figure 10: Mapa de viviendas que cumplen con el caso 1.

Con este filtro, y la figura 10, podemos conseguir viviendas que al menos cumplan parcialmente con los requerimiento, puesto que pueden existir propiedades con mayor número de habitaciones o de baños que puedan relacionarse respecto al precio. Las doce viviendas presentadas se ajustan a lo pedido, tomando en cuenta que podría poseer una especificación adicional de entre estas, como por ejemplo, una baño adicional o una habitación de más. Y si quisiéramos ser más laxos, podría suponer una área mayor, considerando solo como límite inferior el tope de área, y podría obtener una mayor cantidad de viviendas, pero con las \(6\) presentadas es suficiente.

Caso 2:

Realizaremos el mismo proceso para el caso 2.

Paso 1:

Primero, hacemos el filtro respectivo.

datos_caso2 <- datos_finales %>%
  filter(tipo == "Apartamento" & zona == "Sur")

head(datos_caso2)
  zona piso        tipo estrato preciom areaconst parqueaderos banios
1  Sur   01 Apartamento       5     240        87            1      3
2  Sur   01 Apartamento       5     310       137            2      3
3  Sur   03 Apartamento       4     320       108            2      3
4  Sur   05 Apartamento       6     620       135            2      3
5  Sur   06 Apartamento       5     170        56            1      1
6  Sur   03 Apartamento       4     220        75            1      2
  habitaciones  longitud latitud                 barrio comuna
1            3 -76.51700 3.36971                   Lili     17
2            4 -76.53105 3.38296             El Ingenio     17
3            3 -76.53638 3.40770 Santa Anita - La Selva     17
4            4 -76.54252 3.34299    Parcelaciones Pance     22
5            2 -76.51735 3.37251                   Lili     17
6            3 -76.51539 3.36866                   Lili     17

Y verificamos que sean viviendas de la zona sur.

mapabase %>% 
  addCircleMarkers(
    data = datos_caso2,
    lng = ~longitud,
    lat = ~latitud,
    radius = 8,
    color = "green", 
    fillOpacity = 0.8
  )

Figure 11: Filtro de viviendas para el caso 2.

Las viviendas, figura 11, pertenecen al sector requerido.

Paso 2:

Ahora vemos la correlación entre las variables.

datos_analisis_estrato2 <- datos_caso2 %>%
  select(preciom, areaconst, estrato, banios, habitaciones)

corr2 <- ggpairs(datos_analisis_estrato2,
        aes(color = estrato),
        title = "Análisis de correlación del caso 2")

ggplotly(corr2)
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Warning: Can only have one: highlight
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Warning in subplot(columnList, nrows = p$nrow, margin = 0.01, shareX = TRUE, :
Must have a consistent number of axes per 'subplot' to share them.
Warning in subplot(columnList, nrows = p$nrow, margin = 0.01, shareX = TRUE, :
Can only have one: highlight
Warning in subplot(columnList, nrows = p$nrow, margin = 0.01, shareX = TRUE, :
Can only have one: highlight

Figure 12: Correlación para el caso 2.

Como muestra la figura 12, hay pocas viviendas pertenecientes al estrato \(3\), la mayoría se encuentran ubicadas en los estratos más altos. La distribución de precio y área es bastante parecida al caso anterior, donde los precios son más grandes respecto al área de las viviendas, aunque existen viviendas del estrato \(4\) que poseen áreas grandes y precios bajos. Respecto a las gráficas de densidad, presentan formas similares, pero varían, nuevamente, en los valores para los baños y habitaciones, que fluctúa bastante, pues hay viviendas de estrato altos que poseen pocas habitaciones y pocos baños, y viceversa.

Las correlaciones sí son interesantes. El estrato \(5\) es el que menos contribuye a la correlación total del precio versus área. El estrato \(6\) es el que más influye en la correlación del precio y la cantidad de baños, pero tiene igual valor en la cantidad de baños para la correlación precio vs habitaciones con el estrato \(4\). Y la correlación de la cantidad de baños y habitaciones del estrato \(3\) es negativa y débil, pese a que es regular en general.

Paso 3:

m1_2 <- lm(preciom ~ areaconst + estrato + habitaciones + parqueaderos + banios, data = datos_caso2)

summary(m1_2)

Call:
lm(formula = preciom ~ areaconst + estrato + habitaciones + parqueaderos + 
    banios, data = datos_caso2)

Residuals:
     Min       1Q   Median       3Q      Max 
-1044.86   -37.15    -3.57    38.00   891.75 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -13.49979   12.31723  -1.096    0.273    
areaconst      1.26540    0.05697  22.213  < 2e-16 ***
estrato.L    106.26809   11.29783   9.406  < 2e-16 ***
estrato.Q     62.08623    8.01086   7.750 1.47e-14 ***
estrato.C     20.21750    4.81936   4.195 2.85e-05 ***
habitaciones  -6.65211    4.16408  -1.598    0.110    
parqueaderos  75.68897    4.50322  16.808  < 2e-16 ***
banios        47.70640    3.69892  12.897  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 94.34 on 1931 degrees of freedom
Multiple R-squared:  0.7923,    Adjusted R-squared:  0.7916 
F-statistic:  1053 on 7 and 1931 DF,  p-value: < 2.2e-16

Para este caso, el modelo inicial con todas las variables es diferente. Se mantiene ese valor negativo con las habitaciones. Pero el R2 es muchos mayor, casi un \(80\:\%\), que es muy bueno para un modelo sencillo. Podríamos probar una transformación logarítmica por si mejoraría el valor del R2.

m2_2 <- lm(log(preciom) ~ areaconst + estrato + habitaciones + parqueaderos + banios, data = datos_caso2)

summary(m2_2)

Call:
lm(formula = log(preciom) ~ areaconst + estrato + habitaciones + 
    parqueaderos + banios, data = datos_caso2)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.80897 -0.14930  0.00257  0.15335  0.81906 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  4.7023699  0.0290796 161.707  < 2e-16 ***
areaconst    0.0022470  0.0001345  16.707  < 2e-16 ***
estrato.L    0.5485675  0.0266729  20.566  < 2e-16 ***
estrato.Q    0.0243990  0.0189127   1.290 0.197176    
estrato.C    0.0231790  0.0113780   2.037 0.041768 *  
habitaciones 0.0339573  0.0098309   3.454 0.000564 ***
parqueaderos 0.1627249  0.0106316  15.306  < 2e-16 ***
banios       0.1177911  0.0087327  13.488  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.2227 on 1931 degrees of freedom
Multiple R-squared:  0.8103,    Adjusted R-squared:  0.8096 
F-statistic:  1179 on 7 and 1931 DF,  p-value: < 2.2e-16

Y, efectivamente, el R2 es un poco mayor, aunque solo en un punto porcentual.

Paso 4:

# Creamos una lista para almacenar los datos
lista_supuestos <- list()

# Un vector para representar la cantidad
j = 1:4

for (i in j){
  grafico <- graficos_modelo(m2_2, i)
  
  lista_supuestos[[i]] <- grafico
}

# Mostramos los gráficos en una cuadrícula de 2x2
grid.arrange(grobs = lista_supuestos[j], nrow = 2)
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

El modelo logarítmico para el caso 2 tiene la misma tendencia que el anterior. Los errores intentan normalizarse, los valores ajustados pretenden estabilizarse en la línea horizontal y no hay esa forma de embudo en la dispersión.

Paso 5:

Haremos las predicciones usando el segundo modelo.

vivienda2_estrato5 <- data.frame(
  areaconst = 300,
  parqueaderos = 3,
  banios = 3,
  habitaciones = 5,
  estrato = factor(5, levels = levels(datos_finales$estrato))
)

precio_predicho1_2 <- predict(m2_2, newdata = vivienda2_estrato5)

# Transformamos el precio para obtener el valor real (debido a la transformación)
precio_predicho_millones1_2 <- exp(precio_predicho1_2)

sprintf("El precio predicho para la vivienda es: %.2f millones de pesos",  precio_predicho_millones1_2)
[1] "El precio predicho para la vivienda es: 653.74 millones de pesos"
vivienda2_estrato6 <- data.frame(
  areaconst = 300,
  parqueaderos = 3,
  banios = 3,
  habitaciones = 5,
  estrato = factor(6, levels = levels(datos_finales$estrato))
)

precio_predicho2_2 <- predict(m2_2, newdata = vivienda2_estrato6)

precio_predicho_millones2_2 <- exp(precio_predicho2_2)


sprintf("El precio predicho para la vivienda es: %.2f millones de pesos",  precio_predicho_millones2_2)
[1] "El precio predicho para la vivienda es: 874.07 millones de pesos"

Para el caso 2, la vivienda con estas características y un estrato 6 sobrepasa el crédito preaprobado por \(24\) millones de pesos.

Paso 6:

precio_max2 = 850 

# Eliminamos columnas no requeridas
datos_caso2 <- datos_caso2 %>% select(-zona, -tipo, -piso, -comuna)

caso2_filtrado <- datos_caso2 %>% 
  filter(
    areaconst == 300,
    parqueaderos == 3,
    banios == 3,
    habitaciones == 5,
    estrato %in% c("5", "6"),
    preciom <= precio_max2
  )

knitr::kable(caso2_filtrado)
Table 9: Viviendas con especificaciones del caso 2.
estrato preciom areaconst parqueaderos banios habitaciones longitud latitud barrio

Nuevamente, un filtro estricto no es adecuado, por lo que haremos el mismo tratamiento anterior:

caso2_filtrado_rango <- datos_caso2 %>%
  filter(
    areaconst >= 300 & areaconst <= 340,
    parqueaderos >= 3 & parqueaderos <= 4,
    banios >= 3 & banios <= 4, 
    habitaciones >= 5 & banios <= 6, 
    estrato %in% c("5", "6"),
    preciom <= precio_max2
  )

knitr::kable(head(caso2_filtrado_rango))
Table 10: Viviendas con un filtro menos estricto para el caso 2.
estrato preciom areaconst parqueaderos banios habitaciones longitud latitud barrio

E, incluso con un filtro más laxo, no se pueden ubicar viviendas adecuadas a los requerimientos. Pero, la predicción no dice que no encontraremos viviendas con ese crédito preaprobado para estrato \(6\), entonces, lo quitaremos del filtro.

caso2_filtrado_rango2 <- datos_caso2 %>%
  filter(
    areaconst >= 300,
    parqueaderos >= 3,
    banios >= 3, 
    habitaciones >= 5, 
    estrato == "5",
    preciom <= precio_max2
  )

knitr::kable(head(caso2_filtrado_rango2))
Table 11: Viviendas solo de estrato 5.
estrato preciom areaconst parqueaderos banios habitaciones longitud latitud barrio

Sin embargo, el filtro que no tiene en cuenta límites superiores, pero sí limites inferiores establecidos, tampoco encuentran viviendas que se ajusten a los precios, pese a que el modelo sí pudo predecir su precio. Un último filtro sugerido sería establecer el precio, el área y el estrato como datos importantes, y los demás dejarlos libres.

caso2_filtrado_rango3 <- datos_caso2 %>%
  filter(
    areaconst >= 300,
    estrato %in% c("5", "6"),
    preciom <= precio_max2
  )

knitr::kable(head(caso2_filtrado_rango3))
Table 12: Viviendas para el caso 2 sin tener en cuenta amenidades.
estrato preciom areaconst parqueaderos banios habitaciones longitud latitud barrio
6 620 480 3 5 5 -76.52662 3.36322 Urbanización Rio Lili
5 690 486 2 4 4 -76.53111 3.38292 El Ingenio
5 650 600 2 4 5 -76.53400 3.38100 El Ingenio
5 170 605 1 2 2 -76.54294 3.39992 El Gran Limonar
6 850 325 2 3 2 -76.54200 3.37500 La Playa
6 850 352 4 3 3 -76.53729 3.34265 Parcelaciones Pance

Esto se da, obviamente, al comportamiento de las variables habitaciones, baños y parqueaderos en las viviendas, puesto que se podría encontrar una vivienda con grandes dimensiones, en los estratos requeridos y cumpliendo el precio, pero que no posee estás amenidades. La sugerencia sería presentarle al cliente las viviendas encontradas y podría sugerirse la leve remodelación a las características requeridas si se encuentra una viviendas con un precio inferior al establecido. Por ejemplo, de la tabla 12 se puede remarcar la vivienda en la primera fila, que posee un precio de \(620\) millones de pesos, con un área mucho mayor a la requerida, de \(480\: m^{2}\), con una cantidad de parqueaderos y habitaciones exacta, pero con dos baños adicionales a lo requrido. Con la diferencia de precio al crédito preaprobado se puede sugerir las remodelaciones requeridas.

Finalmente, las viviendas encontradas son:

mapabase %>% 
  addCircleMarkers(
    data = caso2_filtrado_rango3,
    lng = ~longitud,
    lat = ~latitud,
    radius = 8,
    color = "#01ef63", 
    fillOpacity = 0.8
  )

Figure 13: Ubicación viviendas para el caso 2.