Cours de R pour débutant pressé : les régressions !

Bon­jour à tous et soyez les bien­ve­nus dans ce 3^ème cours de R pour débu­tant pres­sé.

Aujourd’hui, nous allons voir rapi­de­ment ce qu’est une régres­sion (linéaire ou qua­dra­tique), à quoi ça sert et ce que ça peut nous apprendre sur nos don­nées.

Ne vous êtes-vous jamais deman­dé com­ment en apprendre plus sur vos don­nées, com­ment savoir quel para­mètre est le plus impor­tant ou plus sim­ple­ment s’il est pos­sible « de faire une belle ligne sur mon gra­phique » ?

Non ? Ah… et bien mer­ci à bien­tôt.

Si ? Ah c’était une blague… très bien, conti­nuons alors.

Com­men­çons par le début et deman­dons-nous ce qu’est une régres­sion et pour­quoi on en fait. Une régres­sion c’est le fait de faire cor­res­pondre un modèle mathé­ma­tique à nos don­nées, de manière à esti­mer les para­mètres de la régres­sion, qui seront des approxi­ma­tions des para­mètres de nos don­nées réelles. Pour expli­quer tous ces mots bar­bares, pre­nons un petit exemple théo­rique :

Disons que nous avons des don­nées de mor­ta­li­té du nombre d’ewoks en fonc­tion du nombre TS-TT pré­sents sur le champ de bataille d’Endor. Nous savons par exemple que pour 1 TS-TT, 10 ewoks meurent, pour 2 TS-TT, 20 ewoks meurent et pour 3 TS-TT, 30 ewoks meurent.

Sous R ça donne ça :

> tstt = c(1,2,3)

> ewoks = c(10,20,30)

> plot(ewoks~tstt)
La relation ici est assez simple, puisque ewoks = 10tstt
Passons maintenant à quelques explications mathématiques, mais promis, rien de bien difficile. Le nombre d’Ewoks (ewoks) est la variable que l’on cherche à expliquer, que l’on appelle en général Y. Le nombre de TS-TT (tstt) est une variable explicative (qui explique la variable à expliquer), que l’on appelle en général X, ou A. Aujourd’hui, nous n’entrerons pas dans des modèles à plusieurs variables explicatives (A, B, C), donc appelons-la X. Le modèle mathématique le plus simple et qui fonctionne parfaitement ici est donc :
Y = 0 + 10X
En effet, en mathématiques la relation la plus simple entre deux variables est :
Y = α + βX
La première question qui doit vous venir à l’esprit si vous vous êtes déjà retrouvés dans l’exercice périlleux de la rédaction d’article, de thèse ou de rapport, c’est : « ok, mais quand tu dis que ça fonctionne parfaitement, comment le prouves-tu ? ». Merci, c’est une très bonne question, qui me permet d’introduire notre meilleur ami dans R : lm.
lm est une commande R du package stats (inclus de base dans R), qui calcule les modèles correspondant à nos données et nous fournit quelques détails intéressants :
> lm(ewoks~tstt)
Call:

lm(formula = ewoks ~ tstt)
Coefficients:

(Intercept)         tstt

-4.102e-15   1.000e+01
Ah ouais, là, ça ne nous dit pas grand chose… en effet, en tapant la commande 
lm(ewoks~tstt)
 , tout ce qu’on fait c’est appeler la fonction. Ce qui nous intéresse c’est plutôt un résumé de cet appel :
> summary(lm(ewoks~tstt))
Call:

lm(formula = ewoks ~ tstt)
Residuals:

1         2         3

7.252e-16 -1.450e-15 7.252e-16
Coefficients:

Estimate Std. Error   t value Pr(>|t|)

(Intercept) -4.102e-15 2.713e-15 -1.512e+00   0.372

tstt         1.000e+01 1.256e-15 7.961e+15   <2e-16 ***

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.776e-15 on 1 degrees of freedom

Multiple R-squared:     1,     Adjusted R-squared:     1

F-statistic: 6.338e+31 on 1 and 1 DF, p-value: < 2.2e-16
Et hop, on obtient beaucoup d’informations intéressantes d’un coup ! Explications :
Call:

lm(formula = ewoks ~ tstt)


Tout d’abord le « Call » nous rappelle quel modèle a été appelé. En effet, si l’on place notre fonction dans une variable, on peut ne plus se souvenir de ce qu’on avait appelé, donc ça nous est rappelé.
Residuals:

1         2         3

7.252e-16 -1.450e-15 7.252e-16


Les « Residuals » sont les différences entre les valeurs données dans notre X et celles estimées par le modèle. Ici, on est dans les environs de 10-16. Je pense que vous conviendrez que ce n’est pas beaucoup.
Coefficients:

Estimate Std. Error   t value Pr(>|t|)

(Intercept) -4.102e-15 2.713e-15 -1.512e+00 0.372

tstt         1.000e+01 1.256e-15 7.961e+15  <2e-16 ***

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Les « Coefficients » sont les valeurs que le modèle estime, leur écart-type, leur t-value et la probabilité de rejeter l’hypothèse H0 pour laquelle le coefficient est égal à 0. Dans cet exemple, il n’y a pas de preuve statistique permettant de dire que l'ordonnée à l'origine (Intercept) est différente de 0, alors que le tstt est différent de 0 avec un test très significatif (***). C’est donc ici que l’on apprend que notre modèle est bien :
ewoks = 0 + 10tstt = 10tstt
A noter : le -4.102e-15 est devenu 0 car le test statistique n'est pas significatif (Pr(>|t|) = 0.372). Le 10 vient du 1.000e+01.
Residual standard error: 1.776e-15 on 1 degrees of freedom

Multiple R-squared:      1,     Adjusted R-squared:     1

F-statistic: 6.338e+31 on 1 and 1 DF, p-value: < 2.2e-16


Enfin, ici nous avons quelques informations pour savoir si notre modèle colle bien avec nos données. La « Residual standard error » est l’écart type des résidus, avec leur nombre de degrés de liberté. Les « Multiple R-squared » et « Adjusted R-squared » sont les coefficients de corrélation, plus ils sont proches de 1 et meilleur est notre modèle. Ici le modèle semble parfait ! Pour finir la « F-statistic » qui est un peu plus compliquée (et dont je ne maîtrise pas toutes les ficelles) mais qui a une chose que tout le monde comprend : une p-value ! Elle devra suivre en gros les codes significatifs habituels.
Donc au final, ce modèle a l’air de bien coller avec nos données ;-)
On peut aussi évaluer visuellement notre modèle avec :
> layout(matrix(1:4, 2, 2))

> plot(lm(ewoks~tstt))


Ceci nous affiche 4 graphiques. Vu qu'on est pressés je ne vais pas entrer dans les détails, juste en dire quelques mots.

Ce qu’il faut savoir, c’est que les graphiques du haut (Residuals vs Fitted et Scale-Location) ne doivent pas donner de tendance claire (ne doivent pas aller tous croissants ou tous décroissants), que le normal Q-Q plot (page wikipedia) doit montrer des points répartis autour de la ligne pointillée et que les points la suivent à peu près, et enfin que le dernier ne doit pas montrer de point qui dépasse 1 en abscisse. Vous trouverez tout le détail nécessaire sur ce PDF, notamment en page 5.
Bon, en sciences on aime bien être bien bien bien sûrs des choses, alors on peut faire un AIC (page wikipedia) qui dit que plus sa valeur est basse, plus le modèle est bon et colle à nos données :
> AIC(lm(ewoks~tstt))

[1] -192.5675


-193 c’est très bon. Cependant, il est vrai qu'on utilise plutôt l'AIC en comparaison, en lui-même il n'est pas exploitable dans un article scientifique...
Introduisons ici le dernier élément dont nous allons nous servir aujourd’hui, l’ANOVA, qui va nous renseigner, quand on ne lui donne qu’un modèle, sur l’aspect significatif des variables explicatives :
> anova(lm(ewoks~tstt))

Analysis of Variance Table
Response: ewoks

Df Sum Sq Mean Sq   F value   Pr(>F)

tstt       1   200     200 6.3383e+31 < 2.2e-16 ***

Residuals 1    0       0

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Message d'avis :

In anova.lm(lm(ewoks ~ tstt)) :

les tests F d'ANOVA sur un ajustement pratiquement parfait ne sont pas fiables
On a encore une fois trois étoiles pour notre tstt, qui explique donc très bien le nombre d’ewoks morts. On a de plus quelque chose de très intéressant, le message d’avis, qui nous dit au final que nos ajustements sont pratiquement parfaits.
On peut s’en apercevoir quand on trace la ligne de régression sur notre plot précédent (je remets le code du plot, au cas où mais il n’est pas nécessaire si vous n’avez pas fermé votre fenêtre) :
> plot(ewoks~tstt)

> abline(lm(ewoks~tstt))



Super !
On complique un peu l’affaire ? Parce que bon, je ne suis pas sûr que vous rencontriez un jour des données aussi faciles. Ni d’ewoks d’ailleurs. En revanche pour les TS-TT j’ai encore de l’espoir.
Alors disons maintenant que nos données ressemblent plus à ça :
> tstt = c(10, 20, 30, 40, 50)

> ewoks = c(1000, 900, 600, 300, 100)
Où ewoks est cette fois le nombre d’ewoks survivants. Plaçons cette fois notre modèle dans une variable :
> lm.linear = lm(ewoks~tstt)
Le graphique est le suivant :
> plot(ewoks~tstt)

> abline(lm.linear)
Le résumé du modèle nous donne ceci :
> summary(lm.linear)
Call:

lm(formula = ewoks ~ tstt)
Residuals:

1         2         3         4         5

-6.000e+01 8.000e+01 2.000e+01 -4.000e+01 2.842e-14
Coefficients:

Estimate Std. Error t value Pr(>|t|)

(Intercept) 1300.00 66.33   19.6 0.00029 ***

tstt         -24.00 2.00   -12.0 0.00125 **

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 63.25 on 3 degrees of freedom

Multiple R-squared: 0.9796,           Adjusted R-squared: 0.9728

F-statistic:   144 on 1 and 3 DF, p-value: 0.001245
On voit cette fois que le modèle est :

ewoks = 1300 - 24tstt

On voit aussi que les coefficients de corrélation sont bons, que la p-value est assez faible. Si vous faites l’AIC (je vous laisse deviner comment faire), vous devez obtenir 59.10551, ce qui ne nous dit pas grand chose mais va être intéressant pour la suite.
Il est toujours intéressant de voir si un modèle plus complexe donne de meilleurs résultats, pour pouvoir justifier que l’on a eu raison de garder un modèle linéaire ou pas. Le niveau supérieur de complexité d’un modèle linéaire est un modèle dit polynomial, dans lequel on va ajouter à nos variables explicatives une puissance. Si vous voulez, c’est comme se dire que les solutions offertes dans le plan des puissances 1 ne sont pas suffisantes, alors on passe aux puissances 2. On peut aussi noter que notre α est en fait αX0, où X0 = 1, donc les solutions du plan de puissance 0. L’équation est la suivante :

Y = α + βX1 +γX2
Pour faire cela dans R, on va utiliser la commande :
> lm.polynomial = lm(ewoks~tstt+I(tstt^2))
Ici, on notera le
I()
 , qui est la fonction d’identité qui permet d’utiliser des termes dans le modèle incluant des symboles mathématiques normaux.
On peut faire un plot de la ligne de régression de ce modèle, mais vu que ce n’est pas une droite, il faut ruser un peu et calculer sa valeur à différents points. Pour cela on utilise les coefficients calculés :
> lm.polynomial
Call:

lm(formula = ewoks ~ tstt + I(tstt^2))
Coefficients:

(Intercept)         tstt     I(tstt^2)

1200.0000          -15.4286 -0.1429
> newtstt = seq(10, 50, 1)

> fit = 1200 + -15.4286 * newtstt + -0.1429 * newtstt ^2

> plot(ewoks~tstt)

> lines(newtstt, fit, lty=1)
Ça a l’air moins bien, n’est-ce pas ?

Les valeurs qu'on entre dans « fit » sont celles données en coefficients.

Je vous laisse regarder le plot du modèle, le résumé nous donne ceci :
> summary(lm.polynomial)
Call:

lm(formula = ewoks ~ tstt + I(tstt^2))
Residuals:

1       2       3       4       5

-31.429 65.714 -8.571 -54.286 28.571
Coefficients:

Estimate Std. Error t value Pr(>|t|)

(Intercept) 1200.0000   145.0123   8.275   0.0143 *

tstt         -15.4286    11.0509  -1.396   0.2975

I(tstt^2)     -0.1429     0.1807  -0.791   0.5120

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 67.61 on 2 degrees of freedom

Multiple R-squared: 0.9845,           Adjusted R-squared: 0.9689

F-statistic: 63.31 on 2 and 2 DF, p-value: 0.01555
Et on peut utiliser l’AIC et l’ANOVA pour comparer les deux modèles :
> AIC(lm.linear, lm.polynomial)

df     AIC

lm.linear     3 59.10551

lm.polynomial 4 59.74584
Les AIC sont très proches, en général on dit qu’il faut 2 de différence d’AIC pour justifier le passage à un modèle plus complexe et que 10 est une différence très forte. Ici il semble donc que le passage au modèle polynomial ne soit pas nécessaire et même qu’il est moins bon que le linéaire.
> anova(lm.linear, lm.polynomial)

Analysis of Variance Table
Model 1: ewoks ~ tstt

Model 2: ewoks ~ tstt + I(tstt^2)
Res.Df     RSS Df Sum of Sq     F Pr(>F)

1     3 12000.0

2     2 9142.9 1   2857.1 0.625 0.512
L'ANOVA nous dit que la différence entre les deux n’est pas significative (on ne peut pas rejeter l’hypothèse H0).
On a donc plusieurs tests qui prouvent qu’il vaut mieux rester au modèle linéaire, qui est d’ailleurs meilleur pour nos données.
Enfin et pour conclure, ces modèles nous permettent de prédire des données, grâce à la commande « predict ». Elle a besoin du modèle pour lequel on veut prédire des données, d’un vecteur des données à prédire et de l’intervalle de confiance désiré :
> newtstt = c(0,1,54,55)

> predict(lm.linear,data.frame(tstt = newtstt), level = 0.95)

1   2   3   4

1300 1276   4 -20
On sait donc qu’on avait au départ 1300 ewoks (pour 0 TS-TT, résultat 1), qu’un seul TS-TT en tue environ 24 (1300 - 1276, résultat 1 - résultat 2), et que c’est le 55ème TS-TT qui aura raison de ces sales petites bêtes (résultat 4) !
Evidemment, dans le monde rien n'est aussi simple et je pense que vos expériences ne le seront pas non plus. Il est toujours possible d'ajouter de la complexité et R est un outil suffisamment puissant pour vous suivre partout. On pourrait ici s'amuser à ajouter l'effet des Jedi ou des Wookies, ou encore de l'humidité de la planète sur les rouages ! ;-)

Je tiens à remercier tous mes relecteurs : m4rsu, Gildas et Estel.
Je tiens aussi à préciser qu'aucun ewok n'a été maltraité durant la rédaction de cet article.
      
          
        Vous avez aimé ? Dites-le nous !      
            
      
                  
                      
                  
                      
                  
                      
                  
                      
                  
                      
              
    
    
          Soumettre        
      Moyenne : 0 / 5. Nb de votes : 0    
    
      Pas encore de vote pour cet article.