Tarea 4
Sebastián Lillo, Fabián González, Santeri Korkee
2022-06-04
Tarea 4
P1
Cargue los paquetes “data.table”, “ggplot2” y “caret”, junto con la base de datos.
rm(list=ls())
library(data.table)
library(ggplot2)
library(caret)
library(jtools)
library(scales)
base <- data.table(fread("base.csv", encoding = "Latin-1"))
base<- na.omit(base)
base<- base[complete.cases(base)]P2
Muestre un histograma de ingreso por comuna. Limite el histograma a las observaciones con ingresos menores a 2 millones de pesos.
P2 <- base[base$IngresoFinal<2000000]
ggplot(P2, aes(x=`comunahg`, y=`IngresoFinal`, fill=`comunahg`))+
geom_col(color="blue",size=1) +
theme_minimal() +labs(x="Comuna",
y="Ingresos", title="Ingresos por comuna", subtitle = "Quinta region", caption="Fuente: SEREMI de Transportes")+ theme(axis.text.x = element_text(angle=90, vjust=0.5)) 
P3
Realice un gráfico de dispersión que muestre la relación entre el ingreso y el tiempo de viaje promedio.
##**Ingrespos hasta 2 millones**
ggplot(P2, aes(x=IngresoFinal, y=minviajeprom))+ geom_point()+geom_smooth() + theme_minimal()+labs(x="Ingresos",
y="Viaje promedio", title="Relacion de viaje e ingresos", subtitle = "Ingresos de 2 millones", caption="Fuente: SEREMI de Transportes")+ theme(axis.text.x = element_text(angle=90, vjust=0.5))
P4
Un miembro de su equipo propone un modelo que calcule la incidencia del sexo y la educación en el tiempo de viaje, sin considerar variables adicionales. Usted quiere mostrarle que aquel modelo está incompleto. En orden de lograr esto, haga el modelo de regresión anteriormente mencionado.
f1 <- formula(minviajeprom ~ sexo + educ)
reg01 <- lm(formula = f1, data = base)
sum01<-summary(reg01)P4.1
Calcule los minutos de viaje promedios predichos para cada observación.
Por educación
lista <- c("No tiene", "Secundarios", "Primarios", "Profesional", "Técnico Profesional", "Prebásicos")
for (i in lista){
mujer <- predict(reg01, data.table(sexo = "Mujer", educ=i))
print(paste0("Estimación para una mujer edc:
"
, i," es ", mujer))
hombre <- predict(reg01, data.table(sexo ="Hombre", educ = i))
print(paste0("Estimación para un hombre edc:
"
, i, " es ", hombre))
Sys.sleep(1)
}## [1] "Estimación para una mujer edc:\nNo tiene es 19.0311206830537"
## [1] "Estimación para un hombre edc:\nNo tiene es 22.1271167261285"
## [1] "Estimación para una mujer edc:\nSecundarios es 28.5557754732689"
## [1] "Estimación para un hombre edc:\nSecundarios es 31.6517715163437"
## [1] "Estimación para una mujer edc:\nPrimarios es 23.651555611224"
## [1] "Estimación para un hombre edc:\nPrimarios es 26.7475516542989"
## [1] "Estimación para una mujer edc:\nProfesional es 30.3187628634423"
## [1] "Estimación para un hombre edc:\nProfesional es 33.4147589065171"
## [1] "Estimación para una mujer edc:\nTécnico Profesional es 29.549034662505"
## [1] "Estimación para un hombre edc:\nTécnico Profesional es 32.6450307055799"
## [1] "Estimación para una mujer edc:\nPrebásicos es 19.2400944619213"
## [1] "Estimación para un hombre edc:\nPrebásicos es 22.3360905049961"Para cada observación
f2 <- formula(minviajeprom ~ sexo + educ)
reg03 <- lm(formula = f2, data = base)
base[,minviajeprom_predict := predict(reg03)]
sum02<-summary(reg03)
sum02##
## Call:
## lm(formula = f2, data = base)
##
## Residuals:
## Min 1Q Median 3Q Max
## -32.41 -13.56 -3.56 8.88 620.85
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 31.1237 2.0463 15.209 < 2e-16 ***
## sexoMujer -3.0960 0.2670 -11.595 < 2e-16 ***
## educNo tiene -8.9966 2.1744 -4.138 3.52e-05 ***
## educPrebásicos -8.7876 2.1798 -4.031 5.57e-05 ***
## educPrimarios -4.3761 2.0665 -2.118 0.0342 *
## educProfesional 2.2911 2.0622 1.111 0.2666
## educSecundarios 0.5281 2.0559 0.257 0.7973
## educTécnico Profesional 1.5214 2.0805 0.731 0.4646
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 19.28 on 20927 degrees of freedom
## Multiple R-squared: 0.03249, Adjusted R-squared: 0.03216
## F-statistic: 100.4 on 7 and 20927 DF, p-value: < 2.2e-16P4.2
Dado el modelo anterior, calcule la predicción de los minutos de viaje promedio para una mujer y para un hombre con educación Profesional.
mujer1 <- predict(reg01, data.table(sexo= "Mujer", educ="Profesional"))
hombre1 <- predict(reg01, data.table(sexo="Hombre", educ="Profesional"))
data.table(mujer1, hombre1)## mujer1 hombre1
## 1: 30.31876 33.41476P5
Teniendo ya el modelo propuesto por su colega, contrarreste usted con un modelo de predicción en base a los datos con los que ya cuenta.
f3 <- formula(minviajeprom ~ sexo+educ+edad+IngresoFinal+comunahg+distprom+Macrozonahg)P5.1
Calcule los minutos de viaje promedios predichos para cada observación.
reg05 <- lm(formula = f3, data = base)
base[,minviajeprom_predict_2 := predict(reg05)]
sum03<-summary(reg05)
sum03##
## Call:
## lm(formula = f3, data = base)
##
## Residuals:
## Min 1Q Median 3Q Max
## -67.44 -8.89 -2.09 6.45 633.53
##
## Coefficients: (4 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.531e+01 1.722e+00 8.888 < 2e-16 ***
## sexoMujer -6.593e-01 2.180e-01 -3.025 0.002492 **
## educNo tiene -7.273e+00 1.727e+00 -4.211 2.56e-05 ***
## educPrebásicos -6.262e+00 1.733e+00 -3.614 0.000302 ***
## educPrimarios -2.765e+00 1.640e+00 -1.687 0.091666 .
## educProfesional -2.356e+00 1.642e+00 -1.435 0.151233
## educSecundarios -1.264e+00 1.634e+00 -0.774 0.439134
## educTécnico Profesional -2.168e+00 1.654e+00 -1.311 0.190022
## edad -4.374e-02 5.557e-03 -7.872 3.67e-15 ***
## IngresoFinal -8.140e-07 3.341e-07 -2.437 0.014833 *
## comunahgQuilpue 6.223e+00 7.972e-01 7.807 6.15e-15 ***
## comunahgValparaiso 8.849e+00 7.819e-01 11.317 < 2e-16 ***
## comunahgVilla Alemana 4.951e+00 6.474e-01 7.647 2.15e-14 ***
## comunahgVina del Mar 1.117e+01 8.198e-01 13.625 < 2e-16 ***
## distprom 1.943e-03 1.828e-05 106.296 < 2e-16 ***
## MacrozonahgConcon Oriente 1.928e+00 2.943e+00 0.655 0.512372
## MacrozonahgConcon Poniente NA NA NA NA
## MacrozonahgEl Belloto -3.344e-01 7.312e-01 -0.457 0.647408
## MacrozonahgEl Belloto Norte -1.989e+00 9.510e-01 -2.091 0.036546 *
## MacrozonahgForestal 8.825e-01 8.092e-01 1.091 0.275484
## MacrozonahgMarga-Marga -5.736e+00 1.049e+00 -5.467 4.63e-08 ***
## MacrozonahgMiraflores -2.512e+00 8.588e-01 -2.925 0.003447 **
## MacrozonahgPenablanca -4.876e+00 9.298e-01 -5.244 1.59e-07 ***
## MacrozonahgPlacilla-Curauma -8.969e+00 7.781e-01 -11.527 < 2e-16 ***
## MacrozonahgPlan Valparaiso -5.413e+00 1.169e+00 -4.630 3.67e-06 ***
## MacrozonahgPlan Vina -3.720e+00 7.911e-01 -4.702 2.59e-06 ***
## MacrozonahgPlaya Ancha -3.874e+00 7.090e-01 -5.465 4.69e-08 ***
## MacrozonahgQuilpue Norte -1.108e+00 9.602e-01 -1.154 0.248479
## MacrozonahgQuilpue Poniente -1.014e+00 7.226e-01 -1.404 0.160414
## MacrozonahgQuilpue Sur NA NA NA NA
## MacrozonahgRecreo -3.531e+00 7.544e-01 -4.680 2.88e-06 ***
## MacrozonahgRenaca -8.043e+00 9.480e-01 -8.484 < 2e-16 ***
## MacrozonahgRodelillo -1.392e+00 7.503e-01 -1.855 0.063664 .
## MacrozonahgSanta Julia -3.296e+00 6.935e-01 -4.752 2.03e-06 ***
## MacrozonahgValparaiso Alto -1.350e+00 6.652e-01 -2.030 0.042384 *
## MacrozonahgVilla Alemana Norte 1.824e+00 6.519e-01 2.798 0.005153 **
## MacrozonahgVilla Alemana Poniente NA NA NA NA
## MacrozonahgVina del Mar Oriente NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 15.28 on 20901 degrees of freedom
## Multiple R-squared: 0.3929, Adjusted R-squared: 0.3919
## F-statistic: 409.9 on 33 and 20901 DF, p-value: < 2.2e-16P6
Contando ya con su modelo, usted puede realizar una comparación directa entre ambos. Calcule los errores de predicción dentro de muestra para ambos modelos.
R^{2}
MAE1 <- sum(abs(base$minviajeprom_predict-base$minviajeprom))/nrow(base)
MAE2 <- sum(abs(base$minviajeprom_predict_2-base$minviajeprom))/nrow(base)
data.table(MAE1, MAE2)## MAE1 MAE2
## 1: 14.32717 10.49857R^{2} ajustado
RMSE1 <- (sum((base$minviajeprom_predict-base$minviajeprom)^2)/nrow(base))^(1/2)
RMSE2 <- (sum((base$minviajeprom_predict_2-base$minviajeprom)^2)/nrow(base))^(1/2)
data.table(RMSE1,RMSE2)## RMSE1 RMSE2
## 1: 19.27886 15.27172Comprobar criterios
data.table(regresion=c('Modelo 01','Modelo 02'),
`R2`=c(sum01$r.squared,sum03$r.squared),
`R2 ajustado`=c(sum01$adj.r.squared,sum03$adj.r.squared),
`Mean Absolute Error` = c(MAE(reg01$fitted.values,base$minviajeprom),
MAE(reg05$fitted.values,base$minviajeprom)),
`Root Mean Squared Error`=c(RMSE(reg01$fitted.values,base$minviajeprom),
RMSE(reg05$fitted.values,base$minviajeprom)))## regresion R2 R2 ajustado Mean Absolute Error Root Mean Squared Error
## 1: Modelo 01 0.03248765 0.03216402 14.32717 19.27886
## 2: Modelo 02 0.39288666 0.39192811 10.49857 15.27172P7
Interprete la diferencia en los errores de predicción entre el Modelo 1 y el Modelo 2. ¿Qué modelo hace una mejor predicción dentro de muestra?
En el modelo 01 el R cuadrado y el R cuadrado ajustado son cerca 0.03. En el modelo 02 las cifras correspondientes son cerca de 0.39. Lo dice que el modelo 02 as mucho más preciso. En el Modelo 01 el MAE es 14.33, en el modelo 02 es 10.5. por lo cual, el Modelo 02 es más presciso. Finalmente, considerando la presicion del modelo 2, podemos decir que presenta una mayor explicacion de los datos que el modelo 1, puesto que estimara con presicion cerca de un 40% de los datos.
P8
Realice validación cruzada (CV) a los modelos de la pregunta anterior por el método K-folds con 5 folds. ¿Se mantienen las conclusiones obtenidas en el análisis dentro de muestra?
K-folds con porcentaje
set.seed(12345)
training.samples <- createDataPartition(base$minviajeprom,p = 0.8, list = FALSE)
train.data <- base[training.samples, ]
test.data <- base[-training.samples, ]
model <- lm(f2, data = train.data)
predictions01 <- predict(model,train.data)
predictions02 <- predict(model,test.data)
data.table(' ' = c('Dentro de muestra','Fuera de muestra'),
R2 = c(R2(predictions01, train.data$minviajeprom),
R2(predictions02, test.data$minviajeprom)),
RMSE = c(RMSE(predictions01, train.data$minviajeprom),
RMSE(predictions02, test.data$minviajeprom)),
MAE = c(MAE(predictions01, train.data$minviajeprom),
MAE(predictions02, test.data$minviajeprom)))## R2 RMSE MAE
## 1: Dentro de muestra 0.03331561 19.21935 14.27693
## 2: Fuera de muestra 0.02898701 19.51771 14.51030Regreción Lineal, por K-folds
set.seed(12345)
setupRKCV <- trainControl(method = "repeatedcv" , number = 5, repeats= 3)
predRKCV<-train(f2,data=base,method="lm",trControl= setupRKCV)
print(predRKCV)## Linear Regression
##
## 20935 samples
## 2 predictor
##
## No pre-processing
## Resampling: Cross-Validated (5 fold, repeated 3 times)
## Summary of sample sizes: 16748, 16748, 16748, 16748, 16748, 16748, ...
## Resampling results:
##
## RMSE Rsquared MAE
## 19.2676 0.03215984 14.33307
##
## Tuning parameter 'intercept' was held constant at a value of TRUEset.seed(12345)
setupRKCV <- trainControl(method = "repeatedcv" , number = 5, repeats= 3)
predRKCV<-train(f3,data=base,method="lm",trControl= setupRKCV)
print(predRKCV)## Linear Regression
##
## 20935 samples
## 7 predictor
##
## No pre-processing
## Resampling: Cross-Validated (5 fold, repeated 3 times)
## Summary of sample sizes: 16748, 16748, 16748, 16748, 16748, 16748, ...
## Resampling results:
##
## RMSE Rsquared MAE
## 15.25732 0.3941653 10.51855
##
## Tuning parameter 'intercept' was held constant at a value of TRUERespuesta
Por lo tanto, considerando que el R cuadrado osease la explicación de nuesto modelo se encuentra en un 39% considerando la regreción de f3, podriamos considerar que es un mejor modelo que el planteado solo con 2 varaiables, por lo que no cambiaria la conclusión anterior.