Taller Regresión Lineal
Katherin Quiñones % Lizeth Suarez
Agosto 2022
Taller de Regresión Lineal
Punto 1.
Se busca realizar la predicción del precio de las acciones de Ecopetrol a partir de la variación del precio del barril de petróleo WTI producido en colombia. Para esto se plantea realizar una regresión lineal.
Se van a considerar los siguietnes datos para su construcción.
library(readr)
datos <- read.csv("C:/Users/kaqu/Documents/Maestria Ciencia de Datos/Segundo_semestre/metodos_estadisticos/DatosTaller.csv",sep=";", dec = ",")
head(datos)## ï..Fecha PrecioAccionesEcopetrol PrecioPetroleoWTI
## 1 dic-14-15 1090 35.62
## 2 dic-15-15 1170 36.31
## 3 dic-16-15 1160 37.35
## 4 dic-18-15 1230 34.95
## 5 dic-21-15 1155 34.53
## 6 dic-22-15 1165 35.81Inicialmente se analizara de manera descriptiva el comportamiento de ambas variables
datos$PrecioAccionesEcopetrol = as.numeric(datos$PrecioAccionesEcopetrol)
datos$PrecioPetroleoWTI = as.numeric(datos$PrecioPetroleoWTI)
m = round(mean(datos[,"PrecioAccionesEcopetrol"]),2)
sd = round(sd(datos[,"PrecioAccionesEcopetrol"]),2)
print(paste("PrecioAccionesEcopetrol", "Media =",m,"y sd =", sd))## [1] "PrecioAccionesEcopetrol Media = 1108.39 y sd = 78.42"Encontramos que en promedio las acciones de Ecopetrol cuestan $1.108 con una desviación estandar de $78.42.
m = round(mean(datos[,"PrecioPetroleoWTI"]),2)
sd = round(sd(datos[,"PrecioPetroleoWTI"]),2)
print(paste("PrecioPetroleoWTI", "Media =",m,"y sd =", sd))## [1] "PrecioPetroleoWTI Media = 35.53 y sd = 2.12"En cuanto al valor medio de un barril de petróleo esta en 35.53 dolares con una desviación estandar de 2.12 dolares.
Con la finalidad de observar la relación de estas dos variables, se grafica un diagrama de dispesión para ver si es posible explicar la variabilidad del precio de las acciones de Ecopetrol mediante la variación del precio del barril de petróleo.
library(ggplot2)
ggplot(datos, aes(x = PrecioPetroleoWTI, y = PrecioAccionesEcopetrol)) +
geom_point() +
geom_smooth()
Se observa, una relación entre esta dos variables, que inicialmente pareciera no ser necesariamente lineal. Sin embargo al analizar el coeficiente de correlación de Pearson este nos da del 70.74% con lo que podriamos considerar que podemos generar la regresión lineal inicialmente.
cor(datos$PrecioAccionesEcopetrol, datos$PrecioPetroleoWTI)## [1] 0.7074373a. Estimación de la Regresión Lineal
Procedemos a estructurar la regresión de la siguiente forma
\[PrecioAccionesEcopetro = \beta_0 + \beta_1 PrecioPetroleoWTI + \epsilon\]
lm = lm(datos$PrecioAccionesEcopetrol~datos$PrecioPetroleoWTI)
coef(lm)## (Intercept) datos$PrecioPetroleoWTI
## 177.76779 26.19213Generando la estimación de la recta de la siguiente forma:
\[PrecioAccionesEcopetrol = 177.77+ 26.19* PrecioPetroleoWTI + \epsilon\]
Esta estimación de la recta presenta un \(R^2 = 50.05\), lo que nos dice que el modelo propuesto solo nos explica el 50% de la variabilidad que presenta el precio de las acciones de ecopetrol.
## [1] 0.5004675b. Significancia de Modelo
Al analizar la significacia del modelo, lo hacemos validando que el \(\hat{\beta}_1\) sea significativo, por tanto la hipotesis a evaluar son las siguientes:
\[H_0: \beta_1 = 0\] \[H_1: \beta_1 \neq 0\]
Al realizar la prueba de significancia observamos que el valor p asociado al coeficiente del precio del barril de petróleo es del 1.02%, lo que nos dice que contrastando este resultado con un nivel de significacia del 5%, no se cuenta con evidencia estadistica suficiente para suponer que el \(\beta_1\) estimado no sea significativo.
summary(lm)##
## Call:
## lm(formula = datos$PrecioAccionesEcopetrol ~ datos$PrecioPetroleoWTI)
##
## Residuals:
## Min 1Q Median 3Q Max
## -59.90 -40.74 -15.94 33.40 136.82
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 177.768 232.828 0.764 0.45627
## datos$PrecioPetroleoWTI 26.192 6.542 4.004 0.00102 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 57.13 on 16 degrees of freedom
## Multiple R-squared: 0.5005, Adjusted R-squared: 0.4692
## F-statistic: 16.03 on 1 and 16 DF, p-value: 0.001024c. Interpretacion de los betas
En el caso estudiado vemos que los betas estimados se pueden interpretar de la siguiente forma:
\(\beta_0\): En el caso dado donde el precio del barril de petróleo llegase a ser gratis, el valor de las acciones de Ecopetrol costarían $ 177.77.
\(\beta_1\): Este nos quiere decir que por cada dolar adicional que sume el valor del barril de petróleo el precio de las acciones de Ecopetrol incrementaran en $26.192.
d. Analisis de los supuestos del modelo
Linealidad
datos$yest = predict(lm)
ggplot(datos,aes(x = datos$PrecioPetroleoWTI, y = yest)) +
geom_point() +
geom_smooth() +
labs(title = paste("Media residuos =", round(mean(lm$residuals),3)))
Independencia
Las hipotensis a contrastar son:
\(H_0:\) Los errores no presentan autocorrelación
\(H_1:\) Los errores presentan autocorrelación
library(lmtest)
dw = dwtest(lm, alternative = "two.sided", iterations = 1000)
plot(lm$residuals, main = paste("dw P-Value = ",round(dw$p.value,3)))
Varianza Constante
Las hipotensis a contrastar son:
\(H_0:\) Los errores tiene homocedasticidad de varianza
\(H_1:\) Los errores no tiene homocedasticidad de varianza
grupos<-numeric()
grupos[datos$yest<1101]<-1
grupos[datos$yest>=1001 & datos$yest<1101]<-2
grupos[datos$yest>=1101]<-3
table(grupos)## grupos
## 1 2 3
## 2 4 12bt = bartlett.test(lm$residuals~grupos)
plot(datos$yest, lm$residuals,
main = paste("dt P-Value = ",round(bt$p.value,3)))
Normalidad
Las hipotesis a contrastar son:
\(H_0:\) Los errores presentan una distribución normal
\(H_1:\) Los errores no presentan una distribución normal
par(mfrow=c(1,2))
qqnorm(lm$residuals)
qqline(lm$residuals)
hist(lm$residuals)
shapiro.test(lm$residuals)##
## Shapiro-Wilk normality test
##
## data: lm$residuals
## W = 0.89259, p-value = 0.04276e. Validez del Modelo
Analizando el modelo propuesto, observamos que usando unicamente el precio del petroleo no podemos llegar a predecir las acciones dado que el modelo propuesto solo explica el 50% de la informacion el precio de las acciones, además el modelo propuesto no cumple con la normalidad, por tanto, se podría explorar una transformación de la variables con la finalidad de poder ver si este supuesto se puede corregir.
Punto 2.
Se busca realizar la predicción del Salario Minimo (SMMLV) a partir de la variación de la inflación en colombia. Para esto se plantea realizar una regresión lineal.
Se van a considerar los siguientes datos para su construcción.
datos <- read.csv("C:/Users/kaqu/Documents/Maestria Ciencia de Datos/Segundo_semestre/metodos_estadisticos/DatosTallerPunto2.csv", sep=";", dec = ",")
head(datos)## ï..ANIO INFLACION SMLM
## 1 1999 9.23 236460
## 2 2000 8.75 260100
## 3 2001 7.65 286000
## 4 2002 6.99 309000
## 5 2003 6.49 332000
## 6 2004 5.50 358000Inicialmente se analizara de manera descriptiva el comportamiento de ambas variables
m = round(mean(datos[,"INFLACION"]),2)
sd = round(sd(datos[,"INFLACION"]),2)
print(paste("INFLACION", "Media =",m,"y sd =", sd))## [1] "INFLACION Media = 5.35 y sd = 2.32"Encontramos que en promedio la inflación en Colombia es de 5.35% entre 1999 y 2015 con una desviación estandar de 2.32%.
m = round(mean(datos[,"SMLM"]),2)
sd = round(sd(datos[,"SMLM"]),2)
print(paste("SMLM", "Media =",m,"y sd =", sd))## [1] "SMLM Media = 437078.65 y sd = 129182.57"En cuanto al Salario Minimo en Colombia entre 1999 y 2015, el valor promedio es de $ 437.078.7 con una desviación estandar de $ 129.183.
Con la finalidad de observar la relación de estas dos variables, se realiza un diagrama de dispesión con el fin de ver si es posible explicar la variabilidad del salario minimo en funcion de la inflación anual.
Se observa una relación inversa entre la inflación y el salario minimo en Colombia, es decir, entre mayor es la inflación menor ha sido el salario minimo.
library(ggplot2)
ggplot(datos, aes(x = INFLACION, y = SMLM)) +
geom_point() +
geom_smooth()
Al analizar el coeficiente de correlación de Pearson este nos da -70.87% con lo que podriamos considerar que podemos generar la regresión lineal inicialmente.
cor(datos$INFLACION, datos$SMLM)## [1] -0.7086581Estimación de la Regresión Lineal
Procedemos a estructurar la regresión de la siguiente forma
\[SMLM = \beta_0 + \beta_1 Inflación + \epsilon\]
lm = lm(datos$SMLM~datos$INFLACION)
coef(lm)## (Intercept) datos$INFLACION
## 648485.93 -39489.33Generando la estimación de la recta de la siguiente forma:
\[SMLM = 648.485,93 - 39.489,33*Inflacion + \epsilon\]
Significancia del Modelo
Al analizar la significacia del modelo, lo hacemos validando que el \(\hat{\beta}_1\) sea significativo, por tanto la hipotesis a evaluar son las siguientes:
\[H_0: \beta_1 = 0\] \[H_1: \beta_1 \neq 0\]
Al realizar la prueba de significancia observamos que el valor p asociado al coeficiente de la inflación es de 1.45%, lo que nos dice que contrastando este resultado con un nivel de significacia del 5%, no se cuenta con evidencia estadistica suficiente para suponer que el \(\beta_1\) estimado no sea significativo.
summary(lm)##
## Call:
## lm(formula = datos$SMLM ~ datos$INFLACION)
##
## Residuals:
## Min 1Q Median 3Q Max
## -75463 -63456 -42854 17623 263207
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 648486 58947 11.00 1.4e-08 ***
## datos$INFLACION -39489 10151 -3.89 0.00145 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 94130 on 15 degrees of freedom
## Multiple R-squared: 0.5022, Adjusted R-squared: 0.469
## F-statistic: 15.13 on 1 and 15 DF, p-value: 0.00145Coeficiente de Correlacion
Esta estimación de la recta presenta un \(R^2 = 50.22\), lo que nos dice que el modelo propuesto solo nos explica el 50% de la variabilidad que presenta el salario minimo.
s = summary(lm)
s$r.squared## [1] 0.5021963### Interpretacion de los betas
En el caso estudiado vemos que los betas estimados se pueden interpretar de la siguiente forma:
\(\beta_0\): En el caso dado donde la inflación llegase a ser 0, el valor del ingreso minimo sería de $ 648.485,93.
\(\beta_1\): Este nos quiere decir que por cada unidad porcentual adicional que incremente la inflación, el salario minimo disminuirá en $39.489,33.
Analisis de los supuestos del modelo
Linealidad
El modelo cumple con el supuesto de linealidad, dado que los estimaciones realizados por este, se encuentran cercanas a los valores reales obervados.
datos$yest = predict(lm)
ggplot(datos,aes(x = datos$INFLACION, y = yest)) +
geom_point() +
geom_smooth() +
labs(title = paste("Media residuos =", round(mean(lm$residuals),3)))
Independencia
Frente a la independecia, pareciera que existiera alguna correalción en el tiempo dado los resultado de los residuales, sin embargo, se podria realizar una prueba formal para validadr esta hipotesis.
plot(lm$residuals)
Varianza Constante
Frente a la homogeneidad de varianza, observamos que si bien tenemos pocos datos parecietan que esto no presentan una tendencia en los datos, por tanto, el modelo cumple con el supuesto de homogeneidad de varainza.
plot(datos$yest, lm$residuals)
Normalidad
Analizando las graficas del comportamiento de los residuales en el QQ plot y en el histograma, podemos observar que el modelo propuesto no cumple con este supuesto. Por tanto, no tenemos certezas sobre la significancia del modelo propuesto, dado que no podemos aplicar prueba de hipotesis sobre los betas estimados.
par(mfrow=c(1,2))
qqnorm(lm$residuals)
qqline(lm$residuals)
hist(lm$residuals)
Convenencia de usar el modelo
Dado los resultado presentados anterior, el modelo propuesto no lo recomendaria para predecir el valor de SMML, dado que este factor depende de muchas otras cosas más, y no unicamente de la inflacion y esto lo podemos evidenciar en los resultados obtendidos donde el \(R^2\) es del alrededor del 50%. Además de la inconsistencia con el comportamiento normal de las variables, dado que se esperaria que a mayor inflación el salario minimo fuera mayor. Quizas se podría tratar de modelar el incremente del salario minimo en lugar del valor neto.
Punto 3
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),el area de la vivienda (m2), entre otras variables de interes.
library(readxl)
datos_vivienda <- read_excel("C:/Users/kaqu/Documents/Maestria Ciencia de Datos/Segundo_semestre/metodos_estadisticos/Datos_Vivienda.xlsx")
attach(datos_vivienda)
names(datos_vivienda)## [1] "Zona" "piso" "Estrato"
## [4] "precio_millon" "Area_contruida" "parqueaderos"
## [7] "Banos" "Habitaciones" "Tipo"
## [10] "Barrio" "cordenada_longitud" "Cordenada_latitud"Con base en los datos de precios de vivienda de la actividad en clase realizar un informe que contenga los siguientes puntos utilizando R y RMarkdown (publicar el informe final en Rpubs presentando código, resultados e interpretaciones).
a. Realice un filtro a la base de datos e incluya solo las ofertas de apartamentos, de la zona norte de la ciudad con precios inferiores a los 500 millones de pesos y áreas menores a 300 mt2. Presente los primeros 3 registros de la base y algunas tablas que comprueben la consulta. (¿Adicional un mapa con los puntos de la base, discutir si todos los puntos se ubican en la zona norte o se presentan valores en otras zonas, por qué?).
Base antes de hacer fitros
head(datos_vivienda)## # A tibble: 6 x 12
## Zona piso Estrato precio_millon Area_contruida parqueaderos Banos
## <chr> <chr> <dbl> <dbl> <dbl> <chr> <dbl>
## 1 Zona Sur 2 6 880 237 2 5
## 2 Zona Oeste 2 4 1200 800 3 6
## 3 Zona Sur 3 5 250 86 NA 2
## 4 Zona Sur NA 6 1280 346 4 6
## 5 Zona Sur 2 6 1300 600 4 7
## 6 Zona Sur 3 6 513 160 2 4
## # ... with 5 more variables: Habitaciones <dbl>, Tipo <chr>, Barrio <chr>,
## # cordenada_longitud <dbl>, Cordenada_latitud <dbl>A continuación se miraran algunas descriptivas relacionadas a las variables que se debe hacer filtros.
|
|
Se puede observar que en el tipo de vivienda predominan los apartamentos, con cerca de un 62% de los registros analizados. Hay presencia de valores faltantes con 3 registros, lo que podría evaluarse si fuese necesario, no contar con ellos. En cuanto a la zona, la zona a la que pertenecen la mayoria de viviendas analizadas es a la zona sur y en menor cantidad la zona centro con el 1.5 %.
|
|
En las dos variables cuantitativas, parece haber presencia de datos atípicos. El área promedio de las viviendas correspondientes a esta muestra de 175 m2, mientras que el 50% de las viviendas tiene un precio hasta 330 millones de pesos.
Filtrando Variables
dim(datos_vivienda)[1]## [1] 8322vars <-c('Tipo', 'Zona', 'Area_contruida', 'precio_millon')
cond <-c('Apartamento', 'Zona Norte', 300, 500)
datos_vivienda_sub <- datos_vivienda %>% filter(
.data[[vars[[1]]]]== cond[[1]],
.data[[vars[[2]]]]== cond[[2]],
.data[[vars[[3]]]]< cond[[3]],
.data[[vars[[4]]]]< cond[[4]]
)
dim(datos_vivienda_sub)[1]## [1] 315Inicialmente se cuenta con 8322 registros, después de realizar los 4 filtros sugeridos, quedamos con 315 registros, a continuación se presenta un encabezado del dataframe resultante:
head(datos_vivienda_sub)## # A tibble: 6 x 12
## Zona piso Estrato precio_millon Area_contruida parqueaderos Banos
## <chr> <chr> <dbl> <dbl> <dbl> <chr> <dbl>
## 1 Zona Norte NA 5 340 106 2 2
## 2 Zona Norte 1 3 135 103 1 3
## 3 Zona Norte NA 5 390 102 2 3
## 4 Zona Norte 11 5 350 137 2 3
## 5 Zona Norte NA 5 430 105 NA 3
## 6 Zona Norte NA 3 110 100 NA 1
## # ... with 5 more variables: Habitaciones <dbl>, Tipo <chr>, Barrio <chr>,
## # cordenada_longitud <dbl>, Cordenada_latitud <dbl>b. Realice un análisis exploratorio de datos enfocado en la correlación entre la variable respuesta (precio del apartamento) en función del área construida, estrato y si tiene parqueadero. Use gráficos interactivos con plotly e interprete los resultados.
Antes de realizar el análiss exploratorio sugerido se plantea hacer algunas recodificaciones necesarias en las variables de interes.
Revisión de datos Faltantes
datos_vivienda_sub$parqueaderos = as.numeric(datos_vivienda_sub$parqueaderos)
apply(is.na(datos_vivienda_sub), 2, sum)## Zona piso Estrato precio_millon
## 0 0 0 0
## Area_contruida parqueaderos Banos Habitaciones
## 0 46 0 0
## Tipo Barrio cordenada_longitud Cordenada_latitud
## 0 0 0 0Podemos observar como la variable de parqueaderos contiene 46 datos faltantes, en este caso se asumira esta falta de información como si no tuviera parqueadero, con el fin, de poder recategorizar esta variable mas adelante, en con parqueadero o sin parqueadero. Cabe aclarar que esto representa casi el 15% sobre los registros resultantes después de los correspondientes filtros.
datos_vivienda_sub$parqueaderos[is.na(datos_vivienda_sub$parqueaderos)] <- 0
apply(is.na(datos_vivienda_sub), 2, sum)## Zona piso Estrato precio_millon
## 0 0 0 0
## Area_contruida parqueaderos Banos Habitaciones
## 0 0 0 0
## Tipo Barrio cordenada_longitud Cordenada_latitud
## 0 0 0 0datos_vivienda_sub$parqueaderos = as.numeric(datos_vivienda_sub$parqueaderos)
datos_vivienda_sub$Estrato = as.factor(datos_vivienda_sub$Estrato)
datos_vivienda_sub$parqueaderos_cat=cut(datos_vivienda_sub$parqueaderos, breaks = c(-1,0,4), include.lowest = F, labels = c("sin_parqueadero","con_parqueadero"))
tabla_parque = tabla_freq(datos_vivienda_sub$parqueaderos_cat)
tabla_estrato = tabla_freq(datos_vivienda_sub$Estrato)
kbl(list(tabla_parque, tabla_estrato),
caption = "<center><strong>Distribución Parqueadero y Estrato</strong></center>") %>%
kable_paper(bootstrap_options = "striped")
|
|
De las anteriores tablas construidas, podemos decir que apróximadamente el 85% de las viviendas analizadas cuentan con parqueadero y el estrato donde pertenecen en su mayoria es el estrato 5 con el 72% de la base depurada.
Tanto el precio de la vivienda como el área construida, evidencian una distribución asimétrica positiva, los gráficos anteriormente mostrados, lo respaldan.
Análisis Bivariado
A continuación se mostrarán algunos gráficos interactivos, entre las variables de interés con el fin de poder establecer alguna relación apriori, es decir de forma descriptiva(visual).
library(plotly)
fig1 <- plot_ly(data = datos_vivienda_sub, x =~Area_contruida , y = ~precio_millon ,
marker = list(size = 10,
color = 'lightblue',
line = list(color = 'blue',
width = 2)))
fig1 <- fig1 %>% layout(title = 'Precio Vivienda vs Area Construida',
yaxis = list(zeroline = FALSE),
xaxis = list(zeroline = FALSE))
fig1fig2 <- plot_ly(datos_vivienda_sub, y = ~precio_millon, color = ~Estrato, type = "box",colors = "Set1")
fig2 <- fig2 %>% layout(title = 'Precio Vivienda vs Estrato',
yaxis = list(zeroline = FALSE),
xaxis = list(zeroline = FALSE))
fig2Se puede observar que aunque parece haber una tendencia de que cuando se aumenta de estrato, el precio de la vivienda aumenta, hay una presencia de valores atipicos en el estrato 6 donde se alcanzan valores de 1300 millones.
datos_vivienda_sub$parqueaderos_cat = as.factor(datos_vivienda_sub$parqueaderos_cat)
fig3 <- plot_ly(datos_vivienda_sub, y = ~precio_millon, color = ~parqueaderos_cat, type = "box",colors = "Set1")
fig3 <- fig3 %>% layout(title = 'Precio Vivienda vs Parqueaderos',
yaxis = list(zeroline = FALSE),
xaxis = list(zeroline = FALSE))
fig3De forma descriptiva, parece no haber una diferencia importante en las medianas del precio de la vivienda, cuando esta tiene o no tiene parqueadero. Es más parece que cuando tiene parqueadero el 50% de los datos indicara un precio menor.
A continuación se llevara a cabo un pre-procesamiento de la base, con el fin de volver dummies las variables categóricas, para ser posteriormente ingresadas al modelo.
library('fastDummies')
datos_vivienda_sub1 <- dummy_cols(datos_vivienda_sub, select_columns = 'parqueaderos_cat')
datos_vivienda_sub1 <- dummy_cols(datos_vivienda_sub1, select_columns = 'Estrato')
head(datos_vivienda_sub1)## # A tibble: 6 x 19
## Zona piso Estrato precio_millon Area_contruida parqueaderos Banos
## <chr> <chr> <fct> <dbl> <dbl> <dbl> <dbl>
## 1 Zona Norte NA 5 340 106 2 2
## 2 Zona Norte 1 3 135 103 1 3
## 3 Zona Norte NA 5 390 102 2 3
## 4 Zona Norte 11 5 350 137 2 3
## 5 Zona Norte NA 5 430 105 0 3
## 6 Zona Norte NA 3 110 100 0 1
## # ... with 12 more variables: Habitaciones <dbl>, Tipo <chr>, Barrio <chr>,
## # cordenada_longitud <dbl>, Cordenada_latitud <dbl>, parqueaderos_cat <fct>,
## # parqueaderos_cat_sin_parqueadero <int>,
## # parqueaderos_cat_con_parqueadero <int>, Estrato_3 <int>, Estrato_4 <int>,
## # Estrato_5 <int>, Estrato_6 <int>### c. Estime un modelo de regresión lineal múltiple con las variables del punto anterior e interprete los coeficientes si son estadísticamente significativos. Las interpretaciones deber están contextualizadas y discutir si los resultados son lógicos. Adicionalmente interprete el coeficiente R2 y discuta el ajuste del modelo e implicaciones (que podrían hacer para mejorarlo).
ml=lm(precio_millon~ Area_contruida+Estrato+parqueaderos_cat,data=datos_vivienda_sub)
summary(ml)##
## Call:
## lm(formula = precio_millon ~ Area_contruida + Estrato + parqueaderos_cat,
## data = datos_vivienda_sub)
##
## Residuals:
## Min 1Q Median 3Q Max
## -329.46 -67.87 -14.19 72.27 638.59
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -17.7088 61.5235 -0.288 0.7737
## Area_contruida 1.6034 0.1958 8.188 7.09e-15 ***
## Estrato4 77.5480 58.8863 1.317 0.1888
## Estrato5 133.2809 56.5569 2.357 0.0191 *
## Estrato6 311.9427 58.6761 5.316 2.03e-07 ***
## parqueaderos_catcon_parqueadero 30.9780 20.4741 1.513 0.1313
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 123 on 309 degrees of freedom
## Multiple R-squared: 0.4167, Adjusted R-squared: 0.4073
## F-statistic: 44.15 on 5 and 309 DF, p-value: < 2.2e-16\(H_0 :\beta_1,.., \beta_6=0\)
\(H_a :\beta_1,.., \beta_6\neq 0\)
Se Puede observar que \(P_{value}\) del estadístico F es \(2.2e ^{-16}\), es menor a \(\alpha=0.05\) entonces se rechaza \(H_0\), y se dice que el modelo podría ser valido.
Observamos que en el modelo inicial reflejan variables que no son “significativas”; no obstante esto puede ser resultante de un problema de multicolinealidad presentado por las variables X; en donde vemos obligados a presentar y esquematizar otro planteamiento que elimine este hallazgo.Un ejemplo es el \(P_{value}\) asociado al beta de estrato 4 y a la variable asociada a la tenencia de parqueadero, eso indicaria, que podriamos encontrar un mejor modelo.
El coeficiente 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 42% de las variaciones del precio de las viviendas a travez del area construida, el estrato y el número de parqueaderos. Es un número pequeño, lo que podría indicar que se deberian hacer ciertas exclusiones o agregar otras variables que si expliquen mejor el precio de la vivienda.
d. Realice la validación de supuestos del modelo e interprete los resultados (no es necesario corregir en caso de presentar problemas solo realizar sugerencias de que se podría hacer).
Interpretación Gráfica
par(mfrow=c(2,2))
plot(ml)
Se puede observar que en el gráfico de residuales versus valores ajustados hay una presencia de valores extremos 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 observar que tiene una tendencia lineal podría pensar en normalidad de no ser por los valores atipicos presentados en el extremo.
Pruebas estadísticas para validación de supuestos
- \(E(e_i )= 0\)
mean(residuals(ml))## [1] 7.63763e-16El valor esperado (promedio) de los residuales es igual a 0.
- Supuesto de Normalidad
\(H_0 : X \sim N(\mu,\sigma^2)\)
\(H_a : X \nsim N(\mu,\sigma^2)\)
shapiro.test(ml$residuals)##
## Shapiro-Wilk normality test
##
## data: ml$residuals
## W = 0.90077, p-value = 1.659e-13Com \(P-value\) es menor a a 0.05 (nivel de significancia escogido), se rechaza \(H_0\), entonces podría pensar que los errores no siguen una distribución normal.
- 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(ml)##
## studentized Breusch-Pagan test
##
## data: ml
## BP = 129.59, df = 5, p-value < 2.2e-16Como $ P-value $ es menor a 0.05 (nivel de significancia escogido), no se rechaza \(H_0\), entonces se podría pensar que los errores cumplen con el supuesto de homocedasticidad.
- 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(ml, alternative = "two.sided")## lag Autocorrelation D-W Statistic p-value
## 1 0.1293537 1.739065 0.024
## Alternative hypothesis: rho != 0Como $ P-value $ es menor a 0.05 (nivel de significancia escogido), n se rechaza \(H_0\), entonces se podría pensar que los errores podrian estar autocorrelacionados.
En este caso se viola el supuesto de normalidad y autocorrelación, por ende se podría pensar en realizar alguna transformación sobre la variable respuesta o inclunsive sobre la variable explicativa. Quizas evaluar el logaritmo o la raiz cuadrada sobre la variable respuesta pudiese ser una buena idea.
e. Con el modelo identificado predecir el precio de un apartamento con 100 mt2, de estrato 4 y con parqueadero. ¿Si este apartamento lo están ofreciendo en 450 millones cual seria su opinión con base en el resultado del modelo considera que es una buena oferta?
A continuación buscamos los coeficientes asociados a los betas generados por los modelos.
ml$coefficients## (Intercept) Area_contruida
## -17.708836 1.603372
## Estrato4 Estrato5
## 77.547984 133.280900
## Estrato6 parqueaderos_catcon_parqueadero
## 311.942728 30.977981Así pues la ecuación nos queda de la siguiente manera:
\(\hat {precio\_apto}=-17.708836+1.603372 X_{area} + 77.547984 X_{estrato4} + 133.280900 X_{estrato5} + 311.942728 X_{estrato6} +30.977981 X_{con_parqueadero}\)
precio_pred=-17.708836+1.603372*100+ 77.547984*1 + 133.280900*0 + 311.942728*0 +30.977981*1
precio_pred## [1] 251.1543Con base en el modelo propuesto, la oferta del apartamento en 450 millones no parece ser una buena oferta, ya que el precio sugerido, esta casi por la mitad del valor solicitado. (251 millones de pesos).
f. Con las predicciones del modelo sugiera potenciales ofertas para una persona interesada en un apartamento en la zona norte con mas de 100 mt2 de área, de estrato 4, que tenga parqueadero y tenga encuentra que la persona tiene un crédito preaprobado de máximo 400 millones de pesos. Realice un análisis y presente en un mapa al menos 5 ofertas potenciales que debe discutir.
Inicialmente se debe depurar la base con las caracteristicas y condiciones solicitadas.
dim(datos_vivienda)[1]## [1] 8322vars <-c('Tipo', 'Zona', 'Area_contruida', 'precio_millon','parqueaderos')
cond <-c('Apartamento', 'Zona Norte', 100, 400,1)
datos_vivienda_sub_pred <- datos_vivienda %>% filter(
.data[[vars[[1]]]]== cond[[1]],
.data[[vars[[2]]]]== cond[[2]],
.data[[vars[[3]]]]> cond[[3]],
.data[[vars[[4]]]]<= cond[[4]],
.data[[vars[[5]]]]>= cond[[5]],
)
ID=1:dim(datos_vivienda_sub_pred)[1]
datos_vivienda_sub_pred=data.frame(ID,datos_vivienda_sub_pred)
dim(datos_vivienda_sub_pred)[1]## [1] 879head(datos_vivienda_sub_pred)## ID Zona piso Estrato precio_millon Area_contruida parqueaderos Banos
## 1 1 Zona Norte 2 3 135 56 1 1
## 2 2 Zona Norte NA 5 340 106 2 2
## 3 3 Zona Norte 1 3 135 103 1 3
## 4 4 Zona Norte NA 4 175 77 1 2
## 5 5 Zona Norte 2 3 110 57 1 2
## 6 6 Zona Norte 3 3 155 62 1 2
## Habitaciones Tipo Barrio cordenada_longitud
## 1 3 Apartamento torres de comfandi -76.46745
## 2 3 Apartamento la flora -76.48200
## 3 4 Apartamento calimio norte -76.48347
## 4 3 Apartamento urbanizaciv=n pacara -76.48395
## 5 3 Apartamento puente del comercio -76.48458
## 6 3 Apartamento villa del prado -76.48647
## Cordenada_latitud
## 1 3.40763
## 2 3.43500
## 3 3.48626
## 4 3.45104
## 5 3.46987
## 6 3.46549datos_vivienda_sub_pred$parqueaderos = as.numeric(datos_vivienda_sub_pred$parqueaderos)
datos_vivienda_sub_pred$Estrato = as.factor(datos_vivienda_sub_pred$Estrato)
datos_vivienda_sub_pred$parqueaderos[is.na(datos_vivienda_sub_pred$parqueaderos)] <- 0
datos_vivienda_sub_pred$parqueaderos_cat=cut(datos_vivienda_sub_pred$parqueaderos, breaks = c(-1,0,4), include.lowest = F, labels = c("sin_parqueadero","con_parqueadero"))
library('fastDummies')
datos_vivienda_sub_pred1 <- dummy_cols(datos_vivienda_sub_pred, select_columns = 'parqueaderos_cat')
datos_vivienda_sub_pred1 <- dummy_cols(datos_vivienda_sub_pred1, select_columns = 'Estrato')
head(datos_vivienda_sub_pred1)## ID Zona piso Estrato precio_millon Area_contruida parqueaderos Banos
## 1 1 Zona Norte 2 3 135 56 1 1
## 2 2 Zona Norte NA 5 340 106 2 2
## 3 3 Zona Norte 1 3 135 103 1 3
## 4 4 Zona Norte NA 4 175 77 1 2
## 5 5 Zona Norte 2 3 110 57 1 2
## 6 6 Zona Norte 3 3 155 62 1 2
## Habitaciones Tipo Barrio cordenada_longitud
## 1 3 Apartamento torres de comfandi -76.46745
## 2 3 Apartamento la flora -76.48200
## 3 4 Apartamento calimio norte -76.48347
## 4 3 Apartamento urbanizaciv=n pacara -76.48395
## 5 3 Apartamento puente del comercio -76.48458
## 6 3 Apartamento villa del prado -76.48647
## Cordenada_latitud parqueaderos_cat parqueaderos_cat_sin_parqueadero
## 1 3.40763 con_parqueadero 0
## 2 3.43500 con_parqueadero 0
## 3 3.48626 con_parqueadero 0
## 4 3.45104 con_parqueadero 0
## 5 3.46987 con_parqueadero 0
## 6 3.46549 con_parqueadero 0
## parqueaderos_cat_con_parqueadero Estrato_3 Estrato_4 Estrato_5 Estrato_6
## 1 1 1 0 0 0
## 2 1 0 0 1 0
## 3 1 1 0 0 0
## 4 1 0 1 0 0
## 5 1 1 0 0 0
## 6 1 1 0 0 0datos_vivienda_sub_pred1=datos_vivienda_sub_pred1[,c("ID","precio_millon","Area_contruida","parqueaderos_cat_con_parqueadero","Estrato_4","Estrato_5","Estrato_6",'cordenada_longitud','Cordenada_latitud','Barrio')]
head(datos_vivienda_sub_pred1)## ID precio_millon Area_contruida parqueaderos_cat_con_parqueadero Estrato_4
## 1 1 135 56 1 0
## 2 2 340 106 1 0
## 3 3 135 103 1 0
## 4 4 175 77 1 1
## 5 5 110 57 1 0
## 6 6 155 62 1 0
## Estrato_5 Estrato_6 cordenada_longitud Cordenada_latitud Barrio
## 1 0 0 -76.46745 3.40763 torres de comfandi
## 2 1 0 -76.48200 3.43500 la flora
## 3 0 0 -76.48347 3.48626 calimio norte
## 4 0 0 -76.48395 3.45104 urbanizaciv=n pacara
## 5 0 0 -76.48458 3.46987 puente del comercio
## 6 0 0 -76.48647 3.46549 villa del pradovars <-c('Estrato_4','parqueaderos_cat_con_parqueadero','Area_contruida')
cond <-c(1,1,100)
datos_vivienda_sub_pred1 <- datos_vivienda_sub_pred1 %>% filter(
.data[[vars[[1]]]]== cond[[1]],
.data[[vars[[2]]]]== cond[[2]],
.data[[vars[[3]]]]>= cond[[3]]
)
print(dim(datos_vivienda_sub_pred1)[1])## [1] 25head(datos_vivienda_sub_pred1)## ID precio_millon Area_contruida parqueaderos_cat_con_parqueadero Estrato_4
## 1 327 380 123 1 1
## 2 515 350 130 1 1
## 3 524 290 108 1 1
## 4 583 185 104 1 1
## 5 602 265 125 1 1
## 6 631 380 126 1 1
## Estrato_5 Estrato_6 cordenada_longitud Cordenada_latitud Barrio
## 1 0 0 -76.51437 3.48618 la flora
## 2 0 0 -76.52100 3.49000 la flora
## 3 0 0 -76.52115 3.48930 la flora
## 4 0 0 -76.52300 3.46400 san vicente
## 5 0 0 -76.52353 3.48157 la flora
## 6 0 0 -76.52432 3.48254 la floraPara encontrar las mejores ofertas u oportunidades, se realiza una diferencia entre el precio predicho y el valor real del apartamento, con el fin de organizar de forma decreciente esta diferencia, y quedarnos con los primeros 5 registros.
datos_vivienda_sub_pred1['prediccion_precio']=round(-17.708836+1.603372*datos_vivienda_sub_pred1$Area_contruida+ 77.547984*datos_vivienda_sub_pred1$Estrato_4 + 133.280900*0 + 311.942728*0 +30.977981*datos_vivienda_sub_pred1$parqueaderos_cat_con_parqueadero)
datos_vivienda_sub_pred1['dif_precio']=datos_vivienda_sub_pred1['prediccion_precio']- datos_vivienda_sub_pred1['precio_millon']
datos_vivienda_sub_pred1=datos_vivienda_sub_pred1[order(datos_vivienda_sub_pred1$dif_precio, decreasing = TRUE), ]
head(datos_vivienda_sub_pred1)## ID precio_millon Area_contruida parqueaderos_cat_con_parqueadero Estrato_4
## 9 678 300 287.00 1 1
## 21 802 280 173.00 1 1
## 4 583 185 104.00 1 1
## 11 685 190 104.12 1 1
## 7 643 270 152.00 1 1
## 13 694 300 163.00 1 1
## Estrato_5 Estrato_6 cordenada_longitud Cordenada_latitud Barrio
## 9 0 0 -76.52673 3.47907 la campiv±a
## 21 0 0 -76.53362 3.46337 santa monica
## 4 0 0 -76.52300 3.46400 san vicente
## 11 0 0 -76.52707 3.46279 san vicente
## 7 0 0 -76.52515 3.46334 versalles
## 13 0 0 -76.52785 3.46701 san vicente
## prediccion_precio dif_precio
## 9 551 251
## 21 368 88
## 4 258 73
## 11 258 68
## 7 335 65
## 13 352 52require(leaflet)datos_vivienda_sub_pred1_top5=datos_vivienda_sub_pred1[1:5,]
leaflet() %>% addCircleMarkers(lng = datos_vivienda_sub_pred1_top5$cordenada_longitud,lat = datos_vivienda_sub_pred1_top5$Cordenada_latitud,radius = 0.3,color = "black",label = datos_vivienda_sub_pred1_top5$ID) %>% addTiles()Según lo visualizado en el mapa se puede observar que tienen una cercania, 4 de las 5 ofertas encontradas. A continuación se muestran los barrios a los cuales estan asociados estas ofertas y algunas caracteriticas.
|
PUnto 4
Con base en los datos de arboles proponga un modelo de regresión lineal múltiple que permita predecir el peso del árbol en función de las covariables que considere importantes y seleccionándolas de acuerdo con un proceso adecuado. Tenga en cuenta realizar una evaluación de la significancia de los parámetros, interpretación y proponga un método de evaluación por medio de validación cruzada. Presente métricas apropiadas como el RMSE y MAE.
library(readxl)
datos_arboles <- read_excel("C:/Users/kaqu/Documents/Maestria Ciencia de Datos/Segundo_semestre/metodos_estadisticos/data_arboles.xlsx",
col_types = c("text", "text", "numeric","numeric", "numeric"))
attach(datos_arboles)
print(names(datos_arboles))## [1] "finca" "mg" "peso" "diametro" "altura"head(datos_arboles)## # A tibble: 6 x 5
## finca mg peso diametro altura
## <chr> <chr> <dbl> <dbl> <dbl>
## 1 FINCA_1 GENOTIPO_1 13.7 4.7 5
## 2 FINCA_1 GENOTIPO_1 14.6 5.3 5.6
## 3 FINCA_1 GENOTIPO_1 15.9 4.8 5.8
## 4 FINCA_1 GENOTIPO_1 8.99 3.2 4.3
## 5 FINCA_1 GENOTIPO_1 6.99 2.2 3.3
## 6 FINCA_1 GENOTIPO_2 19.3 6.3 7.9Se puede observar que en el archivo dispuesto, se encuentran variables como la finca, el tipo de genotipo, estas como variables cualitativas y como variables cuantitativas las variables peso, diametro, y altura de los arboles.. A continuaciòn se hara la construcciòn de dos modelos de regresiòn simple usando las variables continua y concluiremos cual podrìa ser el mejor modelo.
Lo primero que se revisarà es si la base de datos contiene datos faltantes y de ser asì, se buscara hacer un tratamiento de ellos, previo a la construcciòn del ejercicio propuesto.
apply(is.na(datos_arboles), 2, sum)## finca mg peso diametro altura
## 0 0 0 0 0Como se puede observar, no existen datos faltantes en esta base.
|
|
De acuerdo a las caracteristicas cualitivas como finca y tipo de genotipo, existe la misma cantidad de arboles por finca (30) y de la misma forma esta distribuido equitativamente el genotipo (45 por cada uno).
|
|
|
La variable que puede representar en su mayoría datos atípicos mas fuertes puede ser el peso.
Análisis Bivariado
library(plotly)
fig1 <- plot_ly(data = datos_arboles, x =~altura , y = ~peso ,
marker = list(size = 10,
color = 'lightblue',
line = list(color = 'blue',
width = 2)))
fig1 <- fig1 %>% layout(title = 'Peso vs Altura de Àrbol',
yaxis = list(zeroline = FALSE),
xaxis = list(zeroline = FALSE))
fig1Cuando se grafica la altura vs el peso de los arboles, se encuentra una tendencia lineal positiva entre ambos, lo que indicaria es que a mayor altura del arbol, este pesaria mas.
library(plotly)
fig2 <- plot_ly(data = datos_arboles, x =~diametro , y = ~peso ,
marker = list(size = 10,
color = 'lightblue',
line = list(color = 'blue',
width = 2)))
fig2 <- fig2 %>% layout(title = 'Peso vs Diametro de Arbol',
yaxis = list(zeroline = FALSE),
xaxis = list(zeroline = FALSE))
fig2Así mismo se puede evidencia que el diametro y el peso tambien tienen esta tendencia, que podría implicar una relación fuerte y explicativa entre ambas variables.
Construcción Modelos:
Se construiran tres modelos, dos regresiones lineales simple y una múltiple haciendo uso de las variables cuantitativas.
Modelo 1: Altura vs Peso
m1=lm(peso ~altura,data=datos_arboles)
summary(m1)##
## Call:
## lm(formula = peso ~ altura, data = datos_arboles)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.228 -1.969 0.572 2.377 15.106
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -7.0456 1.7046 -4.133 8.14e-05 ***
## altura 3.8906 0.2481 15.684 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.211 on 88 degrees of freedom
## Multiple R-squared: 0.7365, Adjusted R-squared: 0.7335
## F-statistic: 246 on 1 and 88 DF, p-value: < 2.2e-16cor1=cor(datos_arboles$peso,datos_arboles$altura)
print('Correlación Peso vs Altura:',str(cor1))## num 0.858
## [1] "Correlación Peso vs Altura:"rcuad_m1=cor1*cor1
print('Coeficiente de Determinacion Peso vs Diametro:',str(rcuad_m1))## num 0.737
## [1] "Coeficiente de Determinacion Peso vs Diametro:"Modelo 2: Diametro vs Peso
m2=lm(peso ~diametro,data=datos_arboles)
summary(m2)##
## Call:
## lm(formula = peso ~ diametro, data = datos_arboles)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.3775 -2.6594 0.0237 1.8758 11.9876
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -9.0203 1.4129 -6.384 7.86e-09 ***
## diametro 5.1026 0.2508 20.346 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.435 on 88 degrees of freedom
## Multiple R-squared: 0.8247, Adjusted R-squared: 0.8227
## F-statistic: 414 on 1 and 88 DF, p-value: < 2.2e-16cor2=cor(datos_arboles$peso,datos_arboles$diametro)
rcuad_m2=cor2*cor2
print('Correlación Peso vs Diametro:',str(cor2))## num 0.908
## [1] "Correlación Peso vs Diametro:"print('Coeficiente de Determinacion Peso vs Diametro:',str(rcuad_m2))## num 0.825
## [1] "Coeficiente de Determinacion Peso vs Diametro:"Modelo 3: Diametro + Altura vs Peso
m3=lm(peso ~diametro + altura ,data=datos_arboles)
summary(m3)##
## Call:
## lm(formula = peso ~ diametro + altura, data = datos_arboles)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.3083 -2.5121 0.1608 2.0088 11.7446
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -9.1205 1.4305 -6.376 8.44e-09 ***
## diametro 4.7395 0.7128 6.649 2.49e-09 ***
## altura 0.3132 0.5751 0.544 0.587
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.449 on 87 degrees of freedom
## Multiple R-squared: 0.8253, Adjusted R-squared: 0.8213
## F-statistic: 205.5 on 2 and 87 DF, p-value: < 2.2e-16A grandes rasgos, se puede observar que en ambos modelos de regresión lineal simple, los dos modelos cumplen con los supuestos de los betas. En el modelo de regresión lineal múltiple el p valor asociado a la variable altura es mayor a 0.05 lo que significaria, que esta variable es candidata a salir del modelo ya que no significativa.
Validación de Supuestos Modelo 1:
- \(E(e_i )= 0\)
error_medio_m1=mean(residuals(m1))
error_medio_m1## [1] -1.28331e-16El valor esperado (promedio) de los residuales es igual a 0.
- Supuesto de Normalidad
\(H_0 : X \sim N(\mu,\sigma^2)\)
\(H_a : X \nsim N(\mu,\sigma^2)\)
norm_m1=shapiro.test(m1$residuals)
norm_m1##
## Shapiro-Wilk normality test
##
## data: m1$residuals
## W = 0.95439, p-value = 0.003157Com \(P-value\) es menor a a 0.05 (nivel de significancia escogido), se rechaza \(H_0\), entonces podría pensar que los errores no siguen una distribución normal.
- 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\)
homocedasticidad_m1=bptest(m1)
homocedasticidad_m1##
## studentized Breusch-Pagan test
##
## data: m1
## BP = 12.716, df = 1, p-value = 0.0003625Como $ 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.
- 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\)
autocor_m1=dwt(m1, alternative = "two.sided")
autocor_m1## lag Autocorrelation D-W Statistic p-value
## 1 0.537971 0.9021616 0
## Alternative hypothesis: rho != 0Como $ P-value $ es menor a 0.05 (nivel de significancia escogido), se rechaza \(H_0\), entonces se podría pensar que los errores podrian estar autocorrelacionados.
Validación de Supuestos Modelo 2:
- \(E(e_i )= 0\)
error_medio_m2=mean(residuals(m2))
error_medio_m2## [1] 2.553898e-17El valor esperado (promedio) de los residuales es igual a 0.
- Supuesto de Normalidad
\(H_0 : X \sim N(\mu,\sigma^2)\)
\(H_a : X \nsim N(\mu,\sigma^2)\)
norm_m2=shapiro.test(m2$residuals)
norm_m2##
## Shapiro-Wilk normality test
##
## data: m2$residuals
## W = 0.95356, p-value = 0.002793Com \(P-value\) es menor a a 0.05 (nivel de significancia escogido), se rechaza \(H_0\), entonces podría pensar que los errores no siguen una distribución normal.
- 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\)
homocedasticidad_m2=bptest(m2)
homocedasticidad_m2##
## studentized Breusch-Pagan test
##
## data: m2
## BP = 11.76, df = 1, p-value = 0.0006052Como $ 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.
- 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\)
autocor_m2=dwt(m2, alternative = "two.sided")
autocor_m2## lag Autocorrelation D-W Statistic p-value
## 1 0.4533746 1.071861 0
## Alternative hypothesis: rho != 0Como $ P-value $ es menor a 0.05 (nivel de significancia escogido), se rechaza \(H_0\), entonces se podría pensar que los errores podrian estar autocorrelacionados.
Validación de Supuestos Modelo 3:
- \(E(e_i )= 0\)
error_medio_m3=mean(residuals(m3))
error_medio_m3## [1] 1.032498e-16El valor esperado (promedio) de los residuales es igual a 0.
- Supuesto de Normalidad
\(H_0 : X \sim N(\mu,\sigma^2)\)
\(H_a : X \nsim N(\mu,\sigma^2)\)
norm_m3=shapiro.test(m3$residuals)
norm_m3##
## Shapiro-Wilk normality test
##
## data: m3$residuals
## W = 0.95745, p-value = 0.004966Com \(P-value\) es menor a a 0.05 (nivel de significancia escogido), se rechaza \(H_0\), entonces podría pensar que los errores no siguen una distribución normal.
- 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\)
homocedasticidad_m3=bptest(m3)
homocedasticidad_m3##
## studentized Breusch-Pagan test
##
## data: m3
## BP = 14.32, df = 2, p-value = 0.0007772Como $ 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.
- 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\)
autocor_m3=dwt(m3, alternative = "two.sided")
autocor_m3## lag Autocorrelation D-W Statistic p-value
## 1 0.4648495 1.048132 0
## Alternative hypothesis: rho != 0Como $ P-value $ es menor a 0.05 (nivel de significancia escogido), se rechaza \(H_0\), entonces se podría pensar que los errores podrian estar autocorrelacionados.
Revisión Criterio Akaike y Bayesiano
Los indicadores con criterio de información incorporan el compromiso entre la complejidad y la capacidad predictiva de un modelo. Cuanto más complejo sea un modelo peor será su capacidad para predecir en una eventual gama de situaciones. Es decir, cuantas más variables y cuantas más interrelaciones entre los componentes incorpore un modelo más especificas serán sus predicciones pero a su vez serán menos generalizables.
akaike_m1=AIC(m1)
bayesian_m1=BIC(m1)
print('Criterio Akaike m1',str(akaike_m1))## num 518
## [1] "Criterio Akaike m1"print('Criterio Bayesian m1',str(bayesian_m1))## num 526
## [1] "Criterio Bayesian m1"akaike_m2=AIC(m2)
bayesian_m2=BIC(m2)
print('Criterio Akaike m2',str(akaike_m2))## num 482
## [1] "Criterio Akaike m2"print('Criterio Bayesian m2',str(bayesian_m2))## num 489
## [1] "Criterio Bayesian m2"akaike_m3=AIC(m3)
bayesian_m3=BIC(m3)
print('Criterio Akaike m3',str(akaike_m3))## num 483
## [1] "Criterio Akaike m3"print('Criterio Bayesian m3',str(bayesian_m3))## num 493
## [1] "Criterio Bayesian m3"Generacion de valores aleatorios para prediccion
y=runif(90,min(datos_arboles$peso),max(datos_arboles$peso))
x_altura=runif(90,min(datos_arboles$altura),max(datos_arboles$altura))
x_diametro=runif(90,min(datos_arboles$diametro),max(datos_arboles$diametro))
datos_arboles2=data.frame(peso=y,altura=x_altura)
head(datos_arboles2)## peso altura
## 1 11.308581 4.201015
## 2 11.133383 7.732309
## 3 38.323084 5.481903
## 4 29.197223 4.508170
## 5 7.562935 8.025194
## 6 38.353914 9.617337Haciendo Prediccion Modelo 1
prediccion_m1=predict(m1,newdata = datos_arboles2 )
resultados_m1=data.frame(y_real=y,y_predict_m1=prediccion_m1)
head(resultados_m1)## y_real y_predict_m1
## 1 11.308581 9.298719
## 2 11.133383 23.037413
## 3 38.323084 14.282086
## 4 29.197223 10.493721
## 5 7.562935 24.176898
## 6 38.353914 30.371222library(Metrics)
attach(resultados_m1)
SCE_m1=sum((y_predict_m1-y)^2)
print('SCE Modelo 1',str(SCE_m1))## num 22709
## [1] "SCE Modelo 1"MAE_m1=mae(y_real,y_predict_m1)
print('MAE Modelo 1',str(MAE_m1))## num 13.2
## [1] "MAE Modelo 1"RMSE_m1=rmse(y_real,y_predict_m1)
print('RMSE Modelo 1',str(RMSE_m1))## num 15.9
## [1] "RMSE Modelo 1"Haciendo Prediccion Modelo 2
datos_arboles3=data.frame(peso=y,diametro=x_diametro)
head(datos_arboles3)## peso diametro
## 1 11.308581 2.727954
## 2 11.133383 6.766455
## 3 38.323084 2.511269
## 4 29.197223 4.528259
## 5 7.562935 6.495634
## 6 38.353914 3.414398prediccion_m2=predict(m2,newdata = datos_arboles3 )
resultados_m2=data.frame(y_real=y,y_predict_m2=prediccion_m2)
head(resultados_m2)## y_real y_predict_m2
## 1 11.308581 4.899337
## 2 11.133383 25.506104
## 3 38.323084 3.793684
## 4 29.197223 14.085535
## 5 7.562935 24.124218
## 6 38.353914 8.401974attach(resultados_m2)
SCE_m2=sum((y_predict_m2-y)^2)
print('SCE Modelo 2',str(SCE_m2))## num 31587
## [1] "SCE Modelo 2"MAE_m2=mae(resultados_m2$y_real,resultados_m2$y_predict_m2)
print('MAE Modelo 2',str(MAE_m2))## num 15.7
## [1] "MAE Modelo 2"RMSE_m2=rmse(resultados_m2$y_real,resultados_m2$y_predict_m2)
print('RMSE Modelo 2',str(RMSE_m2))## num 18.7
## [1] "RMSE Modelo 2"Haciendo Prediccion Modelo 3
datos_arboles4=data.frame(peso=y,diametro=x_diametro,altura=x_altura)
head(datos_arboles4)## peso diametro altura
## 1 11.308581 2.727954 4.201015
## 2 11.133383 6.766455 7.732309
## 3 38.323084 2.511269 5.481903
## 4 29.197223 4.528259 4.508170
## 5 7.562935 6.495634 8.025194
## 6 38.353914 3.414398 9.617337prediccion_m3=predict(m3,newdata = datos_arboles4)
resultados_m3=data.frame(y_real=datos_arboles4$peso,y_predict_m3=prediccion_m3)
head(resultados_m3)## y_real y_predict_m3
## 1 11.308581 5.124098
## 2 11.133383 25.370267
## 3 38.323084 4.498252
## 4 29.197223 13.752758
## 5 7.562935 24.178442
## 6 38.353914 10.073654attach(resultados_m3)
MAE_m3=mae(resultados_m3$y_real,resultados_m3$y_predict_m3)
print('MAE Modelo 3',str(MAE_m3))## num 15.2
## [1] "MAE Modelo 3"RMSE_m3=rmse(resultados_m3$y_real,resultados_m2$y_predict_m3)
print('RMSE Modelo 3',str(RMSE_m3))## num NaN
## [1] "RMSE Modelo 3"Validación Cruzada Modelo 1
library(rlang)
library(vctrs)
library(caret)
control <- trainControl(method = "cv", number = 5)
model_cv_1 <- train(peso ~ altura, data = datos_arboles, method = "lm", trControl = control)
print(model_cv_1)## Linear Regression
##
## 90 samples
## 1 predictor
##
## No pre-processing
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 71, 73, 73, 71, 72
## Resampling results:
##
## RMSE Rsquared MAE
## 4.091055 0.7522999 3.116309
##
## Tuning parameter 'intercept' was held constant at a value of TRUEValidación Cruzada Modelo 2
control <- trainControl(method = "cv", number = 5)
model_cv_2 <- train(peso ~ diametro, data = datos_arboles, method = "lm", trControl = control)
print(model_cv_2)## Linear Regression
##
## 90 samples
## 1 predictor
##
## No pre-processing
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 71, 73, 72, 73, 71
## Resampling results:
##
## RMSE Rsquared MAE
## 3.470151 0.8443777 2.854688
##
## Tuning parameter 'intercept' was held constant at a value of TRUEValidación Cruzada Modelo 3
control <- trainControl(method = "cv", number = 5)
model_cv_3 <- train(peso ~ diametro + altura, data = datos_arboles, method = "lm", trControl = control)
print(model_cv_3)## Linear Regression
##
## 90 samples
## 2 predictor
##
## No pre-processing
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 72, 72, 73, 72, 71
## Resampling results:
##
## RMSE Rsquared MAE
## 3.667319 0.8366192 2.962462
##
## Tuning parameter 'intercept' was held constant at a value of TRUEConsolidación de Resultados entre modelos propuestos:
Algunas definiciones de métricas utilizadas se presentan a continuación:
RMSE: La raíz del error cuadrático medio. Esto mide la diferencia promedio entre las predicciones hechas por el modelo y las observaciones reales. Cuanto menor sea el RMSE, más fielmente podrá un modelo predecir las observaciones reales.
Rsquared: Esta es una medida de la correlación entre las predicciones hechas por el modelo y las observaciones reales. Cuanto mayor sea el R-cuadrado, más fielmente podrá un modelo predecir las observaciones reales.
MAE: El error absoluto medio. Esta es la diferencia absoluta promedio entre las predicciones hechas por el modelo y las observaciones reales. Cuanto más bajo es el MAE, más cerca puede un modelo predecir las observaciones reales.
| num_modelo | modelo | correlaciones | coef_det | betas | mae | rmse | mae_cv | rmse_cv | normalidad | homocedasticidad | autocorrelacion |
|---|---|---|---|---|---|---|---|---|---|---|---|
| Modelo 1 | Altura vs Peso | 0.8582009 | 0.7365088 | Betas Significativos | 13.20260 | 15.88449 | 3.133 | 4.101 | No Cumple | No Cumple | No Cumple |
| Modelo 2 | Diametro vs Peso | 0.9081230 | 0.8246874 | Betas Significativos | 15.65322 | 18.73409 | 2.769 | 3.377 | No Cumple | No Cumple | No Cumple |
| Modelo 3 | Diametro + Altura vs Peso | NaN | 0.8253000 | B_2 No Significativo | 15.17024 | NaN | 2.854 | 3.531 | No Cumple | No Cumple | No Cumple |
A decir verdad, sin cumplimiento de supuestos no me quedaría con ningún modelo de los propuestos anteriormente. Para solucionar este problema se propone hacer algunas transformaciones a la variable respuesta y/o a las variables explicativas. Ahora si observamos los resultados mas desde el punto descriptivo pensariamos que el segundo modelo (diametro vs peso) podrías ser la mejor opción.