BioDataScience2::learnr_setup() SciViews::R("model", lang = "fr") # datasets reddrum <- read("reddrum", package = "UsingR") reddrum$length <- reddrum$length * 0.0254 reddrum <- labelise(reddrum, label = list(length = "Longueur", age = "Âge"), units = list(length = "m", age = "années")) lm1 <- lm(data = reddrum, length ~ age) lm2 <- lm(data = reddrum, length ~ age + I(age^2)) lm3 <- lm(data = reddrum, length ~ age + I(age^2) + I(age^3)) lm_poly_coef <- tidy(lm3) lm_poly_param <- glance(lm3)
BioDataScience2::learnr_banner()
BioDataScience2::learnr_server(input, output, session)
L'objectif de ce tutoriel est de vous exercer à réaliser une régression polynomiale dans R avec la fonction lm()
.
Le jeu de données reddrum
traite de la croissance de l'ombrine tachetée Sciaenops ocellatus (L., 1766), un poisson da la famille des Scianidae qui vit sur la côte américaine de l'Océan Atlantique, du Massachusetts jusqu'au nord du Mexique. Les données de la longueur du poisson en fonction de son âge sont simulées sur base de mesures réalisées par Porch, Wilson et Nieland (2002). Cette espèce peut atteindre 1,5 mètre et peut vivre plusieurs dizaines d'années jusqu'à 50 ou 60 ans.
Voir Porch, C., C.A. Wilson & D. Nieland (2002). A new growth model for red drum (Scianops ocellatus) that accommodates seasonal and ontogenic changes in growth rates. Fish. Bull., 100:149-152.
SciViews::R("model",lang = "fr") # Configuration du système # Importation des données reddrum <- read("reddrum", package = "UsingR") # Conversion de pouce en mètre reddrum$length <- reddrum$length * 0.0254 # Ajout des labels et unités reddrum <- labelise(reddrum, label = list(length = "Longueur", age = "Âge"), units = list(length = "m", age = "années")) skimr::skim(reddrum)
Vous allez modéliser la longueur (m) de l'ombrine tachetée en fonction de son âge dont voici la représentation graphique.
chart(data = reddrum, length ~ age) + geom_point()
La croissance de ces poissons est particulière. Les juvéniles ont une croissance très rapide et ensuite, la croissance ralentit nettement sans jamais devenir nulle. La forme du nuage de points est telle qu'il est difficile, voire impossible de trouver une transformation capable de le linéariser. Commencez, à titre de référence, par ajuster une régression linéaire dans les données non transformées.
reddrum_lm1 <- lm(data = ___, ___ ~ ___) # Résumé du modèle summary(___) # Graphique du modèle chart(___)
reddrum_lm1 <- lm(data = ___, ___ ~ ___) # Résumé du modèle summary(reddrum_lm1) # Graphique du modèle chart(reddrum_lm1) #### ATTENTION: Hint suivant = solution !####
## Solution ## reddrum_lm1 <- lm(data = reddrum, length ~ age) # Résumé du modèle summary(reddrum_lm1) # Graphique du modèle chart(reddrum_lm1)
grade_code("Voici notre régression linéaire. Le R^2 est de 0.82. On pourrait penser que ce modèle est très intéressant, si on se limite à cette valeur. Cependant, en visualisant notre droite sur le graphique, on observe un problème d'ajustement du modèle pour les jeunes individus.")
Afin de compléter notre analyse, voici les graphiques des résidus (vous devriez être capable de les interpréter par vous-mêmes maintenant).
chart$residuals(lm1)
Réalisez ensuite une régression linéaire polynomiale d'ordre deux avec les mêmes variables.
reddrum_lm2 <- lm(data = ___, ___ ~ ___ + ___(___)) # Résumé du modèle summary(___) # Graphique du modèle chart(___)
reddrum_lm2 <- lm(data = ___, ___ ~ ___ + I(___)) # résumé du modèle summary(reddrum_lm2) # Graphique du modèle chart(reddrum_lm2) #### ATTENTION: Hint suivant = solution !####
## Solution ## reddrum_lm2 <- lm(data = reddrum, length ~ age + I(age^2)) # résumé du modèle summary(reddrum_lm2) # Graphique du modèle chart(reddrum_lm2)
grade_code("Cette fois-ci, nous avons ajusté une parabole dans les données... C'est déjà mieux, mais insuffisant.")
Analyse des résidus pour le modèle polynomial d'ordre deux.
chart$residuals(lm2)
Réalisez à présent une régression linéaire polynomiale d'ordre trois avec les mêmes variables.
reddrum_lm3 <- lm(data = ___, ___ ~ ___ + ___(___) + ___(___)) # Résumé du modèle summary(___) # Graphique du modèle chart(___)
reddrum_lm3 <- lm(data = ___, ___ ~ ___ + I(___) + ___(___)) # résumé du modèle summary(reddrum_lm3) # Graphique du modèle chart(reddrum_lm3) #### ATTENTION: Hint suivant = solution !####
## Solution ## reddrum_lm3 <- lm(data = reddrum, length ~ age + I(age^2) + I(age^3)) # résumé du modèle summary(reddrum_lm3) # Graphique du modèle chart(reddrum_lm3)
grade_code("Observez comme la courbe devient de plus en plus flexible à mesure que l'ordre du polynôme augmente.")
Analyse des résidus pour la régression polynomiale d'ordre trois.
chart$residuals(lm3)
Étudiez les résultats obtenus et répondez aux questions ci-dessous.
quiz( question(text = "Quelle est la valeur du paramètre pour le terme d'ordre deux dans le dernier modèle ?", answer(sprintf("%.4f", lm_poly_coef$estimate[1])), answer(sprintf("%.4f", lm_poly_coef$estimate[2])), answer(sprintf("%.4f", lm_poly_coef$estimate[3]), correct = TRUE), answer(sprintf("%.4f", lm_poly_coef$statistic[2])), answer(sprintf("%.4f", lm_poly_coef$statistic[3])), answer(sprintf("%.4f", lm_poly_param$r.squared[1])), allow_retry = TRUE, random_answer_order = TRUE ), question(text = "Quelle est la part de la variance exprimée par ce modèle ?", answer(sprintf("%.3f", lm_poly_coef$estimate[1])), answer(sprintf("%.3f", lm_poly_coef$estimate[2])), answer(sprintf("%.3f", lm_poly_coef$std.error[1])), answer(sprintf("%.3f", lm_poly_param$r.squared[1])), answer(sprintf("%.3f", lm_poly_coef$statistic[1])), answer(sprintf("%.3f", lm_poly_coef$statistic[2])), answer(sprintf("%.3f", lm_poly_param$adj.r.squared[1]), correct = TRUE), allow_retry = TRUE, random_answer_order = TRUE ) )
Nous avons pu observer que notre modèle s'ajuste de mieux en mieux en augmentant l'ordre de notre polynôme. Une méthode consiste à augmenter l'ordre du polynôme de manière itérative jusqu'à ce que le paramètre relatif au terme d'ordre maximum ne soit plus significativement différent de zéro. Ici, cela implique de tester aussi un modèle polynomial d'ordre 4, 5, 6, 7... Nous ne détaillons pas chaque modèle, mais voici ce que donne le polynôme d'ordre six.
reddrum_lm6 <- lm(data = reddrum, length ~ age + I(age^2) + I(age^3) + I(age^4) + I(age^5) + I(age^6)) summary(reddrum_lm6) chart(reddrum_lm6)
Tous les termes du modèle sont significatifs au seuil $\alpha$ de 5%. La valeur de R^2^ atteint 0.94 alors qu'elle n'était que de 0.82 pour la régression linéaire simple. Les graphiques de l'analyse des résidus sont ci-dessous. Nous avons quelques points avec effet de levier et distance de Cook importants (D), et un petit manque de constance des variances (C), mais pour le reste, ces graphiques ne montrent pas de problèmes particuliers.
chart$residuals(reddrum_lm6)
Voici à présent le résultat pour un polynôme d'ordre sept.
reddrum_lm7 <- lm(data = reddrum, length ~ age + I(age^2) + I(age^3) + I(age^4) + I(age^5) + I(age^6) + I(age^7)) summary(reddrum_lm7) chart(reddrum_lm7)
Nous voyons ici que le paramètre du terme en puissance sept n'est plus significativement différent de zéro au seuil $\alpha$ de 5%. Nous avons franchi la limite. Nous en concluons donc que le meilleur modèle polynomial est ici le polynôme d'ordre six.
Poussons notre réflexion encore un peu plus loin. Si nous calculons une régression polynomiale d'ordre n-1. La courbe va s'adapter parfaitement à la distribution de nos observations. La valeur de R^2^ sera de 1. Ce modèle ne sera d'aucune utilité. Nous serons dans un cas de surajustement. Il ne contient plus aucune information pertinente.
Voici une analyse de la régression linéaire simple à la régression polynomiale d'ordre neuf. Au lieu d'utiliser age + I(age^2) + I(age^3) + ...
, ce qui devient rapidement fastidieux, nous utilisons ici poly(x, n, raw = TRUE)
, avec n
, l'ordre du polynôme (l'explication relative à raw = TRUE
sort du cadre de ce cours, mais raw = FALSE
, la valeur par défaut, conduit à transformer le polynôme de sorte que ses paramètres soient orthogonaux l'un à l'autre, ce qui leur confère plus de stabilité).
# Ajustement des différents modèles au sein d'une liste models <- list( lm1 = lm(data = reddrum, length ~ age), lm2 = lm(data = reddrum, length ~ poly(age, 2, raw = TRUE)), lm3 = lm(data = reddrum, length ~ poly(age, 3, raw = TRUE)), lm4 = lm(data = reddrum, length ~ poly(age, 4, raw = TRUE)), lm5 = lm(data = reddrum, length ~ poly(age, 5, raw = TRUE)), lm6 = lm(data = reddrum, length ~ poly(age, 6, raw = TRUE)), lm7 = lm(data = reddrum, length ~ poly(age, 7, raw = TRUE)), lm8 = lm(data = reddrum, length ~ poly(age, 8, raw = TRUE)), lm9 = lm(data = reddrum, length ~ poly(age, 9, raw = TRUE)) )
Nous pouvons extraire les paramètres du modèle en utilisant glance()
et nous utilisons rmse()
pour calculer une métrique, l'erreur quadratique moyenne, qui quantifie l'ajustement du modèle.
# Extraction en série des paramètres des modèles grâce à glance() mod_params <- purrr::map_dfr(models, glance) # Ajout de RMSE selon le même principe mod_params$rmse <- purrr::map_dbl(models, rmse, data = reddrum) # Ajout du nom des modèles mod_params$model <- names(models) # Visualisation des résultats mod_params[ , c("model", "adj.r.squared", "rmse")]
Le chunk ci-dessus mérite une petite explication. Nous avons neuf modèles. Tous ces modèles sont regroupés dans une liste nommée
models
.purrr::map()
et ses variantesmap_dfr()
,map_dbl()
exécutent une fonction donnée en second argument sur chaque élément de la liste en premier argument, c'est-à-dire, sur chacun des neuf modèles que notre liste contient. Une liste est renvoyée, ou un data frame pourmap_dfr()
, ou un vecteur de valeurs numériques doubles pourmap_dbl()
.
Voici un graphique de la variation du R^2^ ajusté et de l'erreur quadratique moyenne en fonction de l'ordre du polynôme du modèle :
a <- chart(data = mod_params, adj.r.squared ~ model) + geom_line(group = 1) + geom_point() b <- chart(data = mod_params, rmse ~ model) + geom_line(group = 1) + geom_point() combine_charts(list(a, b), ncol = 1)
Nous pouvons observer un gain de performance des modèles jusqu'au modèle polynomial d'ordre cinq ou six. Ensuite, ces gains deviennent nuls ou quasi nuls. Ceci confirme qu'il est inutile d'augmenter l'ordre du polynôme au-delà de six et cela valide aussi la première méthode itérative que nous avons utilisée qui consiste à s'arrêter juste avant que le paramètre du terme de plus forte puissance ne soit plus significativement de zéro au seuil $\alpha$ fixé d'avance.
Pour le polynôme d'ordre deux, il est intéressant de déterminer si nous pouvons laisser tomber le terme de puissance un et simplifier le modèle. Cela signifie alors utiliser $y = \alpha + \beta \cdot x^2$, soit une droite dans les données après avoir élevé $x$ au carré. Pour les polynômes d'ordre plus élevé, il faut considérer le modèle comme un tout. Analysez-le sur base du paramètre de puissance la plus élevée uniquement, comme nous l'avons fait ici, et n'essayez pas de simplifier si des termes intermédiaires apparaissent non significatifs (il y a bien évidemment une forte colinéarité entre les paramètres de ces différents termes, ce qui rend leur interprétation individuelle caduque).
La visualisation du modèle, l'étude du résumé du modèle, l'analyse des résidus et l'étude des indicateurs de performances forment un tout. Vous devez maîtriser l'ensemble des outils pour être capable de valider la pertinence d'un modèle et pour pouvoir comparer des modèles entre eux.
Pour une présentation propre du modèle retenu, nous utilisons tabularise()
:
summary(reddrum_lm6) |> tabularise()
Enfin, il est indispensable de paramétrer le modèle. Nous laissons R le faire grâce à la fonction eq__()
que nous avons déjà employée plusieurs fois :
$$r eq__(reddrum_lm6, use_coefs = TRUE, coef_digits = c(3, 3, 4, 5, 7, 11))
$$
Vous venez de terminer ce tutoriel sur la régression polynomiale. La régression polynomiale d'ordre six nous a permis d'ajuster un modèle dans ces données malgré la forme non linéaire du nuage de points. Ce genre de modèle est utile dans un but de prédiction, mais il n'est d'aucune utilité pour aider à expliquer le mécanisme de la croissance de ce poisson. Pour cela, nous devrons nous tourner vers des modèles non linéaires spécialisés (modèles de croissance) que nous étudierons au module 5. Avec le modèle polynomial, vous devez rester très vigilant à deux aspects : le surajustement, car un polynôme d'ordre suffisamment élevé peut, à la limite, passer par tous les points, mais il inclut alors l'erreur de mesure ce qui est contre-productif, et (2) ne jamais faire des extrapolations lors de prédictions, car la courbe polynomiale s'écarte généralement très rapidement des données aux deux extrémités de plage utilisée lors de l'ajustement.
question_text( "Laissez-nous vos impressions sur ce learnr", answer("", TRUE, message = "Pas de commentaires... C'est bien aussi."), incorrect = "Vos commentaires sont enregistrés.", placeholder = "Entrez vos commentaires ici...", allow_retry = TRUE, submit_button = "Soumettre une réponse", try_again_button = "Resoumettre une réponse" )
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.