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)

Objectifs

L'objectif de ce tutoriel est de vous exercer à réaliser une régression polynomiale dans R avec la fonction lm().

Description des données

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.

Ombrine tachetée, avec son ocelle caractéristique à la base de la nageoire caudale.

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)

Modélisation

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

Choix du meilleur modèle

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 variantes map_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 pour map_dfr(), ou un vecteur de valeurs numériques doubles pour map_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))$$

Conclusion

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


BioDataScience-Course/BioDataScience2 documentation built on April 12, 2025, 9:45 a.m.