Document complémentaire au module 5 du cours SDD II de 2025-2026. Distribué sous licence CC BY-NC-SA 4.0.

Veuillez vous référer au cours en ligne pour les explications et les interprétations de cette analyse.

Installer un environnement R adéquat pour reproduire cette analyse.

Rendement photosynthétique

# Préparation des données de "rapid light curve"
rlc <- as_dtx(tribble(
  ~etr, ~par,
  0.0,  0,
  0.5,  2,
  3.2,  11,
  5.7,  27,
  7.4,  50,
  8.4,  84,
  8.9,  170,
  8.4,  265,
  7.8,  399
))
# Labélisation
rlc <- labelise(rlc,
  label = list(
    etr = "ETR",
    par = "PAR"),
  units = list(
    etr = "µmol électrons/m²/s",
    par = "µmol photons/m²/s")
)
# Graphique général de la RLC
chart(data = rlc, etr ~ par) +
  geom_point()

# Modèle polynomial RLC (inadéquat)
rlc_lm <- lm(data = rlc, etr ~  par + I(par^2) + I(par^3))
chart(rlc_lm)

# Équation du modèle PGH
pgh_model <- function(x, etr_max, alpha, beta)
  etr_max * (1 - exp(-x * alpha/etr_max)) * (exp(-x * beta/etr_max))

# Régression non linéaire en utilisant PGH
rlc_nls <- nls(data = rlc, etr ~ pgh_model(par, etr_max, alpha, beta),
  start = list(etr_max = 9, alpha = 1, beta = 0))
# Note : tabularise() ne peut pas encore traiter les modèles non "self-start"
# (vous verrez ce que cela signifie plus loin) pour les équations. Donc,
# nous demandons le tableau SANS l'équation du modèle
summary_(rlc_nls) |> tabularise(equation = FALSE)

Terme

Valeur estimée

Ecart type

Valeur de tobs.

Valeur de p

signif

etr_max

9.35106

0.26197

35.70

3.22·10-08

***

alpha

0.32738

0.01342

24.40

3.11·10-07

***

beta

0.00398

0.00108

3.67

1.04·10-02

*

0 <= '***' < 0.001 < '**' < 0.01 < '*' < 0.05

Ecart type des résidus : 0.1729 pour 6 degrés de liberté

Nombre d’itérations pour converger : 7
Tolérance atteinte à la convergence : 3.814e-06

# Une autre solution est de fournir l'équation en LaTeX à tabularise()
#summary(rlc_nls) |> tabularise(equation = "ETR = ETR_{max} \\cdot (1 - e^{-PAR \\cdot \\alpha/ETR_{max}}) \\cdot e^{-PAR \\cdot \\beta/ETR_{max}}")

# Graphique de la RLC avec le modèle PGH non linéaire ajusté
chart(data = rlc, etr ~ par) +
  geom_point() +
  stat_function(fun = as.function(rlc_nls), col = "skyblue3", linewidth = 1)

# On peut aussi faire simplement (mais c'est moins flexible) :
#chart(rlc_nls)

# Comparaison des deux modèles par le critère d'Akaike
AIC(rlc_lm, rlc_nls)
##         df       AIC
## rlc_lm   5 31.834047
## rlc_nls  4 -1.699078
# Ajustement du modèle avec trace=TRUE pour voir la progression des étapes
rlc_nls <- nls(data = rlc, etr ~ pgh_model(par, etr_max, alpha, beta),
  start = list(etr_max = 9, alpha = 1, beta = 0),
  trace = TRUE)
## 24.34024    (3.38e+00): par = (9 1 0)
## 4.236559    (2.81e+00): par = (8.478564 0.5324097 -0.0004037391)
## 1.141967    (2.31e+00): par = (8.772552 0.2898063 0.001855493)
## 0.1807006   (8.51e-02): par = (9.361593 0.3255633 0.003954183)
## 0.1793698   (4.48e-03): par = (9.353621 0.3272301 0.003990887)
## 0.1793659   (4.22e-04): par = (9.351297 0.32737 0.003982181)
## 0.1793658   (3.98e-05): par = (9.351085 0.3273833 0.00398142)
## 0.1793658   (3.81e-06): par = (9.351065 0.3273846 0.003981348)

Modèles courants en biologie

# Objects R commençant par "SS" (modèles self-start pour la plupart)
apropos("^SS", ignore.case = FALSE)
##  [1] "SSasymp"     "SSasympOff"  "SSasympOrig" "SSbiexp"     "SSD"         "SSfol"       "SSfpl"       "SSgompertz"  "SSlogis"     "SSmicmen"   
## [11] "SSweibull"
# Fonction exponentielle (pas de modèle self-start)
exponent <- function(x, y0, k)
  y0 * exp(k * x)

# Modèle von Bertalanffy en poids (pas de modèle self-start)
asympOff3 <- function(x, Asym, lrc, c0, m) Asym*(1 - exp(-exp(lrc) * (x - c0)))^3

# Equation de Richards (pas de modèle self-start)
richards <- function(x, Asym, lrc, c0, m) Asym*(1 - exp(-exp(lrc) * (x - c0)))^m

# Equation de Preece Baines 1 de croissance humaine (pas de modèle self-start)
preece_baines1 <- function(x, Asym, Drop, lrc1, lrc2, c0)
  Asym - (2 * (Asym - Drop)) / (exp(exp(lrc1) * (x - c0)) + exp(exp(lrc2) * (x - c0)))

# Équation de Tanaka de croissance indéterminée (pas de modèle self-start)
tanaka <- function(x, a, b, c0, d)
  1 / sqrt(b) * log(abs(2 * b * (x - c0) + 2 * sqrt(b^2 * (x - c0)^2 + a * b))) + d

Choix du modèle

# Configuration de R pour le dialecte SciViews avec modélisation
SciViews::R("model", lang = "fr")
# Chargement et visualisation des données
urchins <- read("urchin_growth", package = "data.io")
chart(data = urchins, diameter ~ age) +
  geom_point()

# Graphique avec points jitterisés
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

# Régression non linéaire avec modèle de Gompertz
urchins_gomp <- nls(data = urchins, diameter ~ SSgompertz(age, Asym, b2, b3))
summary_(urchins_gomp) |> tabularise()

Modèle non linéaire de Gompertz

diameter⁡=Asym⋅e(−b2⋅b3age⁡)+ϵ\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

signif

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

# Graphique du modèle de Gompertz ajusté
urchins_plot +
  stat_function(fun = as.function(urchins_gomp), color = "red", linewidth = 1)

# Régression non linéaire avec modèle logistique
urchins_logis <- nls(data = urchins, diameter ~ SSlogis(age, Asym, xmid, scal))
summary_(urchins_logis) |> tabularise()

diameter⁡=Asym1+e(xmid−age⁡)/scal+ϵ\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

signif

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

# Graphique avec les deux modèles ajustés pour comparaison visuelle
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")

# Comparaison des deux modèles via le critère d'Akaike
AIC(urchins_gomp, urchins_logis)
##               df      AIC
## urchins_gomp   4 43861.73
## urchins_logis  4 44139.13
# Régression non linéaire avec modèle de Weibull
urchins_weib <- nls(data = urchins, diameter ~ SSweibull(age, Asym, Drop, lrc, pwr))
summary_(urchins_weib) |> tabularise()

diameter⁡=Asym−Drop⋅e−elrc⋅age⁡pwr+ϵ\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

signif

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

# Graphique des trois modèles ajusté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) +
  stat_function(fun = as.function(urchins_weib), aes(color = "Weibull"), linewidth = 1) +
  labs(color = "Modèle")

# Comparaison par 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
# Équation de Richards
richards <- function(x, Asym, lrc, c0, m) Asym*(1 - exp(-exp(lrc) * (x - c0)))^m

try({
# Régression non linéaire avec modèle de Richards
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

signif

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.309e-06

# Graphique avec les 4 modèles ajustés
urchins_plot +
  stat_function(fun = as.function(urchins_gomp), aes(color = "Gompertz"), linewidth = 0.7) +
  stat_function(fun = as.function(urchins_logis), aes(color = "logistique"), linewidth = 0.7) +
  stat_function(fun = as.function(urchins_weib), aes(color = "Weibull"), linewidth = 0.7) +
  stat_function(fun = as.function(urchins_rich), aes(color = "Richards"), linewidth = 0.7) +
  labs(color = "Modèle")

# Comparaison par critère d'Akaike des 4 modèles
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