1. Analisis exploratorio univariado

Se cargan las librerías necesarias para el desarrollo del trabajo

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

Inicialmente se cargan los datos que se encuentran en el archivo (vivienda4.RDS) y se hace un resumen de estos datos para entender si es necesario realizar algún tipo de preprocesamiento antes de avanzar con su posterior análisis y modelación.

rm(list=ls())
data(vivienda4)
vivienda <- as.data.frame(vivienda4)
summary(vivienda)
##            zona      estrato     preciom        areaconst     
##  Zona Centro :   8   3:   0   Min.   :207.4   Min.   : 40.00  
##  Zona Norte  : 288   4:1706   1st Qu.:230.7   1st Qu.: 60.00  
##  Zona Oeste  :  60   5:   0   Median :238.8   Median : 75.00  
##  Zona Oriente:   6   6:   0   Mean   :243.7   Mean   : 87.63  
##  Zona Sur    :1344            3rd Qu.:251.5   3rd Qu.: 98.00  
##                               Max.   :309.7   Max.   :200.00  
##           tipo     
##  Apartamento:1363  
##  Casa       : 343  
##                    
##                    
##                    
##

Se logra observar que solo hay viviendas en el estrato 4, que hay 5 zonas donde la que posee la mayor cantidad de datos es la Zona Sur, seguida de la Zona Norte, también hay 343 viviendas bajo la categoría de Casa, esto se debe excluir ya que no hace parte del interés de la inmobiliaria. Se procede entonces con la depuración de los datos

  viviendaf <- vivienda %>%
  dplyr::select(-estrato) %>%
  filter(tipo == "Apartamento") %>%
  filter(areaconst <= 200)
viviendaf$tipo <- droplevels(viviendaf$tipo)

Se realiza nuevamente un resumen ahora sobre los datos de vivienda_filtrados

summary(viviendaf)
##            zona         preciom        areaconst               tipo     
##  Zona Centro :   7   Min.   :207.4   Min.   : 40.00   Apartamento:1363  
##  Zona Norte  : 237   1st Qu.:228.8   1st Qu.: 60.00                     
##  Zona Oeste  :  52   Median :236.1   Median : 70.00                     
##  Zona Oriente:   2   Mean   :237.7   Mean   : 75.48                     
##  Zona Sur    :1065   3rd Qu.:243.6   3rd Qu.: 84.00                     
##                      Max.   :305.2   Max.   :200.00

Se observa que se retira la variable \(estrato\), ya que solo aporta contexto al análisis y como ahora la variable \(tipo\) solo tiene información sobre apartamentos. Quedando finalmente con 1363 registros y 4 variables, donde:

Precio

Precio (\(preciom\)), es una variable continua, con un precio medio de aproximadamente 236.1 millones COP, y el 50% de los apartamentos tienen un precio entre 228.8 y 243.6 millones COP (según el primer y tercer cuartil), por otra parte el precio mínimo de los apartamentos es de 207.4 millones COP hasta un máximo de 305.2 millones COP.

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

Área

Área (\(areaconst\)), es también una variable continua, con una área media construida de aproximadamente 75.48\(m^2\), y el 50% de los apartamentos tienen un área construida entre 60 y 84\(m^2\) (según el primer y tercer cuartil). El área construida de los apartamentos varía desde un mínimo de 40\(m^2\) hasta un máximo de 200\(m^2\).

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

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

2. Analisis exploratorio bivariado

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

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

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

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

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

Por otra parte si no se tuviese normalidad en las variables se puede aplicar el coeficiente de Spearman

cor.test(viviendaf$areaconst,viviendaf$preciom, method = "spearman")
## Warning in cor.test.default(viviendaf$areaconst, viviendaf$preciom, method =
## "spearman"): Cannot compute exact p-value with ties
## 
##  Spearman's rank correlation rho
## 
## data:  viviendaf$areaconst and viviendaf$preciom
## S = 109350121, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.7408906

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

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

3. Modelación

Se va estimar el siguiente modelo

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

modelo <- lm(preciom ~ areaconst, data = viviendaf)
summary(modelo)
## 
## Call:
## lm(formula = preciom ~ areaconst, data = viviendaf)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -26.5139  -5.0886  -0.0031   4.6406  24.3309 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 2.001e+02  6.698e-01  298.67   <2e-16 ***
## areaconst   4.984e-01  8.503e-03   58.62   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.081 on 1361 degrees of freedom
## Multiple R-squared:  0.7163, Adjusted R-squared:  0.7161 
## F-statistic:  3436 on 1 and 1361 DF,  p-value: < 2.2e-16

El coeficiente \(\hat{\beta_0}\) es 200.1, su interpretación sería que con 0 metros cuadrados, el precio promedio de las viviendas sería de 200.1 millones, sin embargo, no es coherente, pues una vivienda de 0 metros no existe.

El coeficiente \(\hat{\beta_1}\) es 0.4984, nos indica cuánto incrementa el precio de las viviendas por cada aumento de un metro de área construída. Es decir que por cada aumento en un metro cuadrado de los apartamentos, el precio incrementa 0.4984 millones. El valor P <2e-16 nos indica que es significativo.

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

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

4. Intervalo de confianza

ic <- confint(modelo, level = 0.95)
ic["areaconst", ]
##     2.5 %    97.5 % 
## 0.4817357 0.5150970

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

Adicional, este resultado se confirma realizando la correspondiente prueba \(t\), sobre el coeficiente \(\hat{\beta_1}\)

r <- summary(modelo)
print(r$coefficients["areaconst", c("t value", "Pr(>|t|)")])
##  t value Pr(>|t|) 
## 58.61576  0.00000

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

5. Indicador de bondad de ajuste

r$r.squared
## [1] 0.7162696

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

6. Predicción

predict(modelo, data.frame(areaconst = 110), interval = "prediction", level = 0.95)
##        fit      lwr      upr
## 1 254.8893 240.9814 268.7971

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

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

predict(modelo, data.frame(areaconst = 110), interval = "confidence", level = 0.95)
##        fit      lwr      upr
## 1 254.8893 254.2014 255.5771

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

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

7. Validación supuestos

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

residuos <- residuals(modelo)

dfr <- data.frame(residuos = residuals(modelo))

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

b

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

d

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

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

Normalidad

shapiro.test(residuos)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuos
## W = 0.99885, p-value = 0.5419

Dado que el valor p es 0.5419, no se rechaza la hipótesis nula y podemos asumir que los residuos del modelo siguen una distribución normal.

Homocedasticidad

bptest <- bptest(modelo)
bptest
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo
## BP = 0.83288, df = 1, p-value = 0.3614

Dado que el valor p es 0.3614, no se rechaza la hipótesis nula y podemos asumir que la varianza de los errores es constante, lo que indica homocedasticidad.

dwtest <- dwtest(modelo)
dwtest
## 
##  Durbin-Watson test
## 
## data:  modelo
## DW = 2.0204, p-value = 0.6435
## alternative hypothesis: true autocorrelation is greater than 0

Autocorrelación

El valor p es 0.6435, no se rechaza la hipótesis nula podemos asumir que no hay autocorrelación en los residuos del modelo.

Por lo anterior se concluye que el modelo cumple todos los supuestos.

8. Transformación

Al evaluar los supuestos del modelo de regresión se ve que no es necesario realizar una trasnformación a los datos debido a que no se viola ninguno de estos. Sin embargo, se va a estimar y contruir un modelo a partir de la metodología de Boxcox con el fin de ver el compartamiento de este

MASS::boxcox(modelo)

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

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

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

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

El nuevo modelo estimado con la transformación es

modelo2 <- lm(preciom_transformado ~ areaconst, data = viviendaf)
summary(modelo2)
## 
## Call:
## lm(formula = preciom_transformado ~ areaconst, data = viviendaf)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -26.5139  -5.0886  -0.0031   4.6406  24.3309 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1.991e+02  6.698e-01  297.18   <2e-16 ***
## areaconst   4.984e-01  8.503e-03   58.62   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.081 on 1361 degrees of freedom
## Multiple R-squared:  0.7163, Adjusted R-squared:  0.7161 
## F-statistic:  3436 on 1 and 1361 DF,  p-value: < 2.2e-16

A este modelo también se decide probar sus supuestos

# Residuales
residuos2 <- residuals(modelo2)

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

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

b2

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

d2

La transformación evidencia una escala similar (entre −20 y 20 ), lo que cambia la interpretación de los coeficientes del modelos. También aparenta mostrar menos cantidad de outliers.

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

Normalidad

shapiro.test(residuos2)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuos2
## W = 0.99885, p-value = 0.5419

Homocedasticidad

bptest(modelo2)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo2
## BP = 0.83288, df = 1, p-value = 0.3614

Autocorrelación

dwtest(modelo2)
## 
##  Durbin-Watson test
## 
## data:  modelo2
## DW = 2.0204, p-value = 0.6435
## alternative hypothesis: true autocorrelation is greater than 0

El supuesto de normalidad, homocedasticidad e independencia de los errores como evidencian las pruebas anteriores siguen cumpliendose.

9. Comparación de modelos

r2 <- summary(modelo2)
r$r.squared
## [1] 0.7162696
r2$r.squared
## [1] 0.7162696

Ya que ambos modelos cumplen los supuestos, el modelo inicial posee un \(R^2\)=0.7162696 , y el modelo con la transformación de box-cox tiene un \(R^2\)=0.7162696. Por lo que reafirmamos no es necesaria una transformación debido a que arroja practicamente los mismos resultados.

10. Estimación de varios modelos

modelo1_1=lm(preciom ~ areaconst, data=viviendaf)           # Lin - Lin
modelo2_1=lm(preciom ~ log(areaconst), data=viviendaf)      # Lin - Log
modelo3_1=lm(log(preciom) ~ areaconst, data=viviendaf)      # Log - Lin
modelo4_1=lm(log(preciom) ~ log(areaconst), data=viviendaf) # Log - Log
modelo5_1=lm(sqrt(preciom) ~ areaconst, data=viviendaf)     #Box-Cox lambda 0.5
modelo6_1=lm((1/preciom) ~ areaconst, data=viviendaf)       #Box-Cox lambda -1


stargazer(modelo1_1, modelo2_1, modelo3_1, modelo4_1, modelo5_1, modelo6_1, type="text", df=FALSE)
## 
## ==================================================================================================
##                                                  Dependent variable:                              
##                     ------------------------------------------------------------------------------
##                              preciom                log(preciom)        sqrt(preciom) (1/preciom) 
##                         (1)          (2)          (3)          (4)           (5)          (6)     
## --------------------------------------------------------------------------------------------------
## areaconst             0.498***                  0.002***                  0.016***    -0.00001*** 
##                       (0.009)                  (0.00004)                  (0.0003)     (0.00000)  
##                                                                                                   
## log(areaconst)                    42.878***                  0.174***                             
##                                    (0.794)                   (0.003)                              
##                                                                                                   
## Constant             200.063***   53.820***     5.318***     4.723***     14.218***     0.005***  
##                       (0.670)      (3.409)      (0.003)      (0.014)       (0.022)     (0.00001)  
##                                                                                                   
## --------------------------------------------------------------------------------------------------
## Observations           1,363        1,363        1,363        1,363         1,363        1,363    
## R2                     0.716        0.682        0.696        0.674         0.706        0.673    
## Adjusted R2            0.716        0.682        0.695        0.674         0.706        0.672    
## Residual Std. Error    7.081        7.496        0.030        0.031         0.230        0.0001   
## F Statistic         3,435.808*** 2,919.088*** 3,110.029*** 2,814.311*** 3,272.232***  2,797.204***
## ==================================================================================================
## Note:                                                                  *p<0.1; **p<0.05; ***p<0.01

Se observa que ninguna de las transformaciones mejora el ajuste del modelo inicial, concluyendo así que este es el mejor modelo estimado en comparación con los otros. Además, cabe recordar que a este modelo ya se le realizo la validación de supuestos de forma gráfica y por medio de pruebas de hipótesis en los puntos anteriores.

11. Conclusiones

Se hace evidente la relación entre las variables \(precio\) y \(área\) de los apartamentos de estrato 4 con menos de 200\(m^2\), y esto lo corrobora la prueba de correlación, sin embargo se hace notable la presencia de outliers en los datos, por lo que resulta de interés para un futuro realizar algún tipo de tratamiento que permita trabajar con estos.

Las transformaciones implementadas no fueron realmente útiles ya que además de que afectan directamente la interpretación de los coeficientes del modelo, no lograron mejorar significativamente el modelo inicial planteado.

Se recomienda la inclusión de mas covariables las cuales aporten mayor información y con esto tener una mejoria en el modelo de regresión.