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é à un objet lm.

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

# ✓ Default language: fr

trees <- read("trees", package = "datasets")
trees_lm <- lm(data = trees, volume ~ diameter)
summary(trees_lm)

# 
# 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’erreur type des résidus, considérant que les degrés de liberté du modèle est le 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.

\[\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é R² de manière générale ou r² 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 R² quantifie la capacité du modèle à prédire la valeur de \(y\) connaissant la valeur \(x\) pour le même 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) peut elle-même être décomposée en une fraction expliquée par notre modèle, et la fraction qui ne l’est pas (les 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 R² (ou r²) se définit comme :

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

La valeur du R² 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 cas, tous les résidus valent zéro). Donc, plus le coefficient R² se rapproche de 1, plus le modèle explique bien les données et aura un bon pouvoir de prédiction.

Dans R, le R² multiple se réfère simplement au R² (ou au r² pour les régressions linéaires simples) calculé de cette façon. L’adjectif multiple indique simplement 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_i\). Les équations indiquées plus haut sont valables lorsque l’ordonnée à l’origine n’est pas figée (\(y = a \ x + b\)). 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 = a \ x\) (avec \(b = 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 R² 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 R² obtenus avec et sans forçage à zéro de l’ordonnée à l’origine !

Adjusted R-squared :

La valeur du coefficient R² 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 R². A la limite, si nous ajoutons une nouvelle variable fortement corrélée avec les précédentes⁴, l’apport en terme d’information nouvelle sera négligeable, mais le R² augmentera malgré tout un tout petit peu. Alors dans quel cas l’ajout d’une nouvelle variable est-il pertinent ? Le R² 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 R² 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, le test de 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é à 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 H₀ 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 associée au 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.

À vous de jouer !

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 aberrant, 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)

# 
# Call:
# lm(formula = volume ~ diameter, data = trees, subset = diameter < 
#     0.5)
# 
# Residuals:
#       Min        1Q    Median        3Q       Max 
# -0.215129 -0.068502 -0.001149  0.070522  0.181398 
# 
# Coefficients:
#             Estimate Std. Error t value Pr(>|t|)    
# (Intercept) -0.94445    0.09309  -10.15 6.98e-11 ***
# diameter     5.31219    0.27540   19.29  < 2e-16 ***
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 
# Residual standard error: 0.1082 on 28 degrees of freedom
# Multiple R-squared:   0.93,   Adjusted R-squared:  0.9275 
# F-statistic: 372.1 on 1 and 28 DF,  p-value: < 2.2e-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 fiablement des 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 !

La corrélation entre les prédicteurs dans un modèle linéaire multiple est un gros problème et doit être évité 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.↩︎