actividad2

#install.packages("devtools") # solo la primera vez
#devtools::install_github("centromagis/paqueteMODELOS", force =TRUE)
# Instalar el paquete "model" desde GitHub (este paquete permite generar tablas ANOVA de manera simplificada)
#remotes::install_github("fhernanb/model")
library(paqueteMODELOS)
data("vivienda")

library(paqueteMODELOS)
data("vivienda")
#install.packages("leaflet")
#install.packages("plotly")   

library(dplyr)
library(tidyr)
library(plotly)
library(ggplot2)  
library(lmtest)
library(leaflet)

# Export to CSV in the working directory
write.csv(vivienda, file = "vivienda.csv", row.names = FALSE)

# Crear el dataframe necesario
df <- as.data.frame(vivienda)

Parte 1: Filtrado de base de datos

# Filtro del dataframe por la zona norte y tipo de vivienda Casa
df1 <- df %>%
  dplyr::filter(tipo == "Casa" & zona == "Zona Norte")

# Gráficar mapa de los datos filtrados
library(leaflet)

map <- leaflet(df1) %>%
  addTiles() %>%
  addMarkers(
    lng = ~longitud,
    lat = ~latitud,
    popup = ~as.character(latitud) # Puedes personalizar el contenido del popup aquí
  )

map 

Al aplicar el filtro a la base de datos para incluir únicamente casas ubicadas en la zona norte, se observa que no todas las ofertas cumplen con este criterio. Algunos registros aparecen fuera de la zona correspondiente, que puede deberse a inconsistencias en la digitación de los datos de las coordenadas geográficas. Estos casos pueden deberse a errores humanos en el ingreso de la información o a imprecisiones en el levantamiento de los datos espaciales.

# Definir una longitud y latitud maxima
lat_min_nor = 3.450051
#lat_max_nor=3.50
#lon_min_nor=-76.55
lon_max_nor=-76.537457

# Nuevo fitro de zonas norte
df_norte <- df1 %>%
  filter(longitud >= lon_max_nor & latitud >= lat_min_nor)

# Gráficar el mapa para validar resultados
map_norte <- leaflet(df_norte) %>%
  addTiles() %>%
  addMarkers(
    lng = ~longitud,
    lat = ~latitud,
    popup = ~as.character(latitud) # Puedes personalizar el contenido del popup aquí
  )

map_norte 

Para mitigar este issue de calidad se aplican filtros por latitud y longitud para crear un perimetro que permita filtrar el dataset correctamente

Parte 2: Análisis exploratorio encofado en correlación

Se realiza un análisis senicillo para la identificación de promedios, máximos, mínimos y los missing values

summary(df_norte[, !(names(df_norte) %in% c("latitud", "longitud", "id","zona","piso"))])

##     estrato        preciom         areaconst     parqueaderos   
##  Min.   :3.00   Min.   : 110.0   Min.   :  30   Min.   : 1.000  
##  1st Qu.:3.00   1st Qu.: 241.5   1st Qu.: 135   1st Qu.: 1.000  
##  Median :4.00   Median : 370.0   Median : 232   Median : 2.000  
##  Mean   :4.16   Mean   : 421.8   Mean   : 254   Mean   : 2.211  
##  3rd Qu.:5.00   3rd Qu.: 530.0   3rd Qu.: 330   3rd Qu.: 3.000  
##  Max.   :6.00   Max.   :1940.0   Max.   :1440   Max.   :10.000  
##                                                 NA's   :190     
##      banios        habitaciones        tipo              barrio         
##  Min.   : 0.000   Min.   : 0.000   Length:555         Length:555        
##  1st Qu.: 2.000   1st Qu.: 3.000   Class :character   Class :character  
##  Median : 3.000   Median : 4.000   Mode  :character   Mode  :character  
##  Mean   : 3.506   Mean   : 4.555                                        
##  3rd Qu.: 4.000   3rd Qu.: 5.000                                        
##  Max.   :10.000   Max.   :10.000                                        
## 

Limpieza del set de datos

# extraer la variable id
df_mod1 <- df_norte[ , !(names(df_norte) %in% "id")]

# cambiar el tipo de variable de piso
df_mod1$piso <- as.numeric(df_mod1$piso)

# identficiación de missing values
df_mod1 %>%
  summarise(across(everything(), ~ sum(is.na(.)))) %>%
  pivot_longer(cols = everything(),
               names_to = "column",
               values_to = "na_count")

## # A tibble: 12 × 2
##    column       na_count
##    <chr>           <int>
##  1 zona                0
##  2 piso              251
##  3 estrato             0
##  4 preciom             0
##  5 areaconst           0
##  6 parqueaderos      190
##  7 banios              0
##  8 habitaciones        0
##  9 tipo                0
## 10 barrio              0
## 11 longitud            0
## 12 latitud             0

df_mod1$piso[is.na(df_mod1$piso)] <- median(df_mod1$piso, na.rm = TRUE)
df_mod1$parqueaderos[is.na(df_mod1$parqueaderos)] <- median(df_mod1$parqueaderos, na.rm = TRUE)

Outliers

Análisis por estratos

ggplot(df_mod1, aes(x = factor(estrato), y = preciom)) +
  geom_boxplot(fill = "skyblue", color = "darkblue") +
  labs(title = "Distribucion de precios por estrato",
       x = "Estrato",
       y = "Precio (preciom)") +
  theme_minimal()

Estrato 3: no es muy frecuente encontrar precios de vivienda por más de 500 millones en estratos 3, sin embargo, dado que el norte es una zona que integra zonas residenciales y comerciales, si es posible encontrar precios un poco más elevados del promedio (230M)

mean_precio_estrato3 <- mean(df_mod1$preciom[df_mod1$estrato == 3], na.rm = TRUE)
mean_precio_estrato3

## [1] 229.1505

Para los demás estratos, hace un poco de ruido encontrar cifras tan elevadas, en especial mayores a los 1000M, por ende se filtrar estos registros y se evidencian que están ubicados en barrios que tienen segmentos lujosos, como lo son los barrios Granada y Santa Mónica, por lo tanto se decide conservar todos los puntos para el modelado

preciom_outlier = 1000

df_mod1_preciom_1000 <- df_mod1 %>%
  filter(preciom >= preciom_outlier) %>%
  select(-longitud, -latitud,-zona,-tipo)

print(df_mod1_preciom_1000)

##    piso estrato preciom areaconst parqueaderos banios habitaciones
## 1     2       4    1800     607.0            2      4            8
## 2     2       5    1100     500.0            4      5            5
## 3     1       5    1500     400.0            2      3            4
## 4     2       6    1250     330.0            6      5            4
## 5     2       4    1650     734.0            2      5           10
## 6     2       5    1940     734.0            3      8           10
## 7     2       6    1300     552.0            7      9            9
## 8     2       5    1000     350.0            2      3            3
## 9     2       5    1270     950.0            4      5           10
## 10    2       6    1500     470.0            3      6            5
## 11    1       5    1600     942.0            4      4           10
## 12    2       5    1200     523.3            2      4            7
## 13    2       5    1200     730.0            4      6            6
## 14    2       5    1400     265.0            2     10           10
##                      barrio
## 1                   granada
## 2                   granada
## 3                     menga
## 4                     menga
## 5               san vicente
## 6               san vicente
## 7  santa mónica residencial
## 8  santa mónica residencial
## 9              santa monica
## 10             santa monica
## 11                versalles
## 12                versalles
## 13                versalles
## 14                   vipasa

Análisis por variables

ggplot(df_mod1, aes(y = preciom)) +
  geom_boxplot(fill = "skyblue", color = "darkblue") +
  labs(title = "Boxplot de preciom",
       y = "Precio (preciom)",
       x = "") +
  theme_minimal()

df_long <- df_mod1 %>%
  select(banios, habitaciones) %>%
  pivot_longer(cols = everything(), names_to = "variable", values_to = "value")

ggplot(df_long, aes(x = variable, y = value, fill = variable)) +
  geom_boxplot(color = "darkblue") +
  labs(title = "Boxplots de Banios y Habitaciones",
       x = "Variable",
       y = "Valor") +
  theme_minimal() +
  theme(legend.position = "none")

Se identifican outliers sin embargo dado que el dataset tiene en cuenta diferentes estratos, si es posible en la ciudad de cali, encontrar viviendas con esas caracteristicas. Lo que si no tiene mucho sentido, es encontrar viviendas con 0 habitaciones o baños, por lo cual se reemplazarán esos valores por su respectiva media

# mean para banios
mean_banios <- round(mean(df_mod1$banios[df_mod1$banios != 0], na.rm = TRUE))
df_mod1$banios[df_mod1$banios == 0] <- mean_banios

# mean para habitaciones
mean_habit <- round(mean(df_mod1$habitaciones[df_mod1$habitaciones != 0], na.rm = TRUE))
df_mod1$habitaciones[df_mod1$habitaciones == 0] <- mean_habit

Correlación para variables Spearman

# correlación de columnas específicas
cor_num <- cor(df_mod1[,c( 'preciom','areaconst', 'banios', 'habitaciones', 'estrato')], use = "complete.obs", method = "spearman")

# Visualización de la matriz 
fig <- plot_ly(
  z = cor_num,
  x = colnames(cor_num),
  y = colnames(cor_num),
  type = "heatmap",
  text = round(cor_num, 2),   # labels
  texttemplate = "%{text}",   # show them
  textfont = list(size = 10, color = "black") # label style
) %>%
  layout(
    width = 500,  # smaller width
    height = 500  # smaller height
  )

fig

Paso 3: Modelo de RLM

# Dado que el estrato es una variable ordinal se aplica la función as factor 
df_mod1$estrato <- as.factor(df_mod1$estrato)

# definición del modelo 1
mod1_original =lm(preciom ~ areaconst + parqueaderos + banios + habitaciones + estrato, data=df_mod1)

summary(mod1_original)

## 
## Call:
## lm(formula = preciom ~ areaconst + parqueaderos + banios + habitaciones + 
##     estrato, data = df_mod1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -863.30  -65.42  -15.24   40.68 1078.87 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -16.55422   21.55974  -0.768 0.442919    
## areaconst      0.72332    0.05341  13.542  < 2e-16 ***
## parqueaderos  20.06497    5.92252   3.388 0.000755 ***
## banios        18.87965    6.55244   2.881 0.004116 ** 
## habitaciones  11.26891    5.26620   2.140 0.032808 *  
## estrato4      92.82873   19.21860   4.830 1.77e-06 ***
## estrato5     145.00066   18.37962   7.889 1.66e-14 ***
## estrato6     304.52783   34.83374   8.742  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 152.2 on 547 degrees of freedom
## Multiple R-squared:  0.6454, Adjusted R-squared:  0.6409 
## F-statistic: 142.3 on 7 and 547 DF,  p-value: < 2.2e-16

Análisis de coeficientes:

De acuerdo con los resultados de mod1, todos los coeficientes tiene un efecto significativo en el precio con significancias del 1% y 10% a excepción de la variable habitaciones que presenta una significancia del 5%

En cuanto al R-squared el modelo está explicando aproximadamente el 64.5% de la variación de los precios de las viviendas. Asimismo, existe un 35.5% sin explicar que puede ser debido a factores que afectan el precio que no están en el modelo como la antigüedad de la vivienda, calidad de acabados entre otros

# Cargar el paquete "model"
library(model)

# Generar la tabla ANOVA para evaluar la significancia de la regresión
anova_table_lm(mod1_original)  # Función del paquete "model" para mostrar la tabla ANOVA del modelo

## Anova Table
##              Sum Sq  Df Mean Sq F value    Pr(>F)    
## Regression 23056440   7 3293777  142.25 < 2.2e-16 ***
## Residuals  12665649 547   23155                      
## Total      35722089 554                              
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

De la tabla ANOVA presentada anteriormente, se observa que el valor−p es 2.2×10−16, lo que es extremadamente pequeño. Esto proporciona evidencia suficiente para decir que las variables tienen un efecto estadísticamente significativo en la explicación con un nivel de significancia del 1%

Aplicando logaritmo natural para mejorar el manejo de las diferencias de rangos que tienen las variables como es el caso del área construida, se obtiene que de esta manera el modelo explicaría aproximadamente el 75.8% de la variación de los precios de las viviendas, mejorando 11.3 puntos

# Definición del modelo log-transformado
mod1 <- lm(log(preciom) ~ areaconst + parqueaderos + banios +
                                 habitaciones + estrato, data = df_mod1)

# Resumen del modelo
summary(mod1)

## 
## Call:
## lm(formula = log(preciom) ~ areaconst + parqueaderos + banios + 
##     habitaciones + estrato, data = df_mod1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.40057 -0.15722 -0.01371  0.15387  1.08376 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.816e+00  3.922e-02 122.790  < 2e-16 ***
## areaconst    1.345e-03  9.716e-05  13.847  < 2e-16 ***
## parqueaderos 2.767e-02  1.077e-02   2.568  0.01049 *  
## banios       5.747e-02  1.192e-02   4.821 1.85e-06 ***
## habitaciones 3.036e-02  9.579e-03   3.170  0.00161 ** 
## estrato4     3.771e-01  3.496e-02  10.787  < 2e-16 ***
## estrato5     5.266e-01  3.343e-02  15.750  < 2e-16 ***
## estrato6     7.469e-01  6.336e-02  11.787  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2768 on 547 degrees of freedom
## Multiple R-squared:  0.7587, Adjusted R-squared:  0.7556 
## F-statistic: 245.7 on 7 and 547 DF,  p-value: < 2.2e-16

Parte 4: Validación de los supuestos

1. Los errores siguen una distribución normal

La figura muestra que los puntos del gráfico de normalidad de los residuales estandarizados presentan desviaciones respecto a la línea de referencia. Este comportamiento sugiere que los errores no se distribuyen de manera aproximadamente normal, lo cual constituye una evidencia en contra del supuesto de normalidad

# Obtener residuales estandarizados
residuales_estandarizados1 <- rstandard(mod1)

# Generar el QQ-Plot para evaluar normalidad
qqnorm(residuales_estandarizados1, main="QQ-Plot de los Residuales Estandarizados", col="blue")
qqline(residuales_estandarizados1, col="red", lwd=2)  # Agregar línea de referencia

Para complementar el análisis gráfico se realiza el siguiente Test de normalidad aplicando la prueba de Shapiro-Wilk

H0: Los errores siguen una distribución normal H1: Los errores no siguen una distribución normal

# Determinar residuos ordinarios y aplicar test de normalidad
ei <- residuals(mod1)
shapiro.test(ei)

## 
##  Shapiro-Wilk normality test
## 
## data:  ei
## W = 0.97873, p-value = 3.188e-07

De la salida anterior se aprecia que el p-value = 3.188x10-07 es menor que un nivel de significancia del 0.001, eso significa que hay fuerte evidencia para rechazar la hipótesis de normalidad de los errores.

2. Los errores tienen varianza constante (homocedasticidad):

Se realiza prueba de Breusch-Pagan para comprobar que los errores tienen varianza constante. Dado P-Value 2.197x10-15 < 0.05 se asume que los residuos NO tienen varianza constante con una significancia del 5% y hay presencia de heterosticidad

H0: Los errores tienen varianza constante H1: Los errores no tienen varianza constante

# Aplicar la prueba de Breusch-Pagan directamente con lmtest
bptest(mod1)

## 
##  studentized Breusch-Pagan test
## 
## data:  mod1
## BP = 83.903, df = 7, p-value = 2.197e-15

# Se obtienen los residuales del modelo
ei <- resid(mod1)  # Residuales ordinarios
y_ajustados <- fitted(mod1)  # Valores ajustados (predichos)


# Visualización de la Homocedasticidad
# Gráfico de Residuales vs Valores Ajustados 
plot1.bt<-ggplot(data = df_mod1, aes(x = y_ajustados, y = ei)) +
  geom_point(color = "blue") +  # Puntos de residuales
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +  # Línea horizontal en cero
  labs(title = "Residuales vs Valores Ajustados",
       x = "Valores Ajustados",
       y = "Residuales") +
  theme_minimal()

#Gráfico de Raíz de Valor Absoluto de Residuales vs Valores Ajustados 
plot2.bt<-ggplot(data = df_mod1, aes(x = y_ajustados, y = sqrt(abs(ei)))) +
  geom_point(color = "blue") +  # Puntos
  geom_smooth(method = "lm", se = FALSE, color = "red") +  # Línea de tendencia
  labs(title = "Raiz Cuadrada del Valor Absoluto de Residuales vs Valores Ajustados",
       x = "Valores Ajustados",
       y = "√|Residuales|") +
  theme_minimal()

print(plot1.bt)

print(plot2.bt)

## `geom_smooth()` using formula = 'y ~ x'

3. Pruebas de independencia

Para verificar el supuesto de independencia de los errores, se aplican pruebas estadísticas diseñadas para detectar la presencia de autocorrelación en los residuos del modelo de regresión.

Las hipótesis de estas pruebas son:

H0: Los errores son independientes, es decir, no presentan autocorrelación H1: Los errores no son independientes, es decir, presentan autocorrelación

# Pruebas de independencia de los errores 

# Test de Durbin-Watson
# Evalúa si los errores presentan autocorrelación de primer orden
dw_test <- dwtest(mod1)
print(dw_test)

## 
##  Durbin-Watson test
## 
## data:  mod1
## DW = 1.5326, p-value = 1.1e-08
## alternative hypothesis: true autocorrelation is greater than 0

Con un p-value tan pequeño (0.000000011), se rechaza la hipótesis de independencia. Es decir, los residuos del modelo presentan autocorrelación significativa

# Visualización de los residuos 
# Se crea un data frame con los residuos y su índice de observación
resid_data1 <- data.frame(
  Observacion = seq_along(residuals(mod1)),  # Índice de cada observación
  Residuales = residuals(mod1)              # Valores de los residuos
)

# Gráfico de residuos usando ggplot2
ggplot(resid_data1, aes(x = Observacion, y = Residuales)) +
  geom_point(color = "deeppink", size = 2) +       # Puntos de los residuos
  geom_line(color = "grey", linetype = "dashed") + # Línea guía de los residuos
  geom_hline(yintercept = 0, color = "black", linetype = "solid") + # Línea en y = 0
  labs(title = "Residuos del Modelo vs Index de Observacion",
       x = "Index de la Observacion", 
       y = "Residuales") +
  theme_minimal()  # Estilo de gráfico limpio

Parte 5: Predecir precios

Intervalo de confianza

# Crear un nuevo conjunto de datos con los valores específicos de Edad para calcular el IC

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

# Calcular el intervalo de confianza al 95% para la media de Resistencia en los valores seleccionados
ic_result1 <- predict(object = mod1, newdata = new_dt1, interval = "confidence", level = 0.95)

# Convertir a data.frame para manipularlo
ic_result1 <- as.data.frame(ic_result1)

# Back-transformar (exp) para obtener resultados en escala original
ic_result1$fit      <- exp(ic_result1$fit)
ic_result1$lwr      <- exp(ic_result1$lwr)
ic_result1$upr      <- exp(ic_result1$upr)

# Mostrar resultados en escala original
ic_result1

##        fit      lwr      upr
## 1 306.6729 288.2714 326.2492
## 2 356.1157 336.3013 377.0975

Para una casa de 200 metros cuadrados estrato 4, que tenga 1 parqueadero, 2 baños y 4 habitaciones en la zona norte de Cali el intervalo de confianza del 95% es (288.2714, 326.2492).

Para una casa de 200 metros cuadrados estrato 5, que tenga 1 parqueadero, 2 baños y 4 habitaciones en la zona norte de Cali el intervalo de confianza del 95% es (336.3013, 377.0975).

De acuerdo con el MRLM el valor de del precio de la vivienda con las características no.1 son:

# Predictions in the log scale
pred_log <- predict(mod1, newdata = new_dt1)

# Back-transform to original scale
pred_original <- exp(pred_log)

pred_original

##        1        2 
## 306.6729 356.1157

Parte 6: sugiera potenciales ofertas

De acuerdo con las predicciones del modelo se sugieren las siguientes viviendas:

# Crear dataframe base de resultados
df_mod1_results <- df_mod1 %>%
  mutate(
    # Back-transform para obtener precios en escala original
    preciom_pred = exp(predict(mod1, newdata = df_mod1))
  )

# Solicitudes de la búsqueda 1
solicitud_1 <- df_mod1_results %>%
  filter(
    preciom_pred <= 350,
    estrato %in% c(4, 5),
    areaconst >= 200, 
    habitaciones >= 4, 
    banios >= 2, 
    parqueaderos >= 1
  )

# Tabla de resultados
kable(
  solicitud_1 %>%
    select(zona, tipo, estrato, areaconst, habitaciones, banios, parqueaderos,
           preciom, preciom_pred, longitud, latitud)
)

zona	tipo	estrato	areaconst	habitaciones	banios	parqueaderos	preciom	preciom_pred	longitud	latitud
Zona Norte	Casa	4	250	4	3	1	485	347.4139	-76.53090	3.48532
Zona Norte	Casa	4	225	4	3	2	430	345.3467	-76.51949	3.47946
Zona Norte	Casa	4	216	5	2	2	370	332.0676	-76.51273	3.48214
Zona Norte	Casa	4	216	4	2	2	360	322.1362	-76.51390	3.48386
Zona Norte	Casa	4	265	4	2	2	550	344.0883	-76.51777	3.48060

# Mapa de los puntos potenciales
map_r1 <- leaflet(solicitud_1) %>%
  addTiles() %>%
  addMarkers(
    lng   = ~longitud,
    lat   = ~latitud,
    popup = ~paste0(
      "<b>Zona:</b> ", zona, "<br>",
      "<b>Tipo:</b> ", tipo, "<br>",
      "<b>Estrato:</b> ", estrato, "<br>",
      "<b>Área:</b> ", areaconst, " m²<br>",
      "<b>Habitaciones:</b> ", habitaciones, "<br>",
      "<b>Baños:</b> ", banios, "<br>",
      "<b>Parqueaderos:</b> ", parqueaderos, "<br>",
      "<b>Precio real:</b> ", round(preciom, 0), "<br>",
      "<b>Precio predicho:</b> ", round(preciom_pred, 0)
    )
  )

map_r1

Parte 7: Modelo para opción 2

Filtro por zona y vivienda

# Filtro del dataframe por la zona norte y tipo de vivienda Casa
df_sur <- df %>%
  dplyr::filter(tipo == "Apartamento" & zona == "Zona Sur")

# Mapa de los puntos filtrados
map_sur <- leaflet(df_sur) %>%
  addTiles() %>%
  addMarkers(
    lng = ~longitud,
    lat = ~latitud,
    popup = ~as.character(latitud) # Puedes personalizar el contenido del popup aquí
  )

map_sur

# Definir una longitud y latitud maxima
lat_min_sur = 3.416638
#lat_max_sur=3.50
#lon_min_sur=-76.55
lon_max_sur=-76.539946

# Nuevo fitro de zonas sur
df_mod2 <- df_sur %>%
  filter(longitud >= lon_max_sur & latitud <= lat_min_sur)

# Gráficar el mapa para validar resultados
map_mod2 <- leaflet(df_mod2) %>%
  addTiles() %>%
  addMarkers(
    lng = ~longitud,
    lat = ~latitud,
    popup = ~as.character(latitud) # Puedes personalizar el contenido del popup aquí
  )

map_mod2

summary(df_mod2[, !(names(df_mod2) %in% c("latitud", "longitud", "id","zona","piso"))])

##     estrato         preciom         areaconst       parqueaderos   
##  Min.   :3.000   Min.   :  78.0   Min.   : 40.00   Min.   : 1.000  
##  1st Qu.:4.000   1st Qu.: 175.0   1st Qu.: 65.00   1st Qu.: 1.000  
##  Median :5.000   Median : 250.0   Median : 82.00   Median : 1.000  
##  Mean   :4.681   Mean   : 292.9   Mean   : 92.81   Mean   : 1.378  
##  3rd Qu.:5.000   3rd Qu.: 330.0   3rd Qu.:105.00   3rd Qu.: 2.000  
##  Max.   :6.000   Max.   :1750.0   Max.   :600.00   Max.   :10.000  
##                                                    NA's   :260     
##      banios       habitaciones       tipo              barrio         
##  Min.   :0.000   Min.   :0.000   Length:1775        Length:1775       
##  1st Qu.:2.000   1st Qu.:3.000   Class :character   Class :character  
##  Median :2.000   Median :3.000   Mode  :character   Mode  :character  
##  Mean   :2.439   Mean   :2.915                                        
##  3rd Qu.:3.000   3rd Qu.:3.000                                        
##  Max.   :7.000   Max.   :5.000                                        
## 

Limpieza del set de datos

# extraer la variable id
df_mod2 <- df_mod2[ , !(names(df_mod2) %in% "id")]

# cambiar el tipo de variable de piso
df_mod2$piso <- as.numeric(df_mod2$piso)

# identficiación de missing values
df_mod2 %>%
  summarise(across(everything(), ~ sum(is.na(.)))) %>%
  pivot_longer(cols = everything(),
               names_to = "column",
               values_to = "na_count")

## # A tibble: 12 × 2
##    column       na_count
##    <chr>           <int>
##  1 zona                0
##  2 piso              397
##  3 estrato             0
##  4 preciom             0
##  5 areaconst           0
##  6 parqueaderos      260
##  7 banios              0
##  8 habitaciones        0
##  9 tipo                0
## 10 barrio              0
## 11 longitud            0
## 12 latitud             0

#Ajuste de valores 0
df_mod2$piso[is.na(df_mod2$piso)] <- median(df_mod2$piso, na.rm = TRUE)
df_mod2$parqueaderos[is.na(df_mod2$parqueaderos)] <- median(df_mod2$parqueaderos, na.rm = TRUE)

Outliers

Análisis por estratos

ggplot(df_mod2, aes(x = factor(estrato), y = preciom)) +
  geom_boxplot(fill = "skyblue", color = "darkblue") +
  labs(title = "Distribucion de precios por estrato",
       x = "Estrato",
       y = "Precio (preciom)") +
  theme_minimal()

Análisis por variables

ggplot(df_mod2, aes(y = preciom)) +
  geom_boxplot(fill = "skyblue", color = "darkblue") +
  labs(title = "Boxplot de preciom",
       y = "Precio (preciom)",
       x = "") +
  theme_minimal()

# mean para banios
mean_banios2 <- round(mean(df_mod2$banios[df_mod2$banios != 0], na.rm = TRUE))
df_mod2$banios[df_mod2$banios == 0] <- mean_banios2

# mean para habitaciones
mean_habit2 <- round(mean(df_mod2$habitaciones[df_mod2$habitaciones != 0], na.rm = TRUE))
df_mod2$habitaciones[df_mod2$habitaciones == 0] <- mean_habit2

df_long2 <- df_mod2 %>%
  select(banios, habitaciones) %>%
  pivot_longer(cols = everything(), names_to = "variable", values_to = "value")

ggplot(df_long2, aes(x = variable, y = value, fill = variable)) +
  geom_boxplot(color = "darkblue") +
  labs(title = "Boxplots de Banios y Habitaciones",
       x = "Variable",
       y = "Valor") +
  theme_minimal() +
  theme(legend.position = "none")

Correlación para variables Spearman

# correlación de columnas específicas
cor_num2 <- cor(df_mod2[,c( 'preciom','areaconst', 'banios', 'habitaciones', 'estrato')], use = "complete.obs", method = "spearman")

# Visualización de la matriz 
fig2 <- plot_ly(
  z = cor_num,
  x = colnames(cor_num),
  y = colnames(cor_num),
  type = "heatmap",
  text = round(cor_num2, 2),   # labels
  texttemplate = "%{text}",   # show them
  textfont = list(size = 10, color = "black") # label style
) %>%
  layout(
    width = 500,  # smaller width
    height = 500  # smaller height
  )

fig2

Se evidencia una muy fuerte correlación entre el preciom y el área construida que podría estar sesgando un poco el modelo

# Dado que el estrato es una variable ordinal se aplica la función as factor 
df_mod2$estrato <- as.factor(df_mod2$estrato)

# definición del modelo 1
mod2_normal=lm(preciom ~ areaconst + parqueaderos + banios + habitaciones + estrato, data=df_mod2)

summary(mod2_normal)

## 
## Call:
## lm(formula = preciom ~ areaconst + parqueaderos + banios + habitaciones + 
##     estrato, data = df_mod2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -737.08  -29.76   -0.43   30.59  706.43 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -41.04137   12.68643  -3.235  0.00124 ** 
## areaconst      2.07346    0.06472  32.037  < 2e-16 ***
## parqueaderos  49.92498    4.03998  12.358  < 2e-16 ***
## banios        28.70436    3.22067   8.913  < 2e-16 ***
## habitaciones -15.10980    3.65172  -4.138 3.67e-05 ***
## estrato4      19.27642    8.96302   2.151  0.03164 *  
## estrato5      44.92571    9.19777   4.884 1.13e-06 ***
## estrato6     142.62903   10.88550  13.103  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 71.93 on 1767 degrees of freedom
## Multiple R-squared:  0.8387, Adjusted R-squared:  0.838 
## F-statistic:  1312 on 7 and 1767 DF,  p-value: < 2.2e-16

Análisis de coeficientes:

De acuerdo con los resultados de mod2_normal, todos los coeficientes tiene un efecto significativo en el precio con significancias del 1% y 10% a excepción de la variable estrato nivel 4 que presenta una significancia del 5%

En cuanto al R-squared el modelo está explicando aproximadamente el 83.87% de la variación de los precios de las viviendas. Asimismo, existe un 16.13% sin explicar que puede ser debido a factores que afectan el precio que no están en el modelo como la antigüedad de la vivienda, calidad de acabados entre otros

# Cargar el paquete "model"
library(model)

# Generar la tabla ANOVA para evaluar la significancia de la regresión
anova_table_lm(mod2_normal)  # Función del paquete "model" para mostrar la tabla ANOVA del modelo

## Anova Table
##              Sum Sq   Df Mean Sq F value    Pr(>F)    
## Regression 47528078    7 6789725  1312.2 < 2.2e-16 ***
## Residuals   9142849 1767    5174                      
## Total      56670927 1774                              
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

De la tabla ANOVA presentada anteriormente, se observa que el valor−p es 2.2×10−16, lo que es extremadamente pequeño. Esto proporciona evidencia suficiente para decir que las variables tienen un efecto estadísticamente significativo en la explicación con un nivel de significancia del 1%

Aplicando logaritmo natural para mejorar el manejo de las diferencias de escalas que tienen las variables como es el caso del área construida, se obtiene que de esta manera el modelo explicaría aproximadamente el 82.63% de la variación de los precios de las viviendas, dismiuyendo 0.55 puntos

# Definición del modelo log-transformado
mod2 <- lm(log(preciom) ~ areaconst + parqueaderos + banios +
                                 habitaciones + estrato, data = df_mod2)

# Resumen del modelo
summary(mod2)

## 
## Call:
## lm(formula = log(preciom) ~ areaconst + parqueaderos + banios + 
##     habitaciones + estrato, data = df_mod2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.44622 -0.13455  0.00374  0.13749  0.97168 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.2970632  0.0357575 120.172  < 2e-16 ***
## areaconst    0.0040025  0.0001824  21.941  < 2e-16 ***
## parqueaderos 0.1114284  0.0113869   9.786  < 2e-16 ***
## banios       0.0870549  0.0090777   9.590  < 2e-16 ***
## habitaciones 0.0270809  0.0102926   2.631  0.00858 ** 
## estrato4     0.2854709  0.0252628  11.300  < 2e-16 ***
## estrato5     0.5181714  0.0259245  19.988  < 2e-16 ***
## estrato6     0.7459445  0.0306814  24.313  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2027 on 1767 degrees of freedom
## Multiple R-squared:  0.8269, Adjusted R-squared:  0.8263 
## F-statistic:  1206 on 7 and 1767 DF,  p-value: < 2.2e-16

Validación de los supuestos

1. Los errores siguen una distribución normal

La figura muestra que los puntos del gráfico de normalidad de los residuales estandarizados presentan desviaciones respecto a la línea de referencia en especial en la parte inicial. Este comportamiento sugiere que los errores no se distribuyen de manera aproximadamente normal

# Obtener residuales estandarizados
residuales_estandarizados2 <- rstandard(mod2)

# Generar el QQ-Plot para evaluar normalidad
qqnorm(residuales_estandarizados2, main="QQ-Plot de los Residuales Estandarizados2", col="orange")
qqline(residuales_estandarizados2, col="black", lwd=2)  # Agregar línea de referencia

Para complementar el análisis gráfico se realiza el siguiente Test de normalidad aplicando la prueba de Shapiro-Wilk

H0: Los errores siguen una distribución normal H1: Los errores no siguen una distribución normal

# Determinar residuos ordinarios y aplicar test de normalidad
ei2 <- residuals(mod2)
shapiro.test(ei2)

## 
##  Shapiro-Wilk normality test
## 
## data:  ei2
## W = 0.98659, p-value = 8.367e-12

De la salida anterior se aprecia que el p-value = 8.367x10-12 es menor que un nivel de significancia del 0.001, eso significa que hay fuerte evidencia para rechazar la hipótesis de normalidad de los errores.

2. Los errores tienen varianza constante (homocedasticidad):

H0: Los errores tienen varianza constante H1: Los errores no tienen varianza constante

# Aplicar la prueba de Breusch-Pagan directamente con lmtest
bptest(mod2)

## 
##  studentized Breusch-Pagan test
## 
## data:  mod2
## BP = 387.65, df = 7, p-value < 2.2e-16

# Se obtienen los residuales del modelo

y_ajustados2 <- fitted(mod2)  # Valores ajustados (predichos)


# Visualización de la Homocedasticidad
# Gráfico de Residuales vs Valores Ajustados 
plot1.bt<-ggplot(data = df_mod2, aes(x = y_ajustados2, y = ei2)) +
  geom_point(color = "orange") +  # Puntos de residuales
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +  # Línea horizontal en cero
  labs(title = "Residuales vs Valores Ajustados",
       x = "Valores Ajustados",
       y = "Residuales") +
  theme_minimal()

#Gráfico de Raíz de Valor Absoluto de Residuales vs Valores Ajustados 
plot2.bt<-ggplot(data = df_mod2, aes(x = y_ajustados2, y = sqrt(abs(ei2)))) +
  geom_point(color = "orange") +  # Puntos
  geom_smooth(method = "lm", se = FALSE, color = "black") +  # Línea de tendencia
  labs(title = "Raiz Cuadrada del Valor Absoluto de Residuales vs Valores Ajustados",
       x = "Valores Ajustados",
       y = "√|Residuales|") +
  theme_minimal()

print(plot1.bt)

print(plot2.bt)

## `geom_smooth()` using formula = 'y ~ x'

3. Pruebas de independencia

# Visualización de los residuos 
# Se crea un data frame con los residuos y su índice de observación
resid_data2 <- data.frame(
  Observacion2 = seq_along(residuals(mod2)),  # Índice de cada observación
  Residuales2 = residuals(mod2)              # Valores de los residuos
)

# Gráfico de residuos usando ggplot2
ggplot(resid_data2, aes(x = Observacion2, y = Residuales2)) +
  geom_point(color = "lightblue", size = 2) +       # Puntos de los residuos
  geom_line(color = "black", linetype = "dashed") + # Línea guía de los residuos
  geom_hline(yintercept = 0, color = "red", linetype = "solid") + # Línea en y = 0
  labs(title = "Residuos del Modelo vs Index de Observacion",
       x = "Index de la Observacion", 
       y = "Residuales") +
  theme_minimal()  # Estilo de gráfico limpio

Para verificar el supuesto de independencia de los errores, se aplican pruebas estadísticas diseñadas para detectar la presencia de autocorrelación en los residuos del modelo de regresión.

Las hipótesis de estas pruebas son:

H0: Los errores son independientes, es decir, no presentan autocorrelación H1: Los errores no son independientes, es decir, presentan autocorrelación

# Pruebas de independencia de los errores 

# Test de Durbin-Watson
# Evalúa si los errores presentan autocorrelación de primer orden
dw_test2 <- dwtest(mod2)
print(dw_test)

## 
##  Durbin-Watson test
## 
## data:  mod1
## DW = 1.5326, p-value = 1.1e-08
## alternative hypothesis: true autocorrelation is greater than 0

Con un p-value tan pequeño 1.1x10-08, se rechaza la hipótesis de independencia. Es decir, los residuos del modelo presentan autocorrelación significativa

Pronóstico de precios

Intervalo de confianza

# Crear un nuevo conjunto de datos con los valores específicos de Edad para calcular el IC

new_dt2 <- data.frame(
  areaconst    = c(300, 300),
  parqueaderos = c(3, 3),
  banios       = c(3, 3),
  habitaciones = c(5, 5),
  estrato      = factor(c(5, 6), levels = levels(df_mod1$estrato))
)

# Calcular el intervalo de confianza al 95% para la media de Resistencia en los valores seleccionados
ic_result2 <- predict(object = mod2, newdata = new_dt2, interval = "confidence", level = 0.95)

# Convertir a data.frame para manipularlo
ic_result2 <- as.data.frame(ic_result2)

# Back-transformar (exp) para obtener resultados en escala original
ic_result2$fit      <- exp(ic_result2$fit)
ic_result2$lwr      <- exp(ic_result2$lwr)
ic_result2$upr      <- exp(ic_result2$upr)


# Mostrar los resultados
ic_result2

##         fit       lwr       upr
## 1  851.3552  795.4263  911.2165
## 2 1069.1322 1000.6858 1142.2603

Para un apartamento de 300 metros cuadrados estrato 5, que tenga 3 parqueadero, 3 baños y 5 habitaciones en la zona sur de Cali el intervalo de confianza del 95% es (795.4263, 911.2165).

Para un apartamento de 300 metros cuadrados estrato 6, que tenga 3 parqueadero, 3 baños y 5 habitaciones en la zona sur de Cali el intervalo de confianza del 95% es (1000.6858, 1142.2603)

De acuerdo con las predicciones del modelo se sugieren los siguientes apartamentos:

# Crear dataframe base de resultados
df_mod2_results <- df_mod2 %>%
  mutate(
    # Back-transform para obtener precios en escala original
    preciom_pred = exp(predict(mod2, newdata = df_mod2))
  )

# Solicitudes de la búsqueda 1
solicitud_2 <- df_mod2_results %>%
  filter(
    preciom_pred <= 850,
    estrato %in% c(5, 6),
    areaconst >= 300, 
    habitaciones >= 5, 
    banios >= 3, 
    parqueaderos >= 3
  )

# Tabla de resultados
kable(
  solicitud_2 %>%
    select(zona, tipo, estrato, areaconst, habitaciones, banios, parqueaderos,
           preciom, preciom_pred, longitud, latitud)
)

zona	tipo	estrato	areaconst	habitaciones	banios	parqueaderos	preciom	preciom_pred	longitud	latitud

# Mapa de los puntos potenciales
map_r2 <- leaflet(solicitud_2) %>%
  addTiles() %>%
  addMarkers(
    lng   = ~longitud,
    lat   = ~latitud,
    popup = ~paste0(
      "<b>Zona:</b> ", zona, "<br>",
      "<b>Tipo:</b> ", tipo, "<br>",
      "<b>Estrato:</b> ", estrato, "<br>",
      "<b>Área:</b> ", areaconst, " m²<br>",
      "<b>Habitaciones:</b> ", habitaciones, "<br>",
      "<b>Baños:</b> ", banios, "<br>",
      "<b>Parqueaderos:</b> ", parqueaderos, "<br>",
      "<b>Precio real:</b> ", round(preciom, 0), "<br>",
      "<b>Precio predicho:</b> ", round(preciom_pred, 0)
    )
  )

#map_r2

El modelo no arrojó opciones bajo las condiciones establecidas, a pesar de haber presentado un R² superior en comparación con el modelo aplicado para la Zona Norte. No obstante, al explorar el dataset original se confirma que tampoco existen apartamentos dentro del presupuesto definido, ya que el único que cumple con las características supera los 1000M.