TP1 - Régression simple et régression multiple

Author

PB

Quarto

Ce document est écrit en quarto, un langage formaté qui permet de mélanger du texte (en markdown) et du code (ici, du R).

Nous utilisons ici Rstudio, un tutoriel de présentation est disponible ici. Vous pouvez utiliser un autre IDE en fonction de votre préférence, mais ce TP part du principe que vous utilisez Rstudio.

Vous pouvez partir de ce document pour créer votre rapport, ou bien en créer un nouveau depuis Rstudio : File > New File > Quarto Document....

R et documentation

Dans ce TP, nous utilisons plusieurs fonctions de R, qui sont toutes signalées dans le texte. Lorsqu’une nouvelle fonction apparaît, je vous invite à lire l’aide associée à la fonction dans R directement, en tapant directement dans la console ? suivit du nom de la fonction. Par exemple, pour avoir de l’aide sur la fonction lm :

?lm

L’aide de R, écrite directement par les personnes ayant développée les fonctions, est en générale bien écrite. Il est utile de savoir naviguer dans ces pages de documentation. Pour les fonctions de base que nous utilisons, votre LLM favori vous fournira bien entendu des exemples et explications détaillées. L’une des forces de R est que c’est un langage modulable, avec de nombreuses librairies écrites par des chercheurs et chercheuses, pour résoudre des problèmes avancés. Ces librairies auront toujours une documentation, mais ne seront pas toujours bien représentées par les LLM, car nouvelles. Dans une perspective de recherche et développement, il est donc utile de savoir lire cette documentation.

De manière similaire, ce TP en entier peut être “résolu” parfaitement par un LLM, les questions étant très classiques. Comme pour les exercices de TD, l’objectif principal est de vous entraîner, et de traduire les notions vues en cours en des exemples concrets d’application. Je vous invite donc à tenter de répondre aux questions sans aide, au moins dans un premier temps.

Données

On utilise les données Boston qui sont présentes dans le package MASS. Ce TP est inspiré de la section 3.6 du livre An Introduction to Statistical Learning, de G. James, D. Witten, T. Hastie and R. Tibshirani.

Le package MASS est standard, et s’installe si besoin avec la commande :

install.packages("MASS")

On peut ensuite charger et voir les données :

library(MASS)
head(Boston)
     crim zn indus chas   nox    rm  age    dis rad tax ptratio  black lstat
1 0.00632 18  2.31    0 0.538 6.575 65.2 4.0900   1 296    15.3 396.90  4.98
2 0.02731  0  7.07    0 0.469 6.421 78.9 4.9671   2 242    17.8 396.90  9.14
3 0.02729  0  7.07    0 0.469 7.185 61.1 4.9671   2 242    17.8 392.83  4.03
4 0.03237  0  2.18    0 0.458 6.998 45.8 6.0622   3 222    18.7 394.63  2.94
5 0.06905  0  2.18    0 0.458 7.147 54.2 6.0622   3 222    18.7 396.90  5.33
6 0.02985  0  2.18    0 0.458 6.430 58.7 6.0622   3 222    18.7 394.12  5.21
  medv
1 24.0
2 21.6
3 34.7
4 33.4
5 36.2
6 28.7

Une description du jeu de données peut être obtenue dans l’aide, en tapant dans la console :

?Boston
  1. Décrire brièvement le jeu de données. Qu’est-ce qu’un “individu” \(i\) dans ces données ? Combien y en a t il ? Combien y a t il de variables ? Décrivez brièvement ces variables, en précisant leur type (discret ou continu).

Régression simple

On s’intéresse d’abord au prix médian des maisons dans une ville (medv) en fonction du pourcentage de foyers avec un statut socio-économique faible (lstat).

  1. Tracer le nuage de point correspondant à ces deux variables. On pourra utiliser la librairie ggplot2 et le code ci-dessous. Que pensez-vous de la relation entre les deux variables ? Les hypothèses classiques du modèle linéaire semblent elles vérifiées ? (Ces points sont à discuter de manière qualitative, nous verrons plus tard dans le cours des manières plus rigoureuses de tester ces hypothèses).
library(ggplot2)
ggplot(Boston, aes(x = lstat, y = medv)) + 
  geom_point()
  1. Faire la régression linéaire simple de medv en fonction de lstat à l’aide de la fonction lm. Interprétez les coefficients. Sont-ils significativement non nuls ? Donnez un intervalle de confiance à \(90%\) des coefficients. On pourra utiliser les fonctions summary, coef et confint, en utilisant le squelette suivant :
fit1 <- lm(medv ~ lstat, data = Boston)
summary(fit1)
confint(fit1)
  1. Tracez la droite de régression sur le nuage de points en utilisant le code suivant. Que pensez-vous de l’ajustement du modèle ? Est-ce cohérent avec le résultat de l’ajustement linéaire ?
ggplot(Boston, aes(x = lstat, y = medv)) + 
  geom_point() + 
  geom_smooth(method = "lm")

  1. Donnez un intervalle de prédiction à \(95%\) pour la valeur de medv lorsque lstat vaut \(2\), \(15\), \(30\) et \(50\) pour cent à l’aide de la fonction predict. Quand est-ce que l’intervalle est le plus grand ? Est-ce surprenant ? Discutez de la validité des prédictions, en particulier pour lstat=50.

  2. Par défaut, ggplot dans le code de la fonction ci-dessus trace une zone d’incertitude en grisé. S’agit-il d’une bande de confiance, ou d’une bande de prédiction ? Quelle est la différence ?

  3. Ajoutez une bande de prédiction à l’aide du code (à compléter) suivant. Interprétez.

data_pred <- data.frame(lstat = seq(min(Boston$lstat), max(Boston$lstat), by = 0.1))
pred_grid <- predict(lm.fit, data_pred, interval = "prediction") # lm.fit is the fit object, result of `lm` call.
pred_grid <- as.data.frame(pred_grid)
pred_grid$lstat <- data_pred$lstat
ggplot(Boston, aes(x = lstat, y = medv)) + 
  geom_point() + 
  geom_smooth(method = "lm") +
  geom_ribbon(data = pred_grid, aes(y = fit, ymin = lwr, ymax = upr), alpha = 0.1)

Régression multiple

En régression multiple, on explique une variable à prédire (ici, medv), en fonction de plusieurs prédicteurs. Nous verrons en cours très prochainement comment cette régression multiple est formulée d’un point de vue mathématique et statistique. D’un point de vue pratique, il s’agit simplement d’ajouter des termes à la formule de prédiction, par exemple :

fit2 <- lm(medv ~ lstat + indus, data = Boston)

représente le modèle linéaire suivant : \[ \mbox{mdev}_i = \beta_0 + \beta_1 \times \mbox{lstat}_i + \beta_2 \times \mbox{indus}_i + \epsilon_i \] où \(\mbox{mdev}_i\), \(\mbox{lstat}_i\) et \(\mbox{indus}_i\) sont les variables mesurées pour l’individu \(i\) (\(1 \leq i \leq n\)), \(\epsilon_i\) l’erreur associée au modèle, qui est supposée i.i.d. gaussienne, centrée et de variance \(\sigma^2\) (\(\epsilon_i \sim \mathcal{N}(0, \sigma^2)\)). La commande lm trouve les estimateurs \(\hat{\beta}_0\), \(\hat{\beta}_0\), \(\hat{\beta}_0\) et \(\hat{\sigma}^2\) associés au modèle (par la méthode des moindres carrés), et les fonctions summary, confint, predict, etc utilisent la distribution de ces estimateurs pour réaliser des tests statistiques ou construire des intervalles de confiance.

  1. Faire la régression linéaire simple de medv en fonction de lstat et indus. Interprétez les coefficients en utilisant la fonction summary. Sont-ils significativement non nuls ? Qu’est-ce que cela signifie pour l’influence des variables lstat et indus sur la variable à prédire medv ? Comparez le \(R^2\) associé à ce modèle à celui associé au modèle de régression simple de la section précédente. Que pouvez-vous en dire ?

  2. Faire la régression de medv contre toutes les autres variables du jeu de données. On pourra utiliser la formule medv ~ ., qui est un raccourci pour éviter d’avoir à écrire toutes les variables à la main. Interprétez les résultats. Y-a-t-il des coefficients non significatifs ? Que vaut le \(R^2\) ?

  3. En utilisant la formula medv ~ . -age - indus, faites la régression de medv contre toutes les variables, sauf age et indus. Interprétez. Quel modèle préférez-vous ?

  4. Nous verrons en cours comment utiliser un test de Fisher pour choisir de manière rigoureuse entre les deux modèles. Ce test peut être fait à l’aide de la fonction anova. En utilisant la documentation, utilisez cette fonction, et interprétez. Est-ce cohérent avec les résultats de la question précédente ?

  5. Faites la régression linéaire de medv contre les seules variables age et indus. Interprétez le résultat.

Croissance du \(R^2\)

  1. Générez une variable aléatoire bruit indépendante des autres variables, avec la loi de votre choix. Ajouter cette variable à la régression avec la formule medv ~ . + bruit. Comment a évolué le \(R^2\) ? Est-ce normal ? Que ce passe-t-il si vous ajoutez plusieurs variables de bruit ?

Nous verrons en cours pourquoi ce comportement est attendu, et plusieurs stratégies pour y remédier.

Transformation de données

  1. Faites la régression de medv contre lstat et log(lstat) à l’aide de la formule medv ~ lstat + log(lstat). Interprétez.

  2. Faites la régression de medv contre lstat et lstat^2 à l’aide de la formule medv ~ lstat + I(lstat^2). Interprétez.

Note : la fonction I() est utile ici pour protéger la fonction carrée, l’opérateur ^ ayant une autre signification dans le cadre d’une formule.

  1. Faites la régression de medv contre lstat, lstat^2, lstat^3, …, lstat^10 à l’aide de la formule medv ~ poly(lstat, 10). Interprétez.

Remarque : la fonction poly retourne les polynomes orthogonalisés, et non les puissances simples, comme détaillé dans l’aide ?poly.

Cette section illustre deux points principaux : premièrement, qu’il est possible de faire des régression linéaires contre des variables préalablement transformer pour rendre compte de relations complexes, et deuxièmement qu’il peut être difficile de choisir la bonne transformation à appliquer.