Dotaki : Exercice 2, Question 1 avec analyse exploratoire et méthodologie détaillée
Antoine Villatte
06/11/2020
1 Chargement des packages
On utilisera les packages suivants :
library(ggplot2)
library(ggrepel)
library(data.table)
library(dplyr)
library(caret)
library(randomForest)
library(gbm)
library(gridExtra)
library(xgboost)
library(e1071)
library(knitr)
library(kableExtra)
library(glmnet)2 Chargement du dataset
df <- read.csv("ab_test_sample.csv", sep = ";")
# conversion en data table :
dt <- data.table(df)3 Partition des données
Le but final de cet exercice est de proposer une méthode de suggestion de la variante d’un site à afficher à un utilisateur.
Pour cela, on créera au final un modèle prédictif. Pour fiabiliser les résultats, on partitionne le dataset fourni dans l’exercice en deux partitions : un training set et un testing set :
# On crée un objet inTrain, un vecteur indiquant quels éléments du dataset initial
# sont à intégrer au training set
set.seed(12)
inTrain <- createDataPartition(dt$converted, p = 0.7, list = FALSE)
# On crée le training set à partir de ce vecteur
training <- dt[inTrain, ]
# Tous les éléments du dataset initial n'ayant pas été intégrés dans le training
# set le sont dans le testing set
testing <- dt[-inTrain, ]3.1 Concaténation des partitions
Pour faciliter la manipulation des données, on concatène les training et testing set, mais le testing set aura des valeurs manquantes NA pour la variable de réponse converted.
Ainsi, tout preprocessing ou feature engineering peut-être fait de la même façon sur les deux partition à la fois, et l’analyse exploratoire se fait également facilement sur les lignes du dataset où les valeurs de converted ne sont pas NA
Avant cela, on fait deux choses : on sauvegarde les valeurs de converted du training set dans un vecteur testing_converted avant de les supprimer, et on vérifie qu’il n’y ait aucune valeur manquante de base dans le dataset fourni;
sum(is.na(dt$converted))## [1] 0Il n’y a aucune valeur manquante pour la variable converted. Il n’y a d’ailleurs aucune valeur manquante dans le dataset en général :
sum(is.na(dt))## [1] 0On sauvegarde alors les valeurs de converted, puis on les remplace par des NA :
testing_converted <- testing$converted
testing$converted <- NAet on concatène les deux sets :
comb <- data.table(rbind(training, testing))4 Analyse exploratoire
4.1 Variables
dim(comb)## [1] 8538 9Notre dataset contient 8538 observations et 9 variables.
Regardons quelles sont ces variables
names(comb)## [1] "user_id_zws" "variationname" "converted" "A"
## [5] "O" "N" "E" "C"
## [9] "segment"- Les variables A, O, N, E, C et converted sont numériques
- Les variables variationname et segment sont des châines de caractère et doivent être converties en variables catégoriques (dites “factor”)
- La variable *user_id_zws" est une variable d’identification que l’on supprime du dataset
# Conversion en variables catégoriques :
comb[, variationname := factor(variationname)]
comb[, segment := factor(segment)]
# Supression de la variable ID :
comb <- select(comb, -user_id_zws)4.2 Distribution des données
Regardons comment sont distribuées les données :
qplot(data = comb[!is.na(comb$converted), ], x = converted, geom = "histogram") +
facet_wrap(variationname ~ .) +
stat_bin(binwidth = 1, geom = "text", aes(label = ..count..), hjust = -0.5) +
ggtitle("Distribution de la variable de réponse 'converted' ")
qplot(data = comb[!is.na(comb$converted), ], x = A, geom = "histogram")+
ggtitle("Distribution de la variable A")
qplot(data = comb[!is.na(comb$converted), ], x = O, geom = "histogram")+
ggtitle("Distribution de la variable O")
qplot(data = comb[!is.na(comb$converted), ], x = N, geom = "histogram")+
ggtitle("Distribution de la variable N")
qplot(data = comb[!is.na(comb$converted), ], x = E, geom = "histogram")+
ggtitle("Distribution de la variable E")
qplot(data = comb[!is.na(comb$converted), ], x = C, geom = "histogram")+
ggtitle("Distribution de la variable C")
ggplot(data = comb[!is.na(comb$converted), ], aes(x = segment)) + geom_bar(stat = "count") +
ggtitle("Distribution de la variable 'segment' ") +
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1))
ggplot(data = comb[!is.na(comb$converted), ], aes(x = variationname)) + geom_bar(stat = "count") +
ggtitle("Distribution de la variable 'variationname' ") 
On remarque que les distributions des données numériques sont assez peu centrées ; on devra les pré-processer pour fiabiliser les résultats des modèles prédictifs.
On remarque également que la prévalence des utilisateurs convertis est faible, environ 10% quelle que soit la variation, avec un taux un peu plus élevées pour les variations 2 et 1 :
comb[!is.na(comb$converted), ] %>% group_by(variationname) %>%
summarise(taux_conv = (sum(converted)/ length(converted)))## # A tibble: 3 x 2
## variationname taux_conv
## <fct> <dbl>
## 1 Original 0.112
## 2 Variation 1 0.116
## 3 Variation 2 0.117On évitera alors de transformer la variable converted en variable catégorique. Si l’on crée un modèle de classification, il aura tendance à classifier tous les utilisateurs comme “non convertis” (converted = 0) pour obtenir une précision de 90%.
On préfèrera donc un modèle de régression, qui proposera des probabilités de conversion. Pour les mêmes raisons que pour le modèle de classification, la faible prévalence des convertis fera que les modèles donneront surement des probabilités de conversion faibles quoiqu’il arrive.
Notre objectif sera donc d’obtenir un modèle qui montre une nette séparation entre les probabilités de conversion des utilisateurs effectivement convertis et celles de ceux qui ne le sont pas.
4.3 Préparation des données
4.4 Séparation entre variables numériques et catégoriques
Les variables numériques et catégoriques seront traitées différemment, on doit donc les séparer :
# Liste des noms des variables numériques :
numcols <- sapply(comb, function(y){(class(y) == "integer") | (class(y) == "numeric")})
# Variables numériques :
num_comb <- comb[, .SD, .SDcols = numcols]
# Variables catégoriques
cat_comb <- comb[, .SD, .SDcols = -numcols]La variable converted doit être retirée des variables numériques puisqu’on ne souhaite pas la modifier :
num_comb <- select(num_comb, -converted)4.5 Normalisation des variables numériques
Les modèles classiques de machine learning préfèrent en général avoir en entrée des variables avec une distribution gaussienne centrée. On pré-process donc les données dans ce sens :
PreObj <- preProcess(num_comb, method = c("center", "scale"))
normnum_comb <- predict(PreObj, num_comb)
# conversion en data.table
normnum_comb <- data.table(normnum_comb)4.6 Encodage des variables catégoriques
Les modèles de machine learning classiques accueillent en général assez mal les variables catégoriques. On les convertit donc en variable numériques dites “dummy variables”
On utilisera ici le one-hot encoding (OHE), plus adapté aux variables catégoriques contenant plus de deux niveaux :
# Création d'un objet "ohe" qui permettra de transformer les variables catégoriques
# en dummies
ohe <- dummyVars(~ ., data = cat_comb)
# Utilisation de l'objet ohe sur le dataset cat_comb
ohecat_comb <- predict(ohe, cat_comb)
# Conversion de la matrice en data.table
ohecat_comb <- data.table(ohecat_comb)4.7 Compilation des jeux de données processés
Nos données étant maintenant pré-processées, on les concatène ensemble :
preproc_comb <- data.table(cbind(normnum_comb, ohecat_comb), converted = comb$converted)On sépare à présent le training set du testing set. Pour rappel, le training set contient toutes les observations où les valeurs de converted ne sont pas NA.
On supprimera également la variable converted du testing set puisque c’est celle que l’on cherche à prédire.
preproc_train <- preproc_comb[!is.na(preproc_comb$converted)]
preproc_test <- preproc_comb[is.na(preproc_comb$converted)] %>% select(- converted)5 Model Building
5.1 XGBM
Notre premier modèle sera un modèle XGBM, en général très performant et assez peu sensible à l’overfitting.
La première étape sera une recherche automatisée du meilleur calibrage des hyper paramètres.
Cette étape demande énormément de calculs, on utilisera donc les packages parallel et doParallel pour faire du parrallel computing.
library(parallel)
library(doParallel)
cluster1 <- makeCluster(detectCores() - 2) # on alloue tous les coeurs sauf 2
# au calcul
registerDoParallel(cluster1) # Démarrage du parallel computing
# l'entraînement du modèle se fera avec 5 cross-validations
trctrl <- trainControl(method = "cv", number = 5, allowParallel = T,
verboseIter = F)
# fourchette de recherche du meilleur hyperparamétrage
xgbGrid <- expand.grid(nrounds = 750,
max_depth = seq(2, 10, by = 2),
colsample_bytree = seq(0.2, 0.6, length.out = 5),
eta = seq(0.005, 0.015, length.out = 3),
gamma=0,
min_child_weight = c(0.8, 1, 1.2),
subsample = 1
)
# Fige la RNG de l'ordinateur pour que l'on puisse retomber sur le même modèle
# si on refait l'entrainement
set.seed(1234)
modfit_xgbm <- caret::train(converted ~ ., data = preproc_train,
method = "xgbTree",
trControl = trctrl,
tuneGrid = xgbGrid, verbose = F)
# Arrêt du parallel processing :
stopCluster(cluster1)
registerDoSEQ()On regarde quels sont les meilleurs hyperparamètres :
modfit_xgbm$bestTune## nrounds max_depth eta gamma colsample_bytree min_child_weight subsample
## 1 750 2 0.005 0 0.2 0.8 1On les enregistre comme ceux par défaut :
xgbparam <- list(objective = "binary:logistic",
booster = "gbtree",
max_depth = 2,
colsample_bytree = 0.2,
eta = 0.005,
gamma=0,
min_child_weight = 0.8,
subsample = 1
)On fait maintenant un deuxième série de cross-validations pour déterminer nombre d’itérations suffisantes pour avoir le meilleur modèle.
Si le modèle ne s’améliore pas après 5 itérations, la recherche s’arrête.
set.seed(123456)
xgb_cv <- xgb.cv(params = xgbparam, data = as.matrix(preproc_train),
nrounds = 10000, nfold = 5, showsd = T, stratified = T,
print_every_n = 15, early_stopping_rounds = 5, maximize = F,
label = preproc_train$converted)## [1] train-error:0.045425+0.055640 test-error:0.046509+0.057051
## Multiple eval metrics are present. Will use test_error for early stopping.
## Will train until test_error hasn't improved in 5 rounds.
##
## Stopping. Best iteration:
## [2] train-error:0.022966+0.045932 test-error:0.022241+0.044482Deux itérations suffisent. C’est anormalement peu, ce qui signifie que le modèle a tendance à vouloir classifier toutes les probabilités comme étant très faible, c’est à dire classifier toutes les valeurs de converted comme étant 0 pour avoir 90% de précision.
On peut d’ores et déjà prédire que le modèle XGBM séparera assez peu les données des convertis par rapport aux non-convertis.
On connait maintenant tous les hyperparamètres pour notre modèle XGBM. On fixera cependant le nrounds à 20 au lieu de 2.
trctrl <- trainControl(method = "cv", number = 5, allowParallel = T,
verboseIter = F)
xgbGrid2 <- expand.grid(nrounds = 20,
max_depth = 2,
colsample_bytree = 0.2,
eta = 0.005,
gamma=0,
min_child_weight = 0.8,
subsample = 1
)
set.seed(1234)
modfit_xgbm2 <- caret::train(converted ~ ., data = preproc_train,
method = "xgbTree",
trControl = trctrl,
tuneGrid = xgbGrid2, verbose = F)## [17:09:12] WARNING: amalgamation/../src/objective/regression_obj.cu:174: reg:linear is now deprecated in favor of reg:squarederror.
## [17:09:12] WARNING: amalgamation/../src/objective/regression_obj.cu:174: reg:linear is now deprecated in favor of reg:squarederror.
## [17:09:12] WARNING: amalgamation/../src/objective/regression_obj.cu:174: reg:linear is now deprecated in favor of reg:squarederror.
## [17:09:12] WARNING: amalgamation/../src/objective/regression_obj.cu:174: reg:linear is now deprecated in favor of reg:squarederror.
## [17:09:12] WARNING: amalgamation/../src/objective/regression_obj.cu:174: reg:linear is now deprecated in favor of reg:squarederror.
## [17:09:12] WARNING: amalgamation/../src/objective/regression_obj.cu:174: reg:linear is now deprecated in favor of reg:squarederror.Notre modèle est prêt. On peut déjà observer quelles sont les variables ayant le plus d’importance dans sa prédiction :
imp_mat <- xgb.importance(feature_names = colnames(preproc_train), model = modfit_xgbm2$finalModel)
xgb.ggplot.importance(importance_matrix = imp_mat, rel_to_first = T)
5.2 Random Forest
Notre deuxième modèle sera un RandomForest.
On utilisera la cross-validation définie pour le modèle XGBM pour définir le meilleur nombre de branches à utiliser.
library(parallel)
library(doParallel)
cluster1 <- makeCluster(detectCores() - 2)
registerDoParallel(cluster1)
set.seed(12345)
modfit_rf <- caret::train(converted ~ ., data = preproc_train,
method = "rf",
trControl = trctrl)
stopCluster(cluster1)
registerDoSEQ()On peut regarder quels sont les hyperparamètres retenus pour le meilleur résultat :
modfit_rf$finalModel##
## Call:
## randomForest(x = x, y = y, mtry = param$mtry)
## Type of random forest: regression
## Number of trees: 500
## No. of variables tried at each split: 2
##
## Mean of squared residuals: 0.1006672
## % Var explained: 0.41Ce modèle nous permet également de visualiser quels sont les variables les plus importantes dans ses prédictions :
imp_rf <- importance(modfit_rf$finalModel)
imp_rf <- data.table(variables = row.names(imp_rf), inc_mse = imp_rf[, 1])
imp_rf <- arrange(imp_rf, desc(inc_mse))
ggplot(data = imp_rf,
aes(x = reorder(variables, inc_mse), y = inc_mse, fill = inc_mse)) +
geom_bar(stat = "identity") +
theme(legend.position = "none")+
labs(x = "Variables",
y = "% d'augmentation du MSE si la variable est réordonnée au hasard")+
coord_flip()
5.3 SVM
Notre dernier modèle sera un SVM, très performant en classification mais qui a une forte tendance à l’overfitting.
library(parallel)
library(doParallel)
cluster1 <- makeCluster(detectCores() - 2)
registerDoParallel(cluster1)
set.seed(12345)
# Fourchette de recherche des meilleurs hyperparamètres :
tuned = tune.svm(factor(converted) ~ ., data = preproc_train,
cost = c(seq(0.01, 1, by = 0.1), seq(1, 10, by = 1), seq(10,100, by=10)),
coef0 = c(0.5, 1, 1.5),
tunecontrol=tune.control(cross = 5))
stopCluster(cluster1)
registerDoSEQ()On regarde quels sont les meilleurs hyperparamètres retenus :
tuned$best.parameters## coef0 cost
## 1 0.5 0.01tuned$best.model##
## Call:
## best.svm(x = factor(converted) ~ ., data = preproc_train, coef0 = c(0.5,
## 1, 1.5), cost = c(seq(0.01, 1, by = 0.1), seq(1, 10, by = 1),
## seq(10, 100, by = 10)), tunecontrol = tune.control(cross = 5))
##
##
## Parameters:
## SVM-Type: C-classification
## SVM-Kernel: radial
## cost: 0.01
##
## Number of Support Vectors: 1458Et on construit un modèle SVM avec ces hyperparamètres :
set.seed(123456)
modfit_svm <- svm(converted ~ ., data = preproc_train,
type = "eps-regression",
kernel = "radial",
cost = 0.01,
coef0 = 0.5,
epsilon = 0.1)6 Predictions sur le training set
Avant de les tester sur le testing set, regardons comment se comportent ces modèles sur les données sur lesquelles on les a entraînés.
On regroupe dans une data table toutes les prédictions ainsi que les valeurs réelles de converted
tp_xgbm <- predict(modfit_xgbm2, preproc_train)
tp_rf <- predict(modfit_rf, preproc_train)
tp_svm <- predict(modfit_svm, preproc_train)
tp_dt <- data.table(xgbm = tp_xgbm,
rf = tp_rf,
svm = tp_svm,
reference = factor(preproc_train$converted))On trace ensuite les valeurs prédites par chacun des modèles pour quand la valeur réelle de converted est 0 ou 1 :
qplot(data = tp_dt, x = reference, y = xgbm, geom = "boxplot") +
coord_cartesian(ylim = c(0.46,0.465)) + ggtitle("XGBM") +
labs(x = "Valeur réelle", y = "Valeur prédite")
qplot(data = tp_dt, x = reference, y = rf, geom = "boxplot") +
ggtitle("Random Forest") + labs(x = "Valeur réelle", y = "Valeur prédite")
qplot(data = tp_dt, x = reference, y = svm, geom = "boxplot") +
coord_cartesian(ylim = c(0.0313,0.0320)) + ggtitle("SVM") +
labs(x = "Valeur réelle", y = "Valeur prédite")
Les prédictions faites avec les modèles XGBM et SVM séparent assez peu les probabilités de conversions attribuées aux profils convertis par rapport à celles de ceux qui ne le sont pas, tandis que le modèle Random Forest affiche une nette séparation.
C’est donc ce modèle que l’on utilisera par la suite. Regardons comment évolue cette séparation lorsque le modèle est appliqué au training set :
pred_rf <- predict(modfit_rf, preproc_test)
tempdt <- data.table(reference = testing_converted, predicted = pred_rf)
qplot(data = tempdt, x = factor(reference), y = predicted, geom = "boxplot") +
ggtitle("RandomForest") +
labs(x = "Valeur réelle", y = "Valeur prédite")
Cette séparation est beaucoup moins nette ; le modèle a donc tendance à faire du overfitting sur son training set, mais on garde néanmoins des probabilités de conversion légèrement supérieures pour les profils effectivement convertis.
7 Déterminer la meilleure variante pour un utilisateur
Pour un utilisateur donné, en fonction de ses paramètres A,O,N,E,C et Segment, on lui applique le modèle prédictif pour chacune des 3 variantes, et on retient celle qui donne la plus haute probabilité de conversion.
7.1 Exemple
Prenons le premier utilisateur de notre testing set :
preproc_test[1, ]## A O N E C variationname.Original
## 1: -1.134807 -0.4784724 -1.257618 1.154172 0.2066722 0
## variationname.Variation 1 variationname.Variation 2 segment.Altruists
## 1: 0 1 0
## segment.Analyticals segment.Attentives segment.Audacious
## 1: 0 0 1
## segment.Conventionals segment.Emotives segment.Hedonists segment.Homebodies
## 1: 0 0 0 0
## segment.Straightforwards segment.Trend Setters
## 1: 0 0On crée un dataset avec 3 fois ce même utilisateur :
- Une fois avec variationname.Original = 1 et les autres 0
- Une fois avec variationname.Variation1 = 1 et les autres 0
- Une fois avec variationname.Variation2 = 1 et les autres 0
method_test <- data.table(rbind(preproc_test[1, ], preproc_test[1, ], preproc_test[1, ]))
method_test[1, ]$`variationname.Variation 2` <- 0
method_test[1, ]$variationname.Original <- 1
method_test[2, ]$`variationname.Variation 2` <- 0
method_test[2, ]$`variationname.Variation 1` <- 1On prédit le taux de conversion pour chaque configuration :
method_test$prediction <- predict(modfit_rf, method_test)
method_test$variation <- c("Original", "Variation 1", "Variation 2")
kable(method_test[, .SD, .SDcols = c("variation", "prediction")]) %>% kable_styling(bootstrap_options = c("striped", "hover"), position = "left") | variation | prediction |
|---|---|
| Original | 0.1090783 |
| Variation 1 | 0.1253192 |
| Variation 2 | 0.1223705 |
La Variation 1 donne la probabilité de conversion la plus haute ; c’est donc celle-ci que l’on proposera à l’utilisateur.