Modélisation théorique et pratique (cours 2/4) : construction de modèle et sélection de variables

Après avoir vu en introduction les enjeux de la modélisation et les différents effets que l'on pouvait introduire dans un modèle, reprenons maintenant l'écriture du modèle et les choix de modèles. Donc reprenons l'écriture du modèle de régression linéaire multiple. En rouge apparaissent les paramètres du modèle qui doivent être estimés à partir des données. Le nombre de paramètres estimés est égal au nombre de variables explicatives+ 1 car il faut estimer la constante. Pour estimer tous ces paramètres, nous réécrivons le modèle sous une forme matricielle Y= X× β+ e avec Y le vecteur des valeurs observées, X la matrice qui contient toutes les valeurs des variables explicatives, plus précisément une première colonne avec uniquement des 1 qui permettra d'estimer la constante bêta 0, puis une colonne pour chaque variable explicative. Beta est le vecteur qui regroupe tous les coefficients du modèle à estimer et e est le vecteur des résidus. En plus de l'équation Y= X bêta+ e, on retrouve les hypothèses sur les résidus, à savoir qu'ils sont sans biais, ce qui est donné par l'espérance de e qui vaut 0, de même variance sigma carré et indépendants 2 à 2, ce qui est donné par la variance de e qui est une matrice de variance-covariance qui a sur la diagonale les variances et hors diagonale les covariances entre 2 résidus. Comme on a sigma carré fois la matrice identité, on retrouve bien sur la diagonale que tous les résidus ont la même variance sigma carré et hors diagonale on retrouve que les résidus sont indépendants 2 à 2 puisque la covariance est égale à 0. Si on lit le modèle écrit sous la forme matricielle, on lit pour la première ligne que y1 est égal à 1 fois bêta0 plus x11 fois bêta1 plus etc. plus x1p fois bêta p plus epsilon1. En lisant cela, on retrouve l'équation donnée par l'écriture du modèle sous forme indicée pour i 1. C'est un peu plus difficile à voir mais les modèles d'analyse de variance et d'analyse de covariance peuvent aussi s'écrire sous la forme y X bêta E. Et c'est parce que c'est la même écriture qu'on parle de modèle linéaire quelle que soit la nature des variables explicatives et que les fonctions des logiciels permettant d'estimer les coefficients des modèles sont les mêmes. Pour les variables qualitatives, le codage nécessite d'utiliser plusieurs colonnes pour une même variable. Plus précisément autant de colonnes que de modalités. Et pour les interactions entre variables qualitatives, on aura besoin d'autant de colonnes que de couples de modalités. Nous verrons plus tard dans la prochaine vidéo que des contraintes seront nécessaires pour bien estimer tous les paramètres et qu'il est possible d'utiliser I-1 colonnes pour une variable à I modalités et colonnes pour les interactions de 2 variables à I et J modalités. Pour estimer les paramètres, nous utilisons le critère des moindres carrés qui consiste à minimiser la somme des des carrés des écarts entre les valeurs observées pour la variable Y et les prévisions du modèle. Pour donner l'idée, on cherche à estimer le vecteur de paramètres bêta tel que X bêta soit le plus proche possible de Y. Si on multiplie à gauche par la transposée de X, on a X prime X bêta est le plus proche possible de X prime Y, avec X prime la transposée de X. Et si la matrice X'X est inversible, alors on peut dire qu'on estime bêta par bêta chapeau égal X'X moins 1 X'Y. Cet estimateur est sans biais et la variance de cet estimateur est connue et dépend de X'X et donc simplement de la matrice X, c'est-à-dire du choix des expériences que l'on effectue. Ceci est très important, notamment si on peut planifier les expériences. Dans ce cas, on choisira des expériences telles que X' prenne des valeurs les plus faibles possibles. La variance de bêta chapeau dépend aussi de sigma carré et donc de la variance des résidus, c'est-à-dire la variance de ce que le modèle n'explique pas. Comment estimer cette variance des résidus justement? Eh bien en calculant la somme des carrés des résidus et en divisant par les degrés de liberté associés. Les 2 degrés de liberté de la résiduelle se calculent toujours de la façon suivante: c'est le nombre de données moins le nombre de paramètres estimés à partir des données. Donc pour la régression, ce sera n- car n est le nombre de données et p+ 1 est le nombre de paramètres que l'on estime. Une fois le modèle construit, une question légitime est de se demander à quel point le modèle explique la variable que l'on cherche à étudier. Les valeurs des Y observées varient et cette variabilité des valeurs de Y peut se décomposer en une somme de la variabilité expliquée par le modèle et de la variabilité non expliquée par le modèle, c'est-à-dire la variabilité résiduelle. Ainsi, la variabilité totale de Y, qui est égale à la variance de Y multipliée par n- 1, est bien expliquée par le modèle si la variabilité résiduelle est petite. Et donc si les prédictions du modèle sont proches des valeurs observées. On calcule alors le R², appelé le coefficient de détermination, qui est la part de variabilité totale de la réponse Y qui est expliquée par le modèle. Ce coefficient R² varie entre 0, le modèle n'est pas bon du tout, à 1, le modèle prévoit parfaitement les valeurs observées. Si la part de variabilité expliquée par le modèle est grande, le modèle s'ajuste bien aux données observées. Attention toutefois, un modèle peut bien s'ajuster à des données, c'est-à-dire coller aux données observées, mais ne pas être un bon modèle pour prédire de nouvelles données. Cela arrive notamment quand le modèle contient beaucoup de variables et nécessite donc l'estimation de nombreux paramètres. Plus le nombre de paramètres estimés augmente, plus il est facile de s'ajuster aux valeurs observées. Mais pour autant le modèle n'est pas nécessairement bon pour prédire de nouvelles valeurs. Les modèles qui prédisent bien de nouvelles valeurs sont des modèles qui s'ajustent bien et qui ont nécessité l'estimation d'un nombre raisonnable de paramètres. Cette variabilité expliquée par le modèle peut être décomposée en parts de variabilité expliquées par chacune des variables explicatives et interactions. Cependant il y a une difficulté ici. Car 2 variables explicatives, ou plus de 2 d'ailleurs, peuvent expliquer à peu près la même chose et donc à peu près la même part de variabilité de Y. C'est le cas si 2 variables quantitatives sont très corrélées par exemple. À quelle variable attribuer alors cette variabilité sachant que, lorsque le modèle contient les 2 variables, la variabilité n'est expliquée qu'une seule fois? Il y a alors 2 façons de décomposer la variabilité expliquée par le modèle. La première façon est de calculer la part de variabilité expliquée par la première variable du modèle, puis la part de variabilité expliquée par la deuxième variable et qui n'a pas encore été expliquée par la première, puis la part expliquée par la troisième variable non encore expliquée par les deux premières, etc. L'avantage de cette décomposition est que la variabilité expliquée par l'ensemble du modèle est égale à la somme des variabilités de chaque variable. L'inconvénient majeur est que l'ordre d'apparition des variables dans le modèle modifie la part de variabilité expliquée par une variable. Une autre façon de décomposer la variabilité du modèle est de calculer la part de variabilité expliquée exclusivement par une variable. Ainsi, l'ordre d'apparition des variables dans le modèle ne modifie pas la variabilité expliquée par la variable. En revanche, si une part de variabilité est expliquée par plusieurs variables, s'il y a une variabilité commune expliquée par plusieurs variables, alors cette variabilité commune n'apparaît pas dans la décomposition. La somme des variabilités de toutes les variables sera alors inférieure à la variabilité du modèle. Une variabilité commune est expliquée par plusieurs variables quand les variables explicatives ne sont pas indépendantes. Quand les données sont équilibrées, c'est-à-dire quand les variables explicatives sont indépendantes 2 à 2, La variabilité du modèle se décompose parfaitement et la somme des variabilités expliquées par chaque variable est égale à la variabilité expliquée par le modèle. L'interprétation des résultats est alors beaucoup plus facile. Pour construire le modèle et obtenir la variabilité expliquée par chaque variable, on utilise la fonction linearModel du package FactoMineR. On écrit alors la formule du modèle avec la variable à expliquer, le maximum d'ozone, puis tilde, et enfin les variables explicatives séparées par un plus. Donc on écrit par exemple maxO3~ T9+ T12+ T15+ vent+ pluie en mettant toutes les variables quantitatives et qualitatives. Quand le nombre de variables explicatives est grand, on peut écrire~. Car cela indique que l'on utilise toutes les variables du jeu de données comme variables explicatives, sauf bien entendu le maximum d'ozone qui est la variable réponse. Les sorties du modèle donnent alors l'écart-type résiduel qui vaut 14.51 et a 97 degrés de liberté. Le R² du modèle vaut 0.7686, ce qui signifie que 76,86% de la variabilité du maximum d'ozone est expliqué par ce modèle. Ensuite, un test est construit pour tester la significativité du R². Nous reviendrons sur ce test plus tard. Et ensuite nous avons un tableau avec en ligne chaque variable et dans la colonne SS, pour sum of squares, les sommes de carrés, ce qui correspond à la variabilité expliquée exclusivement par chaque variable. Là encore, nous avons des tests que nous allons décrire ultérieurement. Le modèle ainsi construit contient toutes les variables que nous avons choisies comme variables explicatives et toutes les interactions que nous avons pu considérer. On peut alors se demander si tous ces effets sont utiles ou si certains sont superflus. On peut aussi se demander si aucune variable n'est utile, auquel cas les variables explicatives du modèle seraient très mal choisies. Pour cela, des tests existent et ils sont tous construits selon le même principe. La question que l'on se pose est alors: Est-ce que cet ensemble de V variables apporte des informations complémentaires intéressantes sachant que les autres variables sont dans le modèle? Notez bien ici que l'on teste l'intérêt d'un ensemble de variables sachant que d'autres variables sont déjà présentes dans le modèle, ce qui revient à tester si les variables de l'ensemble V apportent un plus. Mathématiquement, on va tester une hypothèse H0 contre une hypothèse alternative H1. L'hypothèse H0 est tous les coefficients associés aux variables de l'ensemble V sont nuls et l'hypothèse alternative est que l'un au moins des coefficients n'est pas nul. La statistique de test utilisée pour ce test est le ratio entre le carré moyen de l'ensemble V et le carré moyen de la résiduelle. Le carré moyen d'un effet ou d'un ensemble d'effets est simplement la variabilité expliquée par cet effet ou par cet ensemble d'effets divisé par les degrés de liberté de cet effet ou de cet ensemble d'effets. Dans cette statistique de test, le carré moyen résiduel sert de comparaison. Si la variabilité expliquée par l'ensemble V de variables est du même ordre de grandeur que la variabilité résiduelle, alors cela veut dire que l'ensemble V de variables explique peu de choses, à peu près autant que ce que la résiduelle explique, c'est-à-dire autant que la variabilité que l'on n'arrive pas à expliquer. Dans ce cas, on aura envie de dire que l'ensemble V de variables n'apporte pas une information significativement intéressante pour mieux expliquer Y. Et au contraire, si l'ensemble V apporte beaucoup plus d'information que la résiduelle, c'est-à-dire si le F observé est grand, c'est-à-dire si le F obs est grand, Alors on aura envie de dire que l'ensemble V apporte une information significativement intéressante pour comprendre Y. La loi de la statistique de test sous l'hypothèse H0 est une loi de Fisher avec les degrés de liberté du numérateur et les degrés de liberté du dénominateur. Avec ce test, retenir l'hypothèse H0 revient à considérer que le sous-modèle, sans les variables de l'ensemble V, est aussi bon que le modèle complet. Le test conduit à choisir le sous-modèle si la variabilité supplémentaire expliquée par l'ensemble V de variables est du même ordre de grandeur que la variabilité résiduelle, c'est-à-dire la statistique de Fisher n'est pas significativement différente de 1. Le plus souvent, l'ensemble V est restreint à une seule variable, ce qui revient à tester si cette variable apporte une information supplémentaire intéressante sachant que les autres variables sont déjà présentes dans le modèle. L'autre situation classique est au contraire d'avoir l'ensemble V qui contient toutes les variables. L'hypothèse H0 correspond au modèle nul, c'est-à-dire au modèle sans variable explicative. Ce modèle nul n'est absolument pas intéressant pour prédire et indique, s'il est préféré au modèle complet, que le modèle complet n'apporte pas d'information complémentaire intéressante par rapport à un modèle qui n'explique rien. Et par conséquent le modèle complet n'est pas intéressant et donc les variables explicatives sont mal choisies. Si l'ensemble V contient toutes les variables, le test revient à tester si le R² est significativement différent de 0, c'est-à-dire si toutes les variables explicatives ont un coefficient égal à 0 et sont donc inutiles, ou au contraire si au moins une variable est utile pour prédire Y. Pour construire le test, on calcule la variabilité expliquée par l'ensemble V de variables sachant que toutes les autres variables sont déjà dans le modèle. On est donc bien en train de voir si l'ensemble V de variables apporte un plus par rapport aux autres variables présentes dans le modèle. Et on divise cette variabilité par les degrés de liberté associés à cet ensemble de variables. Pour calculer les degrés de liberté, il suffit de sommer les degrés de liberté de chaque variable et interaction présente dans l'ensemble V sachant qu'une variable quantitative a 1 degré de liberté, Une variable qualitative a I-1 degrés de liberté et une interaction a comme degré de liberté le produit des degrés de liberté de chaque facteur. Il peut être intéressant comme exercice d'écrire le test dans les cas particuliers suivants: tester l'effet d'une variable qualitative, tester l'effet d'une interaction, tester l'effet de toutes les variables. Pour l'instant, nous avons construit un modèle utilisant toutes les variables explicatives qui nous semblent pouvoir influer sur Y et nous avons également ajouté des interactions qui peuvent également avoir un effet sur Y. Cependant, ce modèle peut nécessiter l'estimation de nombreux paramètres et il est essentiel de sélectionner les variables utiles qui seront conservées dans le modèle. Cela permet d'une part de simplifier la compréhension du phénomène que l'on étudie en supprimant quelques variables ou interactions. Mais cela permet aussi de construire un modèle plus parcimonieux pour lequel on aura estimé moins de paramètres et qui sera plus performant d'un point de vue prédictif. Se pose alors la question de comment sélectionner un bon modèle. Et qu'est-ce qu'un bon modèle? Un bon modèle est un modèle qui permet de bien comprendre les relations entre les variables explicatives et la variable réponse et qui permet également de bien prédire la variable réponse pour de nouvelles observations. Pour sélectionner un bon modèle, on peut alors utiliser 2 grandes stratégies. La première est de sélectionner le sous-modèle pour lequel la probabilité critique associée au test de significativité du R² est la plus petite. On choisit alors le sous-modèle pour lequel on rejette le plus fortement l'hypothèse que le sous-modèle n'est pas intéressant, ce qui indique que le sous-modèle est intéressant. Une autre façon de faire est de sélectionner le modèle à partir d'un critère comme le critère d'Akaïke noté AIC ou encore le critère d'information bayésien noté BIC pour Bayesian Information Criterion. Le sous-modèle qui a l'AIC ou le BIC le plus petit est le meilleur. Sélectionner selon l'AIC et le BIC peut conduire à des sous-modèles différents, le critère BIC retenant des sous-modèles avec moins de variables. On utilisera de préférence le critère AIC quand on veut bien prédire et le critère BIC quand on veut sélectionner les variables explicatives et comprendre le phénomène. Ces 2 critères fonctionnent selon le même principe qui est de retenir un modèle avec la variabilité résiduelle minimum, en toute rigueur l'opposé du logarithme de la vraisemblance minimum, mais tout en ayant peu de paramètres. Les critères sont d'autant plus pénalisés que le nombre de variables du modèle est grand. Quand les variables explicatives sont quantitatives et nombreuses, on recommande d'utiliser le critère BIC et quand les variables explicatives sont qualitatives, on recommande plutôt l'AIC. Comment maintenant trouver le sous-modèle qui a la p-value la plus petite pour le test du R² ou qui a l'AIC le plus petit ou le BIC le plus petit? Une stratégie très basique consiste à construire tous les sous-modèles possibles, à calculer le critère qui nous intéresse pour chaque sous-modèle et à retenir le sous-modèle qui minimise notre critère. Cependant, le nombre de sous-modèles possibles devient très vite extrêmement grand et les temps de calcul prohibitifs dès que le nombre de variables explicatives, variables explicatives ou interactions, du modèle est grand. Pour contourner ce problème, plusieurs stratégies sont possibles. La première, appelée méthode descendante, consiste à construire le modèle complet puis à supprimer la variable ou l'interaction la moins intéressante, donc qui diminue très peu le critère qui nous intéresse. Rappelons-nous que, quel que soit le critère que l'on utilise, on cherche à le minimiser. Et on itère jusqu'à ce que, quelle que soit la variable que l'on supprime, le critère ne diminue plus. On a alors trouvé un sous-modèle qui minimise notre critère parmi tous les sous-modèles qui ont été explorés. Cela ne garantit pas que nous avons trouvé le sous-modèle avec le critère minimal, mais cela reste un bon sous-modèle. Une autre stratégie consiste à faire l'inverse, à partir du modèle nul sans variable et à ajouter la variable qui améliore, c'est-à-dire qui diminue, le plus possible notre critère d'intérêt. Ensuite, on garde ce sous-modèle et on itère en ajoutant la variable ou l'interaction qui améliore le plus le critère. Et on continue jusqu'à ce que, quelle que soit la variable ou l'effet ajouté, le critère se mette à augmenter. Cette méthode est appelée méthode ascendante. On peut ensuite avoir une méthode nommée stepwise qui combine la méthode descendante et la méthode ascendante, ce qui permet de rajouter une variable qui aurait été supprimée lors des premières étapes ou à enlever une variable qui aurait été rajoutée et qui n'est plus intéressante une fois que d'autres variables ont été ajoutées. Cette méthode stepwise est rapide et permet souvent de trouver un sous-modèle proche du meilleur modèle qu'on ne peut obtenir avec certitude qu'avec la méthode de construction exhaustive de tous les sous-modèles. Reprenons notre exemple sur les données ozone et utilisons la fonction linear_model du package FactoMineR pour choisir un sous-modèle avec le critère BIC. La fonction linear_model retourne les résultats pour le modèle complet et pour le sous-modèle sélectionné. On voit par exemple que le R² du sous-modèle est légèrement plus petit, 0.7622 contre 0.7686 pour le modèle complet. Ainsi, le pourcentage de variabilité du maximum d'ozone expliqué par le sous-modèle est quasiment le même qu'avec le modèle complet. Alors qu'on utilise beaucoup moins de variables et donc qu'on estime moins de paramètres. Les valeurs de la probabilité critique, de l'AIC et du BIC du sous-modèle sélectionné sont plus petites que celles du modèle complet. Le sous-modèle sélectionné est écrit ici et contient seulement 4 variables: T12, NE9, VX9 et MAXO3V. On peut ensuite regarder le tableau qui donne pour chaque variable explicative et interaction les résultats du test de la significativité de cet effet sachant que les autres variables sont présentes dans le modèle. Dans le tableau des tests F donc, on peut lire que toutes les p-values sont inférieures à 5% et donc chacune des 4 variables apporte une information significative supplémentaire utile sachant que les 3 autres variables sont dans le modèle. Dans certains cas, on peut avoir des résultats avec certaines variables ou interactions qui ont une probabilité critique supérieure à 5%. Le test F nous indiquerait que l'on peut supprimer la variable ou l'interaction qui a la p-value la plus élevée alors que le critère BIC a conservé cette variable dans le sous-modèle. Ces situations peuvent arriver et conduisent à une légère contradiction au niveau de l'interprétation sur l'utilité de garder ou non cette variable dans le modèle. Mais cette variable aura généralement peu d'effet sur la prédiction de nouvelles valeurs et la garder ou non dans le sous-modèle ne modifient pas beaucoup les prédictions. Nous verrons ensuite comment interpréter le 2ème tableau avec les tests t de Student dans la prochaine vidéo de cours. Quelle démarche adopter alors pour la modélisation? La première étape est de lister toutes les variables qui entrent en jeu pour expliquer ou prédire la variable réponse. Cette étape est essentielle et il ne faut surtout pas oublier une variable qui peut influer sur Y. Prendre une variable en trop n'est pas grave car les tests statistiques permettront d'éliminer cette variable. 2ème étape, il faut visualiser les données et notamment les liaisons avec la variable réponse pour s'assurer qu'il n'y a pas de valeurs aberrantes et pour avoir quelques idées sur les liaisons entre variables. On peut également faire une ACP en mettant en supplémentaire la variable réponse et en actif toutes les variables explicatives. Ensuite on écrit puis construit le modèle en choisissant à la fois les effets et les interactions qui peuvent expliquer la variable réponse. On ne mettra pas systématiquement et sans réfléchir toutes les interactions possibles dans le modèle. On peut souvent écarter a priori certaines interactions quand on connaît le phénomène que l'on étudie. Et surtout on évite de mettre des interactions d'ordre supérieur à 2. Ensuite, à partir du modèle complet, on sélectionne le sous-modèle qui minimise l'AIC ou le BIC, ou bien on fait de la sélection à la main, en supprimant les interactions et effets non utiles, puis on construit ce sous-modèle. Ensuite, et c'est le plus important, on interprétera les résultats, les effets significatifs comme les non significatifs. On interprétera les coefficients du modèle en faisant bien attention aux possibles confusions notamment pour les variables quantitatives. Et enfin, on utilise le sous-modèle pour prédire de nouvelles valeurs si on a un objectif de prévision. Nous verrons comment interpréter et prédire dans la prochaine vidéo de cours.