2.1 Outils de diagnostic (suite)

La régression linéaire est une matière complexe et de nombreux outils existent pour vous aider à déterminer si le modèle que vous ajustez tient la route ou non. Il est très important de le vérifier avant d’utiliser un modèle. Ajuster un modèle quelconque dans des données est à la portée de tout le monde, mais choisir un modèle pertinent et pouvoir expliquer pourquoi est nettement plus difficile !

2.1.1 Résumé avec summary()(suite)

Reprenons la sortie renvoyée par summary() appliquée à un objet lm (notez ici que nous analysons la sortie de summary() “brute”, mais vous savez comment la transformer en un tableau plus présentable dans un rapport avec tabularise()).

SciViews::R("model", lang = "fr") # Configuration du système

trees <- read("trees", package = "datasets") # Importation des données
trees_lm <- lm(data = trees, volume ~ diameter) # Régression
summary(trees_lm) # Résumé du modèle
# 
# Call:
# lm(formula = volume ~ diameter, data = trees)
# 
# Residuals:
#       Min        1Q    Median        3Q       Max 
# -0.231211 -0.087021  0.003533  0.100594  0.271725 
# 
# Coefficients:
#             Estimate Std. Error t value Pr(>|t|)    
# (Intercept) -1.04748    0.09553  -10.96 7.85e-12 ***
# diameter     5.65154    0.27649   20.44  < 2e-16 ***
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 
# Residual standard error: 0.1206 on 29 degrees of freedom
# Multiple R-squared:  0.9351,  Adjusted R-squared:  0.9329 
# F-statistic: 417.8 on 1 and 29 DF,  p-value: < 2.2e-16

Nous n’avons pas encore étudié la signification des trois dernières lignes de ce résumé. Voici de quoi il s’agit.

  • Residual standard error :

Il s’agit de l’écart type des résidus, considérant que les degrés de liberté du modèle correspond au nombre d’observations \(n\) (ici 31) soustrait du nombre de paramètres à estimer (ici 2, la pente et l’ordonnée à l’origine de la droite). C’est donc une mesure globale de l’importance (c’est-à-dire de l’étendue) des résidus de manière générale et qui a les mêmes unités que y.

\[\sqrt{\frac{\sum(y_i - \hat{y_i})^2}{n-2}}\]

  • Multiple R-squared :

Il s’agit de la valeur du coefficient de détermination du modèle noté R2 de manière générale ou r2 dans le cas d’une régression linéaire simple. Il exprime la fraction de variance exprimée par le modèle. Autrement dit, le R2 quantifie la capacité du modèle à prédire la valeur de \(y\) connaissant la valeur \(x\) pour un individu. C’est donc une indication du pouvoir prédictif de notre modèle autant que de sa qualité d’ajustement (goodness-of-fit en anglais).

Souvenons-nous que la variance totale respecte la propriété d’additivité. La variance est composée au numérateur d’une somme de carrés, et au dénominateur de degrés de liberté. La somme des carrés totaux (de la variance, \(SC(total)\)) peut elle-même être décomposée en une fraction expliquée par notre modèle \(SC(rég)\), et la fraction qui ne l’est pas (les résidus \(SC(résidus)\)) :

\[SC(total) = SC(rég) + SC(résidus)\]

avec :

\[SC(total) = \sum_{i=0}^n(y_i - \bar{y_i})^2\]

\[SC(rég) = \sum_{i=0}^n(\hat{y_i} - \bar{y_i})^2\]

\[SC(résidus) = \sum_{i=0}^n(y_i - \hat{y_i})^2\]

À partir de la décomposition de ces sommes de carrés, le coefficient R2 (ou r2) se définit comme :

\[R^2 = \frac{SC(rég)}{SC(total)} = 1 - \frac{SC(résidus)}{SC(total)}\]

La valeur du R2 est comprise entre 0 (lorsque le modèle est très mauvais et n’explique rien) et 1 (lorsque le modèle est parfait et “capture” toute la variance des données ; dans ce dernier cas, tous les résidus valent zéro). Donc, plus le coefficient R2 se rapproche de 1, plus le modèle explique bien les données et aura un bon pouvoir de prédiction pour autant que l’on ne cherche pas à extrapoler des valeurs en dehors du domaine étudié.

Dans R, le R2 multiple se réfère simplement au R2 (ou au r2 pour les régressions linéaires simples) calculé de cette façon. L’adjectif multiple indique que le calcul est valable pour une régression multiple telle que nous verrons plus loin.

Par contre, le terme au dénominateur considère en fait la somme des carrés totale par rapport à un modèle de référence lorsque la variable dépendante \(y\) ne dépend pas de la ou des variables indépendantes \(x\). Les équations indiquées plus haut sont valables lorsque l’ordonnée à l’origine n’est pas figée (\(y = \alpha + \beta \ x\)). Dans ce cas, la valeur de référence pour \(y\) est bien sa moyenne, \(\bar y\).

D’un autre côté, si l’ordonnée à l’origine est fixée à zéro dans le modèle simplifié \(y = \beta \ x\) (avec \(\alpha = 0\) obtenu en indiquant la formule y ~ x + 0 ou y ~ x - 1), alors le zéro sur l’axe \(y\) est considéré comme une valeur appartenant d’office au modèle et devient valeur de référence. Ainsi, dans les équations ci-dessus il faut remplacer \(\bar y\) par 0 partout. Le R2 est alors calculé différemment, et sa valeur peut brusquement augmenter si le nuage de points est très éloigné du zéro sur l’axe Y. Ne comparez donc jamais les R2 obtenus avec et sans forçage à zéro de l’ordonnée à l’origine !

  • Adjusted R-squared :

La valeur du coefficient R2 ajusté, noté \(\bar{R^2}\) n’est pas utile dans le cadre de la régression linéaire simple, mais est indispensable avec la régression multiple. En effet, à chaque fois que vous rendez votre modèle plus complexe en ajoutant une ou plusieurs variables indépendantes, le modèle s’ajustera de mieux en mieux dans les données, même par pur hasard. C’est un phénomène que l’on appelle l’inflation du R2. À la limite, si nous ajoutons une nouvelle variable fortement corrélée avec les précédentes5, l’apport en matière d’information nouvelle sera négligeable, mais le R2 augmentera malgré tout un tout petit peu. Alors dans quel cas l’ajout d’une nouvelle variable est-il pertinent ? Le R2 ajusté apporte l’information désirée ici. Sa valeur n’augmentera pour l’ajout d’un nouveau prédicteur que si l’ajustement est meilleur que ce que l’on obtiendrait par le pur hasard. Le R2 ajusté se calcule comme suit (il n’est pas nécessaire de retenir cette formule, mais juste de constater que l’ajustement fait intervenir p, le nombre de paramètres du modèle et n, la taille de l’échantillon) :

\[ \bar{R^2} = 1 - (1 - R^2) \frac{n - 1}{n - p - 1} \]

  • F-statistic :

Tout comme pour l’ANOVA, la distribution F de Fisher permet de vérifier la significativité de la régression, car \(MS(rég)/MS(résidus)\) suit une distribution F à respectivement 1 et \(n-2\) degré de liberté, avec \(MS\) les carrés moyens, c’est-à-dire les sommes des carrés \(SC\) divisés par leurs degrés de liberté respectifs.

  • p-value :

Il s’agit de la valeur p associée à la statistique de F, donc à l’ANOVA de la régression linéaire. Pour cette ANOVA particulière, l’hypothèse nulle est que la droite n’apporte pas plus d’explication des valeurs de y à partir des valeurs de x que la valeur moyenne de y (ou zéro, dans le cas particulier d’un modèle dont l’ordonnée à l’origine est forcée à zéro). L’hypothèse alternative est donc que le modèle est significatif au seuil \(\alpha\) considéré. Donc, notre objectif est de rejeter H0 pour ce test ANOVA pour que le modèle ait un sens (valeur p plus petite que le seuil \(\alpha\) choisi).

Le tableau complet de l’ANOVA pour le modèle peut être obtenu à l’aide de anova() :

anova(trees_lm)
# Analysis of Variance Table
# 
# Response: volume
#           Df Sum Sq Mean Sq F value    Pr(>F)    
# diameter   1 6.0762  6.0762   417.8 < 2.2e-16 ***
# Residuals 29 0.4218  0.0145                      
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Nous y retrouvons les mêmes informations, fortement résumées en une ligne à la fin de la sortie de summary(). Ici elles sont présentées sous une forme plus classique de tableau de l’analyse de la variance. N’oubliez pas tabularise() si vous avez besoin d’un tableau propre pour un rapport, aussi bien pour anova() :

anova(trees_lm) |>
  tabularise()

Analyse de la variance
Réponse : volume

Terme

Ddl

Somme des carrés

Carrés moyens

Valeur de Fobs.

Valeur de p

diameter

1

6.076

6.0762

418

< 2·10-16

***

Résidus

29

0.422

0.0145

0 <= '***' < 0.001 < '**' < 0.01 < '*' < 0.05

… que pour summary() :

summary(trees_lm) |>
  tabularise()

Modèle linéaire

Volume de bois [m3]⁡=α+β(Diameˋtre aˋ 1,4m [m]⁡)+ϵ\operatorname{Volume\ de\ bois \ [m^{3}]} = \alpha + \beta_{}(\operatorname{Diamètre\ à\ 1,4m \ [m]}) + \epsilon

Terme

Valeur estimée

Ecart type

Valeur de t

Valeur de p

α\alpha

-1.05

0.0955

-11.0

7.85·10-12

***

β\beta_{}

5.65

0.2765

20.4

< 2·10-16

***

0 <= '***' < 0.001 < '**' < 0.01 < '*' < 0.05

Etendue des résidus : [-0.2312, 0.2717]

Ecart type des résidus : 0.1206 pour 29 degrés de liberté

R2 multiple : 0.9351 - R2 ajusté : 0.9329

Statistique F : 417.8 pour 1 et 29 ddl - valeur de p : < 2.22e-16

À vous de jouer !
h5p

2.1.2 Comparaison de régressions

Vous pouvez à présent comparer ces résultats avec une régression réalisée sans la valeur supérieure à 0,5m de diamètre. Attention ! On ne peut supprimer une valeur sans raison valable. La suppression de points aberrants doit être justifiée (par exemple, vous pouvez montrer qu’il y a eu une erreur de mesure, ou que l’individu mesuré est anormal ou malade, ou qu’il y a un autre problème dans le jeu de données, problème explicable). La raison de la suppression de ce point ici est liée au fait qu’il n’y a qu’un seul arbre mesuré dont le diamètre soit supérieur à 0.5m. Dès lors, nous pourrions souhaiter réduire l’étendue du jeu de données à des arbres plus petits (mais par la même occasion, nous réduisons aussi la portée du modèle, bien sûr). Le nettoyage du jeu de donnée peut se faire avant l’appel de lm(), ou plus simple, nous pouvons utiliser l’argument subset= de la fonction.

trees_lm2 <- lm(data = trees, volume ~ diameter, subset = diameter < 0.5)
chart(trees_lm2)

Le résumé de la régression donne ceci :

summary(trees_lm2) |>
  tabularise()

Modèle linéaire

volume⁡=α+β(diameter⁡)+ϵ\operatorname{volume} = \alpha + \beta_{}(\operatorname{diameter}) + \epsilon

Terme

Valeur estimée

Ecart type

Valeur de t

Valeur de p

α\alpha

-0.944

0.0931

-10.1

6.98·10-11

***

β\beta_{}

5.312

0.2754

19.3

< 2·10-16

***

0 <= '***' < 0.001 < '**' < 0.01 < '*' < 0.05

Etendue des résidus : [-0.2151, 0.1814]

Ecart type des résidus : 0.1082 pour 28 degrés de liberté

R2 multiple : 0.93 - R2 ajusté : 0.9275

Statistique F : 372.1 pour 1 et 28 ddl - valeur de p : < 2.22e-16

Sans oublier l’analyse des résidus.

chart$residuals(trees_lm2)

Comparez ces résultats avec ceux du modèle qui utilise tous les points. Lequel conserver ? Pas facile de décider. Il nous faut des métriques qui permettent de comparer de manière fiable différents modèles pour nous aider. Cela apparaîtra d’ailleurs d’autant plus utile que la situation va passablement se complexifier (dans le bon sens) avec l’introduction de la régression multiple et polynomiale.

À vous de jouer !
h5p

  1. La corrélation entre les prédicteurs dans un modèle linéaire multiple est un gros problème et doit être évitée le plus possible. Cela s’appelle la colinéarité ou encore multicolinéarité. Ainsi, il est toujours préférable de choisir un ensemble de variables indépendantes peu corrélées entre elles dans un même modèle, mais ce n’est pas toujours possible.↩︎