library(readxl)
datos <- read_excel("datos_vivienda.xlsx")
attach(datos)
summary(datos$precio_millon)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 240.0 251.2 305.0 332.1 395.0 480.0
Desviacion_Estandar_Precio <-sd(datos$precio_millon)
Desviacion_Estandar_Precio
## [1] 82.14423
Según estos datos podemos afirmar que la variable Precio de Vivienda (millones de pesos), cuenta con valores entre 240.0 siendo su mínimo y 480.0 siendo su máximo, además cuenta con una media de 332.1.0 y una desviación estándar de 82.14423 millones de pesos.
plot(datos$precio_millon,
main = "Distribución de Precio de Vivienda",
xlab = "Índice del dato",
ylab = "Millones de pesos"
)
hist(datos$precio_millon,
main = "Frecuencia de Precio de Vivienda",
xlab = "Millones de pesos",
ylab = "Frecuencia"
)
summary(datos$Area_contruida)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 80.0 86.0 97.0 115.7 130.0 195.0
Desviacion_Estandar_Area <- sd(datos$Area_contruida)
Desviacion_Estandar_Area
## [1] 35.54332
Según estos datos podemos afirmar que la variable Área Construida (metros cuadrados), cuenta con valores entre 80.0 siendo su mínimo y 195.0 siendo su máximo, además cuenta con una media de 115.7 y una desviación estándar de 35.54332 metros cuadrados.
plot(datos$Area_contruida,
main = "Distribución de Área Construida",
xlab = "Índice del dato",
ylab = "Metros cuadrados"
)
hist(datos$Area_contruida,
main = "Frecuencia de Área Construida",
xlab = "Metros cuadrados",
ylab = "Frecuencia"
)
summary(datos)
## 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
cor(datos$Area_contruida,datos$precio_millon)
## [1] 0.9190295
Según estos datos, podemos afirmar que nuestro modelo puede usar los datos en (80.0, 195.0), y obtener una respuesta en (240.0, 480.0), intentar observar o predecir valores fuera de estos rangos sería una extrapolación para la que el modelo no aporta datos y por ende no es correcto. El Coeficiente de Correlación de Pearson nos arroja un valor de 0.9190295, este coeficiente cuenta con valores en el rango (-1, 1), mientras más positivo, las variables están más relacionadas de manera directamente proporcional, por ende, podemos concluir que al tener un Coeficiente de Correlación de 0.9190295, estas dos variables tienen una relación directamente proporcional, a más aumente una, más aumenta la otra.
plot(datos$Area_contruida,datos$precio_millon
, main = "Precio de Vivienda vs. Área Construida",
xlab = "Metros cuadrados",
ylab = "Millones de pesos"
)
En el diagrama se denota una relación creciente, pero con una aparente curvatura leve que podría indicar que el modelo no estaría correctamente adaptado para hacer uso de la Regresión Lineal Simple y podría necesitar una transformación, esto se corroborará con la validación de supuestos más adelante.
#x y
x <- Area_contruida
y <- precio_millon
Nuestro modelo de regresión lineal simple del precio en función de una función de área y la suma del error es \(precio= 86.234 + 2.124\cdot area +e\).
v <- summary(datos)
mod=lm(y~x)
summary(mod)
##
## Call:
## lm(formula = y ~ x)
##
## 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 ***
## x 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
Si bien \(\small \beta _{0}\) es significativa no tiene sentido interpretar un apartamento de 0 metros cuadrados, \(\small \beta _{1}\) por su parte también es significativa y podemos interpretarle como que por cada aumento unitario en el \(area\) el \(precio\) aumenta en promedio 2.124 millones, lo que quiere decir que su relación es positiva.
b <- summary(mod)
n <- b$df[2]+2
bs <- data.frame(b$coefficients)
b0 <- data.frame(b$coefficients)$Estimate[1]
b1 <- data.frame(b$coefficients)$Estimate[2]
b0d <- c(c(bs$Std..Error[1]),c(bs$t.value[1]),c(bs$Pr...t..[1]))
b1d <- c(c(bs$Std..Error[2]),c(bs$t.value[2]),c(bs$Pr...t..[2]))
sigb0 <- function(a){
ta <- qt(a/2,n-2,lower.tail = F)
if(abs(b0d[2])>ta)return('significativa')
if(abs(b0d[2])<=ta)return('no significativa')
}
sigb1 <- function(a){
ta <- qt(a/2,n-2,lower.tail = F)
if(abs(b1d[2])>ta)return('significativa')
if(abs(b1d[2])<=ta)return('no significativa')
}
#nos regresa si los valores $\beta _{0}$ y $\beta _{1}$ son significativos
sigb0(0.05)
## [1] "significativa"
sigb1(0.05)
## [1] "significativa"
un intervalo de confianza al 95% para el coeficiente \(\small\beta _{1}\) es \(\small(1.740, 2.508)\) esto quiere decir que el verdadero \(\small\beta _{1}\) esta entre esos valores con un 95% de confianza, 0 no está dentro de este intervalo por eso podemos decir con confianza que el valor es diferente de 0 y por lo tanto significativa.
icb1 <- function(a){
ta <- qt(a/2,n-2,lower.tail = F)
ib1 <- c(b1-(ta*b1d[1]),b1+(ta*b1d[1]))
return(ib1)
}
icb1(0.05)# devuelve un intervalo con confianza a para b1
## [1] 1.740170 2.507771
Haciendo una prueba de significancia con \(\small t\) para nuestro \(\small\beta _{1}\) nuestra hipótesis inicial es \(\small\beta _{1}\) es igual a 0 y la alternativa es que \(\small\beta _{1}\) es diferente de 0 y tenemos que nuestro \(\small t _{0}\) es 11.42 y nuestro t es 2.06, ahora bien nuestra hipótesis inicial es rechazada, así que nuestra \(\small\beta _{1}\) es diferente de 0 y por lo tanto significativa.
b1d[2]
## [1] 11.4217
qt(0.05/2,n-2,lower.tail = F)
## [1] 2.063899
sigb1 <- function(a){
ta <- qt(a/2,n-2,lower.tail = F)
if(abs(b1d[2])>ta)return('significativa')
if(abs(b1d[2])<=ta)return('no significativa')
}
sigb1(0.05)
## [1] "significativa"
el \(\small R^{2}\) es igual a 0.845 lo que significa que el modelo explica el 84.5% de la variabilidad del precio de las viviendas , y que modelo tiene un ajuste relativamente bueno a la recta de regresión, sin embargo no significa que cumpla con los supuestos de linealidad, aunque podría mejorar.
b$r.squared
## [1] 0.8446152
Para poder hacer el calculo del precio promedio estimado para este apartamento primero tenemos que saber sí la \(x\) es interpolada, y ya que esta dentro del rango de las \(x\), si podemos proceder a hacer cálculos.
xo <- 110
c(xo>min(x),xo<max(x))
## [1] TRUE TRUE
El precio promedio estimado para un apartamento de 110 metros cuadrados esta en un intervalo de 306.3133 millones a 333.4279 millones, por lo tanto, un apartamento con un precio de 200 millones es una gran oferta en esta zona, una consideración adicional a tener en en cuenta es que la razón por la que tenga tan buen precio puede estar relacionada con el estado del apartamento, su poca deseabilidad en el mercado, o que hayan sucedido cosas terribles en el lugar.
ta <- qt(0.05/2,n-2,lower.tail = F)
yi <- predict(mod, newdata = data.frame(x=c(xo)), se.fit = T)
yo <- yi$fit
ys <- yi$se.fit
c(yo-(ta*ys),yo+(ta*ys))
## 1 1
## 306.3133 333.4279
mod = lm(datos$precio_millon~datos$Area_contruida)
summary(mod)
##
## Call:
## lm(formula = datos$precio_millon ~ datos$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 ***
## datos$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
ei = mod$residuals
yi_mod = mod$fitted.values
round(mean(ei),3)
## [1] 0
Gracias a esto validamos que la media es efectivamente 0, se valida el supuesto de Media Cero.
plot(yi_mod,ei,main="Residuales vs. Valores Ajustados", xlab = "Valores Ajustados", ylab = "Residuales")
abline(h=0,lty=2,col=2)
require(ggplot2)
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 4.0.4
data_val=data.frame(yi_mod,ei)
ggplot(data_val,aes(x=yi_mod,y=ei))+geom_point()+theme_bw()+geom_smooth()
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
Aquí podemos observar que no se estaría cumpliendo gráficamente los requerimientos para validar una varianza constante, el gráfico sólo nos dice que sospechemos la no linealidad del modelo, no se valida el supuesto de Varianza constante.
ei
## 1 2 3 4 5 6 7
## -18.895336 48.137608 32.649962 -51.672533 -28.895336 25.617018 32.691142
## 8 9 10 11 12 13 14
## -22.135041 -26.771366 2.691142 -31.019307 47.245540 21.104664 -26.771366
## 15 16 17 18 19 20 21
## -18.135041 -18.895336 -18.895336 32.649962 14.154080 -1.151513 67.649962
## 22 23 24 25 26
## -11.019307 22.649962 14.732752 -50.408120 -37.308858
hist(ei,col="gray", main = "Histograma de Residuales", xlab ="Residuales", ylab = "Frecuencia")
qqnorm(ei)
qqline(ei,col="red")
shapiro.test(ei)
##
## Shapiro-Wilk normality test
##
## data: ei
## W = 0.95489, p-value = 0.3009
alpha = 0.05
Se observa que el gráfico de normalidad presenta curvaturas y una tendencia a alejarse de la recta teórica, la prueba de Shapiro-Wilk nos arroja un valor P de 0.3009, mayor que el valor alfa de 0.05, entre la prueba gráfica y el test de Shapiro-Wilk prima la validación grafica, que en este caso no se puede otorgar por inconsistencia con los requerimientos de validación, no se valida el supuesto de Normalidad.
Dado que todos estos registros no corresponden a datos en el tiempo no se tiene un orden para realizar la validación del supuesto de Independencia de variables. Se valida por definición del tipo de datos de corte transversal.
mod.lof<-lm(datos$precio_millon~as.factor(datos$Area_contruida))
anova(mod,mod.lof)
## Analysis of Variance Table
##
## Model 1: datos$precio_millon ~ datos$Area_contruida
## Model 2: datos$precio_millon ~ as.factor(datos$Area_contruida)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 24 26212.2
## 2 12 5363.4 12 20849 3.8872 0.01305 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Para un nivel de confianza del 95%, se rechaza la Hipótesis de que el modelo sea lineal, por ende, se le debe hacer alguna transformación para que valide supuestos y linealidad.
Vemos que nuestro modelo actual incumple varios supuestos asi que hemos decidido cambiar de modelo.
area_transformada = 1 / datos$Area_contruida
modelo_transformado <- lm(datos$precio_millon~area_transformada)
summary(modelo_transformado)
##
## Call:
## lm(formula = datos$precio_millon ~ area_transformada)
##
## Residuals:
## Min 1Q Median 3Q Max
## -36.987 -16.743 -5.023 18.547 44.379
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 635.35 18.27 34.77 < 2e-16 ***
## area_transformada -32464.72 1895.32 -17.13 5.84e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 23.05 on 24 degrees of freedom
## Multiple R-squared: 0.9244, Adjusted R-squared: 0.9212
## F-statistic: 293.4 on 1 and 24 DF, p-value: 5.839e-15
El supuesto de la media cero se cumple por defecto.
En el gráfico de residuales vs predichos podemos observar que se cumple el supuesto de varianza constante y que mejora con respecto al modelo orginal.
plot(fitted(modelo_transformado), residuals(modelo_transformado), main = "Residuales vs valores ajustados")
abline(h = 0, lty = 2,col= 2)
No hay evidencia fuerte en contra del supuesto de que los residuales distribuyen normal, ademas se observa que el gráfico de probabilidad normal mejora con respecto al módelo original.
plot(modelo_transformado, 2)
Dado que el valor p = 0.5871 del test shapiro es suficientemente grande confirmamos que se cumple la hipotesis de normalidad.
test_shapiro = shapiro.test(residuals(modelo_transformado))
test_shapiro
##
## Shapiro-Wilk normality test
##
## data: residuals(modelo_transformado)
## W = 0.96859, p-value = 0.5871
Dado que estos registros no corresponden a datos en el tiempo no se tiene un orden temporal para realizar la validación de este supuesto. Se valida por definición del tipo de datos de corte
De la prueba de falta de ajuste podemos observar que el valor F ajuste(1.38) es inferior a f k.2,n-k(0.37) al 95% por tanto no hay evidencia que rechace la linealidad del modelo.
mod.aux <- lm(datos$precio_millon~as.factor(area_transformada))
anova(modelo_transformado,mod.aux)
## Analysis of Variance Table
##
## Model 1: datos$precio_millon ~ area_transformada
## Model 2: datos$precio_millon ~ as.factor(area_transformada)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 24 12755.6
## 2 12 5363.4 12 7392.2 1.3783 0.2935
e <- anova(modelo_transformado,mod.aux)
k<- e$Df[2]+2
fb <- e$F[2]
rlt <- function(a){
fa <- qf(a,k-2,n-k,lower.tail = F)
if(fb>fa) return('no lineal')
if(fb<fa) return('lineal')
}
rlt(0.05)
## [1] "lineal"
Decidimos entonces usar el modelo \(\small y= \beta _{0} + \beta _{1}\cdot\frac{1}{x} +e\) porqué dado que tiene un \(\small R^{2}\) más alto tiene un mejor ajuste en su recta de regresión y alto puede hacer predicciones más útiles, además de que cumple todos los supuestos de una regresión lineal.
#x y
x <- 1/Area_contruida
y <- precio_millon
Veamos la correlación entre el precio en millones de pesos y el reciproco del area en metros cuadrados. Se observa que a mayor cantidad en el area transformada el precio de las viviendas disminuye y su relación es fuerte de acuerdo con el coeficiente de correlación de Pearson (-0,96) mejorando la relación con respecto modelo original con coeficiente de correlación de Pearson(0.92).
area_transformada = 1 / datos$Area_contruida
plot(datos$precio_millon, area_transformada, ylab = "Precio en Millones",
xlab = "1/Area" , main = "Gráfico de dispersión Y = β0 + β1(-1/X) + ε")
cor(datos$precio_millon, area_transformada)
## [1] -0.9614495
Nuestro modelo de regresión lineal simple del precio en función de una función de área y la suma del error es \(precio= 635.3498 - 32464.72\cdot \frac{1}{area} +e\).
x <- 1/Area_contruida
y <- precio_millon
v <- summary(datos)
mod=lm(y~x)
summary(mod)
##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -36.987 -16.743 -5.023 18.547 44.379
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 635.35 18.27 34.77 < 2e-16 ***
## x -32464.72 1895.32 -17.13 5.84e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 23.05 on 24 degrees of freedom
## Multiple R-squared: 0.9244, Adjusted R-squared: 0.9212
## F-statistic: 293.4 on 1 and 24 DF, p-value: 5.839e-15
Si bien \(\small\beta _{0}\) es significativa no tiene sentido interpretar un apartamento de infinitos metros cuadrados negativos, \(\small \beta _{1}\) por su parte también es significativa y podemos interpretarlo como que por cada aumento unitario en \(\small \frac{1}{área}\) el precio disminuye en promedio 32464.72 millones, lo que quiere decir que si el área aumenta el precio también lo hace de un modo reciproco.
b <- summary(mod)
n <- b$df[2]+2
bs <- data.frame(b$coefficients)
b0 <- data.frame(b$coefficients)$Estimate[1]
b1 <- data.frame(b$coefficients)$Estimate[2]
b0d <- c(c(bs$Std..Error[1]),c(bs$t.value[1]),c(bs$Pr...t..[1]))
b1d <- c(c(bs$Std..Error[2]),c(bs$t.value[2]),c(bs$Pr...t..[2]))
sigb0 <- function(a){
ta <- qt(a/2,n-2,lower.tail = F)
if(abs(b0d[2])>ta)return('significativa')
if(abs(b0d[2])<=ta)return('no significativa')
}
sigb1 <- function(a){
ta <- qt(a/2,n-2,lower.tail = F)
if(abs(b1d[2])>ta)return('significativa')
if(abs(b1d[2])<=ta)return('no significativa')
}
#nos regresa si los valores $\beta _{0}$ y $\beta _{1}$ son significativos
sigb0(0.05)
## [1] "significativa"
sigb1(0.05)
## [1] "significativa"
un intervalo de confianza al 95% para el coeficiente \(\small\beta _{1}\) es \(\small (-36376.47, -28552.97)\) esto quiere decir que el verdadero \(\small\beta _{1}\) esta entre esos valores con un 95% de confianza, 0 no esta dentro de este intervalo por eso podemos decir con confianza que el valor es diferente de 0 y por lo tanto significativa.
icb1 <- function(a){
ta <- qt(a/2,n-2,lower.tail = F)
ib1 <- c(b1-(ta*b1d[1]),b1+(ta*b1d[1]))
return(ib1)
}
icb1(0.05)# devuelve un intervalo con confianza a para b1
## [1] -36376.47 -28552.97
Haciendo una prueba de significancia con t para nuestro \(\small\beta _{1}\) nuestra hipótesis inicial es \(\small \beta _{1}\) es igual a 0 y la alternativa es que \(\small \beta _{1}\) es diferente de 0 y tenemos que nuestro \(\small t _{0}\) es -17.13 y nuestro t es 2.06, ahora bien nuestra hipótesis inicial es rechazada, así que nuestra \(\small \beta _{1}\) es diferente de 0 y por lo tanto significativa.
b1d[2]
## [1] -17.12887
qt(0.05/2,n-2,lower.tail = F)
## [1] 2.063899
sigb1 <- function(a){
ta <- qt(a/2,n-2,lower.tail = F)
if(abs(b1d[2])>ta)return('significativa')
if(abs(b1d[2])<=ta)return('no significativa')
}
sigb1(0.05)
## [1] "significativa"
el \(\small R^{2}\) es igual a 0.924 lo que significa que el modelo explica el 92.44% de la variabilidad del precio de las viviendas, además de que el modelo tiene un muy buen ajuste a su recta de regreción y alto puede hacer predicciones útiles.
b$r.squared
## [1] 0.9243852
Para poder hacer el cálculo del precio promedio estimado para este apartamento primero tenemos que saber sí la \(x\) es interpolada, y ya que está dentro del rango de las \(x\), si podemos proceder a hacer cálculos, no hay que olvidar que nuestro modelo cambio y por lo tanto la \(x\) también.
xo <- 1/110
c(xo>min(x),xo<max(x))
## [1] TRUE TRUE
El precio promedio estimado para un apartamento de 110 metros cuadrados esta en un intervalo de 330.8332 millones a 349.5988 millones, por lo tanto un apartamento con un precio de 200 millones es una gran oferta en esta zona, una consideración adicional a tener en en cuenta es que la razón por la que tenga tan buen precio puede estar relacionada con el estado del apartamento, su poca deseabilidad en el mercado, o que hayan sucedido cosas terribles en el lugar.
ta <- qt(0.05/2,n-2,lower.tail = F)
yi <- predict(mod, newdata = data.frame(x=c(xo)), se.fit = T)
yo <- yi$fit
ys <- yi$se.fit
c(yo-(ta*ys),yo+(ta*ys))
## 1 1
## 330.8332 349.5988
Una función para calcular el intervalo de confianza de un \(\small \beta _{1}\) en un \(\small 100(1-a)%\) dados dos vectores X y Y es:
icb1 <- function(x,y,a){
mod=lm(y~x)
b <- summary(mod)
n <- b$df[2]+2
b1 <- data.frame(b$coefficients)$Estimate[2]
b1d <- c(c(bs$Std..Error[2]),c(bs$t.value[2]),c(bs$Pr...t..[2]))
ta <- qt(a/2,n-2,lower.tail = F)
ib1 <- c(b1-(ta*b1d[1]),b1+(ta*b1d[1]))
return(ib1)
}
icb1(x,y,0.05)
## [1] -36376.47 -28552.97