Cargamos las librerias necesarias para trabajar.

library(readr)
library(ggplot2)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(lmtest)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(nortest)

Cargamos nuestro dataframe, el cual fue sacado de https://www.kaggle.com/datasets/danieleduardofajardo/colombia-house-prediction/data

df <- read_csv("/Users/usermac/Documents/especializacion/metodos de regresion/Sales_prediction_Colombia.csv")
## New names:
## Rows: 145552 Columns: 45
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: "," chr
## (31): antiguedad_original, balcon, banoservicio, conjuntocerrado, cuarto... dbl
## (14): ...1, area, areabalcon, areaconstruida, areaterraza, banos, estrat...
## ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
## Specify the column types or set `show_col_types = FALSE` to quiet this message.
## • `` -> `...1`
print(df)
## # A tibble: 145,552 × 45
##     ...1 antiguedad_original  area areabalcon areaconstruida areaterraza balcon 
##    <dbl> <chr>               <dbl>      <dbl>          <dbl>       <dbl> <chr>  
##  1     0 Entre 5 y 10 años   145           10          145            10 Terraza
##  2     1 Entre 0 y 5 años    114           NA          114            NA <NA>   
##  3     2 Entre 5 y 10 años   170           30          170            30 Terraza
##  4     3 Entre 0 y 5 años     61           NA           61            NA Balcón 
##  5     4 Más de 20 años      120.          NA          120.           NA <NA>   
##  6     5 Más de 20 años       56           NA           56            NA <NA>   
##  7     6 Entre 0 y 5 años     58           NA           58            NA <NA>   
##  8     7 Entre 10 y 20 años  211           NA          211            NA <NA>   
##  9     8 Entre 0 y 5 años     57.8         NA           57.8          NA <NA>   
## 10     9 Entre 0 y 5 años    225           NA          225            NA Terraza
## # ℹ 145,542 more rows
## # ℹ 38 more variables: banos <dbl>, banoservicio <chr>, conjuntocerrado <chr>,
## #   cuarto_de_escoltas <chr>, cuartodeservicio <chr>,
## #   depositoocuartoutil <chr>, depositos <chr>, estrato <dbl>,
## #   estudioobiblioteca <chr>, garajecubierto <chr>, garajes <dbl>,
## #   gimnasio <chr>, habitaciones <dbl>, halldealcobasoestar <chr>,
## #   instalaciondegas <chr>, jacuzzi <chr>, jardin <chr>, latitud <dbl>, …

1. Análisis descriptivo de los datos

Incluir: descripción de las variables, hipótesis a probar, gráfica, medidas de asociación entre variables. ### 1.1 Descripción de las variables. A continuación, se presenta la descripción completa de las variables del conjunto de datos sobre vivienda en Colombia. #### Características Generales del Inmueble

Distribución y Espacios Internos

  • bano: número de baño del inmueble. Número total de baños en la vivienda.
  • banosservicio: número de baños de servicio. Número de baños destinados al personal de servicio.
  • cuartodeservicio: Habitación destinada para el personal de servicio doméstico.
  • depositoocuartoutil: Espacio adicional, generalmente ubicado en el sótano o en un área común, destinado al almacenamiento.
  • depositos: Número de depósitos o cuartos útiles que le corresponden al inmueble.
  • estudioobiblioteca: Si, si tiene estudio o biblioteca. Indica con “Si” si la vivienda cuenta con un espacio dedicado a estudio o biblioteca.
  • habitaciones: número de habitaciones. Corresponde al número total de alcobas o dormitorios.
  • halldealcobasoestar: Un espacio de distribución que conecta las habitaciones o que funciona como una pequeña sala de estar.
  • instalaciondegas: Indica si el inmueble cuenta con la infraestructura para el suministro de gas.
  • jacuzzi: Indica la presencia de una bañera de hidromasaje.
  • sauna_yo_turco: Indica si la propiedad cuenta con sauna (calor seco) o baño turco (vapor).
  • zonadelavanderia: Área designada para las labores de lavandería.

Características del Edificio y Conjunto Residencial

  • conjuntocerrado: Si, si pertenence a un conjunto cerrado. Indica con “Si” si la vivienda forma parte de un complejo residencial privado con acceso restringido.
  • cuarto_de_escoltas: Espacio destinado a personal de seguridad o escoltas.
  • gimnasio: Indica la existencia de un gimnasio como parte de las áreas comunes.
  • numeroascensores: Cantidad de ascensores disponibles en el edificio.
  • parqueaderovisitantes: Zona de estacionamiento destinada a los visitantes del conjunto o edificio.
  • piscina: Indica la existencia de una piscina en las áreas comunes.
  • plantaelectrica: Indica si el edificio cuenta con un generador eléctrico de respaldo.
  • porteriaovigilancia: Presencia de un puesto de control en la entrada con personal de seguridad.
  • saloncomunal: Espacio de uso común para eventos y reuniones sociales.
  • terraza: Puede referirse a una terraza privada del inmueble o a una terraza común del edificio.
  • vigilancia: Tipo y nivel del servicio de seguridad en la propiedad.
  • zona_de_bbq: Área común equipada con parrilla para hacer asados.
  • zonasverdes: Espacios al aire libre con césped, plantas o jardines dentro del conjunto residencial.

Estacionamiento

  • garajecubierto: Número de estacionamientos que están bajo techo.
  • garajes: Número total de espacios de estacionamiento asignados al inmueble.
  • tipodegaraje: Describe el tipo de estacionamiento (ej. en línea, en servidumbre, independiente).

Ubicación y Otros

  • estrato: Clasificación socioeconómica del inmueble (en Colombia va del 1 al 6), que afecta las tarifas de los servicios públicos.
  • latitud: Coordenada geográfica de la ubicación del inmueble.
  • longitud: Coordenada geográfica de la ubicación del inmueble.
  • vista: Describe la panorámica o el tipo de paisaje que se observa desde el inmueble (ej. exterior, a un parque, a la ciudad).

1.2 Hipotesis a probar.

Queremos ver la relación entre el metraje de un inmuble y su precio final.

1.3 Graficas.

Primero realizamos una grafica para verificar la relación entre \(area\) y \(valorventa\)

ggplot(data = df, aes(x = area, y = valor)) +
  geom_point(alpha = 0.7) +
  labs(title = "Relación entre area y valorventa",
       x = "area",
       y = "valor"
       ) +
  theme_minimal()
## Warning: Removed 24 rows containing missing values or values outside the scale range
## (`geom_point()`).

Esta relación no parece clara, y parece ser debido a outliers. Ahora hacemos una grafica de distribución de la variable dependiente \(valorventa\)

ggplot(data=df, aes(x=valorventa)) + 
  geom_density(
    fill = "skyblue",
    alpha = 0.6, 
    color = "darkblue"
  ) +
  labs(
    title = "Density of valor venta",
    x = "valor venta",
    y = "density"
  )

Vemos que existen ciertos outliers para esta variable. Utiliznando el rango intercuartil.

outliers_iqr = boxplot.stats(df$valor)$out
head(outliers_iqr, n=20)
##  [1] 1750000000 2925000000 1900000000 1750000000 2500000000 2380241095
##  [7] 1580000000 1650000000 2800000000 1750000000 1650000000 1650000000
## [13] 1600000000 3200000000 1600000000 1900000000 1750000000 2100000000
## [19] 1600000000 3100000000

Filtramos estos valores de nuestro dataframe

df = df[!(df$valor %in%outliers_iqr), ]

Y revisamos de nuevo como se ven estos.

ggplot(data=df, aes(x=valorventa)) + 
  geom_density(
    fill = "skyblue",
    alpha = 0.6, 
    color = "darkblue"
  ) +
  labs(
    title = "Density of valor venta",
    x = "valor venta",
    y = "density"
  )

Hacemos lo mismo con la variable explicativa (o independiente) \(area\), y vemos que tambien tiene valores atipocos.

ggplot(data=df, aes(x=area)) + 
  geom_density(
    fill = "skyblue",
    alpha = 0.6, 
    color = "darkblue"
  ) +
  labs(
    title = "Density of valor venta",
    x = "valor venta",
    y = "density"
  )
## Warning: Removed 24 rows containing non-finite outside the scale range
## (`stat_density()`).

outliers_iqr_a = boxplot.stats(df$area)$out
head(outliers_iqr_a, n=20)
##  [1]  224.00  232.00  251.00 6130.00  220.00  245.00  264.00  274.00  398.00
## [10]  270.00  222.00  307.88  240.00  218.00  247.00  730.00  233.00  227.00
## [19]  381.00  287.00

filtramos estos valores de nuestro dataframe.

df = df[!(df$area %in% outliers_iqr_a), ]
ggplot(data=df, aes(x=area)) + 
  geom_density(
    fill = "skyblue",
    alpha = 0.6, 
    color = "darkblue"
  ) +
  labs(
    title = "Density of valor venta",
    x = "valor venta",
    y = "density"
  )
## Warning: Removed 24 rows containing non-finite outside the scale range
## (`stat_density()`).

Ahora graficamos la relación entre \(area\) y \(valorventa\)

ggplot(data = df, aes(x = area, y = valorventa)) +
  geom_point(alpha = 0.7) +
  labs(title = "Relación entre area y valorventa",
       x = "area",
       y = "valor"
       ) +
  theme_minimal()
## Warning: Removed 24 rows containing missing values or values outside the scale range
## (`geom_point()`).

Vemos que existen valores vacios estas variables.

df=df[!is.na(df$valor), ]
df=df[!is.na(df$area), ]

Y tambien tenemos variables que no tienen sentido de negocio, tales como areas menores de 10 metros cuadrados 0 viviendas de 10 millones de pesos.

df = df[df$area>10 & df$valor>10000000,]
ggplot(data = df, aes(x = area, y = valor)) +
  geom_point(alpha = 0.7) +
  labs(title = "Relación entre area y valorventa",
       x = "area",
       y = "valor"
       ) +
  theme_minimal()

2. Ajuste del modelo.

Presentar el modelo ajustado, interpretación (no trivial) de los coeficientes.

2.1 Regresión lineal simple

Ajustamos una regresión lineal simple.

model1 = lm(formula = valor ~ area, data=df)
summary_model1 = summary(model1)
print(summary_model1)
## 
## Call:
## lm(formula = valor ~ area, data = df)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -1.025e+09 -1.007e+08 -2.387e+07  8.304e+07  1.137e+09 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -77818367    1141964  -68.14   <2e-16 ***
## area          6293868      10977  573.36   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 160400000 on 128941 degrees of freedom
## Multiple R-squared:  0.7183, Adjusted R-squared:  0.7183 
## F-statistic: 3.287e+05 on 1 and 128941 DF,  p-value: < 2.2e-16
model2 = lm(formula = log(valor) ~ area, data=df)
summary_model2 = summary(model2)
print(summary_model2)
## 
## Call:
## lm(formula = log(valor) ~ area, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.3951 -0.2051  0.0082  0.2200  1.7682 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1.878e+01  2.345e-03    8008   <2e-16 ***
## area        1.184e-02  2.255e-05     525   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3294 on 128941 degrees of freedom
## Multiple R-squared:  0.6813, Adjusted R-squared:  0.6813 
## F-statistic: 2.756e+05 on 1 and 128941 DF,  p-value: < 2.2e-16

2.2 Interpretación de coeficientes (no trivial).

para el modelo 1 ambos coeficientes son representativos. El intrecepto es -77.818.367, aunque en este contexto no tiene sentido interpretarlo, ya que no existen inmuebles de area 0. El coeficiente de area es de 6.293.868, y es que por cada metro cuadrado (m^2) adicional en el área del inmueble, el precio de venta esperado aumenta en promedio en $6,293,868 COP.

res = residuals(model1)
res2 = residuals(model2)
(exp(1.184e-02)-1)*100
## [1] 1.191037

Ahora, para el modelo 2, ambos coeficientes representativos. El intercepto es de 1.878e+01, pero no tiene interpretación para este caso ya que no tenemos areas de 0 metros cuadrados. Ahora el coeficiente para el area es de 1.184e-02, es decir que por cada metro cuadrado (m^2) adicional en el area del inmueble, el precio de venta aumenta 1.19% ## 3. Validación del modelo. ### 3.1 Bondad de ajuste.

r_cuadrado <- summary_model1$r.squared
print(paste("R-cuadrado:", r_cuadrado))
## [1] "R-cuadrado: 0.718273287233378"
r_cuadrado_2 <- summary_model2$r.squared
print(paste("R-cuadrado model 2:", r_cuadrado_2))
## [1] "R-cuadrado model 2: 0.68129100220572"

Es decir que nuestra variable \(area\) explica el 71% de la variabilidad de la variable \(valorventa\) Y para el modelo lin-lon, la variable \(area\) explica el 68% de la variabilidad de la variable \(valorventa\)

cor(df$area, df$valorventa)
## [1] 0.8475101

Si vemos la correlación, esta es positiva y fuerte. Es decir que cuando area creace, el precio por metro cuadrado crece. ### 3.2 Analisis de residuos.

3.2.1 Los residuos tienen media cero

mean(res)
## [1] 1.746934e-08
mean(res2)
## [1] 2.39756e-17

Tiene media de los residuos muy cercana a cero. #### 3.2.2. Homoceasticidad

plot(model1, which = 1)

plot(model2, which = 1)

Ahora realizamos una prueba formal. Las hipótesis son: \(H_{0}\): Todas las muestras provienen de poblaciones con varianzas iguales \(H_{1}:\) Al menos dos muestras provienen de poblaciones con varianzas diferentes. Dividiremos el el dataframe en dos muestras, una donde el area es mayor al promedio y otra donde el area es menor al promedio.

area_promedio = mean(df$area)
residuals_over = res[df$area>=area_promedio]
residuals_under = res[df$area<area_promedio] 
bartlett.test(list(bajos=residuals_under, altos=residuals_over))
## 
##  Bartlett test of homogeneity of variances
## 
## data:  list(bajos = residuals_under, altos = residuals_over)
## Bartlett's K-squared = 20739, df = 1, p-value < 2.2e-16

Rechazamos la hipotesis nula, por lo tanto aceptamos la verdadera, la dos muestras tienen varianzas diferentes. Ahora, para el modelo lin-lon

area_promedio = mean(df$area)
residuals_over = res2[df$area>=area_promedio]
residuals_under = res2[df$area<area_promedio] 
bartlett.test(list(bajos=residuals_under, altos=residuals_over))
## 
##  Bartlett test of homogeneity of variances
## 
## data:  list(bajos = residuals_under, altos = residuals_over)
## Bartlett's K-squared = 1730, df = 1, p-value < 2.2e-16

Tambien rechazamos la hipotesis nula, y aceptamos la alternativa de que tienen varianzas diferentes.

Podemos seguir iterando mas pruebas estadisticas, Por ejemplo de Breusch-Pagan test \(H_{0}:\) El modelo es homocedástico (la varianza de los errores es constante) \(H_{1}:\) El modelo es heterocedástico (la varianza de los errroes no es constante. )

bptest(model1)
## 
##  studentized Breusch-Pagan test
## 
## data:  model1
## BP = 11891, df = 1, p-value < 2.2e-16

Dado es p-valor, rechazamos la hipotesis nula y aceptamos la alternativda, es decir el modelo es heterocedástico Por lo tanto tenemos suficiente evidencia par decir que el modelo es heterocedástico.

bptest(model2)
## 
##  studentized Breusch-Pagan test
## 
## data:  model2
## BP = 1389.9, df = 1, p-value < 2.2e-16

Seguimos rechazando la hipotesis nula

3.2.3 Independecia de los errores

Esta la podemos asusmir cierta desde que la información no tiene relación entre filas.

3.2.4 Los residuos tienen distribución normal.

qqnorm(res)
qqline(res)

lillie.test(res)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  res
## D = 0.063448, p-value < 2.2e-16
residuos_muestra <- sample(res, 5000)
shapiro.test(residuos_muestra)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuos_muestra
## W = 0.95669, p-value < 2.2e-16

Los errores tampoco siguen una distribución normal. Ahora para el modelo lin-lon

qqnorm(res2)
qqline(res2)

lillie.test(res2)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  res2
## D = 0.020159, p-value < 2.2e-16
residuos_muestra2 <- sample(res2, 5000)
shapiro.test(residuos_muestra2)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuos_muestra2
## W = 0.99551, p-value = 2.736e-11

Los residuos tampoco tienen una distribución normal para el modelo lin-log ## 4. Uso del modelo. Aunque se ha detectado que el modelo viola los supuestos de homocedasticidad y normalidad, lo que afecta la fiabilidad de los intervalos de confianza, a continuación se presenta una estimación puntual y los intervalos calculados por el modelo. Para obtener predicciones más robustas, sería necesario primero corregir estas violaciones, por ejemplo, mediante una transformación de la variable dependiente

nuevos_datos <- data.frame(area = c(45, 100, 300))
intervalo_confianza <- predict(model1, newdata = nuevos_datos, interval = "confidence")

print(intervalo_confianza)
##          fit        lwr        upr
## 1  205405670  204006281  206805058
## 2  551568381  550688173  552448589
## 3 1810341878 1805860926 1814822830
intervalo_prediccion <- predict(model1, newdata = nuevos_datos, interval = "prediction")

print(intervalo_prediccion)
##          fit        lwr        upr
## 1  205405670 -108952517  519763857
## 2  551568381  237212077  865924686
## 3 1810341878 1495954870 2124728885
intervalo_confianza2 <- predict(model2, newdata = nuevos_datos, interval = "confidence")

print(intervalo_confianza2)
##        fit      lwr      upr
## 1 19.31496 19.31209 19.31784
## 2 19.96597 19.96416 19.96778
## 3 22.33327 22.32406 22.34247
intervalo_prediccion2 <- predict(model2, newdata = nuevos_datos, interval = "prediction")
print(intervalo_prediccion)
##          fit        lwr        upr
## 1  205405670 -108952517  519763857
## 2  551568381  237212077  865924686
## 3 1810341878 1495954870 2124728885