Instalando librerias

library(paqueteMETODOS)
library(summarytools)
library(corrplot)
library(PerformanceAnalytics)

Adquirendo datos

Con base en los datos de ofertas de vivienda descargadas del portal Fincaraiz para apartamento de estrato 4 con área construida menor a 200 m2 (vivienda4.RDS) la inmobiliaria A&C requiere el apoyo de un cientifico de datos en la construcción de un modelo que lo oriente sobre los precios de inmuebles.

data(vivienda4)
df = subset(vivienda4, vivienda4$tipo=="Apartamento" & vivienda4$estrato==4 & vivienda4$areaconst<=200)

Analisis

Exploratorio

  1. Realice un análisis exploratorio de las variables precio de vivienda (millones de pesos COP) y área de la vivienda (metros cuadrados) - incluir gráficos e indicadores apropiados interpretados.
# Antes de iniciar el analisis verificamos que el dataset contenga solo la informacion requerida (Apartamentos con menos de 200 m2 y estrato 4)

head(df,10)
## # A tibble: 10 × 5
##    zona       estrato preciom areaconst tipo       
##    <fct>      <fct>     <dbl>     <dbl> <fct>      
##  1 Zona Norte 4           220        52 Apartamento
##  2 Zona Norte 4           320       108 Apartamento
##  3 Zona Sur   4           290        96 Apartamento
##  4 Zona Norte 4           220        82 Apartamento
##  5 Zona Norte 4           220        75 Apartamento
##  6 Zona Norte 4           162        60 Apartamento
##  7 Zona Norte 4           225        84 Apartamento
##  8 Zona Norte 4           370       117 Apartamento
##  9 Zona Norte 4           155        60 Apartamento
## 10 Zona Norte 4           240        75 Apartamento
summary(df$estrato)
##    3    4    5    6 
##    0 1363    0    0
summary(df$areaconst)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   40.00   60.00   70.00   75.48   84.00  200.00
summary(df$tipo)
## Apartamento        Casa 
##        1363           0
precio_zona = aggregate(preciom ~ zona, data=df, FUN = mean)
print(precio_zona)
##           zona  preciom
## 1  Zona Centro 236.0000
## 2   Zona Norte 204.5907
## 3   Zona Oeste 200.7692
## 4 Zona Oriente 262.5000
## 5     Zona Sur 201.7061
# Preciom 
hist(df$preciom, freq = F, ylim = c(0,0.01))
lines(density(df$preciom))

shapiro.test(df$preciom)                       # Prueba de normalidad
## 
##  Shapiro-Wilk normality test
## 
## data:  df$preciom
## W = 0.89881, p-value < 2.2e-16
boxplot(df$preciom, main="Distribucion Precio m2", horizontal = T)

summarytools::descr(df$preciom)                # Resumen estadisticos Precio
## Descriptive Statistics  
## df$preciom  
## N: 1363  
## 
##                     preciom
## ----------------- ---------
##              Mean    202.44
##           Std.Dev     65.29
##               Min     78.00
##                Q1    153.00
##            Median    185.00
##                Q3    240.00
##               Max    645.00
##               MAD     59.30
##               IQR     86.50
##                CV      0.32
##          Skewness      1.44
##       SE.Skewness      0.07
##          Kurtosis      3.83
##           N.Valid   1363.00
##         Pct.Valid    100.00

Observando el histograma se puede notar una leve tendencia a la normalidad pero debe ser confirmado con las pruebas pertinentes, al realizar la prueba de normalidad tenemos un resultado de p-value < 2.2e-16 por lo que debemos rechazar la hipotesis nula para normalidad de la variable. Analizando la distribucion de quartiles se infiere una mayor concentracion de apartamentos hacia precios por encima de la media (el rango IQR es de solo 87 mientras que la distancia entre Q3 y el valor maximo es de 405), esto es confirmado al observar el boxplot

# Area Construida 
hist(df$areaconst, freq = F, ylim = c(0,0.035))
lines(density(df$areaconst))

shapiro.test(df$areaconst)                     # Prueba de normalidad
## 
##  Shapiro-Wilk normality test
## 
## data:  df$areaconst
## W = 0.82347, p-value < 2.2e-16
boxplot(df$preciom, main="Distribucion Area Construida", horizontal = T)

summarytools::descr(df$areaconst)              # Resumen estadisticos Area
## Descriptive Statistics  
## df$areaconst  
## N: 1363  
## 
##                     areaconst
## ----------------- -----------
##              Mean       75.48
##           Std.Dev       22.56
##               Min       40.00
##                Q1       60.00
##            Median       70.00
##                Q3       84.00
##               Max      200.00
##               MAD       14.83
##               IQR       24.00
##                CV        0.30
##          Skewness        2.08
##       SE.Skewness        0.07
##          Kurtosis        6.32
##           N.Valid     1363.00
##         Pct.Valid      100.00
table(df$zona)                                 # Distribucion Zona
## 
##  Zona Centro   Zona Norte   Zona Oeste Zona Oriente     Zona Sur 
##            7          237           52            2         1065
round(prop.table(table(df$zona)) * 100, 2)     # Porcentajes Zona
## 
##  Zona Centro   Zona Norte   Zona Oeste Zona Oriente     Zona Sur 
##         0.51        17.39         3.82         0.15        78.14

Al igual que con el precio, para el area construida es dificil confirmar si tiene una distribucion normal o no solo con el histograma, por lo que debe ser sometida a pruebas de normalidad, en este caso obtenemos un resultado de p-value < 2.2e-16 por lo que rechazamos la hipotesis nula y la variable area construida no tiene normalidad. Con el area construida tambien hay una distribucion similar al precio, lo que sugiere una posible correlacion (el IQR es de 24, con una diferencia Max-Q3 = 116)

Ambas variables tienen un coeficiente de variacion similar (al rededor del 30%)

Bivariado

  1. Realice un análisis exploratorio bivariado de datos, enfocado en la relación entre la variable respuesta (precio) en función de la variable predictora (area construida) - incluir gráficos e indicadores apropiados interpretados.
x = df$areaconst      # Variable predictora
y = df$preciom        # Variable respuesta

plot(x, y, main="Grafico de dispersion entre X y Y")

cor.test(x,y)         # Correlacion de Pearson
## 
##  Pearson's product-moment correlation
## 
## data:  x and y
## t = 41.595, df = 1361, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.7237938 0.7706237
## sample estimates:
##       cor 
## 0.7481389

El resultado del coeficiente de relacion de Pearson confirma lo que se insinua en el grafico de dispersion, hay una correlacion fuerte (0.75)

Modelado

Estimacion de modelo

  1. Estime el modelo de regresión lineal simple entre precio=f(area)+ε. Interprete los coeficientes del modelo β0, β1 en caso de ser correcto.
ml = lm(y~x)
summary(ml)
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -225.404  -23.902   -4.754   25.763  209.021 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 39.04679    4.09977   9.524   <2e-16 ***
## x            2.16473    0.05204  41.595   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 43.34 on 1361 degrees of freedom
## Multiple R-squared:  0.5597, Adjusted R-squared:  0.5594 
## F-statistic:  1730 on 1 and 1361 DF,  p-value: < 2.2e-16

De acuerdo a este modelo se tiene \(Y=β0+β1X+ε\) donde:

  • X es la variable predictora o independiente, para este ejercicio, el area construida
  • Y es la variable respuesta o dependiente, para este ejercicio, el precio
  • β0 = El intercepto es 39.04, esto quiere decir que si el area es 0, el precio seria de 39.04, aunque cabe resaltar que no es posible que el area construida sea igual a 0.
  • β1 = Es la pendiente, en este caso 2.16. Se puede interpretar como la tasa de crecimiento de β0 por cada incremento en en el Area

Intervalo de confianza

  1. Construir un intervalo de confianza (95%) para el coeficiente β1, interpretar y concluir si el coeficiente es igual a cero o no. Compare este resultado con una prueba de hipótesis t.
confint(ml, level = 0.95)
##                2.5 %    97.5 %
## (Intercept) 31.00423 47.089340
## x            2.06264  2.266826

Como podemos ver aqui, el valor de β1 puede estar dentro del rango de 2.06 y 2.26 millones por cada metro cuadrado

En ninguno de los casos el coeficiente fue igual a 0, por lo que podemos confirmar que el precio y el area construida estan correlacionados.

indicador de bondad \(R^2\).

  1. Calcule e interprete el indicador de bondad \(R^2\)
cat("El R2 es igual a ", summary(ml)$r.squared)
## El R2 es igual a  0.5597117

El coeficiente de determinacion o \(R^2\), es una medida estadística que representa la proporción de la variabilidad de la variable dependiente que es explicada por las variables independientes en un modelo de regresión. Es una medida crucial para evaluar la bondad de ajuste de un modelo de regresión. en este case se tiene un \(R^2\) de 0.5597, lo que indica que la variable area explica un 55.97% de la variabilidad en el precio por metro cuadrado para los apartamentos estudiados

Estimacion de precio

  1. ¿Cuál sería el precio promedio estimado para un apartamento de 110 metros cuadrados? Considera entonces con este resultado que un apartamento en la misma zona con 110 metros cuadrados en un precio de 200 millones sería una atractiva esta oferta? ¿Qué consideraciones adicionales se deben tener?.
predict( ml, data.frame(x=110), interval="confidence", level = 0.95 )
##        fit      lwr      upr
## 1 277.1674 272.9573 281.3775
predict( ml, data.frame(x=110), interval="prediction", level = 0.95 )
##        fit      lwr      upr
## 1 277.1674 192.0449 362.2899

Como se puede observar en los resultados, el precio para un apartamento de 110 metros oscila entre 272.9 y 281.37 con un 95% de confianza, y el precio esperado con el modelo actual es de 277.1 millones de pesos por lo que una oferta de 200 millones es extremadamente atractiva.

El factor mas influyente a tener en cuenta es que como se vio en el analisis exploratorio, la zona tiene un gran impacto en el precio de los apartamentos, generando variaciones promedio de hasta 60 millones de pesos.

Supuestos

  1. Realice la validación de los supuestos del modelo por medio de gráficos apropiados, interpretarlos y sugerir posibles soluciones si se violan algunos de ellos. Utilice las pruebas de hipótesis para la validación de supuestos y compare los resultados con lo observado en los gráficos asociados.
par(mfrow=c(2,2))
plot(ml)

Analizando los graficos se determina que los residuales tienen:

shapiro.test(residuals(ml))        # Prueba de normalidad sobre residuos del modelo de regresion
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(ml)
## W = 0.96486, p-value < 2.2e-16

Se obtienen un p-value < 2.2e-16 por lo que se puede determinar que no hay normalidad, el supuesto no se cumple.

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

Transformacion

  1. De ser necesario realice una transformación apropiada para mejorar el ajuste y supuestos del modelo.
m1 = ml                                          # Lin - Lin
m2 = lm(formula = y ~ log(x))                    # Lin - Log
m3 = lm(formula = log(y) ~ x)                    # Log - Lin    
m4 = lm(formula = log10(y) ~ log10(x))           # Log - Log
  1. De ser necesario compare el ajuste y supuestos del modelo inicial y el transformado.

  2. Estime varios modelos y compare los resultados obtenidos. En el mejor de los modelos, ¿se cumplen los supuestos sobre los errores?

# Se aplica Criterio de Informacion de Akaike para evaluar ajuste de cada modelo (menor es mejor)

aic_m1 = AIC(m1)
aic_m2 = AIC(m2)
aic_m3 = AIC(m3)
aic_m4 = AIC(m4)

AIC_values <- c(aic_m1, aic_m2, aic_m3, aic_m4)

nombres_modelos <- c("Lin - Lin", "Lin - Log", "Log - Lin", "Log - Log")
df_AIC <- data.frame(Modelo = nombres_modelos[order(AIC_values)], AIC = sort(AIC_values))
df_AIC
##      Modelo        AIC
## 1 Log - Log -2909.1990
## 2 Log - Lin  -446.9303
## 3 Lin - Log 14059.7466
## 4 Lin - Lin 14146.4576

Por AIC, la transformacion Log-Log parece tener un mejor ajuste, se procede a evaluar los supuestos

lista_modelos <- list("Lin - Lin" = m1, "Lin - Log" = m2, "Log - Lin" = m3 , "Log - Log" = m4)
lista_residuos <- list()

#Se extraen los residuos para cada modelo
for (i in 1:length(lista_modelos)) {
  modelo <- lista_modelos[[i]]
  residuos <- residuals(modelo)
  lista_residuos[[i]] <- residuos
}

# Realiza el análisis de residuos para cada modelo
for (i in 1:length(lista_residuos)) {
  residuos <- lista_residuos[[i]]
  
  # Gráfico de residuos vs. valores predichos
  plot(predict(lista_modelos[[i]]), residuos, xlab = "Valores Predichos", ylab = "Residuos",
       main = paste("Análisis de Residuos - Modelo", names(lista_modelos)[i]))
  
   # Prueba de Breusch-Pagan para homocedasticidad
  breusch_pagan_test <- lmtest::bptest(lista_modelos[[i]])
  print(paste("Prueba de Breusch-Pagan - Modelo", names(lista_modelos)[i]))
  print(breusch_pagan_test)

  # Prueba de Shapiro-Wilk para normalidad
  shapiro_test <- shapiro.test(residuos)
  print(paste("Prueba de Shapiro-Wilk - Modelo", names(lista_modelos)[i]))
  print(shapiro_test)
  
}

## [1] "Prueba de Breusch-Pagan - Modelo Lin - Lin"
## 
##  studentized Breusch-Pagan test
## 
## data:  lista_modelos[[i]]
## BP = 292.99, df = 1, p-value < 2.2e-16
## 
## [1] "Prueba de Shapiro-Wilk - Modelo Lin - Lin"
## 
##  Shapiro-Wilk normality test
## 
## data:  residuos
## W = 0.96486, p-value < 2.2e-16

## [1] "Prueba de Breusch-Pagan - Modelo Lin - Log"
## 
##  studentized Breusch-Pagan test
## 
## data:  lista_modelos[[i]]
## BP = 214.66, df = 1, p-value < 2.2e-16
## 
## [1] "Prueba de Shapiro-Wilk - Modelo Lin - Log"
## 
##  Shapiro-Wilk normality test
## 
## data:  residuos
## W = 0.95826, p-value < 2.2e-16

## [1] "Prueba de Breusch-Pagan - Modelo Log - Lin"
## 
##  studentized Breusch-Pagan test
## 
## data:  lista_modelos[[i]]
## BP = 150.38, df = 1, p-value < 2.2e-16
## 
## [1] "Prueba de Shapiro-Wilk - Modelo Log - Lin"
## 
##  Shapiro-Wilk normality test
## 
## data:  residuos
## W = 0.99051, p-value = 1.035e-07

## [1] "Prueba de Breusch-Pagan - Modelo Log - Log"
## 
##  studentized Breusch-Pagan test
## 
## data:  lista_modelos[[i]]
## BP = 92.877, df = 1, p-value < 2.2e-16
## 
## [1] "Prueba de Shapiro-Wilk - Modelo Log - Log"
## 
##  Shapiro-Wilk normality test
## 
## data:  residuos
## W = 0.98958, p-value = 2.857e-08

A pesar de que parece visualmente parece cumplirse la homoscedasticidad, la prueba de Breusch-Pagan arroja resultados muy inferiores a 0.05, lo mismo para Shapiro-wilk, por lo que ninguno de los suspuestos se cumple en ninguna transformacion

Informe

  1. Con los resultados obtenidos construya un informe para los directivos de la inmobiliaria, indicando el modelo apropiado y sus principales características. A este informe se deben añadir los anexos como evidencia de la realización de los pasos anteriores.

Resultado de analisis de modelo de regresion lineal para prediccion de precios

Como se puede observar en la comparacion de residuales y supuestos realizada en el punto inmediatamente anterior, ninguno de los modelos o transformaciones son validos para predecir valores de viviendas al no cumplir los supuestos de normalidad y homoscedasticidad.

Basado en los resultados de las pruebas realizadas a lo largo de este reporte, se recomienda utilizar metodos de analisis descriptivo tradicionales para estimar precios de viviendas, o utilizar modelos mas complejos que no dependan de normalidad de valores. El modelo utilizado en este ejercicio, y por ende sus transformaciones son dependientes en varios supuestos que simplemente no se cumplen en el dataset que se tiene disponible, por lo que los resultados que arroje, aunque a simple vista, cercanos, son arbitrarios y no deberian usarse como parte de un analisis o estrategia de mercado.