En esta ocasión nos piden filtrar las variables precio en millones y área construida de la base de datos “vivienda”, recordemos que cuenta con 8319 registros y 12 variables de diferentes viviendas en Cali. Adicionalmente, como nuestro barrio de interés es el ingenio lo filtramos de la siguiente manera:
Datos<-Vivienda[Vivienda$Barrio=="el ingenio",c(4,5,9)]
summary(Datos)
## precio_millon Area_contruida Tipo
## Min. : 175.0 Min. : 56.00 Length:205
## 1st Qu.: 290.0 1st Qu.: 96.25 Class :character
## Median : 352.5 Median : 128.50 Mode :character
## Mean : 430.6 Mean : 193.88
## 3rd Qu.: 488.8 3rd Qu.: 203.75
## Max. :1450.0 Max. :1200.00
## NA's :3 NA's :3
Las viviendas del barrio ingenio son 205, de estas viviendas el 63% son apartamentos y el 36% son casas, el precio de las viviendas ubicadas en este barrio está entre 175 y 1450 millones, por otro lado el área está entre 56 y 1200 metros cuadrados.
Lo primero que realizamos es un modelo de regresión simple para predecir el precio mediante el área construida.
mod<-lm(precio_millon~Area_contruida,data=Datos)
summary(mod)
##
## Call:
## lm(formula = precio_millon ~ Area_contruida, data = Datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -282.09 -61.15 -8.32 45.89 477.33
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 195.71496 11.63854 16.82 <2e-16 ***
## Area_contruida 1.21159 0.04632 26.16 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 105.2 on 200 degrees of freedom
## (3 observations deleted due to missingness)
## Multiple R-squared: 0.7738, Adjusted R-squared: 0.7727
## F-statistic: 684.2 on 1 and 200 DF, p-value: < 2.2e-16
Al ajustar un modelo de regresión lineal simple obtenemos un \(R^2\) de 0.77 , esto quiere decir que el \(77\)% de su variabilidad del precio está siendo explicada por el modelo ajustado en función del área construida, en casas del barrio ingenio.También es necesario resaltar que el primer coeficiente resultó ser significativo y tiene un valor de \(B_{1}=1.21159\), lo que quiere decir que por cada metro cuadrado en el área construida el precio de la vivienda vivienda cambia en 1.2 millones.
A continuación realizamos el análisis de los supuestos de normalidad y homocedasticidad.
En el siguiente gráfico se observa que los cuantiles muestrales no se ajustan a los cuantiles teóricos que esperarían bajo normalidad. Además, al realizar la prueba de Shapiro-Wilks se obtiene un p-valor muy pequeño que nos lleva a rechazar la hipótesis de normalidad.
car::qqPlot(mod)
## [1] 75 125
shapiro.test(mod$residuals)
##
## Shapiro-Wilk normality test
##
## data: mod$residuals
## W = 0.94495, p-value = 5.66e-07
En la gráfica se nota que la varianza de los residuos tiende a crecer a medida que lo hacen los valores ajustados. Esto es una señal de heterocedasticidad. Tras realizar la prueba de Breusch-Pagan obtenemos un p-valor pequeño que permite concluir que no se cumple el supuesto de homocedasticidad.
res.stud<-studres(mod)
mod.fit<-mod$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)
bptest(mod)
##
## studentized Breusch-Pagan test
##
## data: mod
## BP = 46.496, df = 1, p-value = 9.183e-12
Para corregir los supuestos realizamos la transformación por BoxCox. Obtenemos un valor de 0.02 que maximiza la logverosimilitud. Este es el valor al cual elevaremos la variable respuesta.
boxcox = MASS::boxcox(mod,lambda=seq(-3,3,length.out = 1000),
ylab='log-verosimilitud')
boxcox$x[boxcox$y ==max(boxcox$y)]
## [1] 0.2012012
modBC<-lm(precio_millon^0.02~Area_contruida,data=Datos)
summary(modBC)
##
## Call:
## lm(formula = precio_millon^0.02 ~ Area_contruida, data = Datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0202422 -0.0035888 -0.0001349 0.0033511 0.0138245
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.117e+00 5.961e-04 1873.99 <2e-16 ***
## Area_contruida 4.990e-05 2.372e-06 21.03 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.005389 on 200 degrees of freedom
## (3 observations deleted due to missingness)
## Multiple R-squared: 0.6887, Adjusted R-squared: 0.6871
## F-statistic: 442.4 on 1 and 200 DF, p-value: < 2.2e-16
Al realizar la transformación por boxcox y ajustar el modelo vemos que la estimación del error se reduce considerablemente ya que antes tenía un valor de 105.2 y ahora 0.005389, el \(R^{2}=0.68\) y el área sigue siendo significativa.
Para ajustar el modelo de regresión categórica de las viviendas del barrio ingenio necesitamos hacer una transformación en la variable Tipo y agregarle un uno si la vivienda es una casa y cero en el caso contrario.Eso lo hacemos con la siguiente función.
Newdata<-dummy_cols(
Vivienda,
select_columns = "Tipo"
)
Newdata<-Newdata[Newdata$Barrio=="el ingenio",c(4,5,9,13,14)]
View(Newdata)
Ajustamos un modelo de la variable precio en millones en función del área construida y el tipo de vivienda casa.
modC= lm(precio_millon~Area_contruida+Tipo_Casa, data=Newdata)
summary(modC)
##
## Call:
## lm(formula = precio_millon ~ Area_contruida + Tipo_Casa, data = Newdata)
##
## Residuals:
## Min 1Q Median 3Q Max
## -277.45 -48.20 -9.33 53.92 491.45
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 196.11964 10.91567 17.967 < 2e-16 ***
## Area_contruida 1.02197 0.05616 18.196 < 2e-16 ***
## Tipo_Casa 99.24986 18.63125 5.327 2.69e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 98.68 on 199 degrees of freedom
## (3 observations deleted due to missingness)
## Multiple R-squared: 0.802, Adjusted R-squared: 0.8
## F-statistic: 403.1 on 2 and 199 DF, p-value: < 2.2e-16
El primer coeficiente \(B_{1}\) nos dice que si el área aumenta en un metro cuadrado el precio aumenta en 1.02197 millones y el segundo coeficiente nos dice que el precio de las casas en promedio es 99.24986 mayor que el de los apartamentos,ambas variables son significativas sin embargo el error estimado es alto. Por otro lado tenemos el modelo ajustado al tipo de vivienda apartamento.
modA=lm(precio_millon~Area_contruida+Tipo_Apartamento, data=Newdata)
summary(modA)
##
## Call:
## lm(formula = precio_millon ~ Area_contruida + Tipo_Apartamento,
## data = Newdata)
##
## Residuals:
## Min 1Q Median 3Q Max
## -277.45 -48.20 -9.33 53.92 491.45
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 295.36951 21.65885 13.637 < 2e-16 ***
## Area_contruida 1.02197 0.05616 18.196 < 2e-16 ***
## Tipo_Apartamento -99.24986 18.63125 -5.327 2.69e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 98.68 on 199 degrees of freedom
## (3 observations deleted due to missingness)
## Multiple R-squared: 0.802, Adjusted R-squared: 0.8
## F-statistic: 403.1 on 2 and 199 DF, p-value: < 2.2e-16
El primer coeficiente \(B_{1}\) nos dice que si el área aumenta en un metro cuadrado el precio aumenta en 1.02197 millones y el segundo coeficiente nos dice que el precio de los apartamentos en promedio es -99.24986 menor que el de las casas, ambas variables son significativas sin embargo el error estimado es alto.
De lo anterior podemos deducir modelo ajustado en función de la variables área construida, y tipo de vivienda arroja un \(R^{2}=0.80\) cuándo tipo de vivienda es casa o apartamento, que en comparación con el modelo de regresión lineal simple es mayor y también que el precio de la vivienda aumenta cuando el tipo de vivienda es una casa.
A continuación realizamos el análisis de los supuestos de normalidad y homocedasticidad.
En el siguiente gráfico se observa que los cuantiles muestrales no se ajustan a los cuantiles teóricos que esperarían bajo normalidad. Además, al realizar la prueba de Shapiro-Wilks se obtiene un p-valor muy pequeño que nos lleva a rechazar la hipótesis de normalidad.
car::qqPlot(modC)
## [1] 75 125
shapiro.test(modC$residuals)
##
## Shapiro-Wilk normality test
##
## data: modC$residuals
## W = 0.93208, p-value = 4.378e-08
En la gráfica se nota que los la varianza de los residuos tiende a crecer a medida que lo hacen los valores ajustados. Esto es una señal de heterocedasticidad. Tras realizar la prueba de Breusch-Pagan obtenemos un p-valor pequeño que permite concluir que no se cumple el supuesto de homocedasticidad.
res.studC<-studres(modC)
mod.fitC<-modC$fitted.values
par(mfrow=c(1,2))
plot(mod.fitC,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.studC),ylab="|Residuos estudentizados|",xlab="Valores ajustados")
lines(lowess(abs(res.studC)~mod.fitC), col = 2)
bptest(modC)
##
## studentized Breusch-Pagan test
##
## data: modC
## BP = 38.517, df = 2, p-value = 4.326e-09
Para corregir los supuestos realizamos la transformación por BoxCox. Obtenemos un valor de 0.05 que maximiza la logverosimilitud. Este es el valor al cual elevaremos la variable respuesta.
boxcoxC = MASS::boxcox(modC,lambda=seq(-3,3,length.out = 1000),
ylab='log-verosimilitud')
boxcoxC$x[boxcoxC$y ==max(boxcoxC$y)]
## [1] 0.05705706
modBCa<-lm(precio_millon^0.05~Area_contruida+Tipo_Casa,data=Newdata)
summary(modBCa)
##
## Call:
## lm(formula = precio_millon^0.05 ~ Area_contruida + Tipo_Casa,
## data = Newdata)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.034866 -0.008077 -0.001100 0.010970 0.041591
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.319e+00 1.540e-03 856.516 < 2e-16 ***
## Area_contruida 1.089e-04 7.922e-06 13.747 < 2e-16 ***
## Tipo_Casa 2.171e-02 2.628e-03 8.263 1.98e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01392 on 199 degrees of freedom
## (3 observations deleted due to missingness)
## Multiple R-squared: 0.7712, Adjusted R-squared: 0.7689
## F-statistic: 335.3 on 2 and 199 DF, p-value: < 2.2e-16
El modelo de regresión categórica con la transformación disminuyo significativamente la estimación del error sin implicar mucha pérdida de la variabilidad explicada del precio en función del área construida y el tipo de vivienda, el \(R^{2}= 0.77\) y mucho mayor que el modelo transformado de regresión simple el cual era \(R^{2}= 0.66\)
Lo anterior nos da un indicio de que incluir el tipo de vivienda como variable explicativa sería lo adecuado, sin embargo lo probaremos con un anova.
anova(mod,modC)
## Analysis of Variance Table
##
## Model 1: precio_millon ~ Area_contruida
## Model 2: precio_millon ~ Area_contruida + Tipo_Casa
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 200 2214303
## 2 199 1937949 1 276354 28.378 2.687e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
De acuerdo al anterior anova, podríamos decir que el mejor modelo es el de regresión categórica, ya que el valor p es pequeño y se rechaza la hipótesis nula de que el coeficiente asociado a la variable tipo es cero.
A continuación realizamos un modelo polinomial para los datos de masa y edad de un grupo de personas.La base de datos cuenta con 60 registros de las dos variables mencionadas anteriormente.La msa toma valores entre 32 y 119, y la edad entre 41 y 78 años.
Para empezar, se realiza la centralización de los datos de la
variable predictora (Edad) de la siguiente manera:
\(X_i-\overline{X}\)
Tras realizar la transformación, se ajusta un modelo cuadrático.
EdadX<-scale(datosME$Edad,center = T,scale = F)
modSq2<-lm(datosME$Masa~EdadX+I(EdadX^2))
Esta es la curva obtenida por el modelo cuadrático y sus intervalos de confianza.
summary(modSq2)
##
## Call:
## lm(formula = datosME$Masa ~ EdadX + I(EdadX^2))
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.086 -6.154 -1.088 6.220 20.578
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 82.935749 1.543146 53.745 <2e-16 ***
## EdadX -1.183958 0.088633 -13.358 <2e-16 ***
## I(EdadX^2) 0.014840 0.008357 1.776 0.0811 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.026 on 57 degrees of freedom
## Multiple R-squared: 0.7632, Adjusted R-squared: 0.7549
## F-statistic: 91.84 on 2 and 57 DF, p-value: < 2.2e-16
Encontramos que el modelo cuadrático logra explicar el 76.32% de la variable respuesta. Por otro lado, teniendo una significancia de 0.05, el coeficiente \(B_1\) es significativo, sin embargo el coeficiente \(B_2\) no logra un aporte significativo al modelo.
Para la interpretación de los \(Betas\), dado que \(B_2\) toma un valor muy pequeño y \(B_1\) toma un valor negativo y teniendo en cuenta que un valor negativo en la covariable \(EdadX\) son aquellos que están por debajo de la media, podemos decir que la relación entre la edad y la masa es negativa.
Anteriormente se obtuvo que el coeficiente \(B_2\) no tiene un aporte significativo, para corroborar esto realizaremos un ANOVA comparando con un modelo sin el valor cuadrático. Se espera que si el valor cuadrático tiene un aporte significativo, el modelo cuadrático reduzca significativamente la suma de cuadrados de los residuos. A continuación se presentan los resultados:
modSimple<-lm(datosME$Masa~EdadX)
anova(modSimple,modSq2)
## Analysis of Variance Table
##
## Model 1: datosME$Masa ~ EdadX
## Model 2: datosME$Masa ~ EdadX + I(EdadX^2)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 58 3874.4
## 2 57 3671.3 1 203.13 3.1538 0.08109 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Tras obtener un p-valor superior a 0.05 concluimos que podemos eliminar el término cuadrático, ya que no tiene un aporte significativo.
Verificaremos si realizar la transformación \(X_i-\overline{X}\) logra o no solucionar los problemas de multicolinealidad. Para ello creamos un modelo cuadrático sin realizar la transformación y compararemos los resutados.
modSq<-lm(Masa~Edad+I(Edad^2),data=datosME)
vif(modSq); vif(modSq2)
## Edad I(Edad^2)
## 128.256 128.256
## EdadX I(EdadX^2)
## 1.001473 1.001473
dat<-as.data.frame(cbind(datosME$Masa,EdadX,EdadX^2))
dat2<-as.data.frame(cbind(datosME,datosME$Edad^2))
cor(dat)
## V1 V2 V3
## V1 1.0000000 -0.86606398 0.14760734
## V2 -0.8660640 1.00000000 -0.03835694
## V3 0.1476073 -0.03835694 1.00000000
cor(dat2)
## Masa Edad datosME$Edad^2
## Masa 1.0000000 -0.8660640 -0.8525732
## Edad -0.8660640 1.0000000 0.9960939
## datosME$Edad^2 -0.8525732 0.9960939 1.0000000
Los resultados muestran que la transformación sí es efectiva para solucionar la multicolinealidad. La correlación entre la \(Edad\) y \(Edad^2\) es muy baja, produciendo valores del VIF cercanos a 1.