Problema

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.

Con este propósito el equipo de asesores a diseñado los siguientes pasos para obtener un modelo y así poder a futuro determinar los precios de los inmuebles a negociar.

#install.packages('devtools')
#devtools::install_github('dgonxalex80/paqueteMETODOS')
library(paqueteMETODOS)
data(vivienda4)

Paso 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.

precio <- vivienda4$preciom
area <- vivienda4$areaconst

Precio

Resumen numérico

describe(precio)
##    vars    n   mean    sd median trimmed   mad min max range skew kurtosis   se
## X1    1 1706 225.37 85.89    210   214.4 75.61  78 760   682 1.49     3.57 2.08

Esta base de datos contiene 1706 casas, el precio promedio de estas casas es de 225.37 millones de pesos, con una desviación estándar de 85.89 millones de pesos, lo que sugiere que los precios varian considerablemente alrededor de la media. El 50% de las 1706 casas tienen un precio de 210 millones de pesos o menos, mientras que el otro 50% tiene un precio de más de 210 millones de pesos. La casa con el menor y mayor precio cuestan 78 y 760 millones de pesos, respectivamente; obteniendo así un rango de precios de 682 millones de pesos. El coeficiente de asimetría de 1.49 nos indica una asimetría positiva y que los precios se acumulan más a la izquierda, es decir, en precios bajos; y, el coeficiente de curtosis de 3.57 indica leptocurtosis, lo que quiere decir que la distribución de los precios de las casas es en punta y tiene colas poco anchas.

Histograma y gráfico de caja

par(mfrow=c(1,2))
hist(precio, main = 'Precio de viviendas (millones de pesos)', xlab='Precio', ylab='Frecuencia', col = 'lightblue')
boxplot(precio, main='Precio de viviendas (millones de pesos)', ylab='Precio', col = 'lightgreen')

El histograma confirma lo mencionado sobre la forma de la distribución de los precios y a eso le agregamos que la distribución presentan un sesgo a la derecha. Por otro lado, se evidencia en el gráfico de cajas que el precio presenta valores atípicos, es decir, valores de precios muy extremos y que solo se dan en unas pocas viviendas, esto es lo que causa el sesgo de la distribución.

Área

Resumen numérico

describe(area)
##    vars    n  mean    sd median trimmed   mad min max range skew kurtosis   se
## X1    1 1706 87.63 36.35     75   80.97 22.24  40 200   160 1.53     1.68 0.88

El área promedio de estas casas es de 87.63 metros cuadrados, con una desviación estándar de 36.35 metros cuadrados, lo que sugiere que las áreas varian considerablemente alrededor de la media. El 50% de las 1706 casas tienen un área de 75 metros cuadrados o menos, mientras que el otro 50% tiene un área de más de 75 metros cuadrados. La casa con el menor y mayor área cuentan con 40 y 200 metros cuadrados, respectivamente; obteniendo así un rango de área de 160 metros cuadrados. El coeficiente de asimetría de 1.53 nos indica una asimetría positiva y que las áreas se acumulan más a la izquierda, es decir, en áreas bajas; y, el coeficiente de curtosis de 1.67 indica leptocurtosis, lo que quiere decir que la distribución de las áreas de las casas es en punta y tiene colas poco anchas.

Histograma y gráfico de caja

par(mfrow=c(1,2))
hist(area, main = 'Área de viviendas (m^2)', xlab='Área', ylab='Frecuencia', col = 'lightblue')
boxplot(area, main='Área de viviendas (m^2)', ylab='Área', col = 'lightgreen')

Al igual que con la variable de precio, el histograma de la variable de área confirma lo mencionado sobre la forma de la distribución de las áreas y a eso le agregamos que la distribución presentan un sesgo a la derecha, pero esta vez con frecuencias mayores en el sesgo. Por otro lado, se evidencia en el gráfico de cajas que el precio presenta valores atípicos, es decir, valores de precios muy extremos y que solo se dan en unas pocas viviendas, esto es lo que causa el sesgo de la distribución.

Paso 2

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.

Coeficiente de correlación de Pearson

cor(area, precio)
## [1] 0.7630166

El coeficiente de correlación de Pearson de \(0.763\) indica que la relación entre el área de las viviendas (x) y el precio de estas (y) es alta y es directamente proporcional, es decir, entre mayor sea el área de una vivienda, tiende a se mayor el precio. Veamos esto de manera gráfica.

Gráfico de dispersión

plot(area, precio, main='Área (m^2) vs Precio (Millones de pesos)', xlab='Área', ylab='Precio', col='#88aaff')

Podemos confirmar con el gráfico de dispersión lo mencionado sobre la correlación entre el área y el precio de las viviendas; entre mayor es el área, tiende a ser mayor el precio.

Paso 3

Estime el modelo de regresión lineal simple entre \(precio = f(area) + \varepsilon\). Interprete los coeficientes del modelo \(\beta_0\), \(\beta_1\) en caso de ser correcto.

Estimar modelo de regresión

modelo1 <- lm(precio ~ area)
summary(modelo1)
## 
## Call:
## lm(formula = precio ~ area)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -195.86  -31.95   -8.95   27.87  431.17 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   67.381      3.510   19.20   <2e-16 ***
## area           1.803      0.037   48.73   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 55.53 on 1704 degrees of freedom
## Multiple R-squared:  0.5822, Adjusted R-squared:  0.5819 
## F-statistic:  2374 on 1 and 1704 DF,  p-value: < 2.2e-16

Intercepto (\(\beta_0\)): El precio que se espera tener cuando una vivienda tiene área cero es de 67.381 millones de pesos. Lo anterior no tiene sentido o no es significativo en este caso, ya que no ha viviendas de área igual a cero.

Pendiente (\(\beta_1\)): Por cada unidad que aumente el área de una vivienda, el precio aumenta en promedio 1.803 millones de pesos.

Paso 4

Construir un intervalo de confianza (\(95\ \%\)) para el coeficiente \(\beta_1\), interpretar y concluir si el coeficiente es igual a cero o no. Compare este resultado con una prueba de hipótesis t.

Intervalo de confianza del \(95\ \%\) para \(\beta_1\)

confint(modelo1, level = 0.95)
##                 2.5 %    97.5 %
## (Intercept) 60.496208 74.265055
## area         1.730404  1.875547

El intervalo de confianza para \(\beta_1\) no contiene el cero. Por tanto, con una confianza del 95%, no hay evidencias para afirmar que el coeficiente \(\beta_1\) es igual a cero y se concluye que el área de las viviendas sí tiene un efecto significativo sobre el precio de estas.

Prueba t para \(\beta_1\)

\[H_0:\beta_1 = 0\] \[H_1:\beta_1 \neq 0\]

El resumen del modelo de regresión nos proporciona los estadísticos \(t\) y los \(p-valores\) de cada coefiente de la regresión.

summary(modelo1)$coefficients
##              Estimate Std. Error  t value     Pr(>|t|)
## (Intercept) 67.380631 3.51003056 19.19659 1.648523e-74
## area         1.802976 0.03700055 48.72834 0.000000e+00

Como el p-valor de el coeficiente \(\beta_1\) (área) es menor que el nivel de significancia de 0.05, se rechaza \(H_0\), por lo que no hay evidencias para afirmar que el coeficiente \(\beta_1\) es igual a cero y se concluye que el área de las viviendas sí tiene un efecto significativo sobre el precio de estas.

Tanto el intervalo de confianza del \(95\ \%\) como la prueba \(t\) nos llevaron a la misma conlcusión: \(\beta_1\) sí es significativo.

Paso 5

Calcule e interprete el indicador de bondad \(R^2\).

summary(modelo1)$r.squared
## [1] 0.5821944

El modelo de regresión se ajusta un \(58.22\ \%\) a nuestros datos, es decir, \(58.22\ \%\)% de la variabilidad del precio está siendo explicada por el modelo de regresión. Es un valor regular, ya que entre más cercano a 1 mejor. Por ahora, esto puede atribuirse a los valores atípicos y al incumplimiento de un supuesto, más adelante revisaremos esto.

Paso 6

¿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?.

\[Modelo \ : Precio=67.380631+1.802976(Área)\] \[Modelo \ : \hat{y}=67.380631+1.802976(x)\] \[Si \ Área(x)=110\ m^2\] \[ Precio = 67.380631+1.802976(110\ m^2)\]

\[ Precio = 265.708\ millones\ de\ pesos\] Un apartamento en la misma zona con 110 metros cuadrados en un precio de 200 millones sí sería una atractiva esta oferta. Hay que tener en cuenta que el modelo no es tan bueno y que las predicciones hechas con este pueden estar muy lejos de la realidad, se deben validar los supuestos del modelo y mejorar su rendimiento para tener seguridad de que es factible hacer predicciones con este.

Paso 7

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.

Obtener residuales y valores ajustados

ajustados <- modelo1$fitted.values
residuos <- modelo1$residuals

Cargar librerías para pruebas Estadísticas:

#install.packages("lmtest")
library(lmtest)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
#install.packages("car")  # Instala el paquete si aún no lo has hecho
library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
## The following object is masked from 'package:purrr':
## 
##     some
## The following object is masked from 'package:psych':
## 
##     logit

Linealidad

plot(area, precio, main='Área (m^2) vs Precio (Millones de pesos)', xlab='Área', ylab='Precio', col='#88aaff')
abline(modelo1, col = "red")

Como ya se había visto en el análisis exploratorio, si existe una relación lineal entre el área y el precio de las viviendas.

Normalidad

\[H_0 : Los\ errores\ siguen\ una\ distribución\ normal\] \[H_1 : Los\ errores\ no\ siguen\ una\ distribución\ normal\] Prueba gráfica:

par(mfrow=c(1,2))
hist(residuos, main="Histograma de residuales", xlab="Residuos", ylab="Frecuencia", col="#3344aa")
qqnorm(residuos, main="QQplot de residuales", xlab="Cuantiles teóricos", ylab="Cuantile de muestra", col="#44cc1188")
qqline(residuos, col = "red", lwd="1")

Prueba estadística Shapiro-Wilk:

shapiro.test(residuos)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuos
## W = 0.92671, p-value < 2.2e-16

El p-valor es menor al nivel de significancia de 0.05, por tanto, se rechaza \(H_0\) y se concluye que hay evidencias para afirmar que lo residuos no siguen una distribución normal. Además, el histograma muestra una cola derecha bastante ancha, lo que indica una fuerte asimetría positiva; y, en el qqplot los puntos no se ajustan a la línea recta, lo cual apoya la no normalidad de los residuales. Esto puede deberse a la presencia de valores atípicos tanto en la variable predictora como en la variable objetivo, imputar o eliminar los outliers muy extremos o aplicar transformaciones podrían ayudarnos a solucionar el incumplimiento de este supuesto.

Homocedasticidad

\[H_0: Los\ errores\ tiene\ varianza\ constante\ (Homogeneidad)\] \[H_1: Los\ errores\ no\ tienen\ varianza\ constante\ (Heterocedasticidad)\] Prueba gráfica:

plot(ajustados, residuos, main="Ajustados vs Residuales", xlab="Valores ajustados", ylab="Residuales", col="#838bc7")
abline(0,0, col="#d3122f")

Prueba estadística Breusch-Pagan:

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

El gráfico de dispersión de los valores ajustados vs los residuos muestra un patrón de cono o embudo, lo que indica que la dispersión de los residuos cambia a lo largo de los valores ajustados, esto sugiere heterocedasticidad. Además, dado que el p-valor de la prueba Breusch-Pagan es menor al nivel de significancia de 0.05, se rechaza la hipótesis nula, por lo que podemos afirmar que hay evidencias de heterocedasticidad. Los outliers también podrían estar afectando la varianza de los errores. Como el modelo es lineal, la recta se encuentra lejos de precios muy extremos y esto hace más grande el error.

Independencia

\[H_0:No\ existe\ correlación\ serial\ entre\ los\ residuales\] \[H_1:Existe\ correlación\ serial\ entre\ los\ residuales\] Prueba gráfica:

indices <- row.names(vivienda4)
plot(vivienda4$preciom, residuos, xlab="Precio", ylab="Residuos", main = "Precio vs Residuales", col="#670d12")

Prueba estadística Durbin-Watson:

durbinWatsonTest(modelo1)
##  lag Autocorrelation D-W Statistic p-value
##    1       0.1635281      1.671255       0
##  Alternative hypothesis: rho != 0

Cómo el p-valor en la prueba Durbin-Watson es menor al nivel de significancia de 0.05, se rechaza la hipótesis nula y hay evidencias para afirmar que los residuos presentan autocorrelación. Además, en el gráfico de dispersión de los residuos vs la variable objetivo, se observa una dependencia lineal positiva. Igualmente, los valores muy extremos también podrían estar afectando la independencia de los errores. Aunque la relación entre precio y área es lineal, la varianza del precio se hace más grande a medida que el área crece, por tanto, los errores también crecen.

El modelo de regresión lineal simple no cumple con los supuestos de normalidad, homocedasticidad e independencia. Por lo que se cuestiona la validez y precisión del modelo para predecir el precio de las viviendas con el área, lo que puede llevar a resultados sesgados y conclusiones poco confiables. Por tanto, antes de explorar una alternativa que requiera tratar los valores atípicos, se explora la alternativa de aplicar una transformación logarítmica a la variable independiente.

Paso 8

De ser necesario realice una transformación apropiada para mejorar el ajuste y supuestos del modelo.

Transformación logarítmica a la variable predictora o independiente (área).

Modelo con log(área):

area_log <- log(area)
modelo2 <- lm(precio ~ area_log)
summary(modelo2)
## 
## Call:
## lm(formula = precio ~ area_log)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -192.72  -27.27   -3.56   23.58  419.52 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -610.083     16.085  -37.93   <2e-16 ***
## area_log     189.708      3.641   52.11   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 53.35 on 1704 degrees of freedom
## Multiple R-squared:  0.6144, Adjusted R-squared:  0.6142 
## F-statistic:  2715 on 1 and 1704 DF,  p-value: < 2.2e-16

Paso 9

De ser necesario compare el ajuste y supuestos del modelo inicial y el transformado.

Este modelo de regresión se ajusta un \(61.44\ \%\) a nuestros datos, es decir, \(61.44\ \%\) de la variabilidad del precio está siendo explicada por el modelo de regresión. Esto es mejor que el \(58.22\ \%\), pero sigue siendo un valor regular.

Gráfico del modelo con log(area):

ajustados2 <- modelo2$fitted.values
residuos2 <- modelo2$residuals

plot(area_log, precio, main='Log(Área (m^2)) vs Precio (Millones de pesos)', xlab='Log(Área)', ylab='Precio', col='#88aaff')
abline(modelo2, col = "red")

Validación de supuestos

Obervar si el nuevo modelo con log(area) cumple los supuestos:

par(mfrow=c(2,2))

hist(residuos2, main="Histograma de residuales", xlab="Residuos", ylab="Frecuencia", col="#3344aa")
qqnorm(residuos2, main="QQplot de residuales", xlab="Cuantiles teóricos", ylab="Cuantile de muestra", col="#44cc1188")
qqline(residuos2, col = "red", lwd="1")

plot(ajustados2, residuos2, main="Ajustados vs Residuales", xlab="Valores ajustados", ylab="Residuales", col="#838bc7")
abline(0,0, col="#d3122f")

plot(precio, residuos2, xlab="Precio", ylab="Residuos", main = "Precio vs Residuales", col="#670d12")

shapiro.test(residuos2)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuos2
## W = 0.91387, p-value < 2.2e-16
bptest(modelo2)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo2
## BP = 146.29, df = 1, p-value < 2.2e-16
durbinWatsonTest(modelo2)
##  lag Autocorrelation D-W Statistic p-value
##    1        0.153554       1.69083       0
##  Alternative hypothesis: rho != 0

Se evidencia, gráficamente y con por los p-valores menores a 0.05, que el modelo con log(area) tampoco cumple con los supuestos de normalidad, homocedasticidad e independencia de los errores.

Paso 10

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

Se obta por probar otro modelo después de haber tratado los valores atípicos con el método del rango intercuartílico y con el logaritmo de la variable Área.

Excluir valores atípicos:

Q1_area <- quantile(area, 0.25)
Q3_area <- quantile(area, 0.75)

IQR_area <- Q3_area - Q1_area

lim_inf_area <- Q1_area - 2 * IQR_area
lim_sup_area <- Q3_area + 2 * IQR_area


Q1_precio <- quantile(precio, 0.25)
Q3_precio <- quantile(precio, 0.75)

IQR_precio <- Q3_precio - Q1_precio

lim_inf_precio <- Q1_precio - 1.5 * IQR_precio
lim_sup_precio <- Q3_precio + 1.5 * IQR_precio

datos <- vivienda4[,c("preciom", "areaconst")]

datos2 <- datos[datos$areaconst > lim_inf_area & datos$areaconst < lim_sup_area,]
datos_sin_outliers <- datos2[datos2$preciom > lim_inf_precio & datos2$preciom < lim_sup_precio,]

area2 <- datos_sin_outliers$areaconst
precio2 <- datos_sin_outliers$preciom

plot(area2, precio2)

Ajustar nuevo modelo:

modelo3 <- lm(precio2 ~ log(area2))
summary(modelo3)
## 
## Call:
## lm(formula = precio2 ~ log(area2))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -176.799  -24.754   -3.389   23.358  198.773 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -534.616     16.212  -32.98   <2e-16 ***
## log(area2)   171.701      3.722   46.14   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 43.35 on 1586 degrees of freedom
## Multiple R-squared:  0.573,  Adjusted R-squared:  0.5728 
## F-statistic:  2128 on 1 and 1586 DF,  p-value: < 2.2e-16

Graficar modelo:

ajustados3 <- modelo3$fitted.values
residuos3 <- modelo3$residuals

plot(area2, precio2, main='Log(Área (m^2)) vs Precio (Millones de pesos)', xlab='Log(Área)', ylab='Precio', col='#88aaff')
abline(modelo3, col = "red")

Validación de supuestos

Obervar si el nuevo modelo sín aoutliers y con log(area) cumple los supuestos:

par(mfrow=c(2,2))

hist(residuos3, main="Histograma de residuales", xlab="Residuos", ylab="Frecuencia", col="#3344aa")
qqnorm(residuos3, main="QQplot de residuales", xlab="Cuantiles teóricos", ylab="Cuantile de muestra", col="#44cc1188")
qqline(residuos3, col = "red", lwd="1")

plot(ajustados3, residuos3, main="Ajustados vs Residuales", xlab="Valores ajustados", ylab="Residuales", col="#838bc7")
abline(0,0, col="#d3122f")

plot(precio2, residuos3, xlab="Precio", ylab="Residuos", main = "Precio vs Residuales", col="#670d12")

shapiro.test(residuos3)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuos3
## W = 0.97859, p-value = 1.202e-14
bptest(modelo3)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo3
## BP = 148.67, df = 1, p-value < 2.2e-16
durbinWatsonTest(modelo3)
##  lag Autocorrelation D-W Statistic p-value
##    1        0.209104      1.578422       0
##  Alternative hypothesis: rho != 0

Se evidencia, gráficamente y con por los p-valores menores a 0.05, que el modelo sin outliers y con log(area) tampoco cumple con los supuestos de normalidad, homocedasticidad e independencia de los errores.

Conclusiones

  • El sesgo que presentan tanto el precio como el área, no permite ajustar un modelo lineal confiable.

  • Ninguno de los modelos ajustado cumple con los supuestos de normalidad, homocedasticidad e independencia de los errores.

  • El mejor modelo es el modelo2, sin eliminación de outliers y la tranformación logarítmica de la variable independiente. Este modelo de regresión se ajusta un 61.44 % a nuestros datos, es decir, 61.44 % de la variabilidad del precio está siendo explicada por el modelo de regresión. Esto es mejor que el 58.22 % y el 57.3 % de los modelos 1 y 3, respectivamente.

  • De los 3 modelos ajustados, se recomienda usar el modelo2, ya que fue el que dio el mayor porcentaje de explicabilidad. Pero, con precausión porque es muy propenso a dar errores altos en el precio con áreas muy grandes.