Sélection de modèle - régression logistique

On considère la sélection de variable et la pénalisation dans le cadre de la régression logistique. Les mêmes impératifs existent que précédemment, soit le risque de surajustement et la nécessité de valider la performance avec des données de réserve (échantillon/entraînement vs validation vs test), et potentiellement de le faire plusieurs fois si nous avons plusieurs critères à optimiser.

Le modèle de régression logistique est estimé par maximum de vraisemblance; il n’y a pas de solution explicite pour les estimations. Cela rend l’estimation plus coûteuse, mais les critères d’information restent disponible. Chaque variable explicative ajoutée au modèle améliorera la vraisemblance des données d’entraînement, donc il est nécessaire de pénaliser ou de considérer la performance sur les données gardées en réserve, potentiellement via validation croisée. Le LASSO se généralise aussi de la même manière.

Si on veut choisir un modèle optimal, on peut prendre celui qui maximise un critère (score propre) sur des données de validation. Pour des données \(\boldsymbol{y} = (y_1, \ldots, y_n)\) avec \(y_i \in \{0,1\}\), c’est le cas du score de Brier (ici positivement orienté) \[S(\boldsymbol{y}, \boldsymbol{p}) = -\sum_{i=1}^n (y_i - p_i)^2\] ou du score logarithmique \[S(\boldsymbol{y}, \boldsymbol{p}) = \sum_{i=1}^n \{\log p_i y_i + \log(1-p_i)(1-y_i)\}.\] Le code ci-dessous cherche (selon une recherche séquentielle, exhaustive, ou via un algorithme génétique) le modèle avec le plus petit BIC — l’utilisation de la pénalisation cela nous évite de refaire une séparation de l’échantillon. Par la suite, on fait le calcul du point de coupure optimal (par validation croisée), et on obtient la classification.

En principe, il faudrait regarder la performance en faisant une autre validation, mais cela impliquerait nécessairement une validation croisée imbriquée dans une autre. Nous regardons donc la performance directement sur les 100K données tests, en calculant le profit (revenus moins dépenses) pour notre envoi d’échantillon promotionnel.

La première partie du code charge les paquets R, divise les bases de données en échantillons d’apprentissage et de validation et définit la formule du modèle avec toutes les interactions et les termes quadratiques.

library(hecmulti)
library(dplyr, warn.conflicts = FALSE)
data(dbm, package = "hecmulti")
## Code pour la sélection de variables
set.seed(60602)
# 
data(dbm, package = "hecmulti")
valid <- dbm[dbm$test == 1, ] |>
  dplyr::select(!c(ymontant, test))
appr <- dbm[dbm$test == 0, ] |>
  dplyr::select(!c(ymontant, test))
# Formule du modèle avec toutes les interactions
# d'ordre 2 (.^2) et les termes quadratiques I(x^2)
formule <- formula(
  yachat ~
    .^2 +
      I(x2^2) +
      I(x6^2) +
      I(x7^2) +
      I(x8^2) +
      I(x9^2) +
      I(x10^2)
)
# Nouvelles bases de données avec toutes ces variables
# On retire la première colonne (1, ordonnée à l'origine)
appr_c <- data.frame(
  cbind(model.matrix(formule, data = appr)[, -1]),
  y = as.integer(appr$yachat)
)
valid_c <- data.frame(
  cbind(model.matrix(formule, data = valid)[, -1]),
  y = as.integer(valid$yachat)
)
valid_ymontant <- with(dbm, ymontant[test == 1L])

On peut ensuite procéder à l’ajustement des différents modèles: celui avec toutes les variables de bases, l’autres avec les 105 coefficients qui mène à un surajustement, ici lié à la quasi-séparation de variables.

Pour chaque modèle, on répète les étapes suivantes

• ajustement ou sélection du modèle avec pénalisation BIC
• calcul du point de coupure optimal avec notre matrice de gain
• classification des données test
• estimation du profit

# Modèle avec toutes les variables principales
base <- glm(yachat ~ ., data = appr, family = binomial)
# Calcul du point de coupe optimal
#  (par validation croisée)
base_coupe <- hecmulti::select_pcoupe(
  modele = base,
  c00 = 0,
  c01 = 0,
  c10 = -10,
  c11 = 57
)
# Performance sur données de validation
base_pred <-
  predict(object = base, newdata = valid, type = "response") > base_coupe$optim
base_perfo <-
  -10 * sum(base_pred) + sum(valid_ymontant[base_pred], na.rm = TRUE)

# Modèle avec toutes les variables + interactions
# Ajustement
complet <- glm(formula = formule, data = appr, family = binomial)
# Sélection du point de coupure
complet_coupe <- hecmulti::select_pcoupe(
  modele = complet,
  c00 = 0,
  c01 = 0,
  c10 = -10,
  c11 = 57
)

Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1

# Performance sur données de validation
complet_pred <-
  predict(object = complet, newdata = valid, type = "response") >
    complet_coupe$optim
# Revenu
complet_perfo <-
  -10 * sum(complet_pred) + sum(valid_ymontant[complet_pred], na.rm = TRUE)

# Sélection de modèle avec algorithme glouton
# Recherche séquentielle (AIC)
seqAIC <- step(
  object = complet,
  direction = "both", # séquentielle
  k = 2, # AIC
  trace = 0
)

Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1

seqAIC_coupe <-
  hecmulti::select_pcoupe(
    modele = seqAIC,
    c00 = 0,
    c01 = 0,
    c10 = -10,
    c11 = 57
  )
seqAIC_pred <-
  predict.glm(object = seqAIC, newdata = valid, type = "response") >
    seqAIC_coupe$optim
seqAIC_perfo <-
  -10 * sum(seqAIC_pred) + sum(valid_ymontant[seqAIC_pred], na.rm = TRUE)
# Recherche séquentielle (BIC)
seqBIC <- step(
  object = complet,
  direction = "both", # séquentielle
  k = log(nobs(complet)), #BIC
  trace = 0
)

Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1

seqBIC_coupe <- hecmulti::select_pcoupe(
  modele = seqBIC,
  c00 = 0,
  c01 = 0,
  c10 = -10,
  c11 = 57
)
seqBIC_pred <-
  predict.glm(object = seqBIC, newdata = valid, type = "response") >
    seqBIC_coupe$optim
seqBIC_perfo <-
  -10 * sum(seqBIC_pred) + sum(valid_ymontant[seqBIC_pred], na.rm = TRUE)

# Recherche exhaustive par algorithme génétique
# avec moins de variables
appr_r <- data.frame(
  cbind(model.matrix(seqAIC)[, -1], y = appr$yachat)
)
valid_r <- data.frame(
  model.matrix(formula(seqAIC), data = valid)[, -1]
)
library(glmulti)

Le chargement a nécessité le package : rJava

Le chargement a nécessité le package : leaps

exgen <- glmulti::glmulti(
  y = y ~ .,
  #nombre de variables limitées
  data = appr_r,
  level = 1, # sans interaction
  method = "g", # recherche génétique
  crit = "bic", # critère (AIC, BIC, ...)
  confsetsize = 1, # meilleur modèle uniquement
  plotty = FALSE,
  report = FALSE, # sans graphique ou rapport
  fitfunction = "glm"
)

TASK: Genetic algorithm in the candidate set.
Initialization...
Algorithm started...
Improvements in best and average IC have bebingo en below the specified goals.
Algorithm is declared to have converged.
Completed.

# Redéfinir le modèle via "glm"
exgen_modele <-
  glm(exgen@objects[[1]]$formula, data = appr_r, family = binomial)
exgen_coupe <-
  hecmulti::select_pcoupe(
    modele = exgen_modele,
    c00 = 0,
    c01 = 0,
    c10 = -10,
    c11 = 57
  )
exgen_pred <-
  predict(exgen_modele, newdata = valid_r, type = "response") >
    exgen_coupe$optim
exgen_perfo <-
  -10 * sum(exgen_pred) + sum(valid_ymontant[exgen_pred], na.rm = TRUE)

Pour le LASSO, la fonction hecmulti::select_pcoupe n’est pas adaptée. On démontre plutôt comment calculer le profit avec les prédictions individuelles des montants d’achats, pondérés par la probabilité d’achat suite à l’envoi d’un échantillon par publipostage.

# Calcul avec LASSO

formule <-
  formula(
    ~ (x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10)^2 +
      I(x2^2) +
      I(x6^2) +
      I(x7^2) +
      I(x8^2) +
      I(x9^2) +
      I(x10^2)
  )

xmat_train_lm <- model.matrix(
  formule,
  data = dbm |>
    dplyr::filter(test == 0, !is.na(ymontant))
)[, -1]
y_train_lm <- with(dbm, na.omit(ymontant[test == 0]))

xmat_train_log <- model.matrix(
  formule,
  data = dbm |>
    dplyr::filter(test == 0)
)[, -1]
y_train_log <- with(dbm, yachat[test == 0])


xmat_test <- model.matrix(
  formule,
  data = dbm |>
    dplyr::filter(test == 1)
)[, -1]
# Créer une grille de pénalités
lambda_seq <- seq(from = 0.01, to = 5, by = 0.01)
cv_lasso_lm <- glmnet::cv.glmnet(
  # Attention: la fonction `glmnet` prend une matrice de modèle
  # et un vecteur de réponses
  x = xmat_train_lm,
  y = y_train_lm,
  alpha = 1, # garder cette valeur fixe à un pour le lasso
  lambda = lambda_seq
)
# Pénalité qui minimise l'EQM de validation croisée
lambopt_lm <- cv_lasso_lm$lambda.min #ou cv_output$lambda.1se
lasso_lm <-
  glmnet::glmnet(
    x = xmat_train_lm,
    y = y_train_lm,
    alpha = 1,
    lambda = lambopt_lm
  )

pred_ymontant_lasso <- predict(
  lasso_lm, # modèle duquel prédire
  s = lambopt_lm, # valeur de la pénalité optimale
  newx = xmat_test
)
# Validation croisée pour trouver le critère
cvfit_log <- glmnet::cv.glmnet(
  x = xmat_train_log,
  y = y_train_log,
  family = "binomial"
)
lambda_opt_log <- cvfit_log$lambda.min
# Le critère par défaut est la déviance (-2ll)
# Ajuster modèle avec pénalisation
lasso_log <- glmnet::glmnet(
  x = xmat_train_log,
  y = y_train_log,
  family = "binomial",
  lambda = lambda_opt_log
)
# Prédire sur l'échantillon test (probabilités)
probs_lasso <- predict(
  lasso_log,
  newx = xmat_test,
  s = lambda_opt_log,
  type = "resp"
)
# Trouver les individus avec un profit moyen supérieur à 10$
decision <- probs_lasso * pred_ymontant_lasso > 10
# Envoyer échantillon si profit > 10
# Calcul du profit (revenus moins coûts)
lasso_perfo <- sum(dbm$ymontant[dbm$test == 1] * decision, na.rm = TRUE) -
  10 * sum(decision)
lasso_perfo

[1] 991063

C’est parmi les stratégies optimales. On pourrait combiner n’importe lequel des autres modèles précédemment ajustés et faire le même genre de calcul.

modèle	no. variables	pt. coupure	sensibilité	taux bonne classif.	profit
(a)			0.23200	1.0000000	601212
(b)	14	0.12	0.70938	0.4376047	940569
(c)	104	0.10	0.76148	0.4914570	925263
(d)	28	0.16	0.79298	0.5339846	963682
(e)	8	0.16	0.78206	0.5178908	984326
(f)	10	0.19	0.79795	0.5422227	950083