SDD II module 4 : Modèle linéaire généralisé

Document complémentaire au module 4 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.

GLM poisson : ray-grass dans les dunes

# Configuration de R en dialecte SciViews::R pour la modélisation
SciViews::R("model", lang = "fr")
# Importation et préparation des données
dune <- sbind_cols(
  read("dune", package = "vegan") |> sselect(Lolipere),
  read("dune.env", package = "vegan") |> sselect(A1, Moisture, Management)
) %>.%
  smutate(., Management = case_when(
    Management == "NM" ~ "conservation",
    .default = "culture") |> factor())
skimr::skim(dune)

Data summary
Name	dune
Number of rows	20
Number of columns	4
_______________________
Column type frequency:
factor	2
numeric	2
________________________
Group variables	None

Variable type: factor

skim_variable	n_missing	complete_rate	ordered	n_unique	top_counts
Moisture	0	1	TRUE	4	1: 7, 5: 7, 2: 4, 4: 2
Management	0	1	FALSE	2	cul: 14, con: 6

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
Lolipere	0	1	2.90	2.83	0.0	0.0	2.0	6.00	7.0	▇▃▁▂▆
A1	0	1	4.85	2.18	2.8	3.5	4.2	5.73	11.5	▇▂▁▁▁

# Abondance de L. perenne
chart(data = dune, ~ Lolipere) +
  geom_bar() +
  labs(x = "Abondance [nbr de plants/quadrat]",
       y = "Nombre de sites")

# Table de contingence humidité du sol versus mode de gestion des dunes
table(Humidité = dune$Moisture, Gestion = dune$Management) |>
  tabularise()

Humidité		Gestion
Humidité		conservation	culture	Total
1	Count	1 (5.0%)	6 (30.0%)	7 (35.0%)
1	Mar. pct (1)	16.7% ; 14.3%	42.9% ; 85.7%
2	Count	1 (5.0%)	3 (15.0%)	4 (20.0%)
2	Mar. pct	16.7% ; 25.0%	21.4% ; 75.0%
4	Count	0 (0.0%)	2 (10.0%)	2 (10.0%)
4	Mar. pct	0.0% ; 0.0%	14.3% ; 100.0%
5	Count	4 (20.0%)	3 (15.0%)	7 (35.0%)
5	Mar. pct	66.7% ; 57.1%	21.4% ; 42.9%
Total	Count	6 (30.0%)	14 (70.0%)	20 (100.0%)
(1) Columns and rows percentages

# Abondance de L. perenne en fonction de l'épaisseur de la couche A1 du sol
chart(data = dune, Lolipere ~ A1) +
  geom_point() +
  labs(x = "Épaisseur de la couche A1 [cm]",
       y = "Abondance [nbr de plants/quadrat]")

# GLM de famille poisson
dune_glm <- glm(data = dune, Lolipere ~ A1, family = poisson)
summary_(dune_glm) |> tabularise()


$\log ({ E( \operatorname{Lolipere} ) }) = \alpha + \beta_{}(\operatorname{A1})$
Terme	Valeur estimée	Ecart type	Valeur de z	Valeur de p
$\alpha$	3.569	0.591	6.04	1.53·10-09
$\beta_{}$	-0.607	0.154	-3.94	8.11·10-05
0 <= '*' < 0.001 < '' < 0.01 < '*' < 0.05
(Paramètre de dispersion pour une Poisson family: 1) Déviance totale : 68.58 pour 19 degrés de liberté Déviance résiduelle : 41.44 pour 18 degrés de liberté AIC: 85.75 - Nombre d’itérations de la fonction de score de Fisher : 5

# Même modèle, mais en GLM quasi-poisson
dune_glm_quasi <- glm(data = dune, Lolipere ~ A1, family = quasipoisson)
summary_(dune_glm_quasi) |> tabularise()


$\log ({ E( \operatorname{Lolipere} ) }) = \alpha + \beta_{}(\operatorname{A1})$
Terme	Valeur estimée	Ecart type	Valeur de t	Valeur de p
$\alpha$	3.569	0.771	4.63	0.000209
$\beta_{}$	-0.607	0.201	-3.02	0.007367
0 <= '*' < 0.001 < '' < 0.01 < '*' < 0.05
(Paramètre de dispersion pour une Quasi-Poisson family: 1.704) Déviance totale : 68.58 pour 19 degrés de liberté Déviance résiduelle : 41.44 pour 18 degrés de liberté AIC: NA - Nombre d’itérations de la fonction de score de Fisher : 5

# Modèle GLM poisson plus complet (avec plus de variables indépendantes)
dune_glm2 <- glm(data = dune, Lolipere ~ A1 +  Management + Moisture,
  family = poisson)
summary_(dune_glm2) |> tabularise()


$\log ({ E( \operatorname{Lolipere} ) }) = \alpha + \beta_{1}(\operatorname{A1}) + \beta_{2}(\operatorname{Management}_{\operatorname{culture}}) + \beta_{3}(\operatorname{Moisture}_{\operatorname{.L}}) + \beta_{4}(\operatorname{Moisture}_{\operatorname{.Q}}) + \beta_{5}(\operatorname{Moisture}_{\operatorname{.C}})$
Terme	Valeur estimée	Ecart type	Valeur de z	Valeur de p
$\alpha$	0.454	1.002	0.453	0.6508
$\beta_{1}$	-0.344	0.159	-2.167	0.0302
$\beta_{2}$	1.921	0.734	2.616	0.0089
$\beta_{3}$	-1.066	0.431	-2.475	0.0133
$\beta_{4}$	0.153	0.460	0.333	0.7394
$\beta_{5}$	0.732	0.517	1.416	0.1567
0 <= '*' < 0.001 < '' < 0.01 < '*' < 0.05
(Paramètre de dispersion pour une Poisson family: 1) Déviance totale : 68.58 pour 19 degrés de liberté Déviance résiduelle : 14.17 pour 14 degrés de liberté AIC: 66.48 - Nombre d’itérations de la fonction de score de Fisher : 6

# Comparaison des deux modèles imbriqués par test Chi-carré
anova(dune_glm, dune_glm2, test = "Chisq") |> tabularise()


Modèle	Ddl des résidus	Déviance résiduelle	Ddl	Déviance	Valeur de p
Lolipere ~ A1	18	41.4
Lolipere ~ A1 + Management + Moisture	14	14.2	4	27.3	1.75·10-05	***
0 <= '*' < 0.001 < '' < 0.01 < '*' < 0.05

# Exemple de prédictions en utilisant un modèle GLM
new_sites <- dtbl_rows(
 ~A1, ~Moisture, ~Management,
   3,         1, "conservation",
   3,         1, "culture",
   5,         2, "conservation",
   5,         2, "culture"
) %>.%
  mutate(., Moisture   = ordered(Moisture, levels = c(1, 2, 4, 5)),
            Management = factor(Management))
# Probabilité d'abondance de L. perenne en ces sites
new_sites$pred <- predict(dune_glm2,
  newdata = new_sites, type = "response")
new_sites

## # A tibble: 4 × 4
##      A1 Moisture Management    pred
##   <dbl> <ord>    <fct>        <dbl>
## 1     3 1        conservation 1.05 
## 2     3 1        culture      7.18 
## 3     5 2        conservation 0.543
## 4     5 2        culture      3.71

# Graphique des observations avec distinction par gestion et humidité
chart(data = dune, Lolipere ~ A1 %size=% Moisture %col=% Management) +
  geom_point(alpha = 0.7) +
  labs(x = "Épaisseur de la couche A1 [cm]",
       y = "Abondance [nbr de plants/quadrat]")

GLM binomiale avec proportions : maturation d’ovocytes

# Préparation des données de maturation d'ovocytes par dose d'hypoxanthine
ovo <- dtbl_rows(
  ~hypo, ~mat, ~tot,
      4,    0,   32,
      3,    3,   23,
      2,   12,   24,
      1,   24,   32,
    0.5,   26,   29,
   0.25,   28,   30,
      0,   35,   35
) %>.%
  mutate(., prop = mat/tot)
# Graphique général
chart(data = ovo, prop ~ hypo) +
  geom_point() +
  labs(x = "Hypoxanthine [µM]",
       y = "Fraction d'ovocytes matures")

# GLM binomiale avec proportions ; notez bien l'utilisation de weights=... ici !
ovo_glm <- glm(data = ovo, prop ~ hypo,
  family = binomial, weights = ovo$tot)
summary_(ovo_glm) |> tabularise()

Modèle linéaire généralisé
$\log\left[ \frac { P( \operatorname{prop} = \operatorname{1} ) }{ 1 - P( \operatorname{prop} = \operatorname{1} ) } \right] = \alpha + \beta_{}(\operatorname{hypo})$
Terme	Valeur estimée	Ecart type	Valeur de z	Valeur de p
$\alpha$	3.27	0.407	8.03	9.54·10-16
$\beta_{}$	-1.78	0.228	-7.83	4.80·10-15
0 <= '*' < 0.001 < '' < 0.01 < '*' < 0.05
(Paramètre de dispersion pour une Binomial family: 1) Déviance totale : 150.3 pour 6 degrés de liberté Déviance résiduelle : 5.576 pour 5 degrés de liberté AIC: 25.16 - Nombre d’itérations de la fonction de score de Fisher : 4

# Autre formulation de GLM binomiale ne nécessitant pas d'indiquer weights=...
ovo_glm_bis <- glm(data = ovo, cbind(mat, tot) ~ hypo,
  family = binomial)
summary_(ovo_glm) |> tabularise()


$\log\left[ \frac { P( \operatorname{prop} = \operatorname{1} ) }{ 1 - P( \operatorname{prop} = \operatorname{1} ) } \right] = \alpha + \beta_{}(\operatorname{hypo})$
Terme	Valeur estimée	Ecart type	Valeur de z	Valeur de p
$\alpha$	3.27	0.407	8.03	9.54·10-16
$\beta_{}$	-1.78	0.228	-7.83	4.80·10-15
0 <= '*' < 0.001 < '' < 0.01 < '*' < 0.05
(Paramètre de dispersion pour une Binomial family: 1) Déviance totale : 150.3 pour 6 degrés de liberté Déviance résiduelle : 5.576 pour 5 degrés de liberté AIC: 25.16 - Nombre d’itérations de la fonction de score de Fisher : 4

# Graphique de notre modèle GLM binomiale avec proportions
chart(data = ovo, prop ~ hypo) +
  geom_point() +
  stat_smooth(method = "glm", method.args = list(family = binomial),
    formula = y ~ x, se = FALSE) +
  labs(x = "Hypoxanthine [µM]",
       y = "Fraction d'ovocytes matures")

GLM binomiale avec variable binaire : acariens

# Préparation des données sur les acariens
mite <- sbind_cols(
  read("mite", package = "vegan") |> sselect(Oribatl1),
  read("mite.env", package = "vegan") |> sselect(WatrCont, Topo)
) %>.%
  smutate(., Oribatl1 = case_when(
    Oribatl1 == 0 ~ "absent",
    .default = "present") |> factor())
skimr::skim(mite)

Data summary
Name	mite
Number of rows	70
Number of columns	3
_______________________
Column type frequency:
factor	2
numeric	1
________________________
Group variables	None