Faites une analyse exploratoire des variables explicatives:
Quelles variables devraient êtres exclues de la modélisation? Justifiez votre réponse.
Comparez la variable réponse avec les autres variables: y a-t-il des transformations qui améliorerait l’ensemble de variables candidates: interactions, création de variables dychotomiques, transformations (racines carrée, transformation logarithmique, etc.)?
Vérifiez s’il y a des variables catégorielles encodées comme des variables numériques.
On remarque en consultant la documentation de la base de données à l’aide de ?college que le nombre de demandes d’admissions (napplications), le nombre d’admissions parmi ces applications nadmission et le nombre d’offres converties par les candidat.e.s, ninscrits, sont reliées et ne peuvent être employées.
Plusieurs variables seront fortement corrélées parce qu’elles dépendent de la capacité d’accueil de l’établissement d’enseignement. Ainsi, on pourrait créer une variable qui représente le pourcentage de temps partiels pour les premiers cycles (plutôt que le décompte). Les variables pourcentdoctorat et pourcentterminal, sont fortement corrélées puisque la plupart des diplômes terminaux (incluant les titres professionnels, les doctorats de premier cycle en médecine, etc.) sont des PhD: le modèle choisira la variable la plus adéquate. Il n’y a pas de variables catégorielles hors variables binaires. Côté vérifications d’usage, on note que tauxdiplom et pourcentdoctorat sont supérieurs à 100%.
On devra retirer la variable catégorielle nom, qui a une modalité différente pour chaque observation.
Scindez la base de données en échantillon avec données d’entraînement (environ 2/3 des données) et échantillon de validation; utilisez le germe aléatoire 60602 via set.seed(60602).
Sélectionnez un modèle à l’aide d’une des méthodes couvertes, mais en basant votre choix sur l’erreur moyenne quadratique évaluée sur l’échantillon de validation.
On peut désormais considérer une séparation en tiers. Pour ce faire, je vais échantillonner des variables logiques vrais et faux avec une cote de 2 pour 1 et ensuite sélectionner les lignes qui correspondent.
Pour l’estimation, on pourrait simplement calculer l’erreur quadratique moyenne de validation pour tous les cinquante premiers modèles en considérant toutes les interactions d’ordre 2.
Code
# Modèle avec interactions de toutes les variablesformule <-formula(napplications ~ .^2)# Sélection séquentielle ascendanterec_seq <- leaps::regsubsets(x = formule, data = db_a,method ="seqrep",nvmax =50)
Si regsubsets permet de recouvrer les différents modèles, il n’y a pas d’utilitaire pour obtenir le résumé et les prédictions. Le paquet hecmulti inclut une fonction pour faire les prédictions d’un modèle donné. À des fins d’illustration, on considère celui qui a le plus petit critère BIC intra-échantillon d’apprentissage.
Code
# Calculer les BIC bic_mod <-summary(rec_seq)$bic# Prédire du modèle avec id variables# Ici, celui avec le plus petit BIC du lotpredict(rec_seq, id =which.min(bic_mod))
Ici, on s’intéresse uniquement au calcul de l’erreur quadratique moyenne pour les modèles choisis par la procédure de recherche séquentielle. La fonction eval_EQM_regsubsets fait le calcul et nous retourne un vecteur avec toutes les mesures. On choisit le modèle qui minimise l’erreur.
Code
# Modèle avec la plus petite erreur # moyenne quadratique de validation externereqm_seq_ve_list <- hecmulti::eval_EQM_regsubsets(model = rec_seq, select =1:50,formula = formule,data = db_a,newdata = db_v)# La fonction calcule la racine EQM sur # les données de validationmod_seq_ve <-which.min(reqm_seq_ve_list)nvar_seq_ve <- mod_seq_ve +1# nombre variables + ordonnée à l'origine.reqm_seq_ve <- reqm_seq_ve_list[mod_seq_ve]
Une logique similaire s’appliquerait avec le LASSO, même si c’est loin d’être la norme (on utilise d’ordinaire la validation croisée). Ici, on ajusterait le modèle avec plusieurs valeurs de \(\lambda\), puis on calculerait l’erreur quadratique moyenne de validation. On prendra la valeur de la pénalisation \(\lambda\) qui minimise cette dernière. Spécifier la séquence de valeurs de \(\lambda\) à essayer nécessite un peu d’essai/erreur.
Code
# Obtenir X et yXmat_a <-model.matrix(formule, data = db_a)[,-1]y_a <- db_a$napplicationslambda_seq <-exp(seq(from =log(0.1),to =log(50), length.out =100L))lasso <- glmnet::glmnet(x = Xmat_a,y = y_a,alpha =1, lambda = lambda_seq)# Prédictions et calcul de l'EQM# On pourrait remplacer `newx` par # d'autres données (validation externe)eqm_lasso_externe <-rep(0, length(lambda_seq))for(i inseq_along(lambda_seq)){ pred <-predict(object = lasso, s = lambda_seq[i], newx =model.matrix(formule, data = db_v)[,-1]) eqm_lasso_externe[i] <-mean((pred - db_v$napplications)^2)}minl <-which.min(eqm_lasso_externe)lambda_opt <- lambda_seq[minl]# Nombre de coefficients non nulsnvar_lasso_ve <- Matrix::nnzero(lasso$beta[,minl, drop =FALSE])reqm_lasso_ve <-sqrt(eqm_lasso_externe[minl])
Répétez la sélection, cette fois en prenant comme critère pour l’erreur moyenne quadratique évaluée par validation croisée (aléatoire) à cinq plis.
On peut utiliser boot::cv.glm ou caret::train avec les différentes méthodes pour ajuster les modèles. Je fais varier le nombre maximal de variables pour retourner les différentes solutions. Comme le résultat de la validation croisée est aléatoire, on peut répéter cette étape pour réduire l’incertitude de prédiction et obtenir une meilleur valeur de l’écart-type.
Code
# Matrice du modèle complet, moins ordonnée à l'origineXmat <-model.matrix(formule, data = db)[,-1]# Variable réponsey <- db$napplicationsseq_vc <- caret::train(form = formule,data = db,method ='leapSeq',tuneGrid =expand.grid(nvmax =1:50),trControl = caret::trainControl(method ="repeatedcv",number =10,repeats =10))# Graphique de la racine de l'EQM en fonction de# nvmaxplot(seq_vc)
Les notes de cours donnent une approche pas à pas avec glmnet pour ajuster le LASSO en choisissant la valeur de \(\lambda\) par validation croisée, mais on peut changer le type de model dans caret et utiliser les mêmes fonctionnalités.
# Meilleur modèlecoefs_lasso_cv <-coef(object = lasso_vc$finalModel, s = lasso_vc$bestTune$lambda)# Nombre de paramètres non-nuls pour la solutionnvar_lasso_vc <- Matrix::nnzero(coefs_lasso_cv)# Estimation de la racine de l'erreur quadratique moyennereqm_lasso_vc <-min(lasso_vc$results$RMSE)## Calcul des prédictions - comme avec n'importe quel modèle# predict(lasso_vc, newdata = db)
Créez un tableau avec le nombre de coefficients non-nuls de votre modèle final et un estimé de l’erreur moyenne quadratique obtenu par validation externe ou croisée.
Il suffit de colliger notre estimation de l’erreur et le nombre de coefficients.
Commentez sur le meilleur modèle parmi les combinaisons.
Si la taille de la base de données est plus conséquente, la validation croisée reste préférable à la validation externe, surtout si on répète cette méthode plusieurs fois pour réduire l’incertitude.1 Le nombre de paramètres dans le modèle LASSO peut sembler élevé, mais il faut garder en tête que les autres paramètres sont rétrécis vers zéro même s’ils sont non-nuls.
Impossible de savoir quel est le meilleur modèle: le modèle avec la plus petite erreur quadratique moyenne est notre meilleur choix ici à défaut d’autres observations. Rappelez-vous qu’un modèle peut donner une bonne performance, mais être choisi parce qu’il surajuste une valeur aberrante.
Le but de l’exercice était davantage de démontrer le travail que de faire une prédiction convaincante. Puisqu’on essaie de prédire une variable de dénombrement, un modèle plus adapté (modèle de régression binomiale négative, ou un modèle pour le log du nombre d’admission) pourrait probablement être plus adéquat dans la mesure où la variance de la réponse dépend de la taille de l’école.
Notes de bas de page
Rien ne nous empêche de faire de même avec la validation externe, soit dit en passant.↩︎
