Modelación EstadÃstica Actividad 3
Ruben Dario Rodriguez Camargo
2023-10-14
library(lmtest)## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(paqueteMETODOS)## Loading required package: cubature## Loading required package: GGally## Loading required package: ggplot2## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2## Loading required package: MASS## Loading required package: summarytools## Loading required package: psych##
## Attaching package: 'psych'## The following objects are masked from 'package:ggplot2':
##
## %+%, alpha## Loading required package: tidyverse## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.2 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ lubridate 1.9.2 ✔ tibble 3.2.1
## ✔ purrr 1.0.1 ✔ tidyr 1.3.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ psych::%+%() masks ggplot2::%+%()
## ✖ psych::alpha() masks ggplot2::alpha()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ✖ dplyr::select() masks MASS::select()
## ✖ tibble::view() masks summarytools::view()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(nortest)
library(stargazer)##
## Please cite as:
##
## Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary Statistics Tables.
## R package version 5.2.3. https://CRAN.R-project.org/package=stargazerlibrary(foreign)Modelo de regresión lineal
Análisis del mercado inmobiliario de la empresa A&C
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
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.
Analisis incial:
La base de datos analizada cuenta con 1706 registros, identificando la Zona de ubicación, estrato socioeconómico, precio descrito por millones, el área de construida por metro cuadrado y el tipo de vivienda.
Podemos evidenciar que la base de datos no cuenta con campos vacion y una adecuada calidad para el analisis de infromacion
data(vivienda4)
pie(table(vivienda$tipo), main="Distribuciones por tipo de vivienda")
El 79.89% de la inforamcion es de tipo vivienda apartamento, y el 100% de las viviendas se encuentran en estrato 4.
barplot(table(vivienda$zona), main="Distribucion de viviendas por zonas" )
summary(vivienda)## id zona piso estrato
## Min. : 1 Length:8322 Length:8322 Min. :3.000
## 1st Qu.:2080 Class :character Class :character 1st Qu.:4.000
## Median :4160 Mode :character Mode :character Median :5.000
## Mean :4160 Mean :4.634
## 3rd Qu.:6240 3rd Qu.:5.000
## Max. :8319 Max. :6.000
## NA's :3 NA's :3
## preciom areaconst parqueaderos banios
## Min. : 58.0 Min. : 30.0 Min. : 1.000 Min. : 0.000
## 1st Qu.: 220.0 1st Qu.: 80.0 1st Qu.: 1.000 1st Qu.: 2.000
## Median : 330.0 Median : 123.0 Median : 2.000 Median : 3.000
## Mean : 433.9 Mean : 174.9 Mean : 1.835 Mean : 3.111
## 3rd Qu.: 540.0 3rd Qu.: 229.0 3rd Qu.: 2.000 3rd Qu.: 4.000
## Max. :1999.0 Max. :1745.0 Max. :10.000 Max. :10.000
## NA's :2 NA's :3 NA's :1605 NA's :3
## habitaciones tipo barrio longitud
## Min. : 0.000 Length:8322 Length:8322 Min. :-76.59
## 1st Qu.: 3.000 Class :character Class :character 1st Qu.:-76.54
## Median : 3.000 Mode :character Mode :character Median :-76.53
## Mean : 3.605 Mean :-76.53
## 3rd Qu.: 4.000 3rd Qu.:-76.52
## Max. :10.000 Max. :-76.46
## NA's :3 NA's :3
## latitud
## Min. :3.333
## 1st Qu.:3.381
## Median :3.416
## Mean :3.418
## 3rd Qu.:3.452
## Max. :3.498
## NA's :3El 78.78% de las viviendas se encuentran ubicadas en la zona sur de la ciudad, el restante 21.22% se encuentran distribuidas en la zon a norte con un 16.88% y el otro 4.34% se ditribuyen entre la zona oeste, centro, y oriente de Cali.
El precio promedio de la viviendas es de 210 millones, con una area promedio de 87 metros
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.
Apartamento=subset(vivienda4, vivienda4$tipo=="Apartamento")
Casa=subset(vivienda4, vivienda4$tipo=="Casa")
Apartamento|> summarise(media_preciom = mean(preciom),mediana_preciom = median(preciom),varianza_preciom = var(preciom),desvi_preciom = sd(preciom),Q1 = quantile(preciom, probs=0.25),P90 = quantile(preciom, probs=0.90))## # A tibble: 1 × 6
## media_preciom mediana_preciom varianza_preciom desvi_preciom Q1 P90
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 202. 185 4263. 65.3 154. 285Casa|> summarise(media_preciom = mean(preciom),mediana_preciom = median(preciom),varianza_preciom = var(preciom),desvi_preciom = sd(preciom),Q1 = quantile(preciom, probs=0.25),P90 = quantile(preciom, probs=0.90))## # A tibble: 1 × 6
## media_preciom mediana_preciom varianza_preciom desvi_preciom Q1 P90
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 317. 300 9368. 96.8 255 440Al analizar la informacion podemos evidenciar que el precio promedio de los apartamentos es de 202 millones y un desviacion estandar de 65 millones, en contraste a las casas con una mediana en el precio de 316 millones con una desviacion estandar de 96 millones, siendo los apartamentos mas economicos por metro cuadrado.
ggplot(vivienda4, aes(x = areaconst, fill = tipo)) + geom_histogram (bins=30) + labs (title = "Distribucion area construida")
ggplot(vivienda4, aes(x = preciom, fill = tipo)) + geom_histogram (bins=30) + labs (title = "Distribucion de precios")
boxplot (vivienda4$preciom~vivienda4$tipo, main="Precio en millones (COP) segun tipo de vivienda")
En la grafica se evidencia que la distribucion de apartamentos en mayor proporcion por area construida se encuentran en el rango de 50 a 100 metros cuadrados y en contraste de las casa que estan distribuidas desde los 50 a 200 metros cuadrados y la distribucion del valor de los apartamentos se corrobora que que el promedio e de 208 millones y las casas con mayor valor se puede analisar valores por enciama del promedio de 600 millones y mas fuera de los rangos promedios.
ggpairs(vivienda4[,3:4], title="correlacion") 
cor.test(x = vivienda4$areaconst, y = vivienda4$preciom, method = "pearson", digits = 3)##
## Pearson's product-moment correlation
##
## data: vivienda4$areaconst and vivienda4$preciom
## t = 48.728, df = 1704, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.7424432 0.7821521
## sample estimates:
## cor
## 0.7630166La correlacion que se visuliza en el diagrama es de 0.763, siendo una correlacion positiva debil
3. Estime el modelo de regresión lineal simple entre precio=f(area)+ε. Interprete los coeficientes del modelo β0, β1 en caso de ser correcto.
modelocasa=lm(preciom ~ areaconst, data=Casa)
summary(modelocasa)##
## Call:
## lm(formula = preciom ~ areaconst, data = Casa)
##
## Residuals:
## Min 1Q Median 3Q Max
## -158.71 -52.01 -17.93 40.25 414.54
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 153.2252 15.9174 9.626 <2e-16 ***
## areaconst 1.2015 0.1123 10.701 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 83.87 on 341 degrees of freedom
## Multiple R-squared: 0.2514, Adjusted R-squared: 0.2492
## F-statistic: 114.5 on 1 and 341 DF, p-value: < 2.2e-16modeloapartamento=lm(preciom ~ areaconst, data=Apartamento)
summary(modeloapartamento)##
## Call:
## lm(formula = preciom ~ areaconst, data = Apartamento)
##
## 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 ***
## areaconst 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-16Se aplica el modelo de regresion lineal para las casas y apartmentos, observando que el intercepto en las casas es de 153 millones sin el precio bases para una area de 1 metros cuadrado, y para los apartamentos es de 39 millones de pesos, sin embargo por cada metros construido las casas aumentan en 1.201 millones por metro mientras que los aprtamentos aumentan en 2.16 millones por metro construido
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 hipótesis t.
confint(modelocasa, level = 0.95)## 2.5 % 97.5 %
## (Intercept) 121.916592 184.533755
## areaconst 0.980623 1.422293confint(modeloapartamento, level = 0.95)## 2.5 % 97.5 %
## (Intercept) 31.00423 47.089340
## areaconst 2.06264 2.266826El intervalo de conianza para el Intercept para el modelo de casas esta entre 1 y 1.4 millones por cada metro cuadrado construido, y el intervalo para el modelo de los apartamentos se encuentra entre 2.06 y 2.26 millones por metro cuadrado construido.
5. Calcule e interprete el indicador de bondad R2.
En el modelo de las casa el idicador de R2 es de 0.25 por lo cual se interpreta que el area construida influye en el precio en un 25%, y el indicador de R2 para los apartamentos es de 0.55, esto define que el 55% del valor depende del area construida, esto indica que en las casa el valor depende de otros factores en un 75% y en los apafrtmentos en un 45%, evidenciando que el valor de los apartementos depende mas del area que de otros factores.
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?.
area110 <- data.frame(areaconst = 110)
predict (modeloapartamento, area110, interval = "confidence", level = 0.95)## fit lwr upr
## 1 277.1674 272.9573 281.3775Segun la informacion analisada el precio propuesto es de 277 millones de pesos, con un limite inferior de 272 millones y un limite superior de 281 millones, con un indeice de confianza IC del 95%, por lo cual el cobro de 200 millones es una oferta muy atractiva sin embargo se encuentra por debajo del limite inferiro por mas de 70 millones, y seria necesario evaluar el estado que la promocion no genere perdidas.
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.
lillie.test(modeloapartamento$residuals)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: modeloapartamento$residuals
## D = 0.060955, p-value = 6.975e-13Para la prueba de normalidad se aplica el test de test Lilliefors, en la cual la hipotesis nula se rechaza y se acepta la hipotesis alterna.
res.stud <- studres(modeloapartamento)
mod.fit <- modeloapartamento$fitted.values
par(mfrow=c(1,2))
plot(mod.fit,res.stud,ylab="Residuos estudentizados",xlab="Valores ajustados")
abline(h=0,lty=2)
lines(lowess(res.stud~mod.fit), col = 2)
plot(mod.fit,abs(res.stud),ylab="|Residuos estudentizados|",xlab="Valores ajustados")
lines(lowess(abs(res.stud)~mod.fit), col = 2)
Podemos evidenciar los residuos del modelo para los apartamentos, los cuales se encuentran alrededor de la lÃnea de regresión los cuales están distribuidos de forma heterogénea, la dispersión es más significativa a medida que aumentan los valores ajustados, adicionalmente la variancia aumenta junto con los valores por lo cual podemos decir que no cumple con el supuesto de homocedasticidad
par(mfrow=c(1,2))
plot(modeloapartamento)
Para la prueba de linealidad de la distribucion de los residuos, se evidencia hay relacion lineal de la prueba Q-Q no se alinea de manera adecuada especialmente en los extremos.
8. De ser necesario realice una transformación apropiada para mejorar el ajuste y supuestos del modelo.
bc <- boxcox(lm(Apartamento$preciom ~ Apartamento$areaconst), lambda = -1:1)
lambda <- bc$x[which.max(bc$y)]
lambda## [1] -0.1313131nuevo_y = ((Apartamento$preciom ^ lambda) - 1 )/ lambda
modBox=lm(nuevo_y ~ areaconst, data=Apartamento)
summary(modBox)##
## Call:
## lm(formula = nuevo_y ~ areaconst, data = Apartamento)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.48740 -0.06527 -0.00535 0.07957 0.32740
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.4459708 0.0097192 354.55 <2e-16 ***
## areaconst 0.0046628 0.0001234 37.79 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1027 on 1361 degrees of freedom
## Multiple R-squared: 0.5121, Adjusted R-squared: 0.5117
## F-statistic: 1428 on 1 and 1361 DF, p-value: < 2.2e-16El valor de lambda es de -0.13, podemos evidenciar que el valor 0 se encuentra dentro del intervalo de confianza, por lo cual podemos definir que la mejor opcion es la trasformacion logaritmica, con este modelo no aumento el valor de R2, disminuye en un 51% lo cual significa estadisticamente significativo.
9. De ser necesario compare el ajuste y supuestos del modelo inicial y el transformado.
lillie.test(modBox$residuals)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: modBox$residuals
## D = 0.03952, p-value = 3.526e-05Segun los resultados obtenidos del Test de Lilliefors, se rechaza la hipotesis nula y aceptamos la alterna
bptest(modBox)##
## studentized Breusch-Pagan test
##
## data: modBox
## BP = 130.48, df = 1, p-value < 2.2e-16Con el fin de confirmar la hipotesis aplicamos el Test de Breush-Pagan, en el cual se evidencia que la homocedasticidad está presente, y rechazamos la hipotesis nula, aceptando la hipotesis alterna, puesto que el valor BP es menor (130.48), que el modelo de apartamentos, sin embargo en los dos casos el valor P es menor de 0.05, por lo cual se confirma la Heteroscedasticidad para los dos modelos
10. Estime varios modelos y compare los resultados obtenidos. En el mejor de los modelos, ¿se cumplen los supuestos sobre los errores?
Mod1 = lm (preciom ~ areaconst, data = Apartamento)
Mod2 = lm(preciom ~ log(areaconst), data = Apartamento)
Mod3 = lm(log(preciom) ~ areaconst, data = Apartamento)
Mod4 = lm(log(preciom) ~ log(areaconst), data = Apartamento)
stargazer(Mod1, Mod2, Mod3, Mod1, type="text", df=FALSE)##
## =======================================================================
## Dependent variable:
## ---------------------------------------------------
## preciom log(preciom) preciom
## (1) (2) (3) (4)
## -----------------------------------------------------------------------
## areaconst 2.165*** 0.009*** 2.165***
## (0.052) (0.0002) (0.052)
##
## log(areaconst) 195.419***
## (4.445)
##
## Constant 39.047*** -635.532*** 4.551*** 39.047***
## (4.100) (19.092) (0.019) (4.100)
##
## -----------------------------------------------------------------------
## Observations 1,363 1,363 1,363 1,363
## R2 0.560 0.587 0.520 0.560
## Adjusted R2 0.559 0.587 0.519 0.559
## Residual Std. Error 43.339 41.982 0.205 43.339
## F Statistic 1,730.157*** 1,933.199*** 1,473.424*** 1,730.157***
## =======================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01Cuando transformamos el área construida y la variable dependiente precio de los apartamentos, podemos definir que el cuarto modelo cuanta con el R2 que explica el área construida en un 58%.
Normalidad
lillie.test(Mod1$residuals)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: Mod1$residuals
## D = 0.060955, p-value = 6.975e-13lillie.test(Mod2$residuals)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: Mod2$residuals
## D = 0.060407, p-value = 1.22e-12lillie.test(Mod3$residuals)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: Mod3$residuals
## D = 0.033644, p-value = 0.001022lillie.test(Mod4$residuals)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: Mod4$residuals
## D = 0.035893, p-value = 0.0003038HOMOCEDASTICIDAD
res.stud<-studres(Mod1)
mod.fit<-Mod1$fitted.values
par(mfrow=c(1,2))
plot(mod.fit,res.stud,ylab="Residuos estudentizados Modelo 1",xlab="Valores ajustados Modelo 1")
abline(h=0,lty=2)
lines(lowess(res.stud~mod.fit), col = 2)
plot(mod.fit,abs(res.stud),ylab="|Residuos estudentizados Modelo 1|",xlab="Valores ajustados Modelo 1")
lines(lowess(abs(res.stud)~mod.fit), col = 2)
bptest(Mod1)##
## studentized Breusch-Pagan test
##
## data: Mod1
## BP = 292.99, df = 1, p-value < 2.2e-16res.stud<-studres(Mod2)
mod.fit<-Mod2$fitted.values
par(mfrow=c(1,2))
plot(mod.fit,res.stud,ylab="Residuos estudentizados Modelo 1",xlab="Valores ajustados Modelo 1")
abline(h=0,lty=2)
lines(lowess(res.stud~mod.fit), col = 2)
plot(mod.fit,abs(res.stud),ylab="|Residuos estudentizados Modelo 1|",xlab="Valores ajustados Modelo 1")
lines(lowess(abs(res.stud)~mod.fit), col = 2)
bptest(Mod2)##
## studentized Breusch-Pagan test
##
## data: Mod2
## BP = 214.66, df = 1, p-value < 2.2e-16res.stud<-studres(Mod3)
mod.fit<-Mod3$fitted.values
par(mfrow=c(1,2))
plot(mod.fit,res.stud,ylab="Residuos estudentizados Modelo 1",xlab="Valores ajustados Modelo 1")
abline(h=0,lty=2)
lines(lowess(res.stud~mod.fit), col = 2)
plot(mod.fit,abs(res.stud),ylab="|Residuos estudentizados Modelo 1|",xlab="Valores ajustados Modelo 1")
lines(lowess(abs(res.stud)~mod.fit), col = 2)
bptest(Mod3)##
## studentized Breusch-Pagan test
##
## data: Mod3
## BP = 150.38, df = 1, p-value < 2.2e-16res.stud<-studres(Mod4)
mod.fit<-Mod4$fitted.values
par(mfrow=c(1,2))
plot(mod.fit,res.stud,ylab="Residuos estudentizados Modelo 1",xlab="Valores ajustados Modelo 1")
abline(h=0,lty=2)
lines(lowess(res.stud~mod.fit), col = 2)
plot(mod.fit,abs(res.stud),ylab="|Residuos estudentizados Modelo 1|",xlab="Valores ajustados Modelo 1")
lines(lowess(abs(res.stud)~mod.fit), col = 2)
bptest(Mod4)##
## studentized Breusch-Pagan test
##
## data: Mod4
## BP = 92.877, df = 1, p-value < 2.2e-16LINEALIDAD
par(mfrow=c(1,2))
plot(Mod1)
par(mfrow=c(1,2))
plot(Mod2)
par(mfrow=c(1,2))
plot(Mod3)
par(mfrow=c(1,2))
plot(Mod4)
NO AUTOCORRELACION DE ERRORES
lmtest::dwtest(Mod1)##
## Durbin-Watson test
##
## data: Mod1
## DW = 1.443, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0lmtest::dwtest(Mod2)##
## Durbin-Watson test
##
## data: Mod2
## DW = 1.4775, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0lmtest::dwtest(Mod3)##
## Durbin-Watson test
##
## data: Mod3
## DW = 1.3187, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0lmtest::dwtest(Mod4)##
## Durbin-Watson test
##
## data: Mod4
## DW = 1.3214, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0OUTLIERS
ggplot(data = Apartamento, aes(x = seq_along(Mod1$residuals), y = Mod1$residuals)) +
geom_point(aes(color = Mod1$residuals)) + geom_line(linewidth = 0.3) + labs(title = "Distribución de residuos Modelo 1", x = "observaciones aptos", y = "residuo") + geom_hline(yintercept = 0) + theme(plot.title = element_text(hjust = 0.5), legend.position = "none")
ggplot(data = Apartamento, aes(x = seq_along(Mod2$residuals), y = Mod2$residuals)) +
geom_point(aes(color = Mod2$residuals)) + geom_line(linewidth = 0.3) + labs(title = "Distribución de residuos Modelo 2", x = "observaciones aptos", y = "residuo") + geom_hline(yintercept = 0) + theme(plot.title = element_text(hjust = 0.5), legend.position = "none")
ggplot(data = Apartamento, aes(x = seq_along(Mod3$residuals), y = Mod3$residuals)) +
geom_point(aes(color = Mod3$residuals)) + geom_line(linewidth = 0.3) + labs(title = "Distribución de residuos Modelo 3", x = "observaciones aptos", y = "residuo") + geom_hline(yintercept = 0) + theme(plot.title = element_text(hjust = 0.5), legend.position = "none")
ggplot(data = Apartamento, aes(x = seq_along(Mod4$residuals), y = Mod4$residuals)) +
geom_point(aes(color = Mod4$residuals)) + geom_line(linewidth = 0.3) + labs(title = "Distribución de residuos Modelo 4", x = "observaciones aptos", y = "residuo") + geom_hline(yintercept = 0) + theme(plot.title = element_text(hjust = 0.5), legend.position = "none")