5.4 Choix du modèle

Le choix d’un modèle non linéaire fait intervenir des critères identiques à ceux d’un modèle linéaire (qualité d’ajustement évaluée par l’AIC, inspection visuelle de l’ajustement dans le nuage de points), mais il fait aussi intervenir une dimension supplémentaire : le choix de la fonction mathématique à ajuster. Comme nous venons de le voir au travers de quelques modèles courants en biologie, le nombre de fonctions mathématiques résultant en des formes similaire, par exemple de type sigmoïde, est grand. Ainsi, le choix de la meilleure fonction à utiliser dans un cas particulier est rendu plus difficile.

Nous allons illustrer l’ajustement d’une courbe non linéaire par le choix et l’ajustement d’un modèle de croissance dans un jeu de données en modélisant la croissance somatique de l’oursin Paracentrotus lividus (Grosjean 2001). Ce jeu de données est disponible dans le package data.io, sous le nom urchin_growth.

SciViews::R("model", lang = "fr")
urchins <- read("urchin_growth", package = "data.io")
chart(data = urchins, diameter ~ age) +
  geom_point()

Comme vous pouvez le voir, différents oursins ont été mesurés via le diamètre de leur test à son ambitus (zone la plus large du squelette de l’oursin, hors piquants) en mm à différents âges (en années). Les mesures ont été effectuées tous les 3 à 6 mois pendant plus de 10 ans, ce qui donne un bon aperçu de la croissance de cet animal y compris la taille maximale asymptotique qui est atteinte vers les 4 à 5 ans (pour ce genre de modèle, il est très important de continuer à mesurer les animaux afin de bien quantifier cette taille maximale asymptotique). Ainsi, l’examen du graphique nous permet d’emblée de choisir un modèle à croissance finie (pas le modèle de Tanaka, donc), et de forme sigmoïdale. Le modèle logistique, Weibull ou Gompertz pourrait convenir par exemple. Nous pouvons à ce stade, essayer différents modèles et choisir celui qui nous semble le plus adapté.

Le choix du meilleur modèle se fait grâce à deux critères :

Les connaissances théoriques et a priori du modèle que l’on ajuste. En effet, il n’existe qu’un seul modèle linéaire, mais une infinité de modèles curvilinéaires qui peuvent s’ajuster dans les données. Le choix du meilleur modèle se fait en fonction de considérations sur le phénomène sous-jacent qui doivent se refléter dans les propriétés mathématiques de la courbe choisie. Par exemple, si on sait que la croissance est asymptotique vers une taille maximale, nous devrons choisir une courbe mathématique qui présente une asymptote horizontale à son maximum pour représenter au mieux le phénomène étudié.
Le coefficient de détermination R² n’est pas calculé par R pour une régression non linéaire, car sa validité est sujette à discussion entre les statisticiens (d’autres logiciels statistiques le calculent). Nous n’avons donc pas d’estimation de la qualité de l’ajustement par ce biais, comme dans le cas de la régression linéaire. Par contre, il est possible de calculer un autre critère plus fiable que nous avons déjà utilisé : le critère d’Akaike (fonction AIC() dans R). Ce critère tient compte à la fois de la qualité d’ajustement et de la complexité du modèle, exprimée par le nombre de paramètres qu’il faut estimer. Plus le modèle est complexe, plus on peut s’attendre à ce qu’il s’ajuste bien dans les données, car il est plus flexible. Cependant, en ce domaine, la complexité n’est pas forcément un gage de qualité. On recherche plutôt un compromis entre meilleur ajustement et simplicité. Le critère d’information d’Akaike quantifie précisément ce compromis, c’est-à-dire que le modèle qui a un AIC le plus faible est considéré comme le meilleur. Appliquons donc ce concept pour sélectionner le meilleur modèle de croissance pour décrire la croissance somatique de nos oursins après avoir sélectionné les modèles candidats les plus judicieux (modèles sigmoïdaux avec asymptote horizontale au maximum).

Notons toutefois que, comme les animaux sont mesurés aux mêmes âges, tous les 3 à 6 mois, certains points se superposent. Afin d’afficher tous les points, il est utile d’utiliser la fonction geom_jitter() à la place de geom_point() qui décale les points d’une valeur aléatoire pour éviter ces superpositions (l’argument width = indique le décalage maximum à appliquer). Voici ce que cela donne (en ajoutant également un titre avec formatage correct du nom en latin) :

urchins_plot <- chart(data = urchins, diameter ~ age) +
  geom_jitter(width = 0.1, alpha = 0.2, color = "darkgrey") +
  ggtitle(expression(paste("Croissance de l'oursin ", italic("Paracentrotus lividus"))))
urchins_plot

Nous avons ici également représenté les points de manière semi-transparente avec alpha = 0.2 (transparence de 20%) pour encore mieux mettre en évidence les points de mesures qui se superposent.

Ajustons maintenant un modèle de Gompertz (modèle selfStart, notez qu’ici tabularise() est capable de calculer l’équation du modèle tout seul) :

urchins_gomp <- nls(data = urchins, diameter ~ SSgompertz(age, Asym, b2, b3))
summary(urchins_gomp) |> tabularise()

Modèle non linéaire de Gompertz
$\begin{aligned} \operatorname{diameter} = Asym \cdot e^{(-b2 \cdot b3^{\operatorname{age}})} + \epsilon \end{aligned}$
Terme	Valeur estimée	Ecart type	Valeur de tobs.	Valeur de p
Asym	57.404	0.25759	222.9	< 2·10-16	***
b2	3.902	0.04635	84.2	< 2·10-16	***
b3	0.435	0.00385	112.9	< 2·10-16	***
0 <= '*' < 0.001 < '' < 0.01 < '*' < 0.05
Ecart type des résidus : 5.49 pour 7021 degrés de liberté
Nombre d'itérations pour converger : 3
Tolérance atteinte à la convergence : 2.655e-06

Utilisons maintenant l’astuce de stat_function() combinée à as.function.nls() pour ajouter la courbe sur le graphique.

urchins_plot +
  stat_function(fun = as.function(urchins_gomp), color = "red", linewidth = 1)

L’ajustement de cette fonction semble très bon, à l’œil. Voyons ce qu’il en est d’autres modèles. Par exemple, une courbe logistique :

urchins_logis <- nls(data = urchins, diameter ~ SSlogis(age, Asym, xmid, scal))
summary(urchins_logis) |> tabularise()

Modèle non linéaire logistique
$\begin{aligned} \operatorname{diameter} = \frac{Asym}{1 + e^{(xmid - \operatorname{age}) /scal}} + \epsilon \end{aligned}$
Terme	Valeur estimée	Ecart type	Valeur de tobs.	Valeur de p
Asym	54.628	0.20299	269	< 2·10-16	***
xmid	2.055	0.00957	215	< 2·10-16	***
scal	0.765	0.00735	104	< 2·10-16	***
0 <= '*' < 0.001 < '' < 0.01 < '*' < 0.05
Ecart type des résidus : 5.6 pour 7021 degrés de liberté
Nombre d'itérations pour converger : 4
Tolérance atteinte à la convergence : 1.079e-06

Et voici le graphique avec les deux modèles superposés :

urchins_plot +
  stat_function(fun = as.function(urchins_gomp), aes(color = "Gompertz"), linewidth = 1) +
  stat_function(fun = as.function(urchins_logis), aes(color = "logistique"), linewidth = 1) +
  labs(color = "Modèle")

Notez qu’ici, la couleur a été incluse dans le “mapping” (argument mapping =) de stat_function() en l’incluant dans aes(). Cela change fondamentalement la façon dont la couleur est perçue par {ggplot2}. Dans ce cas-ci, la valeur est interprétée non comme une couleur à proprement parler, mais comme un niveau (une couche) à inclure dans le graphique et à reporter via une légende. Ensuite, à l’aide de labs() on change le titre de la légende relatif à la couleur par un nom plus explicite : “Modèle”.

Nous pouvons comparer ces modèles à l’aide du critère d’Akaike.

AIC(urchins_gomp, urchins_logis)

#               df      AIC
# urchins_gomp   4 43861.73
# urchins_logis  4 44139.13

Comme on peut le voir clairement sur le graphique, la courbe logistique donne une autre solution, cette dernière est à peine moins bonne que le modèle de Gompertz, selon le critère d’Akaike. Pourtant, la courbe est assez bien démarquée de celle de Gompertz. Essayons maintenant un modèle de Weibull. Ce modèle est plus complexe, car il a quatre paramètres au lieu de trois pour les deux modèles précédents :

urchins_weib <- nls(data = urchins, diameter ~ SSweibull(age, Asym, Drop, lrc, pwr))
summary(urchins_weib) |> tabularise()

Modèle non linéaire de Weibull
$\begin{aligned} \operatorname{diameter} = Asym - Drop \cdot e^{-e^{lrc} \cdot \operatorname{age}^{pwr}} + \epsilon \end{aligned}$
Terme	Valeur estimée	Ecart type	Valeur de tobs.	Valeur de p
Asym	56.80	0.3270	173.7	< 2·10-16	***
Drop	56.81	0.5939	95.7	< 2·10-16	***
lrc	-1.57	0.0286	-54.9	< 2·10-16	***
pwr	1.68	0.0296	56.7	< 2·10-16	***
0 <= '*' < 0.001 < '' < 0.01 < '*' < 0.05
Ecart type des résidus : 5.47 pour 7020 degrés de liberté
Nombre d'itérations pour converger : 3
Tolérance atteinte à la convergence : 2.532e-06

Ajoutons ce nouveau modèle sur le graphique :

urchins_plot +
  stat_function(fun = as.function(urchins_gomp), aes(color = "Gompertz"), size = 1) +
  stat_function(fun = as.function(urchins_logis), aes(color = "logistique"), size = 1) +
  stat_function(fun = as.function(urchins_weib), aes(color = "Weibull"), size = 1) +
  labs(color = "Modèle")

# Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
# ℹ Please use `linewidth` instead.
# This warning is displayed once every 8 hours.
# Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
# generated.

… et comparons à l’aide du critère d’Akaike :

AIC(urchins_gomp, urchins_logis, urchins_weib)

#               df      AIC
# urchins_gomp   4 43861.73
# urchins_logis  4 44139.13
# urchins_weib   5 43810.72

Ce modèle fait presque jeu égal avec le modèle de Gompertz d’après le critère d’Akaike ; juste un tout petit peu mieux. En fait, les deux courbes sont pratiquement superposées l’une à l’autre, mais le modèle de Weibull à un démarrage de croissance plus lent au début, ce qui reflète mieux ce que l’on observe dans les données. Par contre, il est pénalisé par le fait que c’est un modèle plus complexe qui possède un paramètre de plus. L’un dans l’autre, le critère d’information d’Akaike considère donc les deux modèles sur pratiquement le même plan du point de vue de la qualité de leurs ajustements respectifs.

À ce stade, nous voudrions également essayer un autre modèle flexible à quatre paramètres : le modèle de Richards. Malheureusement, il n’existe pas de fonction selfStart dans R pour ce modèle. Nous sommes donc réduits “à mettre les mains dans le cambouis”, à définir la fonction nous-mêmes, à trouver de bonnes valeurs de départ, etc. Voici comment définir la fonction :

richards <- function(x, Asym, lrc, c0, m) Asym*(1 - exp(-exp(lrc) * (x - c0)))^m

Pour les valeurs de départ, là ce n’est pas facile. Asym est l’ asymptote horizontale à la taille maximum. On voit qu’elle se situe aux environs de 55 mm sur le graphique. Pour les autres paramètres, c’est plus difficile à évaluer. Prenons par exemple 1 comme valeur de départ pour les trois autres paramètres, ce qui donne (les valeurs de départ sont obligatoires ici puisque ce n’est pas un modèle selfStart) :

urchins_rich <- nls(data = urchins, diameter ~ richards(age, Asym, lrc, c0, m),
  start = c(Asym = 55, lrc = 0.1, c0 = 1, m = 1))

# Error in numericDeriv(form[[3L]], names(ind), env, central = nDcentral): Missing value or an infinity produced when evaluating the model

… et voilà ! Un excellent exemple de plantage de l’algorithme de minimisation de la fonction objective suite à un comportement inadéquat de la fonction avec les valeurs testées. Ici, la fonction renvoie l’infini et l’algorithme ne peut donc effectuer la minimisation. La fonction de Richards est effectivement connue pour être difficile à ajuster pour cette raison.

Il nous faut donc soit tester d’autres valeurs de départ, soit utiliser un autre algorithme de minimisation, soit les deux. Après différents essais il apparaît que le changement des valeurs de départ suffit dans le cas présent :

urchins_rich <- nls(data = urchins, diameter ~ richards(age, Asym, lrc, c0, m),
  start = c(Asym = 55, lrc = -0.7, c0 = 0, m = 1))
# N'étant pas une fonction selfStart, tabularise() ne peut pas déterminer
# l'équation du modèle toute seule. Donc, nous ne l'affichons pas.
summary(urchins_rich) |> tabularise(equation = FALSE)

Terme	Valeur estimée	Ecart type	Valeur de tobs.	Valeur de p
Asym	58.143	0.3777	153.93	< 2·10-16	***
lrc	-0.276	0.0315	-8.75	< 2·10-16	***
c0	-0.875	0.2846	-3.08	2.11·10-03	**
m	6.207	1.8234	3.40	6.68·10-04	***
0 <= '*' < 0.001 < '' < 0.01 < '*' < 0.05
Ecart type des résidus : 5.487 pour 7020 degrés de liberté
Nombre d'itérations pour converger : 11
Tolérance atteinte à la convergence : 7.332e-06

Ajoutons ce dernier modèle sur notre graphique (avec des traits un peu plus fins pour mieux distinguer les modèles les uns des autres) :

urchins_plot +
  stat_function(fun = as.function(urchins_gomp), aes(color = "Gompertz"), size = 0.7) +
  stat_function(fun = as.function(urchins_logis), aes(color = "logistique"), size = 0.7) +
  stat_function(fun = as.function(urchins_weib), aes(color = "Weibull"), size = 0.7) +
  stat_function(fun = as.function(urchins_rich), aes(color = "Richards"), size = 0.7) +
  labs(color = "Modèle")

… et comparons à l’aide du critère d’Akaike :

AIC(urchins_gomp, urchins_logis, urchins_weib, urchins_rich)

#               df      AIC
# urchins_gomp   4 43861.73
# urchins_logis  4 44139.13
# urchins_weib   5 43810.72
# urchins_rich   5 43854.53

La courbe est très proche des modèles de Gompertz et Weibull aux jeunes âges, mais l’asymptote maximale est légèrement plus haute que pour les deux autres modèles (58 mm au lieu de 57 mm). Les trois courbes sont très, très proches l’une de l’autre. Le critère d’information d’Akaike est marginalement moins bon pour le modèle de Richards que pour celui de Weibull, mais est tout juste meilleur que celui de Gompertz. En outre, l’écart type pour le paramètre x0 est plus conséquent en comparaison de sa valeur, ce qui démontre une certaine instabilité de la fonction par rapport à ce paramètre, et par conséquent, une incertitude dans son estimation. Pour cette raison, et pour la difficulté à l’ajuster, le modèle de Richards sera écarté dans notre cas au bénéfice du modèle de Weibull, voire de celui de Gompertz plus simple.

Le choix final entre Gompertz ou Weibull dépend de l’usage que l’on veut faire du modèle. Si la simplicité du modèle est primordiale, nous garderons Gompertz. Si la croissance des petits oursins est un aspect important de l’analyse, nous garderons Weibull qui semble mieux s’ajuster aux données à ce niveau.

À vous de jouer !

Effectuez maintenant les exercices du tutoriel B05La_nls (Modèles non linéaires).

BioDataScience2::run("B05La_nls")

Réalisez le travail B05Ia_abies_balsamea.

Travail individuel pour les étudiants inscrits au cours de Science des Données Biologiques II à l’UMONS (Q1 : modélisation) à terminer avant le 2024-12-09 23:59:59.

Initiez votre projet GitHub Classroom

Voyez les explications dans le fichier README.md.

Réalisez en groupe le travail B02Ga_models, partie IV.

Travail en groupe de 4 pour les étudiants inscrits au cours de Science des Données Biologiques II à l’UMONS (Q1 : modélisation) à terminer avant le 2024-12-17 23:59:59.

Initiez votre projet GitHub Classroom

Voyez les explications dans le fichier README.md, partie IV.

Références

Grosjean, Philippe. 2001. “Growth Model of the Reared Sea Urchin Paracentrotus Lividus (Lamarck, 1816).” PhD thesis, ULB. http://go.sciviews.org/thesis_PhG.