Modélisation théorique et pratique (cours 1/4) : introduction à la modélisation

Bonjour. Ce cours de modélisation se divise en 2 parties. Dans une première partie, nous allons voir comment construire et sélectionner un modèle. Dans une seconde partie, nous verrons comment interpréter les résultats d'un modèle et comment utiliser un modèle pour faire de la prévision. Mais tout d'abord, qu'est-ce que la modélisation? On peut résumer succinctement en disant que modéliser, c'est d'une part comprendre et d'autre part prévoir. Et qui a déjà fait de la modélisation? À mon avis, chacun de nous, et même en construisant des modèles complexes, mais le plus souvent sans le savoir. En effet, on va prévoir un temps de trajet en fonction du jour de la semaine, de l'heure de départ, du fait qu'il pleuve ou non, que c'est un jour de départ en vacances, etc. On peut aussi chercher à estimer le temps d'attente à un guichet en fonction du nombre de personnes, du nombre de guichets ouverts, du temps mis par la première personne qu'on a vue au guichet, etc. On peut estimer le prix de location d'un appartement en fonction de sa superficie, du nombre de chambres, de sa localisation, de son exposition, de la présence d'espaces verts et d'autres critères qui conduisent à déterminer le prix d'un appartement. Dans ces 3 premiers exemples, on cherche à estimer une variable quantitative. On peut aussi chercher à prévoir une variable qualitative. Comme lorsqu'on décide comment s'habiller le matin, en prenant ou pas un vêtement de pluie. Pour prédire, on utilise là encore des informations comme la pluviométrie des jours précédents, la présence de nuages dans le ciel, la présence de vent, etc. Cette fois, la variable va-t-il pleuvoir que l'on cherche à prédire est qualitative puisque 2 réponses sont possibles: il pleut ou il ne pleut pas. Ce type de modélisation ne sera pas abordé dans ce cours, même si beaucoup d'idées peuvent être reprises puisqu'on essaie de comprendre l'effet de plusieurs variables, sur une variable réponse, mais ici la variable réponse est qualitative. Et comment faites-vous pour prévoir? Vous listez toutes les variables qui potentiellement peuvent influer sur la variable réponse, c'est-à-dire la variable que vous cherchez à comprendre. Une fois que vous avez l'ensemble des variables, vous ignorez les variables qui sont négligeables et qui n'aident pas à mieux prévoir votre réponse. Et vous essayez ensuite de quantifier au mieux l'effet des variables restantes qui ont été sélectionnées. Et alors, est-ce que cette stratégie est bien raisonnable? Eh bien, en fait, oui, tout à fait. Et c'est exactement la démarche que l'on va adopter en modélisation. Mais alors, on peut se demander à quoi servent les statistiques? Eh bien, elles servent à faire cela avec rigueur, pour comprendre et prédire des variables qui mesurent des phénomènes parfois beaucoup plus complexes, mais l'idée de cette démarche est à garder pour la modélisation. Voici quelques exemples de modèles en agronomie ou alimentation, mais on trouve de tels modèles dans de très nombreux domaines comme l'économie, la démographie, la climatologie, etc. Le premier exemple: on cherche à prédire le temps de cuisson idéal en fonction de la composition et du poids d'un aliment, de la température du four, de l'humidité de l'air, etc. Deuxième exemple: on veut prévoir la production de biogaz en fonction de la quantité de déchets agricoles ou alimentaires, des résidus de culture, d'ordures ménagères, etc. Troisième exemple: on veut prévoir la production électrique d'une éolienne en fonction de la vitesse du vent à 10 m, de la vitesse du vent à 80 m, de la température à 2 m, de la pression atmosphérique, de l'humidité relative à 2 m, de la direction du vent. Et dernier exemple: on veut optimiser le couple durée-température pour maximiser un rendement. On va alors construire un modèle du rendement en fonction du temps et de la température et une fois le modèle connu, on pourra prédire le rendement pour n'importe quel couple et ensuite choisir le couple qui optimise le rendement. Dans tous ces exemples de modélisation, les objectifs sont de comprendre quelles variables influent sur une variable réponse quantitative et également de prévoir les valeurs de la variable réponse pour de nouvelles conditions. Nous allons illustrer ce cours avec un exemple sur la prévision de la qualité de l'air. Plus précisément, nous allons chercher à prédire le maximum d'ozone atteint lors d'une journée, et ce en fonction de diverses variables météorologiques comme la température à 9 heures, à 12 heures, à 15 heures, la nébulosité à 9 heures, 12 heures et 15 heures, la vitesse du vent calculée sur un axe ouest-est. Si le vent vient de l'ouest, la valeur sera négative. Si le vent vient de l'est, la valeur sera positive. Toujours mesurées à 9h, 12h et 15h, la valeur du maximum d'ozone de la veille, la direction du vent avec 4 modalités possibles: nord, est, sud et ouest, ou la pluviométrie considérée comme une variable qualitative à 2 modalités: temps sec ou pluvieux. L'objectif de cette étude est alors de construire un modèle pour comprendre d'une part ce qui peut influer sur le maximum d'ozone et et ensuite de construire un modèle qui permettra de prédire le maximum d'ozone pour de nouvelles données météo. Les lignes du jeu de données correspondent à des jours de mesures, à partir du 1er juin 2001, et nous avons 112 relevés, et ces données ont été relevées à Rennes. Les données sont présentes dans ce fichier disponible sur internet. On peut importer le fichier dans le logiciel R en utilisant la fonction read.table en précisant le lien vers le fichier que la première ligne du jeu de données contient le nom des variables avec header=TRUE et en précisant que les variables textuelles sont en réalité des facteurs, c'est-à-dire des variables qualitatives grâce à l'argument stringAsFactor=TRUE. Avant de se lancer dans la modélisation, il est toujours indispensable de visualiser les données. Visualiser les données permet d'une part de s'assurer que nous n'avons pas de données aberrantes et d'autre part ça permet d'avoir une idée des liaisons entre les variables et notamment entre la variable réponse et les autres variables. Ce graphique représente le maximum d'ozone en fonction de la température à 9 heures. Il permet de voir par exemple si la liaison est plutôt linéaire, c'est-à-dire si les points sont plus ou moins autour d'une droite, ou si la liaison n'est pas linéaire, par exemple une liaison quadratique ou exponentielle, ou encore s'il n'y a pas de liaison entre ces 2 variables. Évidemment, si on visualise une liaison plutôt quadratique, on ne pourra pas bien prédire avec des effets linéaires uniquement. Les lignes de code à côté du graphe permettent de reproduire ce graphique. On utilise le package ggplot2 pour construire le graphique et donc la fonction ggplot appliquée sur le jeu de données ozone. Ensuite, on précise grâce à la fonction aes la variable qui sera sur l'axe des x, ici T9, et celle sur l'axe des y, ici le maximum d'ozone. Geom_point permet de dire qu'on veut dessiner les points et geom_smooth qu'on veut dessiner une courbe et ici comme method= lm on dessine une droite. lm veut dire linear model. Et se= FALSE indique qu'on ne dessine pas la zone de confiance autour de la droite. Et enfin ggtitle permet de préciser le titre du graphique. Il est possible d'écrire un modèle de régression simple entre le maximum d'ozone et la température à 9 heures. Plusieurs écritures sont possibles, mais en toute rigueur, on écrit Yi= bêta 0+ bêta 1 Xi+ epsilon i. L'indice i correspond au numéro de l'individu statistique. Dans notre exemple, Yi sera le maximum d'ozone le jour i. Et on écrit que cette valeur est égale à une constante bêta 0 plus bêta 1 fois la température à 9 heures ce même jour i. Et plus un résidu epsilon_i. On considère alors que la relation est la même pour tous les individus avec les mêmes bêta_0 et bêta_1. Des hypothèses sont faites sur les résidus et ce sont les hypothèses classiques en statistiques, à savoir que les résidus ont une espérance nulle, ils ont tous la même variance sigma carré et ils sont indépendants 2 à 2. Pour généraliser ce modèle qui n'a qu'une seule variable explicative en un modèle de régression avec plusieurs variables explicatives, on va simplement étendre l'écriture en ajoutant les effets de variables W, Z, etc. Chaque variable est multipliée par un nouveau coefficient bêta. On aura donc bêta 0 plus bêta 1 Xi plus bêta 2 Wi plus etc. plus bêta P. Zi, et on ajoute toujours à la fin un résidu epsilon_i. En effet, même si on ajoute beaucoup de variables, il reste un epsilon_i, c'est-à-dire une part de variabilité que l'on n'explique pas. Si le nombre de variables est grand, on va plutôt écrire les variables X en les numérotant xi_1, qui sera la ième valeur de la variable explicative X_1, Xi_2 la ième valeur de la variable explicative X_2, etc., jusqu'à Xi_p la ième valeur de la variable explicative X_p. Pour les résidus epsilon_i, les hypothèses sont à nouveau les hypothèses usuelles, à savoir qu'ils ont une espérance nulle, la même variance sigma carré et sont indépendants 2 à 2. Voilà comment les effets des variables explicatives quantitatives seront pris en compte dans le modèle. Si on s'intéresse maintenant à la liaison entre une variable quantitative, toujours le maximum d'ozone, et cette fois une variable qualitative comme la direction du vent, on représentera les données avec une boîte à moustache pour chaque modalité de la variable qualitative. Ce graphe permet à nouveau de repérer d'éventuelles données aberrantes et donne une idée des différences de distribution qu'il peut y avoir pour la variable réponse lorsque les données proviennent de différentes modalités pour la variable qualitative. Là encore, on utilise le package ggplot2 et la fonction ggplot que l'on applique sur notre jeu de données. On précise ensuite la variable x qualitative, ici le vent, la variable y quantitative, ici le maximum d'ozone, et on précise avec fill qu'on colorie les boîtes en fonction de la variable vent, et avec col on précise qu'on colorie les points en fonction de la variable vent également. geom_boxplot indique que l'on dessine des boîtes à moustaches en supprimant les points qui sont extrêmes avec outlier.shape NA. On supprime les points qui sont extrêmes, ceux qui sont en dehors de la boîte à moustaches, car on dessine ensuite tous les points y compris ceux qui sont en dehors de la boîte à moustaches. Et on colorie la boîte avec une certaine transparence en mettant alpha= 0.4. Il faut mettre pour la transparence une valeur entre 0 et 1 pour avoir plus ou moins de transparence. Plutôt que geom_point, on utilise ici geom_jitter qui permet de bouger légèrement les abscisses des points, ce qui évite que tous les points correspondant à une même direction de vent aient la même abscisse et soient plus ou moins superposés. Enfin, on ajoute un titre au graphique. Alors, comment s'écrit le modèle? Simplement en écrivant que le maximum d'ozone dépend de la direction du vent et plus précisément que le maximum d'ozone du jour j pour une direction de vent i, c'est un maximum d'ozone moyen μ plus un effet alpha_1 si le vent vient de l'est, alpha_2 si le vent vient du nord, alpha_3 si le vent vient de l'ouest et alpha_4 s'il vient du sud. Donc selon la direction du vent, on ajoute la valeur spécifique qui convient. Et enfin, on a toujours un résidu, donc quelque chose que l'on n'explique pas. Mathématiquement, on écrira que la réponse y pour le jème individu qui prend la modalité i pour la variable qualitative est égal à une valeur moyenne mu plus une valeur alpha_i qui dépend donc de la modalité de la variable qualitative. Par exemple, si le vent vient de l'est, le maximum d'ozone sera de +4 par rapport au maximum d'ozone moyen, ou -2 si le vent vient du nord, etc. Et nous avons ensuite le résidu epsilon_ij qui correspond à l'écart entre notre réponse y_ij et ce que prévoit le modèle, à savoir mu plus alpha_i. Et sur ces résidus, on fait à nouveau les hypothèses usuelles, c'est-à-dire que les résidus sont d'espérance nulle, de même variance sigma carré, et les résidus sont indépendants 2 à 2. Si maintenant on étend ce modèle d'indice de variance à 1 facteur à un modèle d'indice de variance à 2 facteurs, il suffit d'ajouter un effet bêta_j pour le 2ème facteur. On utilise alors un indice k supplémentaire pour les répétitions. Et on écrit que la réponse Y pour le k-ième individu qui prend à la fois la modalité i pour le premier facteur et la modalité j pour le deuxième facteur est égale à un effet moyen mu plus un effet alpha_i qui dépend de la modalité i que prend l'individu pour le premier facteur plus un effet bêta_j qui dépend de la modalité j que prend l'individu pour le deuxième facteur et on a toujours un résidu epsilon_ijk. Ici, par exemple, pour le k-ième jour de vent d'est où il pleut, on aura un effet moyen plus un effet alpha est plus un effet bêta pluie plus un résidu. Et à nouveau, on a les hypothèses classiques pour les résidus. Et il n'y a aucune difficulté à avoir plus de variables explicatives. On va simplement ajouter les effets de chaque variable explicative. Donc nous avons vu pour l'instant les effets de variables quantitatives sur la réponse ou les effets de variables qualitatives sur la réponse. Cependant, nous avons pour l'instant considéré que les effets de chaque variable ne font que s'additionner. Voyons maintenant des effets conjoints de 2 variables sur la réponse. Et commençons par voir l'effet conjoint d'une variable quantitative et d'une variable qualitative sur la variable réponse. On va construire un graphique avec la variable réponse en ordonnée et en abscisse on met la variable quantitative, ici la température à 9 heures. On va alors simplement colorier les points en fonction de la variable qualitative, ici la direction du vent. Et on peut dessiner une droite de régression par modalité de la variable vent, c'est-à-dire une droite de régression par sous-nuage de points. L'intérêt est alors de voir si l'effet de la température à 9 heures sur le maximum d'ozone est le même pour toutes les directions du vent, ou bien si l'effet de la température à 9 heures sur le maximum d'ozone dépend de la direction du vent. Si les pentes des 4 droites étaient identiques, donc si les droites étaient parallèles, alors cela signifierait que l'effet de la température à 9 heures sur le maximum d'ozone est le même quelle que soit la direction du vent. Autrement dit, 1 degré de plus à 9 heures conduit à la même augmentation ou la même diminution du maximum d'ozone, et cela quelle que soit la direction du vent. Ça, c'est si les droites sont parallèles. Ici, les droites ne sont pas parfaitement parallèles et on a l'impression que l'effet de la température à 9 heures est moins important quand le vent vient de l'est par rapport aux 3 autres directions du vent. L'augmentation du maximum d'ozone est moindre pour une augmentation de la température de 1 degré à 9 heures puisque la pente est plus faible. Cependant, ce graphique permet de visualiser l'effet conjoint de la température et de la direction du vent sur le maximum d'ozone, mais nous aurons besoin d'un test statistique pour décider si l'effet conjoint est significatif ou non. On parle d'interaction de la température et du vent et on testera si l'effet de l'interaction est significatif ou non. Pour construire le graphique, on précise qu'en x on met la variable de vent, en y le maximum d'ozone, on colorie les individus avec une couleur différente pour chaque direction du vent et on groupe les individus par direction du vent, ce qui permet ensuite de faire des régressions par direction du vent. geom_smooth permet de dessiner la droite de régression, la droite de régression car on met method= lm. On met sce= FALSE pour dire qu'on ne veut pas dessiner de zone de confiance. L'écriture du modèle donne la variable réponse, le maximum d'ozone, et précise que ce maximum d'ozone dépend de la direction du vent, de la température à 9 heures et de l'effet conjoint du vent et de la température à 9 heures. Cela signifie qu'on n'a pas seulement l'effet du vent et l'effet de la température à 9 heures, mais bien un effet conjoint en plus. Plus précisément, le maximum d'ozone s'écrit comme un effet moyen plus un effet lié à la direction du vent, c'est-à-dire une constante spécifique selon la direction du vent alpha 1, alpha 2, alpha 3 ou alpha 4. Ensuite, on ajoute un effet linéaire de la température à 9 heures en précisant qu'il y a un effet moyen bêta de la température à 9 heures plus un effet spécifique gamma_i qui dépend de la direction du vent. Autrement dit, la pente de la droite de régression est égale à bêta en moyenne, mais selon la direction du vent, on ajoute à cette pente une valeur gamma_i qui est positive ou négative. Ainsi, la pente de la droite de régression est plus ou moins forte et donc l'impact de l'augmentation de 1°C de la température à 9 heures sur le maximum d'ozone sera plus ou moins important selon la direction du vent. Mathématiquement, le modèle s'écrit Yij, la réponse pour le jème individu de la modalité i pour la variable qualitative, est égale à un effet moyen mu, un effet alpha_i lié à la direction du vent pour notre exemple, puis une pente moyenne bêta plus un écart de pente gamma_i positif ou négatif selon que la pente est plus ou moins forte que la moyenne. Et cette pente est multipliée par Xij, c'est-à-dire la valeur du jème individu qui prend la modalité i. Pour nous, ce sera la température à 9 heures du jème jour quand la direction du vent est i. Et on a en plus les résidus avec les hypothèses usuelles sur les résidus. Visualisons maintenant l'effet conjoint de 2 variables qualitatives sur la variable réponse. On va parler de l'effet de l'interaction de 2 variables qualitatives sur la variable réponse. Notez bien que, dans le dictionnaire, la définition courante de l'interaction est la réaction réciproque de 2 phénomènes l'un sur l'autre. 2 phénomènes ou variables seulement entrent en jeu et ont chacun un impact sur l'autre. En statistiques, l'interaction met en jeu 3 variables et on s'intéresse à l'effet conjoint de 2 variables sur une 3ème variable, la variable réponse. Donc on va construire un graphique avec la variable réponse en ordonnée, en abscisse une des variables qualitatives et on construit une ligne brisée pour chaque modalité de l'autre variable qualitative. Le choix de la variable à mettre en abscisse et de la variable utilisée pour construire les lignes brisées peut être interverti. On visualise alors de façon un peu différente la même interaction Et on peut donc faire les 2 graphiques et choisir celui qui est le plus facile à interpréter. Souvent, on préférera mettre en abscisse la variable qualitative qui a le plus de modalités pour avoir moins de lignes brisées. Comment sont construites les lignes brisées? On va d'abord positionner les points avec un point qui a pour ordonnée la moyenne des valeurs de Y pour une combinaison des 2 variables qualitatives. Par exemple, Pour le couple vent d'est avec pluie, on a ce point, et pour le couple vent d'est sec, on aura ce point. On positionne alors autant de points qu'il y a de combinaisons de modalités des 2 variables qualitatives, puis on relie les points d'une même modalité de la variable qualitative pour constituer une ligne brisée. Si les lignes brisées sont parallèles, alors cela signifie qu'il n'y a pas d'effet de l'interaction sur la variable réponse. Et au contraire, s'il y a des écarts au parallélisme, cela signifie qu'une interaction est présente. Par exemple ici, l'écart de maximum d'ozone entre avec et sans pluie est beaucoup plus important quand le vent vient de l'est par rapport aux écarts quand le vent vient d'une autre direction. Cela signifie que lorsque le vent vient de l'est, l'effet de la pluie ou du temps sec est exacerbé et le maximum d'ozone est particulièrement élevé quand le temps est sec ou particulièrement faible quand le temps est pluvieux. Pour construire ce graphique avec ggplot2, on utilise le package dplyr qui permet de grouper les données par combinaisons de vent et pluie avec la fonction group_by. Ensuite, on calcule la moyenne du maximum d'ozone pour chaque combinaison avec la fonction summarize Puis, pour construire le graphique, on place en abscisse le vent et en ordonnée le maximum d'ozone, on colorie en fonction de la pluie et on groupe en fonction de la pluie, ce qui permettra d'avoir une ligne brisée par modalité de la variable pluie. On utilise ensuite geom_line pour tracer les lignes et geom_point pour mettre les points. Et on ajoute un titre et les libellés des axes. Pour écrire le modèle, on dit simplement que le maximum d'ozone est expliqué en fonction du vent, de la pluie et de l'interaction vent-pluie symbolisée par les 2 points. Plus dans le détail, on écrit que le maximum d'ozone est égal à un effet moyen mu plus un effet alpha_i selon la direction du vent. Donc si le vent vient de l'est, on ajoute alpha_1, s'il vient du nord alpha_2, etc. C'est donc l'ajout d'une constante spécifique selon la direction du vent. Même chose ensuite avec la pluie, on ajoute bêta_1 s'il pleut et et bêta 2 s'il fait sec. Enfin, pour l'interaction, on ajoute une valeur spécifique pour le couple de modalités vent/pluie qui est présent. Si le vent vient de l'est et qu'il pleut, on ajoute alpha bêta 1 1. Et on aura une valeur différente alpha bêta ij pour chaque couple. Et bien entendu, il reste un résidu. Quand on dit qu'on ajoute alpha bêta 1 1, on ajoute alpha bêta 1 1 en plus de alpha 1 et de bêta 1 qui sont alpha 1 l'effet spécifique quand le vent vient de l'est, bêta 1 l'effet spécifique quand il pleut, mais alpha bêta 1 1 c'est une particularité quand il pleut et que le vent vient de l'est. Donc plus formellement, le modèle s'écrit y i j k la réponse pour le k-ième individu qui prend la modalité i pour la première variable qualitative et la modalité j pour la deuxième variable qualitative. Ce Yij est égal à un effet moyen mu, un effet alpha_i lié à la direction du vent pour notre exemple, un effet bêta_j lié à la pluviométrie, puis un effet alpha bêta_ij lié au couple direction du vent i, pluviométrie j. Et on a en plus les résidus avec les hypothèses usuelles pour les résidus. Nous avons vu 4 types d'effets possibles sur une variable réponse. Ces 4 types d'effets sont les plus couramment utilisés et permettent de construire beaucoup de modèles. Nous retrouvons donc l'effet linéaire d'une variable quantitative sur la variable réponse, l'effet d'une variable qualitative à plusieurs modalités, l'effet joint d'une variable quantitative et d'une variable qualitative sur la variable réponse, Cet effet joint permet de différencier l'effet linéaire de la variable quantitative sur la réponse selon les modalités de la variable qualitative. Et enfin l'effet de l'interaction de 2 variables qualitatives sur la variable réponse. Nous pouvons alors construire des modèles plus ou moins complexes en combinant ces 4 types d'effets selon bien sûr la nature des variables explicatives mais aussi les effets que l'on souhaite inclure dans le modèle. On peut ne pas vouloir mettre certaines interactions dans le modèle si on sait a priori qu'il n'y a pas d'effet de ces interactions. Et si on n'a pas de connaissances a priori, on mettra l'interaction dans le modèle et on prendra une décision sur l'effet ou non de l'interaction grâce au test statistique qui sera mis en place. Ces 4 types d'effets et leur interprétation sont très importants à comprendre Car toute la difficulté réside dans l'interprétation des résultats, les calculs étant faits par l'ordinateur assez simplement. Donc bien comprendre ce que représente chacun de ces effets et comment le visualiser et l'interpréter est vraiment crucial. Voici un tableau récapitulatif avec les noms des modèles selon la nature des variables explicatives. Dans tout le cours, la variable réponse est quantitative. Quand il y a une variable explicative et que celle-ci est quantitative, il s'agit d'une régression linéaire simple. Si la variable explicative est qualitative, il s'agit d'une analyse de variance à 1 facteur. Dans le cas particulier où il n'y a que 2 modalités, l'analyse de variance à 1 facteur équivaut à faire une comparaison de moyennes avec variances égales. Si maintenant il y a plusieurs variables explicatives et que toutes sont quantitatives, on parle de régression linéaire multiple. Si toutes les variables explicatives sont qualitatives, il s'agit d'une analyse de variance à k facteurs. Enfin, si certaines variables explicatives sont quantitatives et d'autres sont qualitatives, on parle d'analyse de covariance. À chaque fois, la variable réponse est quantitative et tous ces modèles sont des modèles linéaires. Si la variable réponse est qualitative, le plus souvent binaire, malade versus non malade, on peut représenter les mêmes 4 effets pour les variables explicatives selon la nature des variables. Et la présence ou non d'interaction. Mais la modélisation est différente et on parlera de régression logistique, une méthode qui ne fait pas l'objet de ce cours. La régression logistique s'intéresse à modéliser la probabilité de survenue ou non d'un événement, ou plus précisément modélise le logarithme du ratio de la probabilité de survenue sur la probabilité de non-survenue d'un événement. Voici pour une introduction générale sur la modélisation. Nous verrons ensuite comment construire un modèle, comment sélectionner les variables d'un modèle, puis comment interpréter les résultats et comment utiliser le modèle pour prédire de nouvelles valeurs de Y.