Salem Nkunda Nyisingize

Les modèles additifs généralisés : splines, pénalisation et lissage

Comment laisser les données choisir la forme d'une relation — sans surajuster. Animations, équations, et code Python que vous pouvez exécuter directement dans cette page.

Pour recruteurs (lecture rapide)

Un GAM remplace les coefficients d'une régression linéaire par des fonctions lisses estimées, tout en restant interprétable. C'est le compromis entre la régression linéaire (trop rigide) et les méthodes de boîte noire (ininterprétables). Ce que ce sujet démontre :

  • Régularisation : la pénalité de rugosité est une régression ridge dans une base de splines.
  • Compromis biais-variance : contrôlé continûment par un seul paramètre λ, mesurable via les degrés de liberté effectifs.
  • Sélection de modèle : λ choisi par GCV ou REML, sans validation croisée coûteuse.
  • Algèbre linéaire numérique : matrice chapeau, trace, systèmes pénalisés bien conditionnés.
Splines pénalisées Régularisation GCV / REML mgcv Modèles interprétables

Pourquoi les GAM ?

Une régression linéaire impose une hypothèse forte : l'effet de chaque variable est une droite. Or la plupart des phénomènes réels n'ont pas cette politesse. La température influence le rendement d'une culture jusqu'à un optimum, puis le dégrade. La dose d'un médicament agit non linéairement. Le trafic d'un site dépend de l'heure selon une courbe, pas une pente.

On pourrait ajouter des termes \(x^2\), \(x^3\), ou découper en catégories — mais ces choix sont arbitraires et fragiles. Le modèle additif généralisé pose la question autrement : et si on laissait les données décider de la forme de la courbe, en pénalisant seulement son irrégularité ?

01 / 06Ouverture de la série : splines, pénalisation, lissage.

L'idée en une phrase. Un GAM s'écrit \(\;g(\mu) = \beta_0 + \sum_j f_j(x_j)\;\) : chaque \(f_j\) est une fonction lisse estimée, dont la flexibilité est contrôlée par une pénalité — et non fixée d'avance.

Du linéaire à l'additif

Partons du modèle linéaire \(y = \beta_0 + \beta_1 x + \varepsilon\). Ajusté par moindres carrés sur des données courbes, il laisse une structure évidente dans les résidus : ils ne sont pas répartis au hasard autour de la droite, ils suivent un motif. C'est le signe qu'une information systématique n'a pas été captée.

02 / 06 Les résidus de la droite (orange) fondent lorsqu'on remplace \(\beta_1 x\) par une fonction lisse \(f(x)\).

Modèle linéaire
Forme imposée
\[y_i = \beta_0 + \beta_1 x_i + \varepsilon_i\]

Un seul degré de liberté pour la forme : la pente. Biais élevé si la vraie relation est courbe.

Modèle additif
Forme apprise
\[y_i = \beta_0 + f(x_i) + \varepsilon_i\]

\(f\) appartient à un espace de fonctions lisses. Il reste à définir cet espace — et à empêcher \(f\) de passer par tous les points.

Un lisse est une somme de bases

Pour estimer une fonction, on l'écrit comme combinaison linéaire de fonctions élémentaires connues, appelées fonctions de base :

Représentation en base
\[f(x) = \sum_{k=1}^{K} \beta_k\, B_k(x)\]

Ce geste est décisif : l'estimation d'une fonction (objet de dimension infinie) devient l'estimation d'un vecteur de coefficients \(\beta\) — un problème de régression linéaire ordinaire dans la matrice de design \(B\). Les \(B_k\) usuelles sont des B-splines : des bosses polynomiales par morceaux, chacune non nulle seulement sur un intervalle court. Cette localité est ce qui rend l'ajustement stable.

03 / 06 Les bases apparaissent, sont pondérées par leur \(\beta_k\), puis sommées en une seule courbe.

Construire la base, à la main

Le code ci-dessous s'exécute réellement dans votre navigateur (Python compilé en WebAssembly). Modifiez les valeurs et relancez : le premier lancement prend quelques secondes le temps de charger Python.

PYTHON · NUMPY
Pourquoi la « partition de l'unité » compte. Les B-splines somment à 1 en tout point. Une constante est donc exactement représentable, et les coefficients \(\beta_k\) se lisent comme la hauteur locale de la courbe — ce qui rend la pénalité de la section suivante très naturelle.

La pénalisation et le paramètre λ

Avec \(K\) grand, les moindres carrés ordinaires font passer la courbe par presque tous les points : surajustement. La solution n'est pas de réduire \(K\) — choix discret et brutal — mais d'ajouter un terme de pénalité qui coûte cher aux courbes trop sinueuses :

Critère pénalisé
\[\hat{\beta} = \arg\min_{\beta}\; \underbrace{\|y - B\beta\|^2}_{\text{fidélité aux données}} \;+\; \lambda \underbrace{\beta^{\top} S \beta}_{\text{rugosité}}\]

La matrice \(S\) mesure l'irrégularité. Dans l'approche P-spline d'Eilers et Marx, on prend \(S = D^{\top}D\) où \(D\) est l'opérateur de différences secondes : pénaliser \(\beta^\top S\beta\) revient à pénaliser les variations brusques entre coefficients voisins. La solution est explicite :

Solution pénalisée (ridge dans la base)
\[\hat{\beta} = (B^{\top}B + \lambda S)^{-1} B^{\top} y\]

04 / 06 λ balayé de \(10^{-4}\) à \(10^{4}\) : la courbe se raidit et les degrés de liberté effectifs fondent de ≈ 29 à ≈ 2.

Degrés de liberté effectifs

Comment quantifier la complexité réellement utilisée ? Par la trace de la matrice chapeau \(H = B(B^{\top}B + \lambda S)^{-1}B^{\top}\), qui projette \(y\) sur l'ajustement :

Degrés de liberté effectifs
\[\mathrm{edf}(\lambda) = \operatorname{tr}(H)\]

C'est la généralisation continue du « nombre de paramètres ». Quand \(\lambda \to 0\), \(\mathrm{edf} \to K\) ; quand \(\lambda \to \infty\), la pénalité de différences secondes force une droite, et \(\mathrm{edf} \to 2\). Entre les deux, la complexité varie continûment — c'est exactement ce que montre l'animation.

PYTHON · NUMPY

Choisir λ automatiquement (GCV)

Laisser l'utilisateur régler λ à la main serait un aveu d'échec. La validation croisée « leave-one-out » a ici une forme close remarquable : pour un lisseur linéaire, on peut calculer l'erreur de validation croisée sans jamais refaire l'ajustement. Sa version stabilisée est la validation croisée généralisée :

Critère GCV
\[\mathrm{GCV}(\lambda) = \frac{n\,\|y - Hy\|^2}{\left(n - \operatorname{tr}(H)\right)^2}\]

Le numérateur récompense l'ajustement, le dénominateur pénalise la complexité : plus \(\operatorname{tr}(H)\) est grand, plus le dénominateur rétrécit et plus le score se dégrade. On minimise sur une grille logarithmique de λ. En pratique, mgcv utilise par défaut le REML, plus résistant au sous-lissage, mais le GCV reste le plus simple à comprendre et à coder.

PYTHON · NUMPY
Sur ces données, le GCV sélectionne \(\lambda \approx 0{,}93\), soit \(\mathrm{edf} \approx 11{,}5\) sur les 29 coefficients disponibles. Le modèle n'utilise donc qu'un tiers de sa flexibilité : la pénalité a fait le travail de sélection de modèle, en continu.

Structure additive

Tout ce qui précède portait sur une seule variable. La force du modèle vient de l'assemblage : on somme plusieurs lisses, un par covariable, chacun avec sa propre pénalité \(\lambda_j\).

Modèle additif
\[y_i = \beta_0 + f_1(x_{1i}) + f_2(x_{2i}) + \dots + f_p(x_{pi}) + \varepsilon_i\]

05 / 06 Trois effets partiels, lisibles séparément — l'atout majeur des GAM face aux boîtes noires.

C'est là que se joue l'interprétabilité. Un gradient boosting ou un réseau de neurones capterait aussi ces relations, mais sans permettre de tracer « l'effet de \(x_2\), toutes choses égales par ailleurs ». Ici, chaque \(f_j\) est un objet graphique autonome, accompagné de son intervalle de confiance. Le prix à payer : les interactions ne sont pas captées par défaut — il faut les demander explicitement (voir Limites).

Le « G » : fonction de lien

Jusqu'ici la réponse était continue et le bruit gaussien. Le « généralisé » lève cette contrainte, exactement comme les GLM étendent la régression linéaire : on autorise une réponse de la famille exponentielle (Poisson pour des comptages, binomiale pour des proportions, Gamma pour des durées), reliée au prédicteur par une fonction de lien \(g\).

Modèle additif généralisé
\[g(\mu_i) = \beta_0 + \sum_{j} f_j(x_{ji}), \qquad y_i \sim \text{famille exponentielle},\;\; \mathbb{E}[y_i] = \mu_i\]

06 / 06 Identité, logarithme, logit : le lien adapte le modèle à la nature de la réponse.

Nature de la réponseFamilleLien canonique
Continue, symétriqueGaussienne\(g(\mu) = \mu\)
Comptage (0, 1, 2, …)Poisson\(g(\mu) = \log \mu\)
Binaire ou proportionBinomiale\(g(\mu) = \log\frac{\mu}{1-\mu}\)
Positive, asymétriqueGamma\(g(\mu) = 1/\mu\) (ou log)
Comptage surdisperséBinomiale négative\(g(\mu) = \log \mu\)

L'estimation se fait alors par moindres carrés pondérés itératifs pénalisés (P-IRLS) : on itère le schéma linéaire précédent sur une variable de travail, en remettant à jour les poids.

En pratique avec mgcv

Tout ce que nous avons codé à la main est industrialisé dans le paquet R mgcv, écrit par Simon Wood lui-même. Le code ci-dessous n'est pas exécutable dans le navigateur (il requiert R), mais il est le point d'arrivée naturel de cet article :

# --- R / mgcv --- library(mgcv) # un lisse par covariable ; k = dimension de base, bs = type de spline m <- gam(y ~ s(x1, k = 20, bs = "tp") + s(x2) + x3, family = gaussian, method = "REML", data = df) summary(m) # edf par terme, p-values approximatives, deviance expliquee plot(m, pages = 1, shade = TRUE) # les effets partiels, un par panneau gam.check(m) # diagnostics : k suffisant ? residus corrects ? # reponse de comptage, avec interaction 2D et effet cyclique m2 <- gam(n ~ s(temp, jour, bs = "tp") + s(mois, bs = "cc", k = 12), family = nb, method = "REML", data = df)
ArgumentRôleConseil
kDimension de la base (notre \(K\))Prendre généreux ; la pénalité fait le tri. Vérifier avec gam.check.
bs = "tp"Spline plaque minceDéfaut robuste, invariant par rotation, marche en dimension > 1.
bs = "cc"Spline cycliquePour mois, heure, angle — impose \(f(\text{début}) = f(\text{fin})\).
method = "REML"Sélection de λPréférable au GCV : moins sujet au sous-lissage.
select = TRUEPénalité supplémentairePermet à un terme d'être ramené exactement à zéro.

Limites et extensions

Hypothèse ou limiteSi problématique → solution
Effets purement additifsProduits tensoriels : te(x1, x2) pour une interaction lisse à échelles différentes
Observations indépendantesGAMM : effets aléatoires, corrélation temporelle ou spatiale
Seule la moyenne est modéliséeGAMLSS : lisses aussi sur la variance, l'asymétrie, l'aplatissement
Grandes bases de donnéesbam() : ajustement par blocs, parallélisable
p-values des termes lissesApproximatives (l'edf est estimé) — à lire avec prudence
Extrapolation hors du domaineDangereuse : la pénalité ne contraint rien au-delà des données observées
Le piège classique. Choisir \(k\) trop petit brise tout : la pénalité ne peut pas rattraper une base incapable de représenter la vraie fonction. Choisir \(k\) trop grand ne coûte que du temps de calcul, jamais du biais. En cas de doute, augmentez \(k\) et laissez λ décider.

Références

  • Wood, S. N. (2017). Generalized Additive Models: An Introduction with R, 2e éd., Chapman & Hall/CRC — la référence, et la source de cet article.
  • Eilers, P. H. C. & Marx, B. D. (1996). Flexible smoothing with B-splines and penalties, Statistical Science — l'article fondateur des P-splines.
  • Hastie, T. & Tibshirani, R. (1990). Generalized Additive Models, Chapman & Hall — l'origine du modèle.
  • Ressources en ligne : mgcv sur CRANGAMs in R (Noam Ross)Pyodide (moteur Python de cette page)

Les animations de cet article ont été produites avec Manim ; le code source des six scènes est disponible sur GitHub.