Foro docente Asig. 3 (Actividad 2)

Máster en Big Data & Business Intelligence

Andrés Pina García

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.

1) Carga de librerías

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

Carga de datasets y previsualización de datos

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.

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).

Visualización de datos

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í.

Modelos predictivos y entrenamiento del modelo

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.