La base acontinuación trabajada contiene información de ofertas de vivienda descargadas del portal Fincaraiz, con variables como el precio de vivienda(millones de pesos COP) y el area de la vivienda (m2).

library(readxl)
#Leyendo la base
datos <- read_excel("C:/Users/kaqu/Documents/Maestria Ciencia de Datos/Metodos_simulacion_estadistica/Unidad_2/datos_vivienda.xlsx")
attach(datos)
names(datos)
## [1] "Area_contruida" "precio_millon"

1. Realice un análisis exploratorio de las variables precio de vivienda (millones de pesos COP) y area de la vivienda (metros cuadrados) - incluir graficos e indicadores apropiados interpretados.

kbl(summary(datos),caption = "<center><strong>Descriptivas Generales
  </strong></center>")%>%   kable_paper(bootstrap_options = "striped")
Descriptivas Generales
Area_contruida precio_millon
Min. : 80.0 Min. :240.0
1st Qu.: 86.0 1st Qu.:251.2
Median : 97.0 Median :305.0
Mean :115.7 Mean :332.1
3rd Qu.:130.0 3rd Qu.:395.0
Max. :195.0 Max. :480.0

De las anteriores descriptivas podemos concluir que el 50% de las viviendas evaluadas tienen un area construida hasta de 97 m2.Por otro parte, en promedio las viviendas participantes en este estudio tienen una precio promedio 332 millones COP.

par(mfrow = c(2,2))
boxplot(Area_contruida, horizontal = TRUE, main = "Boxplot Area de la Vivienda",col="#A6CEE3" )
hist(Area_contruida,main = "Histograma Area de la Vivienda" , ylab = "Frecuencia" , xlab = "Area de la vivienda (m2)",col="#A6CEE3")
boxplot(precio_millon, horizontal = TRUE, main = "Boxplot Precio de Vivienda",col="#A6CEE3" )
hist(precio_millon,main = "Histograma Precio de Vivienda" , ylab = "Frecuencia" , xlab = "Precio de Vivienda (Millones)",col="#A6CEE3")

En cuánto a lo que nos permite ver el boxplot y el histograma realizado para ambas variables son variables que no se consideran simétricas, al contrario tienen una distribución asimétricas positivas, lo que infica que tanto en el precio de vivienda como en la area construida hay presencia de valores atípicos. La mayoria de viviendas evaluadas tienen un área entre 80 y 100 m2.El 50% de las viviendas tienen un precio alrededor de los 300 millones COP.

2. Realice un análisis exploratorio bivariado de datos enfocado en la relación entre la variable respuesta (y=precio) en función de la variable predictora (x=area) - incluir graficos e indicadores apropiados interpretados.

ggplot(data = datos, aes(x=Area_contruida, y=precio_millon)) +
  #geom_histogram(aes(y = ..count.., fill = ..count..)) +
  geom_point()+
  #geom_smooth(method=lm, se=FALSE)+
  scale_fill_gradient(low = "light blue", high = "light blue") +
  #stat_function(fun = dnorm, colour = "black", args = list(mean = mean(precio_millon), sd = sd(precio_millon))) +
  ggtitle("Diagrama de Dispersión ") + labs(x = "Área Construida (m2)") + labs(y = "Precio Vivienda COP")

Del anterior gráfico se puede observar que existe una relación en algunos casos positiva entre el área construida y el precio de la vivienda, esto se puede deducir al evidenciar la tendencia creciente, de que a medida que aumenta el área construida tambien lo hace el precio de vivienda en la mayoría de los casos.

Para poder cuantificar esa relación observada anteriormente es importante calcular medidas como por ejemplo el coeficiente de correlación. Existen otro tipos de coeficientes de correlación, pero en este caso como se trata de poder evidenciar una “relación lineal” utilizaré el coeficiente de correlación de pearson.

round(cor(Area_contruida,precio_millon,method = "pearson"),2)
## [1] 0.92

El coeficiente de correlación (R) nos da un valor cercano a 1, lo que indica que hay una relación positiva relativamente fuerte entre el precio de las viviendas y el area construida.

3. Estime el modelo de regresión lineal simple entre precio = f(area) +e.Interprete los coeficientes del modelo β0, β1 en caso de ser correcto. 

m1=lm(precio_millon~Area_contruida)
summary(m1)
## 
## Call:
## lm(formula = precio_millon ~ Area_contruida)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -51.673 -25.612  -6.085  24.875  67.650 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      86.234     22.479   3.836 0.000796 ***
## Area_contruida    2.124      0.186  11.422 3.45e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 33.05 on 24 degrees of freedom
## Multiple R-squared:  0.8446, Adjusted R-squared:  0.8381 
## F-statistic: 130.5 on 1 and 24 DF,  p-value: 3.45e-11

Para poder suponer que el modelo es válido, tendriamos que hacer todo el ejercicio de validación de supuestos, pero hace parte de un punto más adelante.

Ahora bien, para hablar de los coeficientes arrojados por el modelo en cuanto si es correcto interpretarlos o no, podriamos observar los p-values asociados a cada uno de ellos. En ambos casos son menores a cualquier nivel de significancia habitualmente trabajado (0.01,0.05 o 0.1), lo que sugiere que son relevantes en el modelo (diferentes de 0).

Para evaluar la validez del modelo (Paramétros en conjunto diferentes a 0), se observa que el p valor asociado al estadístico F es menor a 0.05, por ende los parametros son diferentes de 0.

4. 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 hipotesis t. 

kbl(confint(m1),caption = "<center><strong>Intervalo de Confianza B1
  </strong></center>")%>%   kable_paper(bootstrap_options = "striped")
Intervalo de Confianza B1
2.5 % 97.5 %
(Intercept) 39.839830 132.627917
Area_contruida 1.740169 2.507772

Se puede observar en la anterior tabla que el intervalo de confianza para \(\beta_1\) no contiene al 0, por ende se considera que si esta variable si es relevante dentro del modelo y además por ser valores mayores a 0, hay un impacto positivo del área en el precio de la vivienda.

Otra forma de corroborar lo anteriormente encontrado, procedemos a realizar una prueba t que sigue la siguiente hipotésis:

\(H_0 :\beta_1=0\)

\(H_a :\beta_1\neq 0\)

gl=nrow(datos)-2
t_calculada= 10.36
t_tablas=qt(0.025,gl,lower.tail = FALSE)


resultado=data.frame(t_calculada,t_tablas)
resultado
##   t_calculada t_tablas
## 1       10.36 2.063899

Se Puede observar que \(|t_0| > t_{\alpha/2 ,n-2}\), entonces se rechaza \(H_0\), y se dice que \(\beta_1\) es diferente de 0.

5.Calcule e interprete el indicador de bondad y ajuste R2 

coef_r2=round(cor(Area_contruida,precio_millon)*cor(Area_contruida,precio_millon),2)
coef_r2
## [1] 0.84

El coeficiented de determinación representa la proporción de la variabilidad de Y que es posible explicar a travez de x. En este caso el modelo construido explica el 84% de las variaciones del precio de las viviendas a travez del area construida.

6. Cual seria el precio promedio estimado para un apartamento de 110 metros cuadrados? Considera entonces con este resultado que un apartemento en la misma zona con 110 metros cuadrados en un precio de 200 millones seria una buena oferta? Que consideraciones adicionales se deben tener?. 

round(predict(m1,newdata = list(Area_contruida=110)),2)
##      1 
## 319.87

El precio promedio estimado para un apartamento de 110 m2, seria de casi 320 millones COP. Ahora bien para saber si un apartamento con el mismo área cuadrado es una mejor opción se deben tener en cuenta otros factores adicionales como por ejemplo la zona, si cuenta con opción de parqueadero o no, entre otros factores, que al ser un modelo de regresión lineal simple no se estan teniendo en cuenta.

7. Realice la validación de supuestos del modelo por medio de graficos apropiados, interpretarlos y sugerir posibles soluciones si se violan algunos de ellos.

Interpretación Gráfica
par(mfrow=c(2,2))
plot(m1)

Se puede observar que en el gráfico de residuales versus valores ajustados hay una presencia de valores extremos con el dato 4,21 y 25 lo que no hace que del todo se cumpla que el valor esperado este cerca de 0 en todos los casos. En cuanto al QQ plot se puede obser que tiene una tendencia lineal podría pensar en normalidad, aún existiendo valores extremos, pero debería corroborarse con una prueba estadística.

Pruebas estadísticas para validación de supuestos
  1. \(E(e_i )= 0\)
mean(residuals(m1))
## [1] -3.760347e-16

El valor esperado (promedio) de los residuales es igual a 0.

  1. Supuesto de Normalidad

\(H_0 : X \sim N(\mu,\sigma^2)\)

\(H_a : X \nsim N(\mu,\sigma^2)\)

shapiro.test(m1$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  m1$residuals
## W = 0.95489, p-value = 0.3009

Com $ P-value $ es mayor a a 0.05 (nivel de significancia escogido), no se rechaza \(H_0\), entonces podría pensar que los errores siguen una distribución normal.

  1. Supuesto de Homocedasticidad de Varianza (Breush Pagan)

\(H_0 : Los \, residuales \, se \, distribuyen \, con \, la \, misma \, varianza\)

\(H_a : Los \, residuales \, no \ se \, distribuyen \, con \, la \, misma \, varianza\)

bptest(m1)
## 
##  studentized Breusch-Pagan test
## 
## data:  m1
## BP = 5.8737, df = 1, p-value = 0.01537

Como $ P-value $ es menor a 0.05 (nivel de significancia escogido), se rechaza \(H_0\), entonces se podría pensar que los errores no cumplen con el supuesto de homocedasticidad.

  1. Supuesto de Autocorrelación de los errores (Durbin-Watson)

\(H_0 : No \, existe \, correlación \, entre \, los \, errores\)

\(H_a : Existe \, correlación\, entre\, los\, errores\)

dwt(m1, alternative = "two.sided")
##  lag Autocorrelation D-W Statistic p-value
##    1      0.02511179      1.883052   0.786
##  Alternative hypothesis: rho != 0

Como $ P-value $ es mayor a 0.05 (nivel de significancia escogido), no se rechaza \(H_0\), entonces se podría pensar que los errores no estan autocorrelacionados.

En este caso se viola el supuesto de homocedasticidad, por ende se podría pensar en realizar alguna transformación sobre la variable respuesta o inclunsive sobre la variable explicativa.

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

El supuesto de varianza constante es una premisa básica dentro de un análisis de regresión. Podría hacer alguna transformación sobre la variable respuesta (precio de la vivienda) como por ejemplo el logaritmo, sugiero utilizar un método analitico para seleccionar la transformación adecuada conocida como la transformación de Box-Cox que no es mas que la “potencia ideal” de la variable respuesta, conocido como \(\lambda\), para solucionar problemas de no normalidad y/o de varianza no constante.

\(\log(y) \, si \, \lambda =0\)

bc=boxcox(precio_millon~Area_contruida)

lambda=bc$x[which.max(bc$y)]
lambda
## [1] 1.474747

Ahora bien como \(\lambda\) es diferente de 0, entonces se procede a usar la transformación de Box-Cox.

Modelo Sugerido

m2=lm((((precio_millon^lambda)-1)/lambda)~Area_contruida)
summary(m2)
## 
## Call:
## lm(formula = (((precio_millon^lambda) - 1)/lambda) ~ Area_contruida)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -842.74 -371.84  -84.08  333.51 1086.54 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    -321.366    346.416  -0.928    0.363    
## Area_contruida   34.015      2.866  11.869 1.57e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 509.3 on 24 degrees of freedom
## Multiple R-squared:  0.8544, Adjusted R-squared:  0.8484 
## F-statistic: 140.9 on 1 and 24 DF,  p-value: 1.567e-11

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

Coeficiente de Determinación
coef_r2=round(cor(Area_contruida,(((precio_millon^lambda)-1)/lambda))*cor(Area_contruida,(((precio_millon^lambda)-1)/lambda)),2)
coef_r2
## [1] 0.85
Supuestos Modelo 2: (Transformado)
#define plotting area
op <- par(pty = "s", mfrow = c(1, 2))

#Q-Q plot for original model
qqnorm(m1$residuals,main = "Normal Q-Q Plot Modelo 1")
qqline(m1$residuals)

#Q-Q plot for Box-Cox transformed model
qqnorm(m2$residuals,main = "Normal Q-Q Plot Modelo 2(Transf.)")
qqline(m2$residuals)

#display both Q-Q plots
par(op)

Al parecer hay un mejor ajuste del modelo para los puntos lejanos, en el modelo transfromado, es por esto que se harán pruebas estadísticas que respalden esta hipotésis.

  1. \(E(e_i )= 0\)
mean(residuals(m2))
## [1] -1.613239e-14

El valor esperado (promedio) de los residuales es igual a 0.

  1. Supuesto de Normalidad

\(H_0 : X \sim N(\mu,\sigma^2)\)

\(H_a : X \nsim N(\mu,\sigma^2)\)

shapiro.test(m2$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  m2$residuals
## W = 0.96453, p-value = 0.4885

Com $ P-value $ es mayor a a 0.05 (nivel de significancia escogido), no se rechaza \(H_0\), entonces podría pensar que los errores siguen una distribución normal.

  1. Supuesto de Homocedasticidad de Varianza (Breush Pagan)

\(H_0 : Los \, residuales \, se \, distribuyen \, con \, la \, misma \, varianza\)

\(H_a : Los \, residuales \, no \ se \, distribuyen \, con \, la \, misma \, varianza\)

bptest(m2)
## 
##  studentized Breusch-Pagan test
## 
## data:  m2
## BP = 7.181, df = 1, p-value = 0.007368

Com $ P-value $ es menor a 0.05 (nivel de significancia escogido), se rechaza \(H_0\), entonces se podría pensar que los errores no cumplen con el supuesto de homocedasticidad.

  1. Supuesto de Autocorrelación de los errores (Durbin-Watson)

\(H_0 : No \, existe \, correlación \, entre \, los \, errores\)

\(H_a : Existe \, correlación\, entre\, los\, errores\)

dwt(m2, alternative = "two.sided")
##  lag Autocorrelation D-W Statistic p-value
##    1      0.02666335      1.871755   0.726
##  Alternative hypothesis: rho != 0

Com \(P-value\) es mayor a 0.05 (nivel de significancia escogido), no se rechaza \(H_0\), entonces se podría pensar que los errores no estan autocorrelacionados.

Finalmente se puede observar que la explicabilidad del modelo después de la transformación, mejora en un punto pero sigue habiendo problemas con el tema de los supuestos del modelo en cuanto a la varianza constante, por lo cual recomendaría que se debería depronto hacer alguna modificación a la variable explicativa, para ver si se puede mejorar este acercamiento.

Puede considerarse además que la correlación existente, no justamente debe ser lineal, por lo que debe contemplarse la posibilidad de evaluar un modelo polinomico.