Enunciado del problema

A través de un conjunto de datos (TABLON_ENTRENAMIENTO.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 TABLON_VALIDACION.csv. Para ello, elija la métrica que mejor se adapte a la naturaleza de los datos. Pista: la variable TARGET es cuantitativa continua.

Este reto se debe resolver en RStudio y presentar los resultados a través de un documento (html o pdf) creado con RMarkdown.

Debe indagar sobre diferentes naturalezas de modelos predictivos y mostrar 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.

Sólo se puntuarán las contribuciones que resuelvan (de forma parcial o total) el caso planteado.

Si sólo se participa en el segundo se podrá optar a una calificación máxima de 8 puntos. Si participan en el primero o en ambos podrán optar a una calificación máxima planteada.

Ejercicio

Carga de librería y datos

Cargo todas las librerías vistas hasta este punto. Es posible que muchas no se usen, pero por si acaso.

library(mltools)
library(data.table)
library(ggplot2)
library(dplyr)
library(GGally)
library(rpart)
library(rpart.plot)
library(randomForest)
library(caret)
library(mltools)
library(e1071)
library(corrplot)
library(visreg)
library(pROC)
library(scatterplot3d)
library(ggdendro)
library(cluster)
library(factoextra)
library(devtools)
library(ggbiplot)
library(nnet)
library(readr)
library(dplyr)
library(tidyverse)
library(lubridate)
library(data.table)
library(ggplot2)
library(csv)
library(modeest)
library(moments)
library(plotrix)
library(nortest)
library(vcd)
library(gplots)
library(corrplot)
library(car)
library(rpart)
library(rpart.plot)
library(caret)
library(useful)
library(randomForest)
library(tidytext)
library(naivebayes)
library(tm)

A continuación se cargan los datos, que ya están divididos en el conjunto de muestra y de prueba y creo una data.frame que une a ambos dos anteriores para así analizar/filtrar los datos de manera conjunta.
De este modo, puedo corregir (si hubiese que hacerlo), todos los NAs y datos vacíos de ambos conjuntos a la vez (esto lo haría aquí sólo puesto que es un estudio academico, si no lo haría en cada conjunto por separado). Tmbién voy a fijar una seed en este punto, asi uso la misma para todos los modelos de regresión.

## 'data.frame':    5400 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 ...
## TARGET   VAR1   VAR2   VAR3   VAR4 
##      0      0      0      0      0
## [1] 0.127681 4.619065

Puesto que no hay valores nulos, no hace falta filtrar más. Voy a mirar también ahora que rango de valores toman los los datos:

apply(total,2,function(x) length(unique(x)))
## TARGET   VAR1   VAR2   VAR3   VAR4 
##   5397     50     12      3      3
apply(total[2:5],2,function(x) unique(x))
## $VAR1
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50
## 
## $VAR2
##  [1]  0.50  0.75  1.00  2.50  5.00  7.50 10.00  0.80  2.00  4.00  6.00  8.00
## 
## $VAR3
## [1] 1 3 6
## 
## $VAR4
## [1] 300 400 500

Para ver si existen algun relación más o menos fuerte entre TARGET y las VARi, voy a hacer primero una visualización de los datos.

ggplot(total,aes(VAR1,TARGET)) +  geom_point()

ggplot(total,aes(VAR2,TARGET)) +  geom_point()

ggplot(total,aes(VAR3,TARGET)) +  geom_point()

ggplot(total,aes(VAR4,TARGET)) +  geom_point()

Se puede ver que con respecto a VAR1 hay una relación creciente, los valores medios crecen según aumenta el valor de VAR1 y, por otro lado, también hay bastante variabilidad, la cual parece aumentar hasta más o menos VAR1=20, dónde comienza a mantenerse estable.
Para las otras variables es algo más dificil apreciar su relación con TARGET, aunque se puede decir que, en general, para cada valor de cada una de VARi, TARGET puede oscilar a lo largo de su rango completo de valores (Target Min= 0.127681 Target Max= 4.619065)

Dicha variabilidad tamién la podemos ver en la matriz de covarianzas:

M=cor(total)
M
##             TARGET      VAR1          VAR2          VAR3          VAR4
## TARGET  1.00000000 0.9084015  1.275580e-01 -2.439513e-02 -3.366129e-02
## VAR1    0.90840145 1.0000000  0.000000e+00  0.000000e+00  0.000000e+00
## VAR2    0.12755801 0.0000000  1.000000e+00 -6.848732e-19 -4.376592e-20
## VAR3   -0.02439513 0.0000000 -6.848732e-19  1.000000e+00 -1.960702e-21
## VAR4   -0.03366129 0.0000000 -4.376592e-20 -1.960702e-21  1.000000e+00

A su vez, puesto que las matriz de correlación muestra valores muy bajos de covarianza entre Target y VAR2/3 y 4, se puede decir que no hay una relación creceinte/decreciente entre ellas (podría haber algún otro tipo, por ejemplo, funciones trigonométricas). Además, las covarianzas entre las variables son también prácticamente cero, lo que indicaría independencia de la variables entre sí. Miro a ver si, por casualidad, realizar una reducción de dimensionalidad pudiese mejorar el modelo:

pc <- prcomp(total,center = T,scale. = T)
summary(pc)
## Importance of components:
##                           PC1    PC2    PC3    PC4     PC5
## Standard deviation     1.3850 1.0000 1.0000 1.0000 0.28591
## Proportion of Variance 0.3836 0.2000 0.2000 0.2000 0.01635
## Cumulative Proportion  0.3836 0.5837 0.7836 0.9837 1.00000

Los 4 primeros PC parece que son capaces de explicar bastante bien la variabilidad del modelo, con un 98% de la variabilidad explicada. Además, tienen una SD de 1 o más, con lo que sería lógico usarlos. Por otro lado, en nuestro modelo tenemos 4 variables independientes, por lo que realizar esta dimensionalidad no parecería que pudiese producir una gran mejora.

Modelos de predicción

Regresión lineal

Dado que la variable “Target” es una variable continua, se va comenzar con una regresión linear. Ahora ya se va a volver al las dos data bases sugeridas en el enunciado “train” y “test, que ya se cargaron anteriormente.

Creamos el modelo lineal y lo entrenamos, para ver su métrica:

modelo=lm(TARGET~VAR1+VAR2+VAR3+VAR4,data = train)
modelo
## 
## Call:
## lm(formula = TARGET ~ VAR1 + VAR2 + VAR3 + VAR4, data = train)
## 
## Coefficients:
## (Intercept)         VAR1         VAR2         VAR3         VAR4  
##   0.9167535    0.0868552    0.0574772   -0.0196577   -0.0005762
summary(modelo)
## 
## Call:
## lm(formula = TARGET ~ VAR1 + VAR2 + VAR3 + VAR4, data = train)
## 
## 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

Aqui podemos ver que el modelo se formaría de la siguiente manera: Target=0.0868552VAR1+0.0574772VAR2-0.0196577VAR3-0.0005762VAR4+0.9167535 Este modelo tiene un R2 de 0.8453, es decir, explica un 84.35% de la variabiliad de TARGET. A parte de este analisis, tenemos que ver otra serie de métricas, para saber si es bueno o no. Éstas son:

  • Normaliad de los residuos.

  • Las variables son independientes entre si.

  • Homocedasticidad.

  • No autocorrelación en los residuos.

plot(modelo)

  1. La gráfica de Residuals vs. Fitter se ve que sigue una parabola, lo cual parece indicar que hay una relación no lineal que el modelo no es capaz de prever. Esto podría indicar que no es un buen modelo de predicción.

  2. El Q-Q de los residuales también parece desviarse bastante en los extremos, por lo que tampoco parece ajustarse bien a la normalidad.

  3. La curva Scale-Location tiene una forma muy característica en forma de arcos. Es dificil saber si están repartidos de manera homogenea por encima y por debajo de la lines “quasi” horizontales, pero se puede decir que claramente siguen una curva casi en forma de acueducto.

  4. Residuals v Leverage parece no indicar que haya claros outlaier.

Visto lo anterior, no parece muy correcto querer seguir con este modelo, a pesar de un valor de R2 relativamente bueno.

predicciones_lr=predict(modelo,newdata=test)
summary(predicciones_lr)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.6435  1.9862  3.0757  3.0747  4.1576  5.5268
tabla_lr1=cbind(predicciones_lr, test[1])
LR1=c(sqrt(mean((tabla_lr1[,1] - tabla_lr1[,2])^2)),R2(tabla_lr1[,1],tabla_lr1[,2]))
print(sprintf("Mean of squared errors (MSE): %.3f", LR1[1]))
## [1] "Mean of squared errors (MSE): 0.560"
# R2 score
print(sprintf("R2 score: %.3f", LR1[2]))
## [1] "R2 score: 0.844"

Voy a ver los Box-Plots, a ver como afectan los diferentes valores de las variables al valor de Target, para saber si tal vez la eliminación de alguna varible pueda mejorar el análisis. Antes las tengo que covertir a factor:

total$VAR1=as.factor(total$VAR1)
total$VAR2=as.factor(total$VAR2)
total$VAR3=as.factor(total$VAR3)
total$VAR4=as.factor(total$VAR4)
train$VAR1=as.factor(train$VAR1)
train$VAR2=as.factor(train$VAR2)
train$VAR3=as.factor(train$VAR3)
train$VAR4=as.factor(train$VAR4)
test$VAR1=as.factor(test$VAR1)
test$VAR2=as.factor(test$VAR2)
test$VAR3=as.factor(test$VAR3)
test$VAR4=as.factor(test$VAR4)
ggplot(train, aes(VAR1,TARGET)) +  geom_boxplot()

ggplot(train, aes(VAR2,TARGET)) +  geom_boxplot()

ggplot(train, aes(VAR3,TARGET)) +  geom_boxplot()

ggplot(train, aes(VAR4,TARGET)) +  geom_boxplot()

Los gráficos anteriores muestran muy poca variabilidad de las medias y variabilidad en función de VAR3 y VAR4. Por ello, voy a hacer un modelo lineal sólo con VAR1 y VAR2, a ver si los gráficos mejoran. Vuelvo a cargar los archivos, puesto que los había convertido a factores.

train=read.csv("TABLON_ENTRENAMIENTO.csv", header = T, sep = ";")
test=read.csv("TABLON_VALIDACION.csv", header = T, sep = ";")
modelo2=lm(TARGET~VAR1+VAR2,data = train)
modelo2
## 
## Call:
## lm(formula = TARGET ~ VAR1 + VAR2, data = train)
## 
## Coefficients:
## (Intercept)         VAR1         VAR2  
##     0.62075      0.08686      0.05748
summary(modelo2)
## 
## Call:
## lm(formula = TARGET ~ VAR1 + VAR2, data = train)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.42200 -0.49951  0.08893  0.47908  1.01663 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.6207512  0.0228058   27.22   <2e-16 ***
## VAR1        0.0868552  0.0006803  127.68   <2e-16 ***
## VAR2        0.0574772  0.0028469   20.19   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.551 on 3147 degrees of freedom
## Multiple R-squared:  0.8415, Adjusted R-squared:  0.8414 
## F-statistic:  8355 on 2 and 3147 DF,  p-value: < 2.2e-16
plot(modelo2)

Se puede obersvar que los gráficos anteriores no cambian mucho con respecto al primer modelo, con lo que se puede quedar tal y como está, como modelo 1.

predicciones_lr=predict(modelo2,newdata=test)
summary(predicciones_lr)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.7536  1.9798  3.0795  3.0747  4.1512  5.4233
tabla_lr2=cbind(predicciones_lr, test[1])
LR2=c(sqrt(mean((tabla_lr2[,1] - tabla_lr2[,2])^2)),R2(tabla_lr2[,1],tabla_lr2[,2]))
print(sprintf("Mean of squared errors (MSE): %.3f", LR2[1]))
## [1] "Mean of squared errors (MSE): 0.562"
# R2 score
print(sprintf("R2 score: %.3f", LR2[2]))
## [1] "R2 score: 0.843"

EL R2 sigue igual que el anterior, relativamente alto. Sin embargo, el análisis de las gráficas de residuos me hacen dudar de la validez de este método como predictor.

Regresión KNN

Voy a lanzar otro de los métodos vistos. En este caso, un KNN:

trainControl <- trainControl(method="repeatedcv", number=10, repeats=3)
modelo_knn<- train(TARGET~., data=train, method="knn", metric="RMSE", trControl=trainControl)
modelo_knn
## k-Nearest Neighbors 
## 
## 3150 samples
##    4 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold, repeated 3 times) 
## Summary of sample sizes: 2835, 2835, 2836, 2835, 2834, 2836, ... 
## Resampling results across tuning parameters:
## 
##   k  RMSE       Rsquared   MAE       
##   5  0.1139772  0.9935571  0.08428741
##   7  0.1230879  0.9923868  0.09010853
##   9  0.1163597  0.9930433  0.08153529
## 
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was k = 5.
plot(modelo_knn)

El R2 es muy alto, lo cual podría ser un indicador de overfitting. Lo comprobaré con el modelo de test para lo que lanzo el predictor:

predict_knn <- predict(modelo_knn,test,type = "raw")
summary(predict_knn)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.1505  2.2878  3.6771  3.1242  4.2491  4.5987
tabla_knn=cbind(predict_knn, test[1])
TKNN=c(sqrt(mean((tabla_knn[,1] - tabla_knn[,2])^2)),R2(tabla_knn[,1],tabla_knn[,2]))
print(sprintf("Mean of squared errors (MSE): %.3f", TKNN[1]))
## [1] "Mean of squared errors (MSE): 0.062"
# R2 score
print(sprintf("R2 score: %.3f", TKNN[2]))
## [1] "R2 score: 0.999"

El R2 da mejores resultados que en el modelo de train, con lo que parece que no está haciendo un overfitting.

Regresion con Support Vector Machine (SVM)

Pruebo con un modelo de Support Vector Machines (SVM)

modelo_svm <- train(TARGET~., data=train, method="svmLinear", metric="RMSE", trControl=trainControl)
modelo_svm
## Support Vector Machines with Linear Kernel 
## 
## 3150 samples
##    4 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold, repeated 3 times) 
## Summary of sample sizes: 2837, 2834, 2834, 2835, 2835, 2835, ... 
## Resampling results:
## 
##   RMSE       Rsquared   MAE      
##   0.5708391  0.8418944  0.4680641
## 
## Tuning parameter 'C' was held constant at a value of 1

El valor de R2 es peor que en el caso anterior, al igaul que el RMSE. En cualquier caso, calculo los valores al igual que en modelos anteriores:

predict_svm <- predict(modelo_svm,test,type = "raw")
summary(predict_svm)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.8856  2.2051  3.1858  3.1894  4.1868  5.5162
tabla_svm=cbind(predict_svm, test[1])
TSVM=c(sqrt(mean((tabla_svm[,1] - tabla_svm[,2])^2)),R2(tabla_svm[,1],tabla_svm[,2]))
print(sprintf("Mean of squared errors (MSE): %.3f", TSVM[1]))
## [1] "Mean of squared errors (MSE): 0.573"
# R2 score
print(sprintf("R2 score: %.3f", TSVM[2]))
## [1] "R2 score: 0.842"

Se ve que tanto el modelo “train” como los valores predichos con el modelo para el set “test” tienen un R2 muy similar.

Random Forest

modelo_rf<-randomForest(x=train,
               y=train$TARGET,
               importance=TRUE
               )
modelo_rf
## 
## Call:
##  randomForest(x = train, y = train$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.06572344
##                     % Var explained: 96.57
plot(modelo_rf)

En tornoe a los 75 arboles, el modelo no adquire una mejoría mucho mayor. El porcentaje de la variabilidad explicada es del 96.57%, muy alta. Veremos si está haciendo overfitting:

predict_rf <- predict(modelo_rf,test)
tabla_rf=cbind(predict_rf, test[1])
TRF=c(sqrt(mean((tabla_rf[,1] - tabla_rf[,2])^2)),R2(tabla_rf[,1],tabla_rf[,2]))
print(sprintf("Mean of squared errors (MSE): %.3f", TRF[1]))
## [1] "Mean of squared errors (MSE): 0.247"
# R2 score
print(sprintf("R2 score: %.3f", TRF[2]))
## [1] "R2 score: 0.990"

Parece que no, al estar producieno un R2 también muy alto

Resultados

Para ver la comparativa de todos los modelos, creo una table conjunta:

resumen<-rbind(LR1,LR2,TKNN,TSVM,TRF)
colnames(resumen)=c("MSE","R2")
resumen
##             MSE        R2
## LR1  0.56018241 0.8442649
## LR2  0.56246637 0.8429789
## TKNN 0.06215875 0.9985039
## TSVM 0.57322679 0.8423358
## TRF  0.24734879 0.9901354

Se puede ver que, en este caso, el mejor modelo a la hora de predecir es el KNN, con el mejor R2 y MSE, seguido del Random Forest.

Puesto que hay dos modelos capaces de predecir bastante bien (KNN y RF), no procedo con la opción de reducción de dimensionalidad, que por otro lado y dada la independencia entre variables sugerida por la matriz de covarianzas, podría indicar que no mejoraría mucho el modelo.