A través de un conjunto de datos (TargetTrain.csv) se debe entrenar un modelo que realice la predicción de la variable TARGET teniendo en cuenta los inputs dados por las variables (VAR1, VAR2, VAR3, VAR4).
Se debe realizar el test del modelo utilizando los datos del fichero TargetTest.csv. Para ello, elegiremos la métrica que mejor se adapte a la naturaleza de los datos e indagaremos sobre diferentes naturalezas de modelos predictivos, mostrando cuál de ellos es el que mejor resultado arroja para el caso bajo estudio tras realizar la comparación entre las métricas obtenidas.
Procedemos a cargar algunas de las librerías trabajadas durante el curso y que creemos podrían ser útiles para este análisis:
library(readr)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(data.table)
##
## Attaching package: 'data.table'
## The following objects are masked from 'package:dplyr':
##
## between, first, last
library(ggplot2)
library(gplots)
##
## Attaching package: 'gplots'
## The following object is masked from 'package:stats':
##
## lowess
library(corrplot)
## corrplot 0.92 loaded
library(rpart)
library(randomForest)
## randomForest 4.7-1.1
## Type rfNews() to see new features/changes/bug fixes.
##
## Attaching package: 'randomForest'
## The following object is masked from 'package:ggplot2':
##
## margin
## The following object is masked from 'package:dplyr':
##
## combine
library(caret)
## Loading required package: lattice
library(mltools)
library(corrplot)
library(pROC)
## Type 'citation("pROC")' for a citation.
##
## Attaching package: 'pROC'
## The following objects are masked from 'package:stats':
##
## cov, smooth, var
library(cluster)
library(factoextra)
## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa
library(devtools)
## Loading required package: usethis
library(naivebayes)
## naivebayes 0.9.7 loaded
##
## Attaching package: 'naivebayes'
## The following object is masked from 'package:data.table':
##
## tables
library(nnet)
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ lubridate 1.9.3 ✔ tibble 3.2.1
## ✔ purrr 1.0.2 ✔ tidyr 1.3.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ data.table::between() masks dplyr::between()
## ✖ randomForest::combine() masks dplyr::combine()
## ✖ dplyr::filter() masks stats::filter()
## ✖ data.table::first() masks dplyr::first()
## ✖ lubridate::hour() masks data.table::hour()
## ✖ lubridate::isoweek() masks data.table::isoweek()
## ✖ dplyr::lag() masks stats::lag()
## ✖ data.table::last() masks dplyr::last()
## ✖ purrr::lift() masks caret::lift()
## ✖ randomForest::margin() masks ggplot2::margin()
## ✖ lubridate::mday() masks data.table::mday()
## ✖ lubridate::minute() masks data.table::minute()
## ✖ lubridate::month() masks data.table::month()
## ✖ lubridate::quarter() masks data.table::quarter()
## ✖ tidyr::replace_na() masks mltools::replace_na()
## ✖ lubridate::second() masks data.table::second()
## ✖ purrr::transpose() masks data.table::transpose()
## ✖ lubridate::wday() masks data.table::wday()
## ✖ lubridate::week() masks data.table::week()
## ✖ lubridate::yday() masks data.table::yday()
## ✖ lubridate::year() masks data.table::year()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(moments)
##
## Attaching package: 'moments'
##
## The following object is masked from 'package:mltools':
##
## skewness
library(vcd)
## Loading required package: grid
library(car)
## Loading required package: carData
##
## Attaching package: 'car'
##
## The following object is masked from 'package:purrr':
##
## some
##
## The following object is masked from 'package:dplyr':
##
## recode
library(rpart)
library(rpart.plot)
library(e1071)
##
## Attaching package: 'e1071'
##
## The following objects are masked from 'package:moments':
##
## kurtosis, moment, skewness
##
## The following object is masked from 'package:mltools':
##
## skewness
DataTrain=read.csv("TargetTrain.csv", header=TRUE, sep=";")
DataTest=read.csv("TargetTest.csv", header=TRUE, sep=";")
str(DataTrain)
## 'data.frame': 3150 obs. of 5 variables:
## $ TARGET: num 0.128 0.138 0.167 0.277 0.376 ...
## $ VAR1 : int 1 2 3 4 5 6 7 8 9 10 ...
## $ VAR2 : num 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 ...
## $ VAR3 : int 1 1 1 1 1 1 1 1 1 1 ...
## $ VAR4 : int 300 300 300 300 300 300 300 300 300 300 ...
str(DataTest)
## 'data.frame': 2250 obs. of 5 variables:
## $ TARGET: num 0.153 0.164 0.183 0.291 0.402 ...
## $ VAR1 : int 1 2 3 4 5 6 7 8 9 10 ...
## $ VAR2 : num 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 ...
## $ VAR3 : int 1 1 1 1 1 1 1 1 1 1 ...
## $ VAR4 : int 300 300 300 300 300 300 300 300 300 300 ...
summary(DataTrain)
## TARGET VAR1 VAR2 VAR3 VAR4
## Min. :0.1277 Min. : 1.0 Min. : 0.500 Min. :1.000 Min. :300
## 1st Qu.:2.2131 1st Qu.:13.0 1st Qu.: 0.750 1st Qu.:1.000 1st Qu.:300
## Median :3.5503 Median :25.5 Median : 2.500 Median :3.000 Median :400
## Mean :3.0593 Mean :25.5 Mean : 3.893 Mean :3.333 Mean :400
## 3rd Qu.:4.1773 3rd Qu.:38.0 3rd Qu.: 7.500 3rd Qu.:6.000 3rd Qu.:500
## Max. :4.6191 Max. :50.0 Max. :10.000 Max. :6.000 Max. :500
summary(DataTest)
## TARGET VAR1 VAR2 VAR3 VAR4
## Min. :0.1486 Min. : 1.0 Min. :0.80 Min. :1.000 Min. :300
## 1st Qu.:2.3024 1st Qu.:13.0 1st Qu.:2.00 1st Qu.:1.000 1st Qu.:300
## Median :3.7150 Median :25.5 Median :4.00 Median :3.000 Median :400
## Mean :3.1526 Mean :25.5 Mean :4.16 Mean :3.333 Mean :400
## 3rd Qu.:4.2959 3rd Qu.:38.0 3rd Qu.:6.00 3rd Qu.:6.000 3rd Qu.:500
## Max. :4.6055 Max. :50.0 Max. :8.00 Max. :6.000 Max. :500
Observamos la naturaleza de las VAR1-4, identificando que la variable de estudio Target es numérica continua, al igual que VAR2, mientras que el resto son de tipo entero. No obstante, transformamos el conjunto de datos a variables continuas, de cara a favorecer los pasos posteriores.
DataTrain=DataTrain %>%mutate_if(is.integer, as.numeric)
str(DataTrain)
## 'data.frame': 3150 obs. of 5 variables:
## $ TARGET: num 0.128 0.138 0.167 0.277 0.376 ...
## $ VAR1 : num 1 2 3 4 5 6 7 8 9 10 ...
## $ VAR2 : num 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 ...
## $ VAR3 : num 1 1 1 1 1 1 1 1 1 1 ...
## $ VAR4 : num 300 300 300 300 300 300 300 300 300 300 ...
Nuestro data.frame de entrenamiento contempla 3150 registros y el de validación 2250.Tras aplicar la función “summary” se puede ver que ambos conjuntos de datos no tienen valores NA, por lo que contamos con un dataset limpio para proceder a su análisis exploratorio.
Procedemos a valorar la posible relación entre la variable de estudio Target y el resto de variables que conforman el DataTrain. Comenzamos investigando la forma y simetría de las distribuciones de cada una de las variables en el conjunto de datos DataTrain:
skewness(DataTrain$VAR1)
## [1] 0
kurtosis(DataTrain$VAR1)
## [1] -1.202102
skewness(DataTrain$VAR2)
## [1] 0.6235926
kurtosis(DataTrain$VAR2)
## [1] -1.115935
skewness(DataTrain$VAR3)
## [1] 0.2389493
kurtosis(DataTrain$VAR3)
## [1] -1.500952
skewness(DataTrain$VAR4)
## [1] 0
kurtosis(DataTrain$VAR4)
## [1] -1.500952
Los resultados de asimetría y curtosis para cada una de las variables son:
skewness(DataTrain\(VAR1): Asimetría de VAR1 es 0 (perfectamente simétrica). kurtosis(DataTrain\)VAR1): Curtosis de VAR1 es 1.79904 (algo más puntiaguda que una distribución normal).
skewness(DataTrain\(VAR2): Asimetría de VAR2 es 0.6238896 (casi simétrica). kurtosis(DataTrain\)VAR2): Curtosis de VAR2 es 1.885261 (algo más puntiaguda que una distribución normal).
skewness(DataTrain\(VAR3): Asimetría de VAR3 es 0.2390631 (casi simétrica). kurtosis(DataTrain\)VAR3): Curtosis de VAR3 es 1.5 (algo más puntiaguda que una distribución normal).
skewness(DataTrain\(VAR4): Asimetría de VAR4 es 0 (perfectamente simétrica). kurtosis(DataTrain\)VAR4): Curtosis de VAR4 es 1.5 (algo más puntiaguda que una distribución normal).
Seguimos con la visualización de los datos, graficando la relación entra la variable Target y VAR1-4, empezando por calcular la matriz de correlación y posteriormente dos gráficos: barras y puntos:
correlation_matrix=cor(DataTrain)
correlation_target=correlation_matrix["TARGET", -1]
variables=names(correlation_target)
barplot(correlation_target, names.arg = variables,
main = "Correlaciones entre TARGET y otras variables",
xlab = "Variables", ylab = "Correlación")
plot(correlation_target, main = "Correlaciones entre TARGET y otras variables",
xlab = "Variables", ylab = "Correlación", pch = 19)
Hemos podido observar que la correlación entre las VAR1-4 y la variable
de estudio Target son muy distintas entre sí. Contemplamos que Target
presenta una fuerte correlación con VAR1, mientras que el valor de
correlación para Target vs. resto de variables cae bastante, lo que
indicaría una débil correlación entre estas variables.
No obstante, vamos a generar la matriz de correlación para validar esta reflexión:
M=cor(DataTest)
print(M)
## TARGET VAR1 VAR2 VAR3 VAR4
## TARGET 1.00000000 0.9128401 9.902817e-02 -1.780120e-02 -3.322691e-02
## VAR1 0.91284012 1.0000000 0.000000e+00 0.000000e+00 0.000000e+00
## VAR2 0.09902817 0.0000000 1.000000e+00 -5.125467e-19 0.000000e+00
## VAR3 -0.01780120 0.0000000 -5.125467e-19 1.000000e+00 1.411705e-20
## VAR4 -0.03322691 0.0000000 0.000000e+00 1.411705e-20 1.000000e+00
El valor 0.9128401 es la correlación entre TARGET y VAR1. El valor 0.09902817 es la correlación entre TARGET y VAR2. El valor -0.01780120 es la correlación entre TARGET y VAR3. El valor -0.03322691 es la correlación entre TARGET y VAR4.
Estos números indican la fuerza y la dirección de la relación lineal entre las variables. Los valores cercanos a 1 o -1 indican una relación más fuerte, mientras que los valores cercanos a 0 indican una relación débil o nula. En este caso, la alta correlación positiva (0.91) entre Target y VAR1 sugiere una fuerte relación lineal positiva entre esas dos variables, mientras que las correlaciones más cercanas a cero sugieren relaciones más débiles entre Target y las otras variables (VAR2, VAR3, VAR4), lo que indicaría independencia de la variables entre sí.
Comenzaremos entrenando con nuestro conjunto de datos DataTrain un modelo de regresión lineal:
modeloRL=lm(TARGET~VAR1+VAR2+VAR3+VAR4,data = DataTrain)
print(modeloRL)
##
## Call:
## lm(formula = TARGET ~ VAR1 + VAR2 + VAR3 + VAR4, data = DataTrain)
##
## Coefficients:
## (Intercept) VAR1 VAR2 VAR3 VAR4
## 0.9167535 0.0868552 0.0574772 -0.0196577 -0.0005762
summary(modeloRL)
##
## Call:
## lm(formula = TARGET ~ VAR1 + VAR2 + VAR3 + VAR4, data = DataTrain)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.31196 -0.49367 0.09087 0.47884 1.00021
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.9167535 0.0552213 16.601 < 2e-16 ***
## VAR1 0.0868552 0.0006762 128.455 < 2e-16 ***
## VAR2 0.0574772 0.0028297 20.312 < 2e-16 ***
## VAR3 -0.0196577 0.0047486 -4.140 3.57e-05 ***
## VAR4 -0.0005762 0.0001195 -4.822 1.49e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5476 on 3145 degrees of freedom
## Multiple R-squared: 0.8435, Adjusted R-squared: 0.8433
## F-statistic: 4238 on 4 and 3145 DF, p-value: < 2.2e-16
El modelo resultante tiene un R2 de 0.8453, que explicaría en un 84,35% la variabiliad de nuestra variavble de estudio: Target.
Para validar el modelo con el conjunto de entrenamiento, vamos a usar la función “predict” y haremos la consulta del R-squared para poder comparar con el modelo de regresión lineal proyectado durante el entrenamiento:
prediccionesRL1=predict(modeloRL,newdata=DataTest)
summary(prediccionesRL1)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.6435 1.9862 3.0757 3.0747 4.1576 5.5268
r2Test=1 - sum((DataTest$TARGET - prediccionesRL1)^2) / sum((DataTest$TARGET - mean(DataTest$TARGET))^2)
print(r2Test)
## [1] 0.8408305
Tras calcular el R2, podemos comprobar que los datos resultantes para ambos conjuntos (entrenamiento y validación sugieren que el modelo no ha sobreajustado (overfitting) los datos de entrenamiento y que tendría buena capacidad de generalización. Sin embargo, antes de dar por bueno el modelo, vamos a observar un gráfico de residuos que nos permitirá profundizar mejor en su utilidad:
plot(modeloRL)
Los gráficos anteriormente pintados podrían dar pistas de que nos encontramos ante un modelo predictivo que no es especialmente bueno para nuestro conjunto de datos. Podemos observarlo, por ejemplo, en el primer gráfico de residuos vs. fitted. Aquí, nos encontramops con una parábola, que tumbaría la hipótesis de que existe una relación lineal entre nuestras variables.
Seguimos, por tanto, probando otros modelos predictivos con el objetivo de encontrar aquél que mejor rendimiento nos proporcione.Entrenaremos a continuación un modelo de Support Vector Machines (SVM) y otro de Random Forest (RF):
modeloSVM=train(TARGET~., data=DataTrain, method="svmLinear", metric="RMSE")
modeloSVM
## Support Vector Machines with Linear Kernel
##
## 3150 samples
## 4 predictor
##
## No pre-processing
## Resampling: Bootstrapped (25 reps)
## Summary of sample sizes: 3150, 3150, 3150, 3150, 3150, 3150, ...
## Resampling results:
##
## RMSE Rsquared MAE
## 0.5695276 0.8402771 0.4684363
##
## Tuning parameter 'C' was held constant at a value of 1
summary(modeloSVM)
## Length Class Mode
## 1 ksvm S4
El resultado del R2 para el modelo SVM mínimamente más bajo que el que nos encontramos en el estudio del modelo de LR, por lo que probaremos a entrenar a continuación un modelo RF:
modeloRF=randomForest(x=DataTrain, y=DataTrain$TARGET, importance=TRUE)
print(modeloRF)
##
## Call:
## randomForest(x = DataTrain, y = DataTrain$TARGET, importance = TRUE)
## Type of random forest: regression
## Number of trees: 500
## No. of variables tried at each split: 1
##
## Mean of squared residuals: 0.07042453
## % Var explained: 96.32
Tras entrenar nuestro modelo RF observamos que el % de explicación de la variabilidad de Target ascendería a un 96,8%, lo cual es un resultado mucho más positivo que en los modelos desarrollados previamente. Procedmos ahora a validar con nuestro conjunto de datos DataTest este modelo que cuenta con hasta 500 árboles:
prediccionesRF=predict(modeloRF,DataTest)
summary(prediccionesRF)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.5695 2.1466 3.6515 3.1374 4.0717 4.2939
r2Test3= R2(prediccionesRF, DataTest$TARGET)
print(r2Test3)
## [1] 0.989735
Tras validar modelo RF con el conjunto de DataTest, nos hemos encontrado con un R2 que asciende a 0,99, lo que podría llevarnos al escenario de overfitting. Para ver si nuestro modelo ha sobreajustado los datos, lo que lo invalidaría inmediatamente.
plot(modeloRF)
En este punto del análisis predictivo, podríamos decidir recortar el número de árboles del modelo RF, pues vemos que a partir de 75/80 árboles el error permanence constante y el modelo podría tender al overfitting:
modeloRF2=randomForest(x=DataTrain, y=DataTrain$TARGET, ntree = 50, importance =TRUE)
print(modeloRF2)
##
## Call:
## randomForest(x = DataTrain, y = DataTrain$TARGET, ntree = 50, importance = TRUE)
## Type of random forest: regression
## Number of trees: 50
## No. of variables tried at each split: 1
##
## Mean of squared residuals: 0.04361353
## % Var explained: 97.72
Tras una drástica reducción de árboles, seguimos en un porcentaje de explicación de variabilidad de Target similar al modelo inicial de RF, por lo que, dado que el algoritmo Random Forest parece ser, de entre los probados, el que mejor responde en sus predicciones para la variable Target, apostaríamos por seguir trabajando nuestro conjunto de datos con RF.