How to Perform a Nested ANOVA in R (Step-by-Step)

L’analyse de variance, introduite par le statisticien et généticien Ronald Aylmer Fisher dans les années 1920, constitue l’un des piliers méthodologiques fondamentaux de l’inférence statistique moderne. Dans sa forme canonique unidirectionnelle ou factorielle factorisée, elle postule que les différents niveaux de traitement sont combinés de manière orthogonale ou croisée, garantissant ainsi que chaque modalité d’un facteur puisse interagir de façon symétrique avec l’ensemble des modalités des autres facteurs du modèle. Cependant, la complexité intrinsèque des dispositifs expérimentaux réels — qu’il s’agisse des sciences agronomiques, biomédicales, écologiques ou psychosociales — impose fréquemment des contraintes structurelles où l’orthogonalité parfaite s’avère impossible ou non pertinente sur le plan pratique.

Lorsque les unités expérimentales ou les sous-groupes d’observation sont structurellement subordonnés à un niveau hiérarchique supérieur spécifique, le modèle linéaire classique doit s’adapter pour capturer adéquatement cette dépendance stochastique sous-jacente. L’analyse de variance emboîtée (en anglais, Nested ANOVA ou Hierarchical ANOVA) offre le cadre inférentiel idoine permettant de décomposer rigoureusement les sources de variation lorsque les niveaux d’un facteur secondaire n’ont de sens, d’existence physique ou d’assignation expérimentale qu’à l’intérieur d’une modalité donnée du facteur principal. Ignorer cette hiérarchie conduit inévitablement à ce que Stuart Hurlbert (1984) a formellement théorisé sous le terme de pseudo-réplication, une erreur d’inférence majeure provoquant un gonflement drastique du taux d’erreur de première espèce (α) et des conclusions fallacieuses.

Le présent guide a pour vocation de décortiquer, avec un niveau de rigueur académique sans compromis, l’ensemble des fondements théoriques, mathématiques et informatiques présidant à la mise en œuvre de l’ANOVA emboîtée au sein de l’écosystème R. À travers un continuum didactique allant de l’algèbre sous-jacente des carrés moyens attendus jusqu’aux extensions par les modèles mixtes linéaires (Linear Mixed Models), ce manuel exhaustif permettra aux chercheurs, doctorants et biostatisticiens de concevoir des protocoles robustes, d’exécuter des analyses irréprochables et de formaliser leurs résultats conformément aux standards les plus exigeants de la communication scientifique.

1. Introduction conceptuelle à l’ANOVA emboîtée et aux plans hiérarchiques

Définition théorique du plan factoriel emboîté

Sur le plan de l’architecture expérimentale, un plan factoriel est dit emboîté (ou hiérarchisé) lorsqu’un facteur secondaire, traditionnellement noté B, voit ses différentes modalités distribuées de telle manière que chaque niveau spécifique de B n’apparaît et n’est évalué que dans une et une seule modalité du facteur principal, traditionnellement désigné par A. Contrairement à un plan factoriel croisé où chaque niveau de B traverse l’intégralité des niveaux de A (permettant ainsi l’évaluation d’un produit tensoriel d’effets et l’estimation de termes d’interaction de type A × B), le plan emboîté interdit ontologiquement l’existence d’une telle interaction.

Dans cette configuration géométrique, les modalités du facteur sous-jacent B ne partagent aucune identité commune d’un niveau de A à un autre. Par exemple, si le facteur principal A représente trois protocoles thérapeutiques distincts et que le facteur secondaire B correspond à des hôpitaux assignés à ces protocoles, l’hôpital numéro 1 sous le traitement A n’a aucun lien organique, géographique ou statistique avec l’hôpital numéro 1 sous le traitement B. L’étiquetage des niveaux du facteur inférieur est donc purement arbitraire et conventionnel. On symbolise mathématiquement cette relation de dépendance unilatérale sous la notation formelle B(A), signifiant « le facteur B est emboîté dans le facteur A ».

L’ANOVA emboîtée constitue ainsi l’outil analytique par excellence pour partitionner la variance globale d’une variable réponse continue en plusieurs composantes étanches : la variation attribuable aux dissemblances entre les groupes principaux, la variation attribuable à l’hétérogénéité des sous-unités au sein de ces mêmes groupes, et enfin l’erreur résiduelle pure (ou micro-variation d’échantillonnage) constatée entre les répétitions au sein de chaque sous-unité.

Pertinence et cas d’usage en recherche expérimentale

L’omniprésence des structures hiérarchiques dans la nature et dans les dispositifs d’expérimentation confère à l’ANOVA emboîtée un champ d’application considérable. Dans les sciences biologiques et environnementales, l’échantillonnage s’effectue quasi systématiquement selon des gradients spatiaux ou temporels multiniveaux. Considérons l’évaluation de la toxicité d’un xénobiotique sur la biodiversité benthique : un écologue sélectionnera plusieurs bassins versants (facteur principal), échantillonnera plusieurs cours d’eau au sein de chaque bassin (facteur secondaire emboîté), prélèvera plusieurs carottes de sédiments par cours d’eau (facteur tertiaire emboîté) et mesurera l’abondance d’un taxon indicateur sur plusieurs aliquotes par carotte. Tenter d’analyser de telles données via une ANOVA factorielle classique ou en agrégeant prématurément les données conduirait soit à un effondrement de la puissance statistique, soit à une violation du postulat d’indépendance des observations.

En médecine translationnelle et en pharmacologie clinique, les essais précliniques recourent massivement à ces structures. Dans les tests de criblage de molécules, des portées de rongeurs peuvent être soumises à différentes doses d’un principe actif (niveau 1), plusieurs animaux sont sélectionnés par portée (niveau 2 emboîté), et de multiples coupes histologiques ou biopsies d’organes sont analysées pour chaque individu (niveau 3 d’échantillonnage technique). De façon analogue, en psychométrie et sciences de l’éducation, les apprenants sont regroupés au sein de classes scolaires, elles-mêmes intégrées au sein d’établissements territoriaux soumis à des programmes pédagogiques distincts.

Dans l’ensemble de ces situations, l’ANOVA emboîtée permet non seulement de tester l’effet réel du traitement principal en neutralisant statistiquement le « bruit » ou la variabilité des échelons intermédiaires, mais elle offre également l’opportunité d’estimer les composantes de la variance (variance components). Cela permet d’identifier à quel niveau hiérarchique réside l’instabilité majeure du système expérimental, guidant ainsi l’optimisation des protocoles d’échantillonnage futurs.

Présentation du cas d’étude fil conducteur

Afin d’ancrer de manière pragmatique les développements mathématiques et la programmation sous R, nous établissons un fil conducteur expérimental standardisé issu de la physiologie végétale et de l’agronomie de précision. Imaginons un consortium de recherche agronomique évaluant l’efficacité de trois formulations fertilisantes innovantes (le facteur principal Engrais, comportant 3 modalités fixes : « Contrôle », « Organo-Minéral », « Bio-Stimulant ») sur le rendement en biomasse sèche de plants de soja (Glycine max).

Pour chaque type de fertilisant, l’étude sélectionne aléatoirement 4 exploitations agricoles distinctes (le facteur secondaire emboîté Exploitation, comportant 4 niveaux par modalité d’engrais, soit 12 exploitations distinctes au total). Il est capital de souligner qu’aucune exploitation ne teste plus d’un type d’engrais : l’exploitation numéro 1 recevant le « Bio-Stimulant » n’a aucune parenté avec l’exploitation numéro 1 recevant le « Contrôle ». Au sein de chaque exploitation agricole, 6 parcelles expérimentales homogènes sont délimitées et récoltées, fournissant les observations unitaires de biomasse exprimées en quintaux par hectare (q/ha).

Ce protocole constitue un plan emboîté à deux facteurs hiérarchisés équilibré (balanced two-stage nested design), totalisant :

• 3 niveaux pour le facteur principal (a = 3) ;
• 4 exploitations par niveau de fertilisant (b = 4), soit un total de 3 × 4 = 12 exploitations ;
• 6 observations indépendantes par exploitation (n = 6) ;
• Un échantillon global de N = 3 × 4 × 6 = 72 observations expérimentales.

L’objectif de cette étude sera double : déterminer formellement si les formulations d’engrais modulent significativement le rendement moyen en biomasse, tout en quantifiant la variabilité introduite par l’hétérogénéité intrinsèque entre exploitations agricoles traitées de façon identique.

2. Fondements mathématiques et décomposition de la variance

L’équation formelle du modèle linéaire emboîté

Sur le plan de la théorie statistique linéaire générale, le modèle de l’ANOVA emboîtée à deux facteurs s’exprime sous la forme paramétrique suivante :

Y_ijk = μ + α_i + β_j(i) + ε_ijk

Dans cette formulation, chaque terme est rigoureusement défini pour caractériser les strates du dispositif expérimental :

• Y_ijk représente la valeur quantitative observée de la variable réponse pour la k-ième unité expérimentale (k = 1, 2, …, n) appartenant à la j-ième sous-unité (j = 1, 2, …, b) emboîtée dans la i-ième modalité du traitement principal (i = 1, 2, …, a).
• μ désigne la moyenne générale théorique (ou grand mean) de la population parente sous-jacente.
• α_i formalise l’effet principal différentiel du i-ième niveau du facteur supérieur A. Si le facteur A est à effets fixes, cet effet est soumis à la contrainte d’identifiabilité usuelle ∑ α_i = 0.
• β_j(i) matérialise l’effet propre de la j-ième modalité du facteur secondaire B conditionné à son appartenance exclusive à la i-ième modalité de A. L’indice parentésé (i) stipule l’emboîtement structurel. Si le facteur B est considéré comme un facteur aléatoire (cas le plus fréquent en pratique, où les exploitations sont un échantillon d’un univers plus vaste), on postule que β_j(i) suit une loi normale centrée : β_j(i) ~ N(0, σ²_β).
• ε_ijk constitue le terme d’erreur stochastique résiduelle non explicable par la structure hiérarchique du modèle. Les erreurs résiduelles sont postulées indépendantes et identiquement distribuées selon une loi normale : ε_ijk ~ N(0, σ²).

Partitionnement des sommes des carrés et degrés de liberté

L’essence de l’ANOVA réside dans l’orthogonalité de ses projecteurs vectoriels au sens de l’algèbre hilbertienne. La variabilité totale des observations autour de la moyenne globale, exprimée par la Somme des Carrés Totale (SST), se décompose sans ambiguïté en une somme de trois formes quadratiques indépendantes :

SST = SSA + SSB(A) + SSE

Chacune de ces composantes répond à une définition algébrique exacte :

• Somme des Carrés du facteur principal (SSA) : elle capture l’écart pondéré des moyennes de chaque modalité principale Ȳ_i·· par rapport à la moyenne générale Ȳ_··· :
  SSA = b × n × ∑ (Ȳ_i·· – Ȳ_···)², avec df_A = a – 1 degrés de liberté.
• Somme des Carrés du facteur emboîté (SSB(A)) : elle agrège les déviations des moyennes des sous-groupes Ȳ_ij· vis-à-vis de la moyenne de leur modalité principale respective Ȳ_i·· :
  SSB(A) = n × ∑∑ (Ȳ_ij· – Ȳ_i··)², avec df_B(A) = a(b – 1) degrés de liberté.
• Somme des Carrés Résiduelle (SSE) : elle correspond à la variabilité intrinsèque des répétitions individuelles Y_ijk autour de la moyenne de leur propre sous-unité Ȳ_ij· :
  SSE = ∑∑∑ (Y_ijk – Ȳ_ij·)², avec df_E = a × b × (n – 1) degrés de liberté.

Le nombre total de degrés de liberté vérifie l’identité fondamentale d’additivité arithmétique :

df_Total = (a – 1) + a(b – 1) + ab(n – 1) = abn – 1 = N – 1

La division de chaque somme des carrés par son degré de liberté respectif fournit les Carrés Moyens (Mean Squares : MSA, MSB(A), et MSE), indispensables à l’établissement des statistiques inférentielles de test.

Détermination des espérances des carrés moyens et ratios F

Le choix critique du dénominateur dans le calcul de la statistique F de Fisher constitue la complexité nodale des plans emboîtés. L’identification du bon dénominateur ne relève pas d’une convention arbitraire, mais découle strictement de la théorie des Espérances des Carrés Moyens (Expected Mean Squares ou EMS), popularisée par les règles de Cornfield et Tukey.

Dans un modèle mixte classique à deux niveaux où le facteur principal A est fixe (les engrais spécifiquement choisis) et le facteur emboîté B(A) est aléatoire (les exploitations agricoles représentant un échantillon représentatif tiré au sort d’une population d’exploitations), les espérances mathématiques des carrés moyens se formulent ainsi :

• E(MSE) = σ²
• E(MSB(A)) = σ² + n σ²_β
• E(MSA) = σ² + n σ²_β + [b n / (a – 1)] × ∑ α_i²

Pour tester l’hypothèse nulle d’absence d’effet du facteur principal H₀: α₁ = … = α_a = 0 (ce qui équivaut à ∑ α_i² = 0), nous devons construire un ratio F dont l’espérance sous l’hypothèse nulle tend rigoureusement vers 1. Si nous divisons MSA par l’erreur pure MSE (comme le ferait naïvement une routine d’ANOVA non informée), nous constatons que sous H₀, le ratio E(MSA) / E(MSE) vaut (σ² + n σ²_β) / σ², ce qui excède strictement 1 dès lors que la variance inter-exploitations σ²_β est non nulle. Ce choix erroné gonfle de manière désastreuse le taux d’erreur de première espèce.

Par conséquent, la seule construction valide impose d’évaluer l’effet du facteur principal par rapport à la variance de l’échelon emboîté :

F_A = MSA / MSB(A), distribué selon une loi de Fisher à df₁ = a – 1 et df₂ = a(b – 1) degrés de liberté.

En revanche, pour tester l’hétérogénéité des sous-unités emboîtées (H₀: σ²_β = 0), la comparaison canonique se fait contre l’erreur résiduelle :

F_B(A) = MSB(A) / MSE, distribué selon une loi de Fisher à df₁ = a(b – 1) et df₂ = ab(n – 1) degrés de liberté.

3. Distinction cruciale : Facteurs emboîtés versus Facteurs croisés

Critères géométriques et topologiques de discrimination

L’attribution erronée d’un statut croisé à un facteur emboîté — ou réciproquement — constitue une source majeure d’invalidation statistique dans les manuscrits scientifiques soumis aux revues à comité de lecture. La frontière entre ces deux paradigmes repose sur des propriétés d’invariance et de connectivité topologique dans l’espace des données.

Dans un plan factoriel croisé (crossed design, noté A × B), toute modalité b_j du facteur B conserve une entité intrinsèque et stable à travers toutes les modalités du facteur A. Cela implique que l’on peut interroger l’effet différentiel de b_j indépendamment de A, et réciproquement. Typiquement, si les engrais étaient croisés avec des types d’irrigation (Goutte-à-goutte, Aspersion, Gravitaire), chaque type d’irrigation serait testé sous chaque engrais. La matrice d’incidence factorielle est alors saturée : chaque cellule du croisement (i, j) contient des observations.

À l’inverse, dans un plan emboîté (nested design, noté B(A)), il n’existe aucune permutabilité des niveaux de B entre les niveaux de A. Si nous réorganisons les étiquettes numériques du facteur emboîté en les renommant arbitrairement au sein de chaque groupe principal (par exemple, appeler l’exploitation 1 du groupe Bio-Stimulant « Alpha » et l’exploitation 1 du groupe Contrôle « Bêta »), la structure sous-jacente du modèle et les statistiques inférentielles demeurent strictement invariantes. Dans un plan croisé, une telle substitution détruirait l’estimation de l’interaction et altérerait l’interprétation des effets principaux.

Le critère de décision heuristique peut s’énoncer de la façon suivante : « Le niveau 1 du facteur B observé sous le niveau 1 de A partage-t-il une réalité concrète, biologique ou physique avec le niveau 1 du facteur B observé sous le niveau 2 de A ? » Si la réponse est strictement négative, le facteur B est fonctionnellement et irrévocablement emboîté dans A.

Erreurs fréquentes d’identification dans la littérature empirique

La confusion entre emboîtement et croisement induit des pathologies inférentielles lourdes de conséquences empiriques. La première pathologie, ubiquitaire dans les articles faisant l’objet de rétractations, est la pseudo-réplication sacrificielle ou simple. Elle survient lorsqu’un chercheur traite des observations répétées (les parcelles ou sous-échantillons) comme si elles constituaient des répliques entièrement indépendantes du facteur principal, en ignorant l’échelon emboîté intermédiaire. Dans notre cas d’étude, cela reviendrait à traiter les 24 parcelles associées à chaque engrais comme 24 répétitions aléatoires homogènes, ramenant artificiellement les degrés de liberté de l’erreur du test de A de a(b – 1) = 9 à ab(n – 1) = 60. Le dénominateur utilisé devient MSE au lieu de MSB(A), ce qui engendre des valeurs de p artificiellement microscopiques et des déclarations fallacieuses de significativité.

La seconde erreur d’identification réside dans l’introduction injustifiée d’un terme d’interaction dans le modèle statistique. Certains analystes paramètrent par mégarde leur logiciel en spécifiant une syntaxe factorielle saturée (Y ~ Engrais * Exploitation). Or, comme chaque exploitation n’est observée qu’au sein d’un seul engrais, la grande majorité des cellules de la table de contingence bidirectionnelle contiennent un effectif nul. Les degrés de liberté du terme d’interaction Engrais × Exploitation et de l’effet principal de l’exploitation deviennent mathématiquement indissociables (aliased ou colinéaires). L’algorithme linéaire soustrait alors des degrés de liberté de manière incohérente ou renvoie des matrices de rang non plein avec des coefficients non estimables (marqués NA sous R).

La troisième dérive méthodologique réside dans l’agrégation prématurée par calcul de moyennes : pour contourner la complexité de l’emboîtement, certains chercheurs réduisent l’ensemble des mesures d’une exploitation à une seule moyenne par exploitation, avant de conduire une ANOVA à un facteur à 12 observations. Bien que cette approche protège du risque de première espèce sur l’effet principal, elle sacrifie l’intégralité de l’information statistique relative à la variabilité intra-exploitation et interdit toute estimation formelle des composantes de variance.

4. Conditions d’application et postulats statistiques préalables

Distribution normale des résidus à chaque niveau hiérarchique

L’inférence paramétrique exacte fondée sur les distributions F de Snedecor-Fisher repose de manière indéfectible sur la normalité des distributions sous-jacentes. Dans un schéma hiérarchique emboîté, cette hypothèse de normalité ne doit pas être postulée uniquement sur l’ensemble aggloméré de la variable brute, mais doit être scrupuleusement vérifiée à chaque niveau de l’architecture stochastique.

Au niveau le plus fin, les résidus stochastiques d’échantillonnage ε_ijk doivent suivre une distribution normale gaussienne de moyenne nulle et d’écart-type résiduel constant σ :

ε_ijk ~ i.i.d. N(0, σ²)

Tout écart substantiel vis-à-vis de cette distribution gaussienne (asymétrie marquée ou kurtosis excessif) fragilise l’exactitude des valeurs de p calculées pour le ratio MSB(A) / MSE.

De surcroît, lorsque le facteur emboîté est conceptualisé comme aléatoire, l’effet propre de ces sous-unités β_j(i) doit lui-même être distribué selon une loi normale centrée d’écart-type σ_β :

β_j(i) ~ i.i.d. N(0, σ²_β)

Dans les dispositifs expérimentaux où le nombre de sous-unités b par modalité principale est modeste (par exemple b = 3 ou 4), le contrôle formel de la normalité de ces effets aléatoires par des tests de normalité devient délicat en raison d’un manque de puissance. Il relève alors de la responsabilité méthodologique du chercheur de s’assurer de la plausibilité de cette hypothèse par une inspection minutieuse des distributions empiriques et par la vérification de l’absence de valeurs extrêmes ou aberrantes au niveau agrégé.

Homogénéité des variances (homoscédasticité)

Le postulat d’homoscédasticité impose que la variance des résidus soit rigoureusement stable à travers l’ensemble des strates du protocole. Dans une ANOVA emboîtée, cette condition d’égalité des variances se décline sous une double exigence fondamentale :

Premièrement, la variance résiduelle intra-sous-unités (intra-exploitations dans notre modèle) doit être invariante d’un groupe à un autre :

Var(ε_ijk) = σ² ∀ i, j

Si certaines exploitations agricoles manifestent une dispersion intrinsèque de leurs parcelles considérablement plus large que d’autres (en raison, par exemple, d’une granulométrie pédologique hétérogène), le terme MSE aggloméré cesse d’être un estimateur non biaisé et efficace de l’erreur commune.

Deuxièmement, la variance inter-sous-unités conditionnelle doit être homogène d’un traitement principal à l’autre :

Var(β_j(i)) = σ²_β ∀ i

Cela signifie que l’hétérogénéité constatée entre les exploitations recevant l’engrais Contrôle doit être du même ordre de grandeur que l’hétérogénéité constatée entre les exploitations traitées par le Bio-Stimulant. Si un traitement engendre une instabilité phénoménologique se traduisant par un éclatement de la variance inter-sites, le test F conventionnel régissant l’effet principal (F = MSA / MSB(A)) devient asymétrique et perd son calibrage théorique nominal de première espèce.

Indépendance conditionnelle des observations

L’indépendance statistique des erreurs résiduelles résume l’hypothèse la plus critique et la plus sensible aux distorsions inférentielles de l’analyse de variance. Contrairement aux violations mineures de normalité auxquelles le ratio F est modérément robuste en présence d’effectifs équilibrés (grâce au théorème central limite), la présence d’autocorrélation résiduelle altère dramatiquement les variances d’échantillonnage.

Dans un plan emboîté, les observations Y_ijk ne sont explicitement pas mutuellement indépendantes dans l’absolu, puisqu’il existe une corrélation intra-classe (Intraclass Correlation Coefficient, noté ICC ou ρ) induite par le regroupement au sein d’une même sous-unité :

ρ = σ²_β / (σ²_β + σ²)

Le rôle mathématique de l’ANOVA emboîtée est précisément de modéliser explicitement cette covariance intragroupe. Dès lors, le postulat d’indépendance ne s’applique pas aux données brutes, mais stipule que, conditionnellement à l’inclusion de l’effet hiérarchique β_j(i), les résidus résiduels finaux ε_ijk doivent présenter une covariance strictement nulle :

Cov(ε_ijk, ε_i’j’k’) = 0 dès que (i, j, k) ≠ (i’, j’, k’)

Sur le plan expérimental, cette condition impose une randomisation scrupuleuse des parcelles au sein de chaque exploitation, l’absence de transferts physiques de matières entre parcelles adjacentes (effets de bordure neutralisés par des zones tampons adéquates), et l’indépendance des mesures instrumentales dans le temps.

5. Étape 1 : Préparation de l’environnement R et simulation du jeu de données

Configuration de l’environnement et chargement des bibliothèques

Pour mener à bien une analyse d’ANOVA emboîtée robuste et conforme aux standards de l’analytique moderne, nous mobilisons un ensemble d’outils issus de l’écosystème CRAN. L’environnement doit être configuré pour assurer une traçabilité totale et une isolation fonctionnelle des traitements.

Nous faisons appel aux extensions suivantes :

• tidyverse (incluant dplyr, tidyr, purrr, et ggplot2) pour la manipulation expressive des structures de données (data frames ou tibbles) et la grammaire graphique avancée.
• car (Companion to Applied Regression), développé par John Fox, pour le calcul des sommes des carrés de Type II/III et les tests formels d’homoscédasticité de Levene.
• emmeans (Estimated Marginal Means) pour l’inférence a posteriori, l’estimation des moyennes ajustées par les moindres carrés et la conduite de contrastes orthogonaux et tests post-hoc conservateurs.
• lme4 et lmerTest pour la transition vers les modèles linéaires mixtes et le calcul exact des approximations de Satterthwaite et de Kenward-Roger pour les degrés de liberté.
• performance et see (de la suite easystats) pour l’évaluation diagnostique automatisée et globale des postulats sous-jacents aux modèles statistiques.

Pour charger ces bibliothèques dans votre session de travail, exécutez les instructions suivantes :

Syntaxe R de configuration initiale :

options(contrasts = c("contr.sum", "contr.poly"))
library(tidyverse)
library(car)
library(emmeans)
library(lme4)
library(lmerTest)
library(performance)
set.seed(42)

La commande fixant les contrastes orthogonaux (options(contrasts = ...)) est fondamentale : par défaut, R utilise le codage de traitement (contr.treatment), particulièrement adapté aux régressions simples, mais source d’asymétries problématiques lors de l’estimation de modèles factoriels non orthogonaux ou lors de la décomposition des sommes des carrés de Type III.

Création programmatique du jeu de données expérimental

Afin d’illustrer chaque concept statistique sur un substrat tangible et reproductible, nous simulons un ensemble de données fidèlement calibré sur notre cas fil conducteur agronomique. La simulation algorithmique permet de contrôler exactement les paramètres de la population sous-jacente (les moyennes théoriques des engrais, la variance inter-exploitations et la variance résiduelle intra-parcelles), offrant ainsi une clarté didactique optimale pour évaluer la validité de l’estimation.

Nous modélisons les paramètres suivants :

• Moyenne générale de biomasse : μ = 50.0 q/ha.
• Effets différentiels des engrais : α_Contrôle = -5.0 q/ha, α_{Organo-Minéral} = +1.0 q/ha, α_{Bio-Stimulant} = +4.0 q/ha.
• Variance inter-exploitations : σ²_β = 4.0 (soit un écart-type σ_β = 2.0 q/ha).
• Variance résiduelle pure : σ² = 2.25 (soit un écart-type d’erreur d’échantillonnage σ = 1.5 q/ha).

Script de génération de données :

# Paramètres de simulation
n_engrais <- 3
n_fermes_par_engrais <- 4
n_parcelles_par_ferme <- 6

engrais_niveaux <- c("Controle", "Organo_Mineral", "Bio_Stimulant")
engrais_effets <- c(-5, 1, 4)

# Construction du squelette expérimental
donnees_agronomie <- expand.grid(
  Parcelle = 1:n_parcelles_par_ferme,
  Ferme = 1:n_fermes_par_engrais,
  Engrais = factor(engrais_niveaux, levels = engrais_niveaux)
) %>%
  mutate(
    # Assignation d'un identifiant unique de ferme (emboîtement explicite)
    Ferme_ID = paste(Engrais, Ferme, sep = "_"),
    Ferme_ID = factor(Ferme_ID)
  )

# Simulation des effets aléatoires de ferme et du terme résiduel
effets_fermes <- rnorm(n_engrais * n_fermes_par_engrais, mean = 0, sd = 2.0)
names(effets_fermes) <- levels(donnees_agronomie$Ferme_ID)

donnees_agronomie <- donnees_agronomie %>%
  mutate(
    Effet_Engrais = case_when(
      Engrais == "Controle" ~ engrais_effets[1],
      Engrais == "Organo_Mineral" ~ engrais_effets[2],
      Engrais == "Bio_Stimulant" ~ engrais_effets[3]
    ),
    Effet_Ferme = effets_fermes[Ferme_ID],
    Erreur_Residuelle = rnorm(n(), mean = 0, sd = 1.5),
    Biomasse = 50 + Effet_Engrais + Effet_Ferme + Erreur_Residuelle
  ) %>%
  select(Engrais, Ferme, Ferme_ID, Parcelle, Biomasse)

Formatage et vérification des types de données

Une vigilance impérative doit être accordée au formatage des variables catégorielles. Dans les environnements statistiques comme R, l’inclusion d’une variable numérique dans une formule de modèle linéaire force implicitement son interprétation sous forme de covariable continue (régression polynomiale ou linéaire simple), annihilant complètement la logique partitionnelle de l’ANOVA.

Il est donc essentiel de s’assurer que les facteurs Engrais, Ferme et Ferme_ID sont formellement instanciés comme des objets de classe factor. De surcroît, la question de l’encodage de la variable emboîtée (Ferme versus Ferme_ID) est cruciale. Si nous conservons une variable Ferme contenant uniquement les valeurs 1, 2, 3, 4 dupliquées pour chaque engrais, la structure syntaxique de R pourrait être trompeuse sans une spécification explicite de l’opérateur d’emboîtement.

Nous validons la structure du jeu de données par les commandes d’inspection standardisées :

str(donnees_agronomie)
head(donnees_agronomie, 10)

Le tableau généré possède 72 lignes et 5 colonnes. La variable réponse quantitative Biomasse est de type numeric, tandis que Engrais, Ferme et Ferme_ID sont explicitement déclarés en facteurs catégoriels nominaux. L’environnement expérimental virtuel est désormais calibré pour débuter l’exploration descriptive.

6. Étape 2 : Exploration préliminaire et visualisation descriptive des données

Agrégation statistique descriptive multivariée

Avant d’ajuster le moindre modèle linéaire formel, une exploration méticuleuse de la structure empirique des données s’avère indispensable. Elle permet de détecter d’éventuelles anomalies de saisie, d’apprécier visuellement l’ampleur des écarts entre traitements et d’évaluer la variance relative entre les sous-groupes.

En employant les verbes d’action de dplyr, nous calculons les métriques de tendance centrale (moyenne arithmétique) et de dispersion (écart-type, erreur standard de la moyenne et intervalle de confiance à 95 %) à l’échelle globale de chaque engrais, mais également à l’échelon plus fin de chaque combinaison Engrais × Ferme :

# Statistiques descriptives au niveau des exploitations emboîtées
resume_fermes <- donnees_agronomie %>%
  group_by(Engrais, Ferme_ID) %>%
  summarise(
    n = n(),
    Moyenne = mean(Biomasse),
    Ecart_Type = sd(Biomasse),
    Erreur_Standard = Ecart_Type / sqrt(n),
    IC_95_Inf = Moyenne - qt(0.975, df = n - 1) * Erreur_Standard,
    IC_95_Sup = Moyenne + qt(0.975, df = n - 1) * Erreur_Standard,
    .groups = "drop"
  )

# Statistiques globales au niveau des engrais
resume_engrais <- donnees_agronomie %>%
  group_by(Engrais) %>%
  summarise(
    N_Total = n(),
    N_Fermes = n_distinct(Ferme_ID),
    Moyenne_Globale = mean(Biomasse),
    Ecart_Type_Global = sd(Biomasse),
    .groups = "drop"
  )

L’examen de ces tables synthétiques met en lumière la hiérarchie des variabilités : les moyennes par exploitation oscillent substantiellement autour de la moyenne globale de leur engrais respectif, confirmant la présence tangible d’un effet ferme d’origine stochastique, tandis que la dispersion intra-exploitation (évaluée par l’écart-type résiduel au sein de chaque ferme) présente une stabilité satisfaisante.

Visualisation graphique hiérarchisée avec ggplot2

La traduction visuelle d’un dispositif emboîté requiert un affichage multicouche capable de simultanément révéler la dispersion des observations unitaires, la variance inter-sous-unités et le profil moyen des groupes principaux. Le progiciel ggplot2 offre les primitives géométriques idéales pour orchestrer cette représentation synthétique.

Nous concevons une figure combinant les points individuels de chaque parcelle (avec une dispersion horizontale contrôlée via geom_jitter), des boîtes à moustaches (boxplots) dépeignant l’hétérogénéité des exploitations, et des marqueurs proéminents symbolisant la moyenne marginale globale de chaque traitement d’engrais :

graphique_emboite <- ggplot(donnees_agronomie, aes(x = Ferme_ID, y = Biomasse, fill = Engrais)) +
  geom_boxplot(alpha = 0.5, outlier.shape = NA, color = "grey30") +
  geom_jitter(width = 0.2, alpha = 0.7, size = 1.8, aes(color = Engrais)) +
  facet_wrap(~ Engrais, scales = "free_x", space = "free_x") +
  stat_summary(fun = mean, geom = "point", shape = 23, size = 3.5, fill = "white", color = "black") +
  labs(
    title = "Distribution hiérarchisée des rendements de biomasse",
    subtitle = "Décomposition par type d'engrais et par exploitation agricole emboîtée",
    x = "Identifiant de l'exploitation agricole (emboîtée)",
    y = "Biomasse sèche (q/ha)",
    caption = "Les losanges blancs matérialisent les moyennes par exploitation"
  ) +
  theme_bw(base_size = 12) +
  theme(
    legend.position = "none",
    axis.text.x = element_text(angle = 45, hjust = 1),
    strip.background = element_rect(fill = "grey90", color = "grey50"),
    strip.text = element_text(face = "bold")
  )

print(graphique_emboite)

Ce type de visualisation permet d’appréhender d’un coup d’œil deux dimensions analytiques majeures : la séparation nette des distributions entre les modalités d’engrais (qui semble présager un effet principal significatif) et le décalage systématique des niveaux moyens d’une exploitation à l’autre au sein d’une même facette (qui matérialise visuellement la variance du facteur emboîté).

Détection préliminaire des points atypiques ou influents

Les données aberrantes (outliers) peuvent exercer une influence disproportionnée sur l’estimation des carrés moyens, en particulier au niveau des sous-unités emboîtées où les effectifs sont souvent modérés (ici, n = 6 parcelles). Une observation aberrante située au sein d’une seule parcelle peut fausser la moyenne locale de l’exploitation et gonfler artificiellement la somme des carrés intra-groupe ou inter-groupes.

Pour dépister ces valeurs avec rigueur, nous calculons les scores centrés réduits (Z-scores) ou les distances résiduelles par rapport à la médiane locale, au sein de chaque exploitation, en mobilisant la règle des écarts interquartiles (IQR) formalisée par John Tukey :

valeurs_atypiques <- donnees_agronomie %>%
  group_by(Ferme_ID) %>%
  mutate(
    Q1 = quantile(Biomasse, 0.25),
    Q3 = quantile(Biomasse, 0.75),
    IQR_val = Q3 - Q1,
    Est_Atypique = Biomasse < (Q1 - 1.5 * IQR_val) | Biomasse > (Q3 + 1.5 * IQR_val)
  ) %>%
  filter(Est_Atypique)

print(valeurs_atypiques)

Si aucune observation n’est extraite par ce filtre, le jeu de données présente une régularité distributionnelle optimale. Si des observations atypiques apparaissent, elles doivent être soumises à une inspection qualitative pour écarter toute erreur de retranscription, sans jamais être supprimées aveuglément sans justification agronomique ou métrologique irréfutable.

7. Étape 3 : Spécification et ajustement du modèle d’ANOVA emboîtée dans R

Syntaxe de modélisation avec la fonction aov()

La formulation des modèles hiérarchisés sous R repose sur un ensemble d’opérateurs syntaxiques précis implémentés dans la formule symbolique du modèle. Deux opérateurs fondamentaux coexistent pour stipuler l’emboîtement au sein de la fonction standard aov() :

• L’opérateur de division hiérarchique / : L’expression A / B se traduit algorithmiquement par A + B %in% A. Elle indique formellement que le facteur B est emboîté dans le facteur A.
• L’opérateur explicite d’inclusion %in% : L’écriture A + B %in% A (ou A + A:B) déclare expressément l’effet principal de A combiné à l’effet de B restreint aux niveaux de A.

Considérons l’ajustement naïf suivant :

modele_aov_naif <- aov(Biomasse ~ Engrais / Ferme, data = donnees_agronomie)
summary(modele_aov_naif)

Bien que cette syntaxe produise une décomposition arithmétique exacte des sommes des carrés et des degrés de liberté, elle dissimule un écueil théorique fondamental : par défaut, la fonction summary.aov() teste toutes les sources de variation par rapport au carré moyen résiduel MSE (l’erreur résiduelle pure des parcelles). Comme démontré dans notre section consacrée aux fondements mathématiques, le ratio calculé pour le facteur principal devient alors F = MSA / MSE au lieu du ratio légitime F = MSA / MSB(A). Dès lors, le modèle standard est incomplet pour inférer sur l’effet fixe lorsque les sous-unités sont aléatoires.

Spécification de l’effet aléatoire via le terme Error()

Pour forcer l’algorithme d’ANOVA classique de R à sélectionner le bon dénominateur dans le calcul du test F de Fisher, il est impératif d’intégrer le terme spécialisé Error() au sein de la formule passée à la fonction aov(). Ce terme Error instruit R sur la structure de stratification des erreurs (multistratum error design).

La syntaxe canonique appropriée s’écrit de la manière suivante :

modele_aov_correct <- aov(Biomasse ~ Engrais + Error(Engrais/Ferme), data = donnees_agronomie)
summary(modele_aov_correct)

Alternativement, si l’identifiant des exploitations est déjà codé de manière unique (notre variable Ferme_ID), la structure de l’erreur inter-fermes peut être déclarée directement au sein du terme d’erreur :

modele_aov_correct_v2 <- aov(Biomasse ~ Engrais + Error(Ferme_ID), data = donnees_agronomie)
summary(modele_aov_correct_v2)

En analysant la sortie générée par R, on observe que le tableau est désormais partitionné en strates d’erreur étanches :

• Error: Ferme_ID (Strate intermédiaire inter-exploitations) : Cette strate évalue l’effet du facteur principal Engrais en utilisant explicitement le carré moyen résiduel de la strate Ferme_ID (qui n’est autre que MSB(A)) comme dénominateur de Fisher. Les degrés de liberté du dénominateur sont rigoureusement fixés à 9.
• Error: Within (Strate intra-exploitation / résiduelle) : Cette strate quantifie la variabilité résiduelle des parcelles au sein des exploitations (60 degrés de liberté), servant de dénominateur uniquement pour tester d’éventuels facteurs intra-unités.

Comparaison fonctionnelle entre aov() et lm()

Il est instructif de disséquer la relation sous-jacente qui lie la fonction aov() à la fonction générale de régression linéaire lm(). En réalité, aov() n’est qu’une surcouche (wrapper) syntaxique autour de lm() : elle convertit les contrastes en décompositions d’ANOVA et modifie la classe de l’objet résultant pour orienter les méthodes d’impression génériques vers un tableau de variance.

Toutefois, une divergence fonctionnelle majeure sépare les deux approches dès lors que l’on manipule des structures hiérarchiques :

• La fonction lm() traite l’intégralité des prédicteurs comme des effets fixes déterministes. Elle ajuste l’équation matricielle Y = Xβ + ε par la méthode des moindres carrés ordinaires (OLS). Par conséquent, lm(Biomasse ~ Engrais / Ferme, data = donnees_agronomie) estime des coefficients de régression spécifiques pour chaque exploitation, mais elle est intrinsèquement incapable de scinder la matrice de variance-covariance des erreurs en plusieurs strates aléatoires.
• La fonction aov() couplée à Error() déploie l’algorithme des moindres carrés généralisés multistrates (multistratum analysis). Elle projette le vecteur réponse orthogonalement sur les différents sous-espaces linéaires définis par les termes d’erreur, garantissant une inférence exacte sur les ratios de Fisher, mais au détriment de l’accès direct aux fonctions habituelles d’extraction des résidus individuels ou d’estimations marginales directes.

Cette limitation intrinsèque de aov() constitue l’un des arguments méthodologiques majeurs plaidant pour l’adoption des modèles mixtes contemporains basés sur la maximisation de la vraisemblance restreinte (REML).

8. Étape 4 : Interprétation rigoureuse du tableau d’ANOVA et des tests F

Décryptage des sources de variation du tableau de résultats

L’interprétation de la sortie d’un modèle d’ANOVA emboîtée correctement paramétré exige une lecture méthodique des composantes du tableau de variance. Examinons la structure typique renvoyée par l’instruction summary(modele_aov_correct) :

Strate 1 : Error: Ferme_ID

          Df Sum Sq Mean Sq F value   Pr(>F)    
Engrais    2 1014.2  507.10   18.44 0.000676 ***
Residuals  9  247.5   27.50                     
---

Strate 2 : Error: Within

          Df Sum Sq Mean Sq F value Pr(>F)
Residuals 60  135.0    2.25               

Ce tableau révèle l’architecture inférentielle suivante :

• L’effet principal Engrais : Il possède a – 1 = 3 – 1 = 2 degrés de liberté. La somme des carrés s’élève à 1014.2, ce qui engendre un carré moyen MSA = 507.10.
• L’effet de l’erreur inter-fermes (Residuals de la strate Ferme_ID) : Il cumule a(b – 1) = 3 × (4 – 1) = 9 degrés de liberté. Sa somme des carrés est de 247.5, soit un carré moyen MSB(A) = 27.50.
• Le ratio F de Fisher pour le facteur Engrais : F = MSA / MSB(A) = 507.10 / 27.50 = 18.44. Associé aux degrés de liberté (2, 9), il produit une valeur de p = 0.000676. L’effet de l’engrais est donc hautement significatif au seuil α = 0.001. Si l’on avait commis l’erreur de diviser MSA par le carré moyen résiduel de la strate Within (MSE = 2.25), le ratio F aurait atteint la valeur factice de 225.38, illustrant de façon saisissante le péril de la pseudo-réplication.
• L’erreur résiduelle pure (Residuals de la strate Within) : Elle cumule ab(n – 1) = 3 × 4 × (6 – 1) = 60 degrés de liberté et présente un carré moyen MSE = 2.25.

Le test de l’effet ferme s’obtient en évaluant manuellement le ratio :
F_{Ferme(Engrais)} = MSB(A) / MSE = 27.50 / 2.25 = 12.22. Avec 9 et 60 degrés de liberté, la probabilité critique associée 1 - pf(12.22, 9, 60) est largement inférieure à 10^-6, confirmant une variabilité substantielle et hautement significative entre les exploitations.

Calcul et interprétation des tailles d’effet

L’inférence statistique moderne impose de ne pas restreindre l’interprétation des résultats à la seule dichotomie signifiant/non signifiant dictée par la p-valeur, mais de chiffrer la magnitude des effets observés au moyen de mesures standardisées de taille d’effet (effect sizes). Dans le cadre d’une ANOVA emboîtée, deux métriques principales sont mobilisées : l’êta-carré partiel (η²_p) et l’oméga-carré partiel (ω²_p), ce dernier étant nettement moins biaisé pour les échantillons de taille finie.

L’êta-carré partiel pour le facteur principal A s’évalue par la formule :

η²_p = SSA / (SSA + SSB(A))

Notons avec minutie que le dénominateur comprend la somme des carrés de l’effet emboîté SSB(A) et non la somme des carrés de l’erreur pure SSE, puisque c’est SSB(A) qui constitue l’échelon de comparaison direct de A. En utilisant nos valeurs numériques :

η²_p = 1014.2 / (1014.2 + 247.5) = 1014.2 / 1261.7 ≈ 0.804

Pour obtenir un estimateur non biaisé des paramètres populationnels, on calcule l’oméga-carré partiel (ω²_p) :

ω²_p = [df_A × (MSA – MSB(A))] / [df_A × MSA + (N – df_A) × MSB(A)]

L’exécution programmatique de ces calculs sous R s’opère aisément :

# Extraction des composantes de variance
SS_A <- 1014.2
SS_BA <- 247.5
eta_carre_partiel <- SS_A / (SS_A + SS_BA)
cat("Eta-carré partiel de l'effet Engrais :", round(eta_carre_partiel, 4), "n")

Une valeur d’êta-carré partiel de 0.804 indique que 80.4 % de la variance attribuable à l’échelon des exploitations est directement expliquée par la typologie des formulations fertilisantes, ce qui traduit un impact agronomique particulièrement marqué selon les balises conventionnelles de Cohen.

Gestion des plans non équilibrés (Unbalanced Nested Designs)

Dans la pratique de la recherche empirique, la survenue de données manquantes, la mortalité d’individus ou la perte de parcelles expérimentales aboutit fréquemment à un plan déséquilibré (unbalanced design), dans lequel le nombre de sous-unités par traitement (b_i) ou le nombre de répliques par sous-unité (n_ij) n’est pas strictement constant.

Le déséquilibre introduit une rupture majeure : les sommes des carrés calculées par la décomposition canonique séquentielle (dite de Type I) deviennent dépendantes de l’ordre d’entrée des termes dans le modèle, et l’orthogonalité exacte des projecteurs vectoriels est perdue. De surcroît, les espérances des carrés moyens (EMS) deviennent des combinaisons linéaires complexes où les coefficients scalaires associés aux variances ne sont plus de simples constantes entières n ou b, mais des moyennes harmoniques pondérées.

Dans une telle conjoncture :

• L’usage de aov() avec le terme Error() devient formellement déconseillé, car l’algorithme ne gère pas de manière appropriée les décompositions de sommes des carrés de Type II ou III ajustées pour les déséquilibres hiérarchiques.
• Le recours aux Carrés Moyens Synthétiques de Satterthwaite (pseudo-F tests) ou, idéalement, aux Modèles Mixtes Linéaires (LMM) ajustés par le package lme4 s’impose comme la solution mathématique de référence pour maintenir la justesse de l’inférence.

9. Étape 5 : Diagnostic des résidus et vérification empirique des postulats

Évaluation formelle et graphique de la normalité résiduelle

La validation inférentielle d’une analyse de variance emboîtée requiert un contrôle scrupuleux du comportement stochastique des résidus. Pour extraire adéquatement ces résidus sous R sans distorsion, il est fortement conseillé de modéliser l’architecture complète soit via lm() (pour examiner les résidus de niveau 1), soit via un modèle mixte lmer() (permettant d’extraire à la fois les résidus d’échantillonnage et les effets aléatoires prédits de niveau 2).

Ajustons le modèle linéaire complet pour procéder au diagnostic résiduel :

modele_diagnostic <- lm(Biomasse ~ Engrais / Ferme_ID, data = donnees_agronomie)
residus_bruts <- residuals(modele_diagnostic)
residus_standardises <- rstandard(modele_diagnostic)

L’évaluation graphique de la normalité s’effectue traditionnellement par le biais d’un diagramme quantile-quantile (Q-Q plot) normal enrichi d’une bande de confiance à 95 %, complété par un histogramme de distribution :

# Diagnostic graphique de la normalité
par(mfrow = c(1, 2))
qqnorm(residus_standardises, pch = 19, col = "steelblue", main = "Q-Q Plot des Résidus")
qqline(residus_standardises, col = "firebrick", lwd = 2)
hist(residus_standardises, breaks = 12, col = "grey90", border = "grey30", freq = FALSE,
     main = "Distribution des Résidus", xlab = "Résidus standardisés")
curve(dnorm(x), add = TRUE, col = "firebrick", lwd = 2)
par(mfrow = c(1, 1))

Sur le plan inférentiel formel, l’adéquation à la loi normale peut être corroborée par le test de Shapiro-Wilk :

test_shapiro <- shapiro.test(residus_bruts)
print(test_shapiro)

Une p-valeur supérieure au seuil usuel α = 0.05 atteste que l’on ne peut rejeter l’hypothèse nulle de distribution gaussienne des perturbations résiduelles.

Vérification de l’homoscédasticité résiduelle

Le contrôle de la constance de la variance s’opère par l’analyse conjointe des graphiques de dispersion des résidus et des tests formels de non-constance de variance. Le graphique des résidus standardisés en fonction des valeurs ajustées (Fitted values) constitue la clé de voûte de cette évaluation.

plot(fitted(modele_diagnostic), residus_standardises,
     pch = 19, col = "darkgreen",
     xlab = "Valeurs ajustées (q/ha)",
     ylab = "Résidus standardisés",
     main = "Résidus vs Valeurs Ajustées")
abline(h = 0, col = "grey40", lty = 2)
abline(h = c(-2, 2), col = "firebrick", lty = 3)

Une dispersion symétrique des points formant une bande d’épaisseur uniforme sans motif géométrique caractéristique (forme d’entonnoir, de pavillon de trompette ou de parabole) valide l’hypothèse de variance constante.

Pour étayer cette observation visuelle par un test robuste et non paramétrique, nous appliquons le test de Levene modifié (Brown-Forsythe), calculé à la médiane, disponible dans le package car :

test_levene <- leveneTest(Biomasse ~ Engrais, data = donnees_agronomie)
print(test_levene)

L’obtention d’une statistique de Levene non significative atteste que les variances intra-traitements sont statistiquement indiscernables, consolidant ainsi la validité théorique du dénominateur d’erreur employé dans les tests de Fisher.

Remédiations méthodologiques en cas de violation des postulats

Lorsque les données expérimentales dérogent de manière substantielle aux postulats fondamentaux, plusieurs voies de remédiation doivent être examinées méthodiquement :

• Transformations de stabilisation de la variance : En présence d’une relation hétéroscédastique entre la moyenne et la variance (fréquente pour des données de comptage ou des concentrations biologiques), l’application d’une transformation logarithmique log(Y), racine carrée √Y, ou de la famille paramétrique de Box-Cox permet souvent de rétablir conjointement la normalité et l’homoscédasticité.
• Modèles Linéaires Généralisés Mixtes (GLMM) : Si la variable dépendante est intrinsèquement discrète (par exemple, nombre d’insectes ravageurs par parcelle suivant une loi de Poisson ou binomiale négative), la conversion artificielle par ANOVA emboîtée est caduque. L’emploi d’un GLMM via glmer() modélise explicitement la distribution réelle de la variable.
• Modélisation explicite de l’hétérogénéité des variances : Si l’hétéroscédasticité est inhérente aux groupes et ne peut être résorbée par transformation, l’utilisation de la fonction gls() ou lme() du package nlme permet de spécifier une structure de variance par groupe (par exemple, weights = varIdent(form = ~ 1 | Engrais)), autorisant chaque groupe à posséder sa propre variance d’erreur.

10. Étape 6 : Analyses post-hoc et comparaisons multiples pour facteurs hiérarchisés

Adaptation du test Tukey HSD aux architectures emboîtées

Dès lors que l’ANOVA emboîtée rejette formellement l’hypothèse nulle globale d’égalité des moyennes pour un facteur principal comportant au moins trois niveaux (k ≥ 3), la question se pose d’identifier quelles paires spécifiques de traitements manifestent des écarts statistiquement significatifs. C’est l’apanage des procédures de comparaisons multiples post-hoc.

Néanmoins, l’application naïve de la fonction standard TukeyHSD() de R sur un objet aov emboîté conduit à une erreur d’inférence majeure : TukeyHSD() utilise par défaut les degrés de liberté résiduels globaux (60 df dans notre cas) et le carré moyen résiduel de l’erreur pure MSE pour calibrer la distribution de l’étendue studentisée (Studentized Range Distribution, q). Or, comme établi formellement, l’erreur d’échantillonnage appropriée pour comparer les moyennes des engrais est la variance inter-exploitations MSB(A), disposant uniquement de 9 degrés de liberté.

Ignorer cette distinction dénature le contrôle du taux d’erreur par famille (Family-Wise Error Rate, FWER), produisant des intervalles de confiance excessivement étroits et une avalanche de fausses découvertes.

Mise en œuvre avancée avec le package emmeans

Le package emmeans (Estimated Marginal Means), développé par Russell Lenth, représente l’étalon-or moderne pour conduire des analyses post-hoc rigoureuses sur des modèles hiérarchiques et mixtes. Il recalcule avec exactitude les moyennes marginales estimées, leur erreur standard conditionnelle et les degrés de liberté adéquats.

Pour exécuter ces contrastes sur notre modèle d’ANOVA emboîtée, nous ajustons le modèle mixte correspondant (ou le modèle aov multistrate) et appliquons les commandes suivantes :

# Ajustement via lmer pour interface transparente avec emmeans
modele_mixte_pour_posthoc <- lmer(Biomasse ~ Engrais + (1 | Ferme_ID), data = donnees_agronomie)

# Calcul des moyennes marginales estimées
emm_engrais <- emmeans(modele_mixte_pour_posthoc, ~ Engrais)
print(emm_engrais)

# Comparaisons par paires avec correction de Tukey-Kramer
comparaisons_posthoc <- pairs(emm_engrais, adjust = "tukey")
summary(comparaisons_posthoc)

L’inspection du tableau généré par summary(comparaisons_posthoc) montre que les tests utilisent explicitement les degrés de liberté ajustés selon la méthode de Kenward-Roger ou Satterthwaite (ici df = 9). Les estimations des contrastes (par exemple Bio_Stimulant – Controle) sont accompagnées d’erreurs standard intégrant à la fois la variabilité intra et inter-exploitations, et les valeurs de p résultent de la correction conservative de Tukey-Kramer.

Visualisation graphique des comparaisons post-hoc

La restitution visuelle des contrastes post-hoc bonifie substantiellement la lisibilité des publications scientifiques. emmeans s’intègre nativement à ggplot2 via sa fonction dédiée plot(), permettant de visualiser les intervalles de confiance à 95 % des moyennes ajustées, ainsi que des graphiques de comparaisons par paires :

# Graphique des moyennes marginales et intervalles de confiance
plot(emm_engrais, comparisons = TRUE) +
  theme_bw(base_size = 12) +
  labs(
    title = "Moyennes marginales estimées et contrastes post-hoc",
    subtitle = "Intervalles de confiance à 95 % et flèches de comparaison de Tukey",
    x = "Biomasse estimée (q/ha)",
    y = "Formulation d'engrais"
  )

Sur ce graphique, deux groupes dont les flèches de comparaison (comparison arrows) ne se chevauchent pas présentent une divergence statistiquement significative au seuil ajusté de 5 %. Cette clarté visuelle évite toute méprise quant au positionnement relatif des formulations fertilisantes innovantes par rapport au témoin agronomique.

11. Modèles mixtes linéaires (LMM) : L’alternative moderne avec lme4

Transition conceptuelle vers les modèles à effets mixtes

Bien que l’ANOVA emboîtée classique par moindres carrés offre un cadre intellectuel limpide et historique, la statistique appliquée contemporaine a largement opéré une transition vers le paradigme unifié des Modèles Mixtes Linéaires (Linear Mixed Models, LMM). Ce glissement méthodologique répond à des impératifs mathématiques majeurs que l’ANOVA traditionnelle ne peut satisfaire :

• Robustesse totale face aux plans déséquilibrés : Dans un LMM résolu par le critère de Vraisemblance Restreinte (REML), la présence d’effectifs inégaux par sous-unité n’altère en rien la convergence asymptotique et ne requiert pas d’artifices de calcul de carrés moyens synthétiques.
• Estimation directe des composantes de variance : Le modèle mixte ne se contente pas de tester un ratio d’hypothèse nulle ; il quantifie formellement la variance attribuable à chaque strate (σ²_β et σ²), permettant d’évaluer directement le coefficient de corrélation intra-classe (ρ).
• Flexibilité d’extension : Les modèles mixtes autorisent l’intégration simultanée de covariables continues (ANCOVA emboîtée), de structures de covariance temporelle ou spatiale complexe (autorégressive, spatiale gaussienne) et d’effets emboîtés profonds à trois niveaux ou plus sans complexification de la formulation d’analyse.

Implémentation pratique via lmer() du package lme4

Le package lme4, orchestré par Douglas Bates, Martin Mächler et Ben Bolker, constitue le moteur numérique le plus performant et le plus éprouvé sous R pour ajuster les modèles mixtes linéaires par l’algorithme REML.

Dans notre cas d’étude, la spécification d’un modèle mixte où le facteur principal Engrais est un effet fixe et le facteur emboîté Ferme_ID est un effet aléatoire intercept-seul (random intercept) s’écrit de façon concise :

# Spécification et ajustement du modèle linéaire mixte
modele_lmm <- lmer(Biomasse ~ Engrais + (1 | Ferme_ID), data = donnees_agronomie, REML = TRUE)
summary(modele_lmm)

Analysons scrupuleusement la structure de cette commande :

• Le terme Engrais désigne les prédicteurs à effets fixes modélisant l’espérance mathématique conditionnelle E(Y|X).
• La syntaxe (1 | Ferme_ID) spécifie un intercept aléatoire propre à chaque exploitation. Le « 1 » stipule que c’est la valeur moyenne de base (l’ordonnée à l’origine) qui varie d’une exploitation à l’autre selon une loi normale centrée.
• L’argument REML = TRUE engage la méthode de maximisation de la vraisemblance restreinte, éliminant le biais inhérent à l’estimation classique du maximum de vraisemblance (ML) sur les paramètres de variance lorsque les degrés de liberté fixes sont soustraits.

Inférence statistique et calcul des p-valeurs avec lmerTest

Les utilisateurs familiers de R remarquent immédiatement que le package lme4 d’origine ne fournit volontairement aucune valeur de p dans son tableau récapitulatif summary(). Ce choix d’ingénierie statistique, défendu par Douglas Bates, découle du fait que la définition exacte des degrés de liberté du dénominateur pour les ratios F dans les modèles mixtes n’a pas de solution analytique fermée générale dès lors que le plan n’est pas parfaitement équilibré.

Pour remédier à cette contrainte et satisfaire aux réquisits de publication académique, le package lmerTest vient enrichir la classe de l’objet lmerMod en calculant les degrés de liberté approximatifs par l’algorithme d’approximation de Satterthwaite ou la correction plus conservatrice de Kenward-Roger :

# Obtention du tableau d'ANOVA de Type III avec p-valeurs via Satterthwaite
anova_lmm <- anova(modele_lmm, type = "III")
print(anova_lmm)

Le tableau généré restitue :

• La statistique F de Snedecor (18.44).
• Les degrés de liberté du numérateur (NumDF = 2).
• Les degrés de liberté calculés pour le dénominateur (DenDF = 9.00).
• La probabilité critique exacte (Pr(>F) = 0.000676).

On constate une identité parfaite entre les résultats de l’ANOVA emboîtée multistrate classique via aov() et le modèle mixte lmer() en présence d’un plan parfaitement équilibré. Ce constat démontre sans équivoque que le modèle mixte subsume l’ANOVA emboîtée classique, tout en conférant une adaptabilité structurelle immensément supérieure.

12. Rédaction académique des résultats et directives de publication APA

Normes de présentation des résultats selon le format APA 7e édition

La restitution textuelle des analyses statistiques dans les revues académiques internationales doit impérativement se conformer aux prescriptions typographiques et normatives édictées par l’American Psychological Association (APA 7th Edition). Dans le cas d’une ANOVA emboîtée ou d’un modèle mixte équivalent, l’énoncé doit systématiquement mentionner :

• La formulation explicite du plan hiérarchique et la désignation des unités expérimentales et sous-unités emboîtées.
• La statistique F observée, accompagnée de ses degrés de liberté respectifs entre parenthèses : df du numérateur et df du dénominateur, séparés d’une virgule.
• La valeur exacte de probabilité critique (p-valeur rapportée à trois décimales, sans zéro initial si p < .001).
• Une mesure d’ampleur d’effet standardisée (η²_p ou ω²_p).
• Les résultats des comparaisons post-hoc (différences moyennes estimées, intervalles de confiance à 95 %, correction d’erreur globale appliquée).

Exemple de formulation académique rédigée :

« Une analyse de variance emboîtée a été conduite pour évaluer l’incidence des formulations d’engrais (Contrôle, Organo-Minéral, Bio-Stimulant) sur le rendement en biomasse sèche de plants de soja, en contrôlant la variabilité aléatoire issue des exploitations agricoles emboîtées dans chaque modalité de fertilisant (4 exploitations par engrais, 6 parcelles par exploitation, N = 72). L’analyse révèle un effet principal hautement significatif du type d’engrais, F(2, 9) = 18.44, p < .001, η²_p = .80, démontrant qu’une part prépondérante de la variance inter-exploitations est tributaire du protocole fertilisant appliqué. De surcroît, une hétérogénéité significative a été constatée entre les exploitations recevant un même engrais, F(9, 60) = 12.22, p < .001. Les comparaisons multiples post-hoc ajustées par la méthode de Tukey indiquent que la formulation Bio-Stimulant génère un rendement moyen significativement supérieur (M = 54.21 q/ha, ET = 1.62) par rapport au Contrôle (M = 44.87 q/ha, ET = 1.71, p < .001) et à la formulation Organo-Minérale (M = 50.94 q/ha, ET = 1.55, p = .021). »

Création d’un tableau synthétique d’ANOVA prêt pour publication

La standardisation visuelle des tables de variance est une exigence universelle de la publication scientifique. L’incorporation de sorties brutes de logiciels ou de captures d’écran de terminaux R est formellement proscrite. Le tableau doit présenter des bordures horizontales épurées (style booktabs) et des colonnes rigoureusement labellisées.

Le package rempsyc associé à flextable ou le package gtsummary permettent d’exporter directement des tableaux d’ANOVA au format Word ou LaTeX en parfaite conformité avec les standards APA 7 :

# Exemple de mise en forme programmatique d'un tableau synthétique
tableau_synthese <- data.frame(
  Source = c("Engrais (Fixe)", "Exploitation [Engrais] (Aléatoire)", "Erreur résiduelle (Parcelles)"),
  SS = c("1014.20", "247.50", "135.00"),
  df = c("2", "9", "60"),
  MS = c("507.10", "27.50", "2.25"),
  F_Ratio = c("18.44", "12.22", "—"),
  p_value = c("< .001", "< .001", "—"),
  eta_p2 = c(".80", "—", "—")
)

knitr::kable(tableau_synthese, col.names = c("Source de Variation", "Somme des Carrés", "ddl", "Carré Moyen", "F", "p", "η²p"),
             align = c("l", "r", "r", "r", "r", "r", "r"),
             caption = "Tableau 1. Décomposition de la variance du modèle d'ANOVA emboîtée pour le rendement en biomasse.")

Ce tableau isole avec netteté les strates hiérarchiques et indique explicitement la nature fixe ou aléatoire des effets considérés, assurant une intelligibilité immédiate pour les évaluateurs pairs (peer reviewers).

Recommandations éthiques et reproductibilité de la recherche

La crise de la reproductibilité qui secoue l’édifice des sciences empiriques depuis plus d’une décennie impose des protocoles d’intégrité sans faille. Dans le cadre des analyses hiérarchiques, la tentation de masquer l’emboîtement pour « gagner en degrés de liberté » et faire passer une p-valeur sous le seuil fatidique de .05 relève directement de la fraude méthodologique par p-hacking.

Pour garantir la transparence et l’auditabilité totale de vos démarches scientifiques, appliquez systématiquement les principes directeurs suivants :

• Dépôt des scripts et des données brutes : Déposez l’intégralité du script R nettoyé, commenté et autonome, ainsi que la matrice de données brutes anonymisée, sur des répertoires de science ouverte pérennes dotés d’un DOI (tels que le Center for Open Science (OSF) ou Zenodo).
• Conservation de l’environnement logiciel : Archivez l’information de version de R et de tous les packages chargés en consignant la sortie de l’instruction sessionInfo() à la fin de vos scripts, ou utilisez des outils de gestion d’environnement virtuel tels que renv pour geler les dépendances de calcul.
• Pré-enregistrement des protocoles (Pre-registration) : Si l’étude relève d’une démarche confirmatoire, spécifiez en amont la structure emboîtée du modèle, la règle d’attribution des effets fixes/aléatoires et le plan d’analyse post-hoc sur une plateforme dédiée avant la phase de collecte effective des observations.

Références

• Bates, D., Mächler, M., Bolker, B., & Walker, S. (2015). Fitting linear mixed-effects models using lme4. Journal of Statistical Software, 67(1), 1–48. https://doi.org/10.18637/jss.v067.i01
• Cornfield, J., & Tukey, J. W. (1956). Average values of mean squares in factorials. The Annals of Mathematical Statistics, 27(4), 907–949. https://doi.org/10.1214/aoms/1177728067
• Fisher, R. A. (1925). Statistical Methods for Research Workers. Oliver and Boyd.
• Fox, J., & Weisberg, S. (2019). An R Companion to Applied Regression (3rd ed.). SAGE Publications.
• Hurlbert, S. H. (1984). Pseudoreplication and the design of ecological field experiments. Ecological Monographs, 54(2), 187–211. https://doi.org/10.2307/1942661
• Kenward, M. G., & Roger, J. H. (1997). Small sample inference for fixed effects from restricted maximum likelihood. Biometrics, 53(3), 983–997. https://doi.org/10.2307/2533558
• Kuznetsova, A., Brockhoff, P. B., & Christensen, R. H. B. (2017). lmerTest package: Tests in linear mixed effects models. Journal of Statistical Software, 82(13), 1–26. https://doi.org/10.18637/jss.v082.i13
• Lenth, R. V. (2023). emmeans: Estimated Marginal Means, aka Least-Squares Means (R package version 1.8.5). CRAN. https://cran.r-project.org/package=emmeans
• Satterthwaite, F. E. (1946). An approximate distribution of estimates of variance components. Biometrics Bulletin, 2(6), 110–114. https://doi.org/10.2307/3002019
• Winer, B. J., Brown, D. R., & Michels, K. M. (1991). Statistical Principles in Experimental Design (3rd ed.). McGraw-Hill.