Chapitre 6 Faire de la régression sur R

Les problèmes de régression et de classification supervisée consistent à expliquer et/ou prédire une sortie $y\in\mathcal Y$ avec

$\mathcal Y=\mathbb R$ pour la régression
$\mathcal Y$ de cardinal fini pour la classification supervisée,

par des entrées $x\in\mathbb R^p$ . Il s’agit donc de trouver une fonction $m:\mathbb R^p\to\mathcal Y$ à partir de données $(X_1,Y_1),\dots,(X_n,Y_n)$ .

Ces données sont souvent collectées dans un dataframe df de la forme

$Y$	$X_1$	$X_2$	$\dots$	$X_p$
$y_1$	$x_{1,1}$	$x_{1,2}$	$\dots$	$x_{1,p}$
$\vdots$	$\vdots$	$\vdots$	$\vdots$	$\vdots$
$\vdots$	$\vdots$	$\vdots$	$\vdots$	$\vdots$
$y_n$	$x_{n,1}$	$x_{n,2}$	$\dots$	$x_{n,p}$

Le protocole pour construire un algorithme de régression sur R est toujours le même. Il faut spécifier :

la méthode (ou l’algorithme)
la variable à expliquer
les variables explicatives
le jeu de données
les éventuelles options de la méthode considérée.

Par exemple la commande

method(Y~X1+X3,data=df,...)

ajustera le modèle method pour expliquer $Y$ par $X_1$ et $X_3$ avec les données dans df (les points représentent d’éventuelles options). Voici quelques exemples de méthodes :

fonction R	algorithme	Package	Problème
lm	modèle linéaire		Reg
glm	modèle logistique		Class
lda	analyse discriminante linéaire	MASS	Class
svm	Support Vector Machine	e1071	Class
knn.reg	plus proches voisins	FNN	Reg
knn	plus proches voisins	class	Class
rpart	arbres	rpart	Reg et Class
glmnet	ridge et lasso	glmnet	Reg et Class
gbm	boosting	gbm	Reg et Class
randomForest	forêts aléatoires	randomForest	Reg et Class

Remarque: pour glmnet, on ne peut pas utiliser de formule de la forme Y~.. Il faut spécifier une matrice pour les $X$ et un vecteur pour $Y$ . La fonction model.matrix peut se révéler très utile pour calculer la matrice des $X$ .

Puisqu’il existe un grand nombre d’algorithmes pour répondre à un même problème de régression, il est important de définir des critères de performance afin de les comparer. Ces critères sont généralement inconnus et doivent être estimés à l’aide de procédure de type apprentissage/validation ou validation croisée. On a souvent besoin d’utiliser la fonction predict pour calculer ces critères. Cette fonction est une fonction générique : on peut utiliser predict pour une régression linéaire, logistique, un arbre, une forêt aléatoire… Pour obtenir l’aide de cette fonction pour

la régression linéaire : taper help(predict.lm)
la régression logisitque : taper help(predict.glm)
les régressions pénalisées : taper help(predict.glmnet)
les arbres : taper help(predict.rpart)
les forêts aléatoires : taper help(predict.randomForest)
…

Dans la suite on suppose que $\mathcal Y=\mathbb R$ et on considère le modèle de régression $Y=m(X)+\varepsilon.$ La performance d’un estimateur $\widehat{m}$ de $m$ sera mesurée par son erreur quadratique de prédiction : $E[(Y-\widehat m(X))^2].$

6.1 Modèle linéaire : fonctions lm et predict

On considère le modèle de régression linéaire $Y=\beta_0+\beta_1X_1+\dots+\beta_pX_p+\varepsilon$ où $X_1,\dots,X_p$ sont les variables explicatives, $Y$ la variable à expliquer et $\varepsilon$ le terme d’erreur. On fixe $p=5$ et on considère les données suivantes :

n <- 1000
p <- 5
set.seed(1234)
X.mat <- matrix(rnorm(n*p),ncol=p)
eps <- rnorm(n,mean = 0,sd=0.5)
df <- data.frame(X.mat,eps)
df <- df %>% mutate(Y=X1+X2+X3+X4+X5+eps) %>% select(-eps)

Construire un modèle linaire permettant d’expliquer $Y$ par $X_1,\dots,X_5$ (utiliser la fonction lm) et afficher les estimateurs de $\beta_0,\dots,\beta_5$ (on pourra utiliser les fonctions coef et summary).
On considère le jeu de données test suivant.
```
m <- 500
p <- 5
set.seed(12345)
X.mat <- matrix(rnorm(m*p),ncol=5)
eps <- rnorm(m,mean = 0,sd=0.5)
df.test <- data.frame(X.mat,eps)
df.test <- df.test %>% mutate(Y=X1+X2+X3+X4+X5+eps) %>% select(-eps)
```
Calculer, pour chaque individu de ce nouveau jeu de données, les prédictions faites par le modèle de la question précédente (utiliser la fonction predict avec l’option newdata).
Créer un nouveau dataframe qui contiennent les valeurs prédites $\widehat y_i$ à la question précédente sur une colonne et les valeurs observées $y_i$ du jeu de données df.test sur une autre colonne.
A l’aide du verbe summarize, calculer l’erreur quadratique moyenne (estimée) du modèle linéaire : $\frac{1}{m}\sum_{i\in test}(\widehat y_i-y_i)^2.$

6.2 Sélection de variables

On considère les données suivantes

n <- 1000
p <- 105
set.seed(1234)
X.mat <- matrix(rnorm(n*p),ncol=p)
eps <- rnorm(n,mean = 0,sd=0.5)
df <- data.frame(X.mat,eps)
df <- df %>% mutate(Y=X1+X2+X3+X4+X5+eps) %>% select(-eps)

issues du modèle $Y=\beta_0+\beta_1X_1+\dots+\beta_pX_p+\varepsilon$ avec $p=105$ . On remarquera que seules les variables $X_1,\dots,X_5$ sont explicatives.

Ajuster un modèle linéaire (fonction lm) sur df et afficher les estimateurs de $\beta_0,\dots,\beta_{105}$ .
On propose d’utiliser une procédure de sélection de variables backward à partir du critère BIC. Effectuer cette procédure à l’aide de la fonction step (on pourra utiliser les options direction=“backward” et k=log(n)). On appellera ce modèle mod.step.

On a sélectionné un modèle avec 8 variables : les 5 explicatives et 3 variables de bruit.

Calculer les erreurs quadratiques de prévision $\frac{1}{m}\sum_{i\in test}(\widehat y_i-y_i)^2$ des deux modèles (le modèle complet et le modèle sélectionné) en utilisant le jeu de données test suivant.
```
m <- 300
p <- 105
set.seed(12345)
X.mat <- matrix(rnorm(m*p),ncol=p)
eps <- rnorm(m,mean = 0,sd=0.5)
df.test <- data.frame(X.mat,eps)
df.test <- df.test %>% mutate(Y=X1+X2+X3+X4+X5+eps) %>% select(-eps)
```

6.3 Régression logistique et arbre

On considère le jeu de données spam disponible ici

library(kernlab)
data(spam)

Le problème est d’expliquer la variable type (un email est un spam ou non) par les 57 autres variables.

Séparer les données en un échantillon d’apprentissage dapp de taille 3000 et un échantillon test dtest de taille 1601. On pourra utiliser la fonction sample.
Construire un modèle logistique permettant de répondre au problème en utilisant uniquement les données d’apprentissage. On utilisera la fonction glm avec l’option family="binomial".
A l’aide de la fonction step, effectuer une sélection backward (ça peut prendre quelques minutes).
A l’aide de la fonction rpart du package rpart, construire un arbre de régression (toujours sur les données d’apprentissage) pour répondre au problème. On utilisera les paramètres par défaut de la fonction.
Visualiser l’arbre construit à l’aide des fonctions rpart.plot et visTree des packages rpart.plot et visNetwork
Pour les 3 modèles construits (logistique, backward et arbre) calculer les prédictions de la variable type pour les individus de l’échantillon dtest. On pourra regrouper ces prévisions dans un data-frame à 3 colonnes.
Ajouter au data-frame précédent une colonne où on mettra les valeurs observées de la variable à expliquer.
A l’aide de summarize_at calculer les erreurs de classification des 3 modèles.
Représenter les courbes ROC et calculer les AUC. On pourra consulter les pages 346 et 347 dans Cornillon et al. (2018) pour le tracé de courbes ROC sur R.

Références

Cornillon, P. A., A. Guyader, F. Husson, N. Jégou, J. Josse, N. Klutchnikoff, E. Le Pennec, E. Matzner-Løber, L. Rouvière, and B. Thieurmel. 2018. R Pour La Statistique et La Science Des Données. PUR. https://r-stat-sc-donnees.github.io.