Cargamos paquetes a utilizar

library(data.table)
library(ggplot2)
library(caret)
library(jtools)
library(scales)
library(corrplot)
library(plotly)
library(tidyverse)
library(lubridate)
library(texreg)
library(corrgram)
library(knitr)
library(moderndive)

El formato de respuesta es el siguiente:

Abajo del encabezado, debe estar escrita la pregunta. Agregar los códigos en chunks e incluir el output de ser necesario.

La explicación de las respuestas (si es que hay análisis) debe ir como texto, fuera del chunk.

El título del informe debe ser “Tarea 5” y el nombre de todos los integrantes debe ir en autor. La actividad deberá ser entregado en formato HTML. Agregue el número del grupo.

También debe enviar su archivo en formato R Markdown como respaldo.

Recuerde que puede usar todo el material visto para apoyarse.

Todos los integrantes del grupo deben enviar el archivo de Rmarkdown

Tienen hasta el final de la clase para enviar la actividad

El formato tiene puntaje.(3 ptos)

Suprimir warning y mensajes(3 ptos)

Cargue las librerias y base de datos correspondientes.

cargamos base de datos

path <- 'C:/Users/matir/OneDrive/Escritorio/Mati/Universidad/2021/Segundo Semestre/Data Science/Tarea 5/'
base <- fread(paste0(path,'CarPrice_Assignment.csv'))

Con respecto a los predichos:

Tuvimos una breve reunion con la profe y nos dio permiso de hacer un grafico para realizar los valores predichos. De todas maneras, sabemos cómo hacerlo (creando una nueva variable en la base con los valores predichos)

P1

Muestre un histograma de distribución del precio.(4 ptos)

base <- data.frame(base)

ggplot(base, aes(x=price)) + 
  geom_histogram(bindwidth=0.2, color="black", fill="white") + scale_x_continuous(n.breaks = 10) + scale_y_continuous(n.breaks = 20)+ labs(x="Precio", y="Numero", title= "Distribucion del precio")

P2

Realice un gráfico de correlación entre las variables númericas, estas son, “wheelbase”, “carlength” “carwidth”,“carheight”, “curbweight”, “enginesize”,“boreratio”,“stroke”, “compressionratio”, “horsepower”, “peakrpm”,“citympg”, “highwaympg”,“price” (4 ptos)

correlacion <- cor(base[,c("wheelbase", "carlength", "carwidth","carheight", "curbweight", "enginesize","boreratio","stroke", "compressionratio", "horsepower", "peakrpm","citympg", "highwaympg","price")])
corrplot(correlacion, method = 'color')

P3

Realice dos gráficos de puntos, uno que relacione el precio price, el rendimiento en ciudad citympg y que se diferencie por tipo de tracción drivewheel , el segundo entre precio price, el rendimiento en carretera highwaympgy que se diferencie por tipo de tracción drivewheel (8 ptos)

ggplot(base, aes(x=price, y=citympg,color=drivewheel)) + 
  geom_point() + facet_wrap(facets="drivewheel", scales = "free_y", "free_x") + 
  scale_y_continuous(n.breaks = 20)+ geom_smooth(colour="black", size=0.7)+ 
  scale_color_brewer(palette="Set1") + 
  labs(x="Precio", y="Rendimiento en Ciudad en KM", title= "Precio y Rendmiento en Ciudad", subtitle= "por tipo de Traccion", fill= "Tipo de Causal")

ggplot(base, aes(x=price, y=highwaympg,color=drivewheel)) + 
  geom_point() + facet_wrap(facets="drivewheel", scales = "free_y", "free_x") + 
  scale_y_continuous(n.breaks = 20)+ geom_smooth(colour="black", size=0.7)+ 
  scale_color_brewer(palette="Set1") + 
  labs(x="Precio", y="Rendimiento en Carretera en KM", title= "Precio y Rendmiento en Carretera", subtitle= "por tipo de Traccion", fill= "Tipo de Causal")

Con linea de tendencia

ggplot(base, aes(x=price, y=citympg,color=drivewheel)) + 
  geom_point() + facet_wrap(facets="drivewheel", scales = "free_y", "free_x") + 
  scale_y_continuous(n.breaks = 20)+ geom_smooth(aes(colour= drivewheel, fill = drivewheel), method = "lm", fullrange=TRUE) + scale_color_brewer(palette="Set1") + labs(x="Precio", y="Rendimiento en Ciudad en KM", title= "Precio y Rendmiento en Ciudad", subtitle= "por tipo de Traccion")

ggplot(base, aes(x=price, y=highwaympg,color=drivewheel)) + 
  geom_point() + facet_wrap(facets="drivewheel", scales = "free_y", "free_x") + 
  scale_y_continuous(n.breaks = 20)+ geom_smooth(aes(colour= drivewheel, fill = drivewheel), method = "lm", fullrange=TRUE) + scale_color_brewer(palette="Set1") + labs(x="Precio", y="Rendimiento en Carretera en KM", title= "Precio y Rendmiento en Carretera", subtitle= "por tipo de Traccion")

Hicimos un plotly, por si se quiere saber exactamenteel precio y marca-modelo de cada auto.

grafico1 <- ggplot(base, aes(x=price, y=citympg,color=drivewheel, text=paste("Marca y modelo del auto:",CarName,"\n","Precio:",price))) + 
  geom_point() + facet_wrap(facets="drivewheel") + 
  scale_y_continuous(n.breaks = 20)+ 
  scale_color_brewer(palette="Set1") + 
  labs(x="Precio", y="Rendimiento en Ciudad en KM", title= "Precio y Rendmiento en Ciudad", subtitle= "por tipo de Traccion", fill= "Tipo de Causal")
ggplotly(grafico1, tooltip = "text")
grafico2 <- ggplot(base, aes(x=price, y=highwaympg,color=drivewheel, text=paste("Marca y modelo del auto:",CarName,"\n","Precio:",price))) + 
  geom_point() + facet_wrap(facets="drivewheel") + 
  scale_y_continuous(n.breaks = 20)+ 
  scale_color_brewer(palette="Set1") + 
  labs(x="Precio", y="Rendimiento en Carretera en KM", title= "Precio y Rendmiento en Ciudad", subtitle= "por tipo de Traccion", fill= "Tipo de Causal")
ggplotly(grafico2, tooltip = "text")

P4

Realice una regresión múltiple entre precio con las variables enginetype y fuelsystem. ¿Qué podrian decir de la significacia de las variables?(5 ptos)

Calcule el predicho para cada observación.(2 ptos)

base <- data.table(base)

f <- formula(price ~ enginetype  + fuelsystem)
reg1 <- lm(formula = f, data = base)
sum1 <- summary(reg1)


prediccion.todas1 <-  predict(reg1, newdata = base)

base[,predicciont1:= prediccion.todas1]


gg2 <- ggplot(base, aes(price, predicciont1)) +
  geom_point()  + geom_abline(colour="red") +
  labs(title = "Regresión lineal múltiple del precio", subtitle = "Variables: enginetype y fuelsystem", x = "Precio Actual",
       y = "Predicción del precio") +
  theme(plot.title = element_text(family = "Helvetica", face = "bold", size = (15)), 
        axis.title = element_text(family = "Helvetica", size = (10)))
gg2

Tabla ordenada:

regresion <-  lm(formula = price ~ enginetype + fuelsystem,
               data = base)

sum <- summary(regresion)
sum
## 
## Call:
## lm(formula = price ~ enginetype + fuelsystem, data = base)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11599.4  -2871.7   -399.4   1660.6  25431.8 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       9788.7     2682.8   3.649 0.000340 ***
## enginetypedohcv  13284.1     6391.2   2.078 0.039002 *  
## enginetypel      -2804.7     2602.7  -1.078 0.282566    
## enginetypeohc    -2233.2     1941.6  -1.150 0.251503    
## enginetypeohcf     261.6     2446.3   0.107 0.914966    
## enginetypeohcv    6982.0     2458.2   2.840 0.004996 ** 
## enginetyperotor  -2471.4     6391.2  -0.387 0.699418    
## fuelsystem2bbl    -371.1     2010.6  -0.185 0.853751    
## fuelsystem4bbl    4827.7     7370.9   0.655 0.513276    
## fuelsystemidi     8425.5     2354.8   3.578 0.000439 ***
## fuelsystemmfi     5408.5     6413.6   0.843 0.400124    
## fuelsystemmpfi    8327.7     2013.8   4.135 5.31e-05 ***
## fuelsystemspdi    3434.9     2760.0   1.245 0.214823    
## fuelsystemspfi    3492.5     6413.6   0.545 0.586703    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6141 on 191 degrees of freedom
## Multiple R-squared:  0.4468, Adjusted R-squared:  0.4092 
## F-statistic: 11.87 on 13 and 191 DF,  p-value: < 2.2e-16
get_regression_table(regresion) %>% kable()
term estimate std_error statistic p_value lower_ci upper_ci
intercept 9788.723 2682.816 3.649 0.000 4496.971 15080.475
enginetype: dohcv 13284.083 6391.247 2.078 0.039 677.591 25890.576
enginetype: l -2804.675 2602.677 -1.078 0.283 -7938.358 2329.007
enginetype: ohc -2233.177 1941.569 -1.150 0.252 -6062.849 1596.495
enginetype: ohcf 261.554 2446.301 0.107 0.915 -4563.682 5086.790
enginetype: ohcv 6981.968 2458.172 2.840 0.005 2133.317 11830.619
enginetype: rotor -2471.417 6391.247 -0.387 0.699 -15077.909 10135.076
fuelsystem: 2bbl -371.127 2010.611 -0.185 0.854 -4336.981 3594.728
fuelsystem: 4bbl 4827.694 7370.879 0.655 0.513 -9711.085 19366.472
fuelsystem: idi 8425.479 2354.816 3.578 0.000 3780.695 13070.263
fuelsystem: mfi 5408.455 6413.555 0.843 0.400 -7242.040 18058.949
fuelsystem: mpfi 8327.694 2013.788 4.135 0.000 4355.573 12299.814
fuelsystem: spdi 3434.899 2759.955 1.245 0.215 -2009.007 8878.805
fuelsystem: spfi 3492.455 6413.555 0.545 0.587 -9158.040 16142.949 
screenreg(regresion)
## 
## =============================
##                  Model 1     
## -----------------------------
## (Intercept)       9788.72 ***
##                  (2682.82)   
## enginetypedohcv  13284.08 *  
##                  (6391.25)   
## enginetypel      -2804.68    
##                  (2602.68)   
## enginetypeohc    -2233.18    
##                  (1941.57)   
## enginetypeohcf     261.55    
##                  (2446.30)   
## enginetypeohcv    6981.97 ** 
##                  (2458.17)   
## enginetyperotor  -2471.42    
##                  (6391.25)   
## fuelsystem2bbl    -371.13    
##                  (2010.61)   
## fuelsystem4bbl    4827.69    
##                  (7370.88)   
## fuelsystemidi     8425.48 ***
##                  (2354.82)   
## fuelsystemmfi     5408.45    
##                  (6413.56)   
## fuelsystemmpfi    8327.69 ***
##                  (2013.79)   
## fuelsystemspdi    3434.90    
##                  (2759.95)   
## fuelsystemspfi    3492.45    
##                  (6413.56)   
## -----------------------------
## R^2                  0.45    
## Adj. R^2             0.41    
## Num. obs.          205       
## =============================
## *** p < 0.001; ** p < 0.01; * p < 0.05

En cuanto a la significancia, podemos ver que principalmente si los estadisticos serán significativos si su valor p´se encuentra debajo de los niveles de significancia detallados en la parte “Signif. codes: 0 ‘’ 0.001 ‘’ 0.01 ‘’ 0.05 ‘.’ 0.1 ‘ ’ 1”, entonces si su valor p estar bajo: 0,1 - 0,05 - 0,01 - 0,001 , el estadistico será significativo.

Ahora, en este caso: podemos decir que en cuanto a el tipo de motor –> “dohcv”, “pel”y “ohcv” son significativos. En cuanto al sistema de combustible del coche: “idi” y “mpfi” son significativos.

P5

Dado el modelo anterior, calcule la predicción del precio con un tipo de motor rotor y sistema de gasolina igual a idi (5 puntos)

prediccion <- predict(reg1, data.table(enginetype = "rotor", fuelsystem = "idi"))
prediccion
##        1 
## 15742.79

El precio será de 15.743 dolares.

P6

Realice usted un modelo a su preferencia, este debe predecir de mejor manera de mejor manera(*) (7 ptos)

Calcule el predicho para cada observación.(2 ptos)

f2 <- formula(price ~ enginesize + stroke + cylindernumber + horsepower + enginelocation + compressionratio + peakrpm + enginetype + curbweight + carbody + carwidth + carlength)
reg2 <- lm(formula = f2, data = base)
summary(reg2)
## 
## Call:
## lm(formula = f2, data = base)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6021.1 -1113.2   -80.8  1028.8 10074.5 
## 
## Coefficients: (1 not defined because of singularities)
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          -2.673e+04  1.162e+04  -2.300 0.022568 *  
## enginesize            9.616e+01  1.770e+01   5.433 1.78e-07 ***
## stroke               -4.674e+03  7.697e+02  -6.072 7.30e-09 ***
## cylindernumberfive   -1.071e+04  2.294e+03  -4.670 5.88e-06 ***
## cylindernumberfour   -1.243e+04  2.345e+03  -5.300 3.36e-07 ***
## cylindernumbersix    -8.564e+03  1.843e+03  -4.648 6.47e-06 ***
## cylindernumberthree  -4.548e+03  3.697e+03  -1.230 0.220249    
## cylindernumbertwelve -1.563e+04  3.020e+03  -5.174 6.06e-07 ***
## cylindernumbertwo    -4.562e+03  3.072e+03  -1.485 0.139259    
## horsepower            3.486e+01  1.345e+01   2.591 0.010345 *  
## enginelocationrear    6.660e+03  2.209e+03   3.015 0.002942 ** 
## compressionratio      1.878e+02  5.664e+01   3.316 0.001103 ** 
## peakrpm               1.580e+00  5.033e-01   3.139 0.001979 ** 
## enginetypedohcv      -1.343e+04  3.520e+03  -3.816 0.000186 ***
## enginetypel           1.029e+03  1.269e+03   0.810 0.418752    
## enginetypeohc         3.539e+03  8.464e+02   4.182 4.51e-05 ***
## enginetypeohcf        2.606e+02  1.240e+03   0.210 0.833713    
## enginetypeohcv       -6.168e+03  1.127e+03  -5.472 1.47e-07 ***
## enginetyperotor              NA         NA      NA       NA    
## curbweight            4.436e+00  1.361e+00   3.259 0.001338 ** 
## carbodyhardtop       -2.766e+03  1.290e+03  -2.145 0.033333 *  
## carbodyhatchback     -2.699e+03  1.088e+03  -2.479 0.014080 *  
## carbodysedan         -1.459e+03  1.135e+03  -1.285 0.200398    
## carbodywagon         -2.480e+03  1.232e+03  -2.013 0.045561 *  
## carwidth              6.710e+02  2.129e+02   3.152 0.001899 ** 
## carlength            -8.632e+01  3.983e+01  -2.167 0.031536 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2208 on 180 degrees of freedom
## Multiple R-squared:  0.9326, Adjusted R-squared:  0.9236 
## F-statistic: 103.8 on 24 and 180 DF,  p-value: < 2.2e-16
sum2 <- summary(reg2)


prediccion.todas <-  predict(reg2, newdata = base)

base[,predicciont:=prediccion.todas]


gg2 <- ggplot(base, aes(price, predicciont)) +
  geom_point()  + geom_abline(colour="red") +
  labs(title = "Regresión lineal múltiple del precio", subtitle = "Variables mas significativas", x = "Precio Actual",
       y = "Predicción del precio") +
  theme(plot.title = element_text(family = "Helvetica", face = "bold", size = (15)), 
        axis.title = element_text(family = "Helvetica", size = (10)))
gg2

# si añadieramos la variable "CarName" , no tendría sentido dado que ya el hecho de introducir la marca y modelo del auto se puede obtener el precio exacto. Es como si se quisiera predecir algo en una persona "x" y usaramos como variable el nombre de la persona.

(*)al parecer el enunciado estaba mal redactado, por lo tanto supondremos que se quiere preguntar por el modelo que mejor predice el precio.

En este modelo se incluyeron unicamente variables significativas, por lo tanto segun el R cuadrado ajustado que se logró, que la prediccion se ecerca mucho a la realidad. El modelo se ajusta un 92% app a sus datos.

P7

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? (5 puntos)

base <- data.table(base)

base[,price_predict1 := predict(reg1)]
base[,price_predict2 := predict(reg2)]



MAE1 <- sum(abs(base$price_predict1-base$price))/nrow(base)
MAE2 <- sum(abs(base$price_predict2-base$price))/nrow(base)

data.table(MAE1,MAE2)
##       MAE1    MAE2
## 1: 3936.92 1499.11
RMSE1 <- (sum((base$price_predict1-base$price)^2)/nrow(base))^(1/2)
RMSE2 <- (sum((base$price_predict2-base$price)^2)/nrow(base))^(1/2)


data.table(regresion=c('Modelo 01','Modelo 02'),
           `R2`=c(sum1$r.squared,sum2$r.squared),
           `R2 ajustado`=c(sum1$adj.r.squared,sum2$adj.r.squared),
           `Mean Absolute Error` = c(MAE(reg1$fitted.values,base$price),
                                     MAE(reg2$fitted.values,base$price)),
           `Root Mean Squared Error`=c(RMSE(reg1$fitted.values,base$price),
                                      RMSE(reg2$fitted.values,base$price)))
##    regresion        R2 R2 ajustado Mean Absolute Error Root Mean Squared Error
## 1: Modelo 01 0.4468492   0.4092002             3936.92                5927.128
## 2: Modelo 02 0.9326171   0.9236327             1499.11                2068.699

Claramente, el modelo 02 (el que hicimos nosotros) puede explicar de mejor manera el modelo con creces. Esto se debe principalmente por su R2 ajustado, donde es mas del doble que el modelo 01. Esta diferencia se debe a que incorporamos mas variables pero que tambien son significativas.

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? (9 puntos)

set.seed(12345)
setupKCV <- trainControl(method = "cv" , number = 5)

predkfolds1<-train(f,data=base,method="lm",trControl= setupKCV)
print(predkfolds1)
## Linear Regression 
## 
## 205 samples
##   2 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 164, 165, 164, 164, 163 
## Resampling results:
## 
##   RMSE      Rsquared   MAE     
##   6355.632  0.4070385  4375.004
## 
## Tuning parameter 'intercept' was held constant at a value of TRUE
predkfolds2<-train(f2,data=base,method="lm",trControl= setupKCV)
print(predkfolds2)
## Linear Regression 
## 
## 205 samples
##  12 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 162, 165, 164, 164, 165 
## Resampling results:
## 
##   RMSE      Rsquared   MAE     
##   3086.703  0.8776497  2037.523
## 
## Tuning parameter 'intercept' was held constant at a value of TRUE

Si, a pesar de usar pocos folds, las conclusiones se mantienen (quizas nuestro modelo no tiene un R2 tan alto pero supera nuevamente al del primer modelo con menos variables).