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

Bonjour à tous et soyez les bienvenus dans ce 3ème cours de R pour débutant pressé.

Aujourd’hui, nous allons voir rapidement ce qu’est une régression (linéaire ou quadratique), à quoi ça sert et ce que ça peut nous apprendre sur nos données.

Ne vous êtes-vous jamais demandé comment en apprendre plus sur vos données, comment savoir quel paramètre est le plus important ou plus simplement s’il est possible « de faire une belle ligne sur mon graphique » ?

Non ? Ah… et bien merci à bientôt.

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

Commençons par le début et demandons-nous ce qu’est une régression et pourquoi on en fait. Une régression c’est le fait de faire correspondre un modèle mathématique à nos données, de manière à estimer les paramètres de la régression, qui seront des approximations des paramètres de nos données réelles. Pour expliquer tous ces mots barbares, prenons un petit exemple théorique :

5123321886_d1dba055ef_o

Crédits : leg0fenris - Flickr (CC BY-NC-ND 2.0)

Disons que nous avons des données de mortalité du nombre d’ewoks en fonction 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 :

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 :

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 :

Et hop, on obtient beaucoup d’informations intéressantes d’un coup ! Explications :

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é.

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.

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.

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 :

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.

Capture d’écran 2015-07-15 à 10.30.24

plot(lm(ewoks~tstt)) Crédits : Aurélien C. (CC-BY-SA 3.0)

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 :

-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 :

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) :

Capture d’écran 2015-07-01 à 14.23.44

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 :

Où ewoks est cette fois le nombre d’ewoks survivants. Plaçons cette fois notre modèle dans une variable :

Le graphique est le suivant :

Le résumé du modèle nous donne ceci :

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 :

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 :

Ç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 :

Et on peut utiliser l’AIC et l’ANOVA pour comparer les deux modèles :

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.

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é :

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.

  • À propos de
  • Passionné d'informatique et de biologie, après des études à Orléans, un stage de M2 au CSBC de l'Université de Cambridge (UK), une thèse à l'IRBI de Tours sur l'évolution des virus et un Post-doc à Toronto sur les cancers humains, je suis aujourd'hui de retour à Orléans à l'INRA, où j'essaie de lier phénotype des arbres et expression de leurs gènes. Je suis très intéressé par les approches d'apprentissage machine. Mes autres passions sont le basket-ball, le krav maga, la musique ou encore les jeux vidéos.

Laisser un commentaire