Publications du Laboratoire

de
Statistique et Probabilités

Composantes de la variabilité
Plans d’expériences

Jean-Marc Azaı̈s

Laboratoire de Statistique et Probabilités — UMR CNRS C55830

Université Paul Sabatier — 31062 – Toulouse cedex 4.
2
Chapitre 1

Modèles mixtes

1 Analyse de la variance multivariable

Dans cette section nous présentons l’analyse de la variance multivariable dans sa
présentation classique. En fait ses modèles peuvent se traiter comme un cas particulier
des modèles mixtes que nous allons décrire dans la suite de cet monographie.
L’analyse de la variance, la régression et le modèle linéaire en général, ont pour but
d’expliquer une variable quantitative. Cependant, dans de nombreux cas, l’observation est
multiple : on peut par exemple mesurer la longueur et la largeur d’une feuille, le rendement
d’une réaction chimique en plusieurs instants, etc... Une analyse variable par variable est
toujours possible, elle donne des estimations qui sont parfaitement exploitables. Pour ce
qui est des tests, la situation est plus complexe : chaque analyse uni-variable donne un
test différent pour une hypothèse comme par exemple la nullité de l’effet d’un facteur.
Ces tests sont peu puissants car il n’utilisent que l’information d’une seule variable et il
peuvent même être contradictoires. Il est donc très souhaitable de pouvoir regrouper les
tests sur les différentes variables.
Exemple : Pour apprécier et comparer le degré de pollution de rivières, on dose différents
métaux lourds dans les viscères de truites. On dispose par exemple de 6 observations :
Plomb dans le foie, plomb dans l’intestin, plomb dans le colon
Mercure dans le foie, mercure dans l’intestin, mercure dans le colon.

On suppose que l’on observe plusieurs truites par rivières et on note Y ij l’observation
sur la truite j de la rivière i. On se trouve dans une situation d’analyse de la variance à
un facteur à ceci près que l’observation se trouve dans IR 6 . On pose
Yij = ti + ij ; ti ∈ IR6 ; ij ∼ N (0, Σ) ; i = 1, I ; j = 1, ni ,
où Σ est une matrice 6x6 définie positive.
Comme nous l’avons déjà souligné, l’estimation peut parfaitement se faire variable par
variable, par contre le test de la nullité “l’effet rivière ” par exemple aura 6 réponses en
fonction de la variable utilisée.
Pour obtenir un test unique, nous allons paraphraser les formules de l’analyse de la
variance uni-variable en essayant de les généraliser au cas vectoriel.
-On peut toujours estimer ti par Yi. qui est le vecteur des moyennes. On définit ainsi

3
4 CHAPITRE 1. MODÈLES MIXTES

le vecteur des résidus associé à l’unité i, j.

ˆij = Yij − Yi. .
Comme il s’agit toujours d’un vecteur de taille 6 ( que l’on notera maintenant L pour
généraliser), la somme des carrés est remplacé par une matrice L × L. Si on note en
exposant les coordonnées d’un vecteur ou d’une matrice, on a
SCRk,l = ˆkij ˆlij , k = 1, L ; l = 1, L.
X

i,j

Dans le cas uni variable cette quantité suivait une loi σ 2 χ2 (n − I) (n = I × J est le nombre
d’unités). Dans le même état d’esprit on définit

Définition 1 (loi de Whishart de paramètres d et Σ) Soient d vecteurs indépendants

X1 , ..., Xd de loi N (0, Σ). On définit la matrice M des sommes de carrés et produits par
Xi (Xi )0
X
M=
i=1,d

où le 0 désigne la transposée. Cette matrice est composée de sommes de carrés qui suivent
de lois σ 2 χ2 (d) sur la diagonale et de sommes de produits à l’extérieur de la diagonale.
Par définition elle suit la loi de Whishart de paramètres d et Σ.

Par les mêmes techniques qu’en uni-variable, on montre que

SCR suit un loi de WhishartN − I, Σ.
On définit également la somme de carrés et produits associée au facteur :
(Yi. − Y.. )(Yi. − Y.. )0 .
X
SCF :=
i,j

Sous l’hypothèse H0 d’absence d’effet du facteur (rivière dans notre exemple), on montre
que
SCF suit une loi de WhishartI − 1, Σ.

Il reste à construire le test. Rappelons qu’en uni-variable, ce test est basé sur le rapport
SCF/(I − 1)
F̂ :=
SCR/(N − I)
qui suit la loi de Fisher F(I−1),(N −I) . En multivariable, d’un certain coté, la situation est
analogue puisque SCR/(N − I) est un estimateur de Σ et d’un autre coté, différente dans
le sens qu’il n’est pas possible de faire le quotient de deux matrices.
On pose W := SCR/(N − I) (comme within) matrice des variations intra-facteur
et B := SCF/(I − 1) (comme between) matrice de variations inter-facteur. Les tests en
multivariable sont multiples et basés sur les statistiques suivantes
|W |
– le λ de Wilks : |B+W | qui donne le test du rapport de vraisemblance,
– la trace de Pillai : T r(B((B + W )−1 )),
– la trace de Hotelling-Lawley : T r(BW −1 ),
– la plus grande valeur propre de Roy : PGVP (W −1 B).
Fort heureusement, dans la plupart des cas ces tests donnent des résultats cohérents. Les
lois connues et complexes sous H0 ont de très bonnes approximation par des lois de Fisher
qui sont utilisées le plus souvent par les logiciels pour calculer les niveaux de signification.
2. MODÈLES À EFFETS ALÉATOIRES ET MIXTES 5

Mise en œuvre SAS (cas équirépété)

Proc Glm data=sasuser.pollution;

Class riviere;
model Hgf Hgi Hgc Pbf Pbi Pbc = riviere;
means riviere/Tukey;
Manova h=riviere; run;

Ce modèle d’analyse de la variance multivariable peut apparaı̂tre comme un cas par-

ticulier de modèle mixte que nous allons présenter dans les section suivantes.

2 Modèles à effets aléatoires et mixtes

2.1 Facteurs croisés et hiérarchisés

Quand on considère deux facteurs en situation d’analyse de la variance, deux situations

peuvent se rencontrer :
– Chacun d’eux a un sens indépendamment de l’autre, c’est le cas le plus courant, les
facteurs sont dits croisés. Exemple de facteurs en général croisés :
variété × lieu
machine × opérateur
concentration × température (réaction chimique)
– Un des facteur, par exemple B ne signifie rien de concret tant que l’ on ne connaı̂t
pas l’indice associé de l’autre facteur (A). Dans ce cas le facteur B est hiérarchisé
au facteur A. C’est le plus souvent le cas quand le facteur B correspond à une
numérotation à l’intérieur des niveaux du facteur A. Considérons par exemple plu-
sieurs portées de lapereau, et supposons que dans chaque portée, on numérote chaque
lapereau de 1 au nombre de lapereaux dans la portée. Il ne peut y avoir d’effet propre
du numéro de lapereau car il n’y a aucun rapport entre les lapereaux qui porte un
même numéro.
Exemple de facteurs le plus souvent hiérarchisés

famille génétique / no de descendant

lapine / no de portée /no du lapin dans la portée.

Remarque Une autre alternative consisterait a ne pas numéroter le facteur hiérarchisé

a l’intérieur des niveaux du facteur hiérarchisant : dans l’exemple des lapereau, si le
première lapine donne naissance par exemple à 4 lapins, le premier lapereau de la seconde
portée serait numéroté 5. C’est alors un troisième type de relation entre deux facteurs :
l’information de l’un (lapine) est alors comprise dans l’information de l’autre. En fait,
pour des raisons d’organisation des logiciels statistiques cette second approche ne donne
pas des décompositions satisfaisantes : le faire sur un exemple pour s’en convaincre. En
conséquence on préfère souvent la première approche.
Quand on est en présence de facteurs hiérarchisés, il faut prendre soin de le déclarer
avec la syntaxe propre au logiciel que l’on utilise. En effet, il ne faut pas, dans ce cas, définir
un effet propre du facteur hiérarchisé. La décomposition est différente de celle du modèle
croisé dans le sens ou ce qui serait l’effet principal du facteur hiérarchisé est incorporé à
l’interaction.
6 CHAPITRE 1. MODÈLES MIXTES

Plus précisément on écrit

Yijk = µ + αi + γij + ijk

Où γij est l’effet du facteur hiérarchisé, avec les contraintes

i. αi = 0
P
i=1,I
P
ii. pour tout i j=1,J γij = 0.

2.2 Modèles mixtes équirépétés

Exemple :On désire comparer les aptitudes de 4 firmes à produire des insecticides
efficaces. Chaque firme produit de nombreux insecticides mais on échantillonne exactement
deux produits numérotés (1 et 2) par firme. Pour étudier les insecticides on utilise 24 boites
contenant 400 moustiques chacune . Chaque produit est introduit dans 3 boites prises
au hasard parmi les 24 et on compte le nombre de moustiques survivants au bout d’un
certain temps. On considère que le nombre 400 est suffisamment grand pour que l’aspect
binomial discret du problème puisse être oublié. Cela permet éviter d’utiliser un modèle
linéaire généralisé à effets aléatoires. On reconnaı̂t une structure ou l’effet insecticide est
hiérarchisé au facteur firme. Comme les produits ont été échantillonnés le facteur produit
est aléatoire.
On veut répondre aux questions : Y a t’il une différence entre les firmes ? Quelle
est la variabilité la plus importante : celle relative au choix du produit ou la variabilité
résiduelle ?
Pour l’instant nous allons voir comment répondre à la seconde question. On note
Yij,k l’observation dans la kième boite du jème produit de la firme k, i = 1, I(4) ; j =
1, J(2) ; k = 1, K(3) et on pose le modèle

Yijk = µ + fi + Pij + ijk

où les majuscules désignent des effets aléatoires. On suppose qu’il sont tous indépendants,
à l’intérieur d’un effet et d’un effet à l’autre et que les P ij , ijk sont centrés gaussiens
de variances respectivement σP2 et σ 2 . L’importance relative des effets va s’apprécier par
l’estimation des différentes variances. Ici on est dans un cas équirépété, l’estimation en
sera très simple. Il existe un estimateur qui est uniformément meilleur par des techniques
d’exhaustivité (pour un condition exacte pour qu’il en soit ainsi voir Coursol (1980)). On
définit les sommes de carrés résiduelle SCR et associé à l’effet produit SCP . Par utilisation
de l’orthogonalité du modèle il est facile de vérifier que

(Yijk − Yij. )2
X
SCR =
ijk

(Yij. − Yi.. )2 .
X
SCP =
ijk

On laisse à titre d’exercice de montrer que

E(SCP ) = I(J − 1)KσP2 + I(J − 1)σ 2 (1.1)

2
E(SCR) = IJ(K − 1)σ (1.2)
2. MODÈLES À EFFETS ALÉATOIRES ET MIXTES 7

On obtient un système triangulaire donc inversible qui donne un estimateur de σ P2 , σ 2 en

identifiant les sommes de carrés à leurs espérances, c’est à dire en résolvant le système

SCP = I(J − 1)K σ̂P2 + I(J − 1)σ̂ 2 (1.3)

2
SCR = IJ(K − 1)σ̂ (1.4)

Cela correspond à une méthode de moments. Comme le cas équirépété reste un cas par-
ticulier, nous ne détaillerons pas , mais on montre (Coursol 1980) que dans ce cas les
estimateurs ci dessus sont optimaux parmi les estimateurs sans biais.

2.3 Modèles mixtes généraux

Première remarque : la méthode précédente reste toujours applicable. Si on utilise

les sommes de carrés de type III pour l’identification, elle porte le nom de la Méthode
de Henderson III. Le calcul des coefficient dans (1.1)(1.2) est plus complexe, mais des
expressions matricielles peuvent être établies par les techniques que nous allons développer
plus loin, expressions qui peuvent être calculées par l’ordinateur. Cette méthode est parfois
très acceptable surtout pour les dispositifs peu déséquilibrés, et dans tous les cas elle fournit
de bon point de départ pour les méthode itératives que nous allons développer plus loin.
Au départ un modèle mixte est la donnée d’un vecteur Y qui peut s’écrire

Y = Xθ + Z1 β1 + ... + Zk−1 βk−1 + 

où θ est un vecteur de paramètres inconnus, les matrices X, Z 1 , ..., Zk sont des matrices
connues de dimension convenables et β 1 , ..., βk sont de vecteurs indépendants gaussiens
formés de coordonnées indépendantes et de même variance ( on note γ i la variance com-
mune des coordonnées de βi ). Le modèle est mixte dans le sens où il regroupe de effets
déterministes θ et des effets aléatoires. De la formule ci dessus on déduit facilement que

E(Y ) = Xθ

γi Zi Zi0 .
X
V ar(Y ) = Vγ :=
i=1,k

en posant Vk = Id. C’est cette forme que l’on va poser comme définition.

Définition 2 (modèles mixtes) Un vecteur aléatoire gaussien Y de taille n est dit suivre
un modèle mixte statistique si
E(Y ) = Xθ
X
V ar(Y ) = Vγ := γi Vi .
i=1,k

où X, V1 , ..., Vk sont de matrices connues, θ est un vecteur de paramètres inconnus variant
dans IRp , γ := (γ1 , ..., γk ) varie dans l’ensemble S := {γ : Vγ > 0}.

Ce modèle est noté : Y ∼ M M (X, V1 , ..., Vk ). Les (γ1 , ..., γk ) sont appelées composantes
de la variance. On pose Y = Xθ +  avec Var() = V γ .
8 CHAPITRE 1. MODÈLES MIXTES

2.4 Estimation des effets fixes

Pour simplifier on supposera X de plein rang, on a vu que l’on pouvait s’y ramener le
plus souvent. Si γ est connu, le modèle mixte se ramène a un modèle linéaire ordinaire et
l’estimateur optimal parmi les estimateurs linéaires sans biais est

θ̂ = (X 0 Vγ−1 X)−1 X 0 Vγ−1 Y. (1.5)

Démonstration : Par diagonalisation il existe une matrice T symétrique que l’on notera
−1/2
par la suite Vγ telle que

T T = Vγ−1 ; T Vγ T = Id

On pose donc Ỹ = T Y = T Xθ + T  = X̃θ + ˜. Il est facile de voir que V ar(˜

) = Id et que
donc Ỹ suit un modèle linéaire à variance connue. Dans ce modèle l’estimateur optimal
est
θ̂ = ((X̃)0 X̃)−1 (X̃)Ỹ = (X 0 Vγ−1 X)−1 X 0 Vγ−1 Y.
2
Notons que l’équation (1.5) est numériquement facile à résoudre même pour des modèles
comprenant plusieurs centaines de paramètres.
De la même manière si on considère une hypothèse linéaire sur le vecteur θ le test de
Fisher dans le modèle Ỹ = X̃θ + ˜ sera optimal. Quand γ ne sera plus connu mais estimé,
on estimera encore les effets fixes en utilisant les équations ci- dessus avec γ̂ comme s’il
s’agissait de la vraie valeur. Nous verrons plus loin que cela correspond à la solution du
maximum de vraisemblance gaussien. En conclusion le problème statistique essentiel dans
un modèle mixte est l’estimation des composantes de la variance.

2.5 Estimation par MIVQUE dans un modèle mixte

Supposons que nous voulons estimer une composante γ i de γ ou plus généralement une
combinaison linéaire
h = h0 γ
des composantes. h ∈ IRk donné. On cherche un estimateur
(i) invariant : ĥ(Y ) = ĥ(Y + Xθ), ∀θ ∈ IRp
(ii) quadratique : ĥ(Y ) = Y 0 BY où B est une matrice (il est basé sur des sommes de
carrés).
(iii) sans biais : E(ĥ(Y )) = h
(iv) de faible variance.

Proposition 1 (1) On peut toujours supposer B symétrique et on le fera tout le temps.

(2) L’estimateur est invariant ssi BX = 0 (B est maintenant symétrique)
(3) Si Y est gaussien et si BX = 0 , V ar(Y 0 BY ) = V ar(0 B) = 2T r(BVγ BVγ ).

La difficulté principale vient du fait que la variance de ĥ dépend en général de la valeur

de γ que l’on cherche à estimer. Il est donc sans espoir (sauf dans les cas équirépétés) de
chercher un estimateur optimal. On va donc se donner un a priori γ 0 et minimiser la
variance gaussienne au voisinage de γ 0 : T r(BVγ0 BVγ0 ).
2. MODÈLES À EFFETS ALÉATOIRES ET MIXTES 9

Théorème 1 Soit γ0 un a priori. Un estimateur vérifiant (i) (ii) (iii) et de γ 0 -variance

minimale est appelé MIVQUE (Minimum Variance Quadratic Unbiased Estimator). Il est
donné par les équations
k
0
Q0γ0 Vγ−1 δi Vi Vγ−1
X
ĥ(Y ) = Y 0 0
Q γ0 Y
i=i

avec
- Qγ0 := Id − X(X 0 Vγ−1
0
X)−1 X 0 Vγ−1
0
,
- δ solution de δ 0 Mγ0 = h
- Mγ0 = {T r(Vi Vγ−1
0
Qγ0 Vj Vγ−1
0
Qγ0 ), i = 1, k, j = 1, k}.
En particulier si la matrice Mγ0 est inversible, il existe un MIVQUE de toute compo-
sante.

MIVQUE vectoriel

Supposons Mγ0 inversible, comme toute combinaison est estimable, on obtient un

MIVQUE de chaque coordonnée γi . Pour ce faire, on cherche δ i tel que

M γ0 δ i = f i

où fi est le ième vecteur de la base canonique (avec un 1 en ième position). Cela implique
que
{δli , i = 1, k, j = 1, k} = (Mγ0 )−1
On définit le vecteur des sommes de carrés S = {Y 0 Q0γ0 Vγ−1
0
Vl Vγ−1
0
Qγ0 Y, l = 1, k}. On a
alors
γ̂ = (Mγ0 )−1 S.

2.6 Estimation par maximum de vraisemblance restreinte

On suppose que X est de plein rang, Y est un vecteur gaussien. La vraisemblance vaut

(2π)−n/2 |Vγ |−1/2 exp(−1/2 (Y − Xθ)0 Vγ−1 (Y − Xθ))

On passe à −2 log de cette expression et à une constante près on est conduit à maximiser

L(θ, γ) := log(|Vγ |) + (Y − Xθ)0 Vγ−1 (Y − Xθ)

Pour trouver les équations du maximum de vraisemblance, on remarque

∂Vγ−1 ∂Vγ −1
= −Vγ−1 V
∂γi ∂γi γ

∂log|Vγ | ∂Vγ0 −1
= T r( V ).
∂γi ∂γi γ
On en déduit que
∂L(θ, γ)
= 0 ⇔ (X 0 Vγ−1 X)θ = X 0 Vγ−1 Y
∂θ
10 CHAPITRE 1. MODÈLES MIXTES

Cette équation est connue sous le nom d’équation de Gauss-Markov. Elle permet de re-
trouver l’équation(1.5).

∂L(θ, γ)
= 0 ⇔ T r(Vi Vγ−1 ) = (Y − Xθ)0 Vγ−1 Vi Vγ−1 (Y − Xθ)
∂γi

Cependant, dans beaucoup de cas les simulations montrent que le maximum de vrai-
semblance est biaisé. Comme la vraie difficulté, comme on l’a vu réside dans l’estimation
des composantes de la variance, on va d’une façon “concentrer” la vraisemblance sur cette
estimation.

Définition 3 On appelle vraisemblance restreinte, la vraisemblance de T 0 Y où T est un

matrice quelconque de rang maximal telle que T 0 X = 0. Il est facile de voir que T 0 Y ∼
N (0, T 0 Vγ T ) et donc on pose

LR (γ) := log(|T 0 Vγ T |) + Y 0 T (T 0 Vγ T )−1 T 0 Y

Proposition 2 LR ne dépend de T qu’à travers l’ajout d’une constante qui ne change

rien et de plus

∂LR (γ)
= 0 ⇔ T r(Vi Vγ−1 Qγ ) = Y 0 Vγ−1 Qγ Vi Vγ−1 Qγ Y
∂γi

En conséquence en comparant les équations, on constate qu’un point fixe du MIVQUE (γ̂ =
γ0 ) est une solution des équations du maximum de vraisemblance restreinte. L’itération
du MIVQUE est donc une façon (parmi d’autres) de résoudre ces équations. Pour plus de
détail on pourra consulter Azaı̈s Bardin Dhorne (1993).

2.7 Tests

Pour ce qui concerne les tests d’hypothèses sur les effets fixes, on réalise des tests de
Fisher en supposant que les estimateurs des composantes de la variance sont en fait les
valeurs exactes et en utilisant la méthode de la section 2.4.
Le plus souvent les tests sur les effets aléatoires correspondent à la nullité d’une va-
riance : on peut chercher à tester la nullité d’un effet “famille génétique”, ou “sujet” par
exemple. La première solution consiste à utiliser un test exact de Fisher. En effet la nullité
d’un effet aléatoire correspond strictement à l’absence d’effet, c’est à dire également à la
nullité d’un effet déclaré en effet fixe. En résumé, pour tester la nullité d’une composante
de la variance, on peut déclarer l’effet correspondant en fixe et utiliser le test de Fisher
correspondant. Sauf dans le cas équilibré (voir Coursol 1980) ce test n’est plus optimal.
Mais il est exact dans le sens où son niveau réel est toujours égal au niveau nominal.
L’autre alternative est d’utiliser les tests classiques asymptotiques associés à la méthode
du maximum de vraisemblance : le test du rapport de vraisemblance et le test de Wald.
Pour les présenter, nous utilisons les notations traditionnelles où θ est le paramètre du
modèle statistique. Le θ de cette partie est donc en fait égal au θ, γ du modèles mixte dans
le cas de la vraisemblance classique, il est égal au seul γ dans le cas de la vraisemblance
restreinte.
2. MODÈLES À EFFETS ALÉATOIRES ET MIXTES 11

Dans une expérience qui comprend une grand nombre de répétitions indépendantes et
sous des hypothèses de régularité (modèles de vraisemblance réguliers) qui sont vérifiées
dans notre cas, on sait que l’estimateur du maximum de vraisemblance est asymptoti-
quement sans biais, normal et de variance donnée par l’inverse de l’information de Fisher
(Dacunha-Castelle & Duflo, 1983). L’information de Fisher I(θ) est donnée par
 ∂ 2 log V (θ) 
(I(θ))ij = −E
∂θi ∂θj

où V est la vraisemblance, θ le paramètre du modèle. Comme, I(θ) tend vers l’infini avec
le nombre de répétition, la phrase “ asymptotiquement sans biais, normal et de variance
donnée par l’inverse de l’information de Fisher” doit se comprendre comme

(I(θ))1/2 (θ̂ − θ) → N (0, Id)

quand le nombre de répétitions tend vers l’infini.

Une telle propriété reste vraie pour la vraisemblance restreinte et dans des cadres plus
généraux, mais sous des conditions très techniques (Jiang 1996). Ces conditions sont en
particulier que les composantes de la variance correspondent à des effets aléatoires β 1 ...βk
au sens de notre première définition d’un modèles mixte et que ces composantes soient
toutes positives. Dans ces conditions il est possible de construire un test de Wald de
l’hypothèse
θi = θi,0
dont la région de rejet est
−1/2
(I(θ)−1
ii ) (θ̂i − θi,0 ) > Zα

où Zα est la valeur critique pour la valeur absolue d’une loi normale standard.
Compte tenu des hypothèses de l’article Jiang (1996), ce test n’est pas valide pour
tester la nullité d’une composante. En pratique on peut vérifier que si on l’applique quand
même et et ce surtout pour les petits échantillons, il est peu puissant et conservatif (le
niveau réel est nettement plus important que le niveau nominal).
Des calculs élémentaires quoique longs, montrent que dans le cas de la vraisemblance
restreinte, l’information de Fisher vaut
1
Iγ = { T r(Vi V −1 QVj V −1 Q), i, j = 1, k}
2

Un autre test possible qui est en général plus puissant quoique parfois non-conservatif
est le test du rapport de vraisemblance. Si L g est la log vraisemblance prise au
maximum et Lp est la vraisemblance ou maximum sous l’hypothèse nulle, on montre que
sous les mêmes hypothèses que précédemment
L 1 2 0
Lg − L p → χ (k ),
2
où k 0 est la différence de dimension paramétrique entre les deux modèles : hypothèse
générale, hypothèse nulle.
Ces deux tests : Wald et Rapport de Vraisemblance peuvent être indifféremment utilisés
pour le maximum de vraisemblance comme pour le maximum de vraisemblance restreint.
12 CHAPITRE 1. MODÈLES MIXTES
Chapitre 2

Plans d’expériences randomisés

1 Introduction

Dans cette partie nous allons présenter l’idée forte suivante : la statistique n’a pas
comme seul objet de “ traiter des données” mais également d’en préparer le recueil pour
en améliorer la qualité. D’importants gains sont possibles lors de cette étape. Les méthodes
que nous allons présenter s’appliquent plutôt à des expériences de “labo” dans le sens le
plus général possible, plutôt qu’à des situations où les données sont recueillies “comme
elles viennent” par exemple dans les enquêtes.
Les deux buts de la planification sont : (1) de permettre une interprétation claire en
évitant les confusions, (2) de maximiser la précision de l’expérience.
Pour illustrer le premier point, prenons l’exemple de la scolarisation en maternelle. Des
études incontestables ont montré que les enfants scolarisés en maternelle ont de meilleurs
résultats dans la suite de leur scolarité que les enfants qui ne rejoignent l’école qu’au
primaire. Doit on en déduire qu’il faut rendre la scolarisation en maternelle obligatoire
pour lutter contre l’échec scolaire ? Une réponse directe : oui n’est pas possible. En effet
deux interprétations sont possibles : (a) c’est effectivement la scolarité en maternelle qui
améliore les résultats ; (b) dans la France actuelle les élèves qui ne vont pas en maternelle
sont une minorité qui correspond à des groupes sociaux bien particuliers, ce qui peut
expliquer la différence de réussite scolaire par l’origine sociale.
Dans cet exemple on pourra affiner l’analyse en contrôlant toute les variables indiquant
la situation sociale, avec toujours le risque d’en oublier une. Mais il est clair que sans cette
information complémentaire les données de départ sont sans valeur pour répondre à la
question. Une solution de type planification à ce problème serait de définir un groupe
d’enfants test et d’au hasard les répartir en deux groupes l’un qui serait scolarisé en
maternelle l’autre non. Bien sûr, c’est impossible.
Pour illustrer le second point, considérons le problème de pesée de deux objets A et
B avec une balance sans biais qui donne chaque résultat avec une erreur indépendante de
variance σ 2 . On suppose que la balance est capable de peser les deux objets ensembles.
1ère méthode : On pèse A et B séparément ; Le coût est de deux pesées la précision
est σ 2 .
2ème méthode : On pèse A + B et A − B ; Le coût est de deux pesées, la précision est
σ 2 /2,car les poids de A et B sont maintenant obtenues comme moyennes de deux pesées.

13
 14 CHAPITRE 2. PLANS D’EXPÉRIENCES RANDOMISÉS

2 Nécessité de la randomisation

Comme nous allons essayer de vous en convaincre, la randomisation est la seule méthode
qui évite les confusions (en fait elle en contrôle la probabilité), permet de faire une
expérience équitable, de le prouver et enfin permet d’apprécier la précision des résultats.
Nous voulons comparer l’efficacité de deux médicaments A et B contre la grippe sur
40 malades.
expérience 1 On administre le médicament A aux 20 premiers malades qui se présentent. On note
leur état d’amélioration, ensuite on administre B aux 20 suivants et on note leur
état. A la fin de l’expérience, on calcule les moyennes et on déclare comme meilleur
le médicament qui a la meilleure moyenne.
Cette expérience est désastreuse : (1) durant la durée de l’expérience la maladie,
la température extérieure, la fatigue des personnes qui ont réalisé l’expérience ont
pu évoluer : l’expérience n’est plus équitable dû à la confusion entre le temps et le
médicament ; (2) Certains participants de l’expérience : malades ou médecins qui
connaissent parfaitement les médicaments administrés peuvent fausser le résultat
inconsciemment : c’est l’effet placebo bien connu.
expérience 2 Au fur et à mesure qu’un patient arrive en consultation, on alterne strictement A et
B, le plan est donc ABABABABABAB... Ce plan est certainement meilleur que le
précédent mais souffre encore de deux gros défauts : en premier lieu la systématicité
rend impossible le “double aveugle” le médecin et dans une certaine mesure le malade
sauront toujours la nature du produit administré ; en second lieu on ne dispose pas
de méthode statistique valide pour choisir entre les situations : A meilleur que B ; B
meilleur que A ; A et B équivalents.
expérience 3 On pourrait construire une variante de l’expérience 2 où l’on essayerait de répartir
au mieux les individus entre les deux groupes en fonction de l’âge, du poids, des
antécédents. Ce plan séduisant a exactement les mêmes inconvénients que le précédent.
expérience 4 On tire au hasard 20 personnes parmi les 40 premiers malades auquel on administre
A. Les autres reçoivent B. On pose un modèle conceptuel : on imagine avoir admi-
nistré les deux médicaments à chacun des malades et note R ik la réponse du malade
i au médicament k. On pose
Rik = mi + ak

qui suppose l’additivité entre l’effet du malade a k et l’effet du médicament mi . Pour

se débarrasser d’une indétermination on pose k ak = 0. Maintenant on note Yij
P

la réponse du jème malade parmi ceux qui ont reçu le traitement i. Du modèle
précédent on déduit
Yij = mi + bij

où les bij sont tirés sans remise parmi les ak . Les propriétés du tirage sans remise
montrent que les bij sont d’espérance nulle et que leur matrice de variance (rangés
par exemple en ordre lexicographique) vaut

σ 2 (Id − J/n),

où n est le nombre total de données, J est la matrice n × n composée de 1 et

1 P
σ 2 = n−1 2
k=1,n ak
3. PLANS D’EXPÉRIENCES CLASSIQUES 15

Non seulement, le terme extra-diagonal σ 2 /n est petit mais surtout il est constant :
son influence sur les comparaisons entre traitements est nulle, de sorte que l’on peut
faire une analyse de la variance sur les données présentes, ce qui permet de savoir si
les médicaments sont significativement différents.
Exercice 1 Prouver l’affirmation ci-dessus :

2σ 2
Var(Yi. − Yi0 . =
r

Enfin l’expérience 4 permet le travail en “double aveugle”, ainsi elle garantit et

prouve l’équitabilité de l’expérience.

3 Plans d’expériences classiques

dans la présentation on utilisera les notations classiques suivantes : t, r, b, k, n sont
respectivement les nombres de traitements, répétitions de chaque traitement, de blocs,
d’unités par bloc et enfin de données.

Plan en randomisation totale

C’est le plan de l’expérience 4 précédente. Dans cas général, on a t traitements, rt unités

expérimentales. On tire au hasard et sans remise, r unités pour le premier traitement, puis
r pour le second etc... les données sont analysées par une analyse de la variance à 1 facteur :
le traitement. A part dans les expérience médicales ce plan est en fait peu utilisé.

Plan en blocs complets

[Fisher 1931] Ce plan suit les trois principes de R. Fisher : répétition ; randomisation
et contrôle local. Dans le même cadre que le plan en randomisation totale, on regroupe
les rt unités expérimentales en r blocs homogènes de taille t. Dans un exemple médical
les blocs peuvent par exemple correspondre aux sexe, à l’âge etc. Dans une expérience
agronomique ce sera un ensemble de parcelles contiguës, dans toute expérience de labo,
ce sera les unités traitées le même jour, par la même personne. Le principe du plan est le
suivant : dans chaque bloc on alloue indépendamment une unité exactement à chaque
traitement et ce de façon aléatoire.
Exemple : Une association de consommateurs désire comparer le fonctionnement de
t aspirateurs. Les essais sont conduits par des ménagères membres de l’association. On
suppose qu’il y a b ménagères et que chacune veut bien conduire t essais. Chaque ménagère
a donc à sa disposition t “cases” pour un essai. L’ensemble de ces cases de taille tb constitue
l’ensemble des unités, les ménagères constituent les blocs. Le plan revient à tirer au hasard
l’ordre dans lequel chaque ménagère va expérimenter chaque aspirateur.
Pour être efficace le plan doit maximiser la variabilité inter-bloc et minimiser la varia-
bilité intra-bloc : en d’autres termes il faut rendre les blocs le plus homogènes possible.
Dans l’exemple, on aura intérêt à ce que les t essais d’une même ménagère soient aussi
rapprochés que possible. Il est beaucoup moins gênant par contre qu’au total l’expérience
se déroule sur un long laps de temps.
16 CHAPITRE 2. PLANS D’EXPÉRIENCES RANDOMISÉS

Des calculs analogues à ceux fait pour le plan en randomisation totale mais un peu
plus complexes montrent que l’on peut valider une analyse de la variance à deux facteurs
additifs : traitement et bloc. Cette analyse est équirépétée.
Le plan en blocs complets est quasi-toujours préférable au plan en randomisation totale.
C’est le plan le plus employé. C’est celui que l’on essayera d’utiliser a priori.

Plan en blocs incomplets

Plans équilibrés : Supposons que nous ayons 9 bières à comparer (facteur d’intérêt :
facteur traitement) et 12 dégustateurs (facteur bloc). Il est clair que passé un certain
nombre de dégustations, un dégustateur est incapable de comparer ses sensations. On
supposera donc qu’un dégustateur ne peut comparer que 3 bières. Nous avons donc 36
unités réparties en bloc de 3 et on propose la répartition suivante :

Deg. 1 1 2 3 Deg. 4 1 4 7
Deg. 2 4 5 6 Deg. 5 2 5 8
Deg. 3 7 8 9 Deg. 6 3 6 9
Deg. 7 1 5 9 Deg. 10 1 8 6
Deg. 8 4 8 3 Deg. 11 2 4 9
Deg. 9 7 2 6 Deg. 12 3 5 7

Cette répartition qui n’a rien d’aléatoire, pour l’instant, est équilibrée : chaque traitement
se retrouve une fois et une seule exactement avec chaque autre traitement.
Un plan en blocs incomplets équilibré est un plan possédant la propriété ci- dessus, il est
décrit (partiellement) par :
le nombre de traitements t (9 dans l’exemple)
le nombre de répétitions r (4)
le nombre de blocs b (12)
la longueur d’un bloc k (3)
l’indice de concurrence λ : le nombre de fois ou deux traitements se retrouvent ensemble
(1). En comptant de deux manières différentes le nombre de parcelles et le nombre de
voisins, on obtient :
rt = bk
r(k − 1) = λ(t − 1)
Il n’existe pas des solutions pour toutes les tailles vérifiant ces équations. Il existe des
tables de plans : Raghavarao (1971).

Les lattices : Nous présentons ci dessous une méthode pour construire des plans
qui sont sous certaines conditions équilibrés et qui conservent de toute façons de bonnes
propriétés : les plans lattice.

Définition 4 (lattices) Un (h, p2 ) lattice est un plan pour p2 traitements avec p premier
ou puissance de premier avec h répétitions et des blocs de longueur p. Aux p 2 traitements
sont associés deux facteurs A et B à p niveaux et dont les niveaux sont numérotés par des
éléments du corps Fp à p éléments. Le plan est constitué de h répliques : ensembles de p
blocs qui contiennent une fois et une seule chaque traitement. Dans la réplique 1 on confond
3. PLANS D’EXPÉRIENCES CLASSIQUES 17

A avec les blocs, dans la réplique 2 B et dans les autre successivement A + iB, i ∈ F p \0.
On montre que le lattice (p + 1), p2 est équilibré.

Exercice 2 En utilisant la description de F 4 comme l’ensemble {0, 1, x, 1 + x} muni de

l’addition modulo 2 et de la multiplication modulo 2 et le polynôme irréductible 1 + x + x 2 ,
construire le lattice 5, 42 . Par exemple sur F4 (1 + x)(1 + x) = 1 + 2x + x2 = x

Définition 5 (Plans circulants) Un plan circulant est un plans pour t traitements en

t blocs de tailles t − 1. Pour construire les blocs on élimine tour à tour chacun des traite-
ments.
Un plan circulant est toujours équilibré, son indice de concurrence λ vaut t − 2.

Plans non équilibrés : Quand il existe un plan équilibré, on montre qu’il est optimal
(Pour plus de détail voir Druilhet, 19xx). ma lheureusement cela n’est pas toujours possible
en particulier l’équilibre demande souvent un nombre de répétition élevé. Dans ce cas il
existe des méthodes pour construire des plan conservant certaines propriétés, par exemple
l’équilibre partiel (Coursol, 1980). mais de toutes façons les propriétés suivantes restent
vraies que le plan soit équilibré ou non.
Randomisation : La randomisation se fait en deux étapes,
(1) “ mélange des blocs” : dans l’exemple on affecte un numéro de dégustateur à un
dégustateur réel (M Dupond) au hasard ;
(2)” mélange des traitements par bloc” : dans l’exemple, les trois bières devant être
présentées à un dégustateur, le sont dans un ordre aléatoire.

Analyse : On montre que la randomisation valide un modèle avec des effets traitement
fixes et des effets blocs aléatoires : c’est un modèle mixte. Si on ne dispose pas des moyens
de traiter un tel modèle, on peut toutefois utiliser un modèle avec blocs et traitements
fixes qui correspond à une légère perte d’information.

Carrés latins

Exemple : on veut comparer 4 peintures sur 4 maisons carrées ayant 4 façades de

mêmes expositions N, S, E, O. Le facteur d’intérêt est le facteur peinture, les facteurs
maison et orientation étant des facteurs parasites. Si on veut équilibrer les relations du
premier facteur avec chacun des deux autres, on peut être amené (avant randomisation)
à utiliser la répartition suivante
peinture orientation
maison N S E O
A 1 2 3 4
B 2 3 4 1
C 3 4 1 2
D 4 1 2 3
On appelle carré latin, un plan pour trois facteurs au même nombre n de niveaux,
comprenant n2 unités et tel que le nombre de répétitions d’une paire de niveaux pour
deux facteurs est toujours 1.
18 CHAPITRE 2. PLANS D’EXPÉRIENCES RANDOMISÉS

Randomisation : La répartition du tableau ci dessus n’est pas très aléatoire, on

satisfera à l’exigence de randomisation en faisant dans l’ordre que l’on voudra un “mélange
des lignes” et “un mélange des colonnes”. Cette randomisation valide, sous les mêmes
hypothèses d’additivité que pour le plan en randomisation totale, une analyse de la variance
à trois facteurs additifs.
La méthode de construction “en diagonale” fonctionne pour toute valeur de n. Mais
il existe de nombreuses méthodes de construction. La plus classique est basée sur un
groupe fini. On construit la valeur du troisième facteur comme la table d’addition de deux
permutations du groupe (noté additivement).
Plans en ligne et colonnes non équilibrés : On peut s’intéresser a un plan dont les
unités sont placées sur un réseau de l lignes et de c colonnes mais dont l’allocation des trai-
tements n’a plus la belle propriété du carré latin. On utilise encore la même randomisation,
mais elle valide maintenant une analyse avec un effet ligne et colonne aléatoire.

Exercice 3 Vérifiez informatiquement que pour un carré latin on peut déclarer indifféremment
les effets lignes et colonnes comme fixes ou aléatoires, cela ne change pas les résultats.

Split-plot (parcelle subdivisée)

Exemple : on désire comparer 3 recettes de gâteau au chocolat et 6 températures de

cuisson 180,(10)...,230. On réalise l’expérience suivante : chaque jour, on réalise 3 pâtes
correspondant aux 3 recettes. Chaque pâte est divisée en 6 sous-partie cuites à chacune des
6 températures différentes. Après cuisson, on mesure le mœlleux du gâteau en mesurant
l’angle de rupture α d’une tranche. On a donc pour expliquer cette variable quantitative
trois facteurs : la température (que l’on décide de considérer comme qualitatif), la recette
et le jour (facteur bloc). Ce qui est particulier à cet exemple est qu’il y a deux sources
d’erreurs : une attachée à la mesure de l’angle et l’autre à la confection d’une pâte.

Définition 6 On appelle plan split-plot un plan pour deux facteurs traitements. Le pre-
mier A à t niveaux et le second B à s niveaux. Pour le premier facteur, on construit
un plan en blocs complets à rt “grandes unités” avec sa randomisation. Ensuite, chaque
“grande unité” est divisée en s sous unités auxquelles sont affectées dans un ordre aléatoire
les s valeurs du traitement B

Si i, j, k sont les niveaux de A, B, bloc dans l’ordre, le modèle auquel conduit la rando-
misation est

Yijk = µ + ai + bj + a ∗ bij + blk + Eik + ijk ; i = 1, t ; j = 1, s ; k = 1, r.

où E est l’erreur associée aux “grandes unités”.

Ce modèle peut s’analyser bien évidemment à l’aide d’un programme de modèle mixte,
mais il peut également s’analyser à partir des projections inter et intra “grande unité”. La
projection inter revient à travailler sur Y i.k qui avec des contraintes classiques donne

Yi.k = µ + ai + blk + Eik + i.k .

3. PLANS D’EXPÉRIENCES CLASSIQUES 19

Les deux aléas peuvent être confondus, on retrouve ainsi le modèle du plan en blocs
complets. La projection intra est théoriquement basées sur les Y ijk − Yi.k . On montre
qu’elle est équivalente au modèle complet dans lequel l’effet E jk est supposé fixe, ceci à
condition de se limiter à l’estimation et aux test sur le facteur B et sur l’interaction A ∗ B.
En conclusion on dit que le facteur A est totalement estimable inter-grandes unités et
l’interaction A ∗ B et le facteur B sont totalement estimable intra-grandes unités.
20 CHAPITRE 2. PLANS D’EXPÉRIENCES RANDOMISÉS
Chapitre 3

Plans fractionnaires

1 Introduction
On considère une réaction chimique qui dépend pour simplifier de trois facteurs : le
ph P H, avec une valeur standard de 7 , la température, T avec comme valeur standard
30◦ C, et la dose D, avec comme valeur standard 100. On sait que l’on peut faire varier
chacun de ces facteurs entre deux limites et on cherche à savoir s’il ont une influence sur
la réponse, par exemple le rendement. On va comparer deux expériences
Expérience 1 : On fait varier d’abord le facteur P H
On fait 4 répétitions à P H = 6.5, T = 30 ◦ , D = 100 et et 4 répétitions à P H = 7.5,
T = 30◦ , D = 100 .
On fait varier ensuite le facteur température
On fait 4 répétitions à P H = 7, T = 25 ◦ , D = 100 et et 4 répétitions à P H = 7,
T = 35◦ , D = 100 .
On fait varier enfin le facteur dose
On fait 4 répétitions à P H = 7, T = 30 ◦ , D = 90 et et 4 répétitions à P H = 7,
T = 30◦ , D = 110.
Le coût total de cette expérience est de 24 unités et la puissance expérimentale est la
comparaison de moyennes de 4.
Expérience 2 : On réalise toutes les 8 combinaisons entre les deux valeurs hautes et
basses des trois facteurs. C’est à dire les huit unités

(6.5, 25, 90), (6.5, 25, 110), (6.5, 35, 90), (6.5, 35, 110),

(7.5, 25, 90), (7.5, 25, 110), (7.5, 35, 90), (7.5, 35, 110)

Le coût total est maintenant de 8 unités. Si on fait l’analyse à l’aide d’un modèle additif
à trois facteurs : H, T, D, ce modèle est orthogonal, on comparera donc des moyennes de 4
observations pour tester la significativité d’un facteur. Donc en première approximation,
on a construit une expérience trois fois mois chère qui a la même puissance.
Remarquons bien, qu’en présence d’interactions entre les facteurs, il sera possible de les
détecter dans la seconde expérience : en effet il reste 4 degrés de liberté dans la résiduelle,
il encore possible d’en occuper un ou deux pour un ou deux termes d’interaction. Par

21
22 CHAPITRE 3. PLANS FRACTIONNAIRES

contre comme l’expérience 1 est subdivisées en sous-expériences dans lesquelles on ne fait

varier qu’un facteur à la fois, il ne sera pas possible de détecter la présence d’interaction.
En conclusion, quand un phénomène dépend de plusieurs facteurs, on a toujours intérêt
à les étudier globalement. Le plan le plus simple est celui de l’expérience 2 où l’on fait
toutes les combinaisons des niveaux (dans notre exemple chaque facteur a deux niveaux)
de tous les facteurs. Il est appelé plan factoriel complet. Si il y a p facteurs tous à 2
niveaux, le plan factoriel complet demande 2 p unités. Cela est souvent trop coûteux. On
va donc faire une partie : une fraction du plan complet. Un tel plan est appelé plan
factoriel fractionnaire

2 Cadre général pour des facteurs à deux niveaux

On considère p facteurs à deux niveaux noté −1, +1. On appelle traitement une com-
binaison (i1 , ..., ip ) de ces facteurs. L’ensemble des traitements est de cardinal 2 p et corres-
pond au plan factoriel complet. Soit E l’espace des réponses aux divers traitements. Cet
espace peut être vu comme l’ensemble des fonctions

[−1, +1]p → IR
p
ou comme l’espace vectoriel IR2 .
On se place dans le cas où n := 2p est trop grand pour pouvoir expérimenter toutes
les combinaisons : on va n’en faire qu’une partie, une fraction d’où le nom.
En premier, nous allons définir proprement les interactions multiples et les effets prin-
cipaux. Dans E il existe des éléments particulier : les fonction coordonnées. On définit la
kième fonction coordonnée Ak , (1 ≤ k ≤ p) par
A
(i1 , ..., ip ) →k ik

Soit B ⊂ {1, ..., p}, on définit la fonction A B élément de E par

AB := A11 ...App ; k = 1 si k ∈ B ; k = 0 sinon

Proposition 3 Quand B varie dans l’ensemble des parties de {1, ..., p}, les A B forment
√
une base orthogonale de E de norme n.

On définit les espaces suivants pour B partie de {1, ..., p}.

VB =fonction des seules coordonnées dans B
WB orthogonal dans VB de tous les VD pour D ⊂ B, D 6= B

Exemple 1 Par exemple si p = 2 un élément de E peut s’écrire comme un tableau 2 de

la forme !
a b
.
c d
Un élément de V{1} est de la forme
!
a a
.
b b
2. CADRE GÉNÉRAL POUR DES FACTEURS À DEUX NIVEAUX 23

Un élément de W{2} :
!
a −a
.
a −a

Proposition 4 dim VB = 2|B| , où —B— est le nombre d’éléments de B

si D ⊂ B, AD ∈ VB ,
{AD pour D ⊂ B} engendre VB ,
AB ∈ W B ,
AB engendre WB qui est donc de dimension 1.

WB est défini comme l’espace de l’interaction entre les facteurs de B avec les exception
suivantes : si B est réduit à un facteur , il s’agit alors de l’effet principal de ce facteur, si
B = ∅, W∅ est l’espace de la “moyenne générale”. Il est de dimension 1 quelle que soit la
taille de B.
Les différents résultats de la proposition s’enchaı̂nent aisément de sorte que leur démonstration
est aisée.

Définition de la valeur d’un effet principal ou d’une interaction : La vraie

réponse f au traitement est un élément de E qui se décompose sur la base des générateurs
AB . Nous appellerons “effet de AB ” noté e(AB ) le coefficient de cette décomposition.

e(AB )AB .
X
f=
B⊂{1,...,p}

Exercice 4 construisez les différents espaces : V ∅ , V{1} , V{2} , V{1,2} , W∅ , W{1} , W{2} , W{1,2}
dans le cas de” p = 3 facteurs.

Exemple avec trois facteurs.

Le tableau suivant décrit le plan factoriel complet ainsi que les différents générateurs
avec nos notations

A∅ A{1} A{2} A{3} A{1,2} A{1,3} A{2,3} A{1,2,3} traitement

1 −1 −1 −1 +1 +1 +1 −1 (−1, −1, −1)
1 −1 −1 +1 +1 −1 −1 +1 (−1, −1, +1)
1 −1 +1 −1 −1 +1 −1 +1 (−1, +1, −1)
1 −1 +1 +1 −1 −1 +1 −1 (−1, +1, +1)
1 +1 −1 −1 −1 −1 +1 +1 (+1, −1, −1)
1 +1 −1 +1 −1 +1 −1 −1 (+1, −1, +1)
1 +1 +1 −1 +1 −1 −1 −1 (+1, +1, −1)
1 +1 +1 +1 +1 +1 +1 +1 (+1, +1, +1)

Notons A := A{1} , B := A{2} , C := A{3} . Choisissons les 4 données pour lesquelles

ABC = 1, ce qui définit un “demi plan” ou “fraction”. Même en supposant observer la
réponse sans erreur, on est amené à résoudre le système de quatre équations :

f (−1, −1, +1) = e(1) + e(ABC) − e(A) − e(BC) − e(B) − e(AC) + e(C) + e(AB)
24 CHAPITRE 3. PLANS FRACTIONNAIRES

f (−1, +1, −1) = e(1) + e(ABC) − e(A) − e(BC) + e(B) + e(AC) − e(C) − e(AB)
f (+1, −1, −1) = e(1) + e(ABC) + e(A) + e(BC) − e(B) − e(AC) − e(C) − e(AB)
f (+1, +1, +1) = e(1) + e(ABC) + e(A) + e(BC) + e(B) + e(AC) + e(C) + e(AB)
D’après les relations d’orthogonalité déjà prouvées, ce système a une solution unique en
les inconnues
e(1) + e(ABC)
e(A) + e(BC)
e(B) + e(AC)
e(C) + e(AB).
Dans le tableau restreint correspondant à ABC = 1, certaines relations sont vérifiées
entre les vecteurs. Elle découlent directement de la relation qui a défini la fraction : 1 =
ABC, A = BC, B = AC, C = AB. Elles impliquent directement les confusions observées.
En effet dans l’écriture générale :

e(AB )AB ,
X
f=
B⊂{1,...,p}

certains AB sont confondus sur la fraction, on regroupe donc les coefficients.

3 Méthode des facteurs de base

On appellera plan 2p−q (p > q), un plan comprenant 2p−q unités pour p facteurs à deux
niveaux. Le taux de fraction est donc 2 −q . Une méthode de construction qui donne des
fractions régulières (terme que l’on définira ultérieurement) est la méthode des facteurs
de base : on construit le plan factoriel complet pour les p − q premiers facteurs. On
définit ensuite les valeurs des autres facteurs en fonction des p premiers (sous forme de
produit exclusivement). On appelle relations de définition ces q relations que l’on exprimera
sous la forme canonique 1 = .... Il est facile de vérifier que la fraction ainsi définie est
exactement le sous ensemble du plan complet qui vérifie les relations de définition. On
définit également les relations complètes de définition qui sont les 2 q terme égaux que l’on
obtient en combinant les relations de définition. Ces relations donnent l’alias de 1 c’est
à dire les termes confondus avec 1 . A partir de cet alias et par combinaison, on obtient
l’alias d’un terme quelconque donné qui comprend 2 q − 1 termes.

Exemple 2 Considérons le plan 25−2 pour 5 facteurs A, B, C, D, E avec A, B, C comme

facteurs de base et les relations de définition : D = AB, E = AC. Il s’agit donc d’un
quart de plan à 8 unités. Les relations complètes de définition sont

1 = ADB = ACE = DCBE.

En effet A2 = B 2 = ... = E 2 = 1. Il y a bien 4 effets confondus. L’alias de A par exemple

comprend trois termes :
A = DB = CE = ABCDE.
Il y a donc confusion dans cette fraction entre entre l’effet principal du facteur A et les
interactions D ∗ B et C ∗ E.
4. PLAN POUR L’ÉTUDE DES EFFETS PRINCIPAUX ET DES INTERACTION DOUBLES25

Estimation des effets

On suppose que les observations sont faite avec une variance σ 2 et qu’elle sont non
corrélées. Choisissons une suite de générateurs :A B1 , ..., ABm , m ≤ 2p−q telle qu’il y ait
au plus un représentant de chaque alias. On n’est pas obligé de les prendre tous. Alors le
modèle linéaire m
ẽ(ABi )ABi + i1 ,...ip ,
X
Y (i1 , ...ip ) =
i=1

où ẽ(ABi ) est l’alias de e(ABi ) au sens de la somme de effets de tous les générateurs qui
lui sont confondus, est orthogonal sur l’ensemble des données du plan fractionnaire. La
matrice d’information : X 0 X vaut nId (n est le nombre de données : 2 p−q ) et les estimateurs
Bi
ˆẽ(ABi ) = < Y, A > ,
n
sont non corrélés de variance σ 2 /n.
Un fraction qui a la propriété que deux effets sont orthogonaux ou confondus est appelé
une fraction régulière.

4 Plan pour l’étude des effets principaux et des interaction

doubles
Un plan est dit de résolution ρ entière si
ρ =inf{ nb de symboles des éléments de l’alias de 1}
Exemples :
Résolution III : tous les effets principaux sont non confondus.
Résolution IV : un effet principal ne peut être confondu avec un interaction, mais deux
interactions peuvent être confondues.
Résolution V : on peut poser un modèle avec toutes les interactions et effets principaux
sans confusion.
De manière générale, la résolution V est considérée comme suffisante dans toutes les
situations, alors que la résolution III est considérée comme une propriété minimale.

Définition 7 (Aberration) Cette notion est un raffinement de la notion de résolution.

L’aberration d’un plan pour p facteurs est un p-uple dont l’élément i contient le nombre
de symboles de i lettres dans l’alias de 1 (1 exclu). Par exemple l’aberration du plan de la
partie 3 vaut : (0, 0, 2, 1, 0)

Les aberrations sont munies de l’ordre lexicographique, on parle de plan de minimum

d’aberration.

Proposition 5 Si le plan est de résolution III, tous les effets principaux sont estimables
donc en comptant les degrés de liberté, le nombre p de facteurs vérifie

p≤n−1
26 CHAPITRE 3. PLANS FRACTIONNAIRES

Ce maximum est atteint dans le sens où pour tout r il existe un plan de résolution III avec
2r unités et 2r − 1 facteurs.
Une fraction régulière de résolution IV avec 2 r unités comprend au maximum 2r−1
facteurs. Ce maximum est atteint.
En résolution V, il n’y a pas de résultat général, on connaı̂t simplement le nombre de
facteurs maximaux pour les petites valeurs

r 4 5 6 7 8 9
nombres d’unités = 2 r 16 32 64 128 256 512
nb max de facteurs p 5 6 8 11 17 23
ddl du modèle 16 22 37 67 154 277
ddl du modèle avec 1 Fac. de plus 22 29 46 79 172 301

Compléments
– Quand le plan est de taille trop importante pour pouvoir être conduit de manière
homogène, il doit être découpé en blocs. Les facteurs blocs peuvent être considérés
comme des facteurs ordinaires sauf que les notions de résolution et aberration n’ont
plus la même pertinence. Le mieux est de regarder en détail les confusions en étant
bien conscient que tout effet confondu avec un facteur bloc sera totalement non-
estimable intra-bloc (voir exemple ci-dessous).
– Si un ou plusieurs des facteurs ont un nombre de niveaux égaux à 4, 8 (ou plus
généralement une puissance de 2) on peut se ramener au cas précédent en recodant
les niveaux à l’aide de 2 ou 3 pseudo-facteurs : Si A possède 4 niveaux notés 1, 2, 3, 4
on peut le recoder à l’aide de 2 pseudo-facteurs en utilisant la table suivante :
A A 1 A2
1 1 1
2 1 −1
3 −1 1
4 −1 −1
– Exemple. On veut expérimenter 5 facteurs : A, B, C, D, E à l’aide de 32 unités ce qui
correspond à un plan complet. Malheureusement on considère qu’il est impossible de
réaliser plus de 8 unités de manière homogène de sorte qu’il faut introduire un facteur
bloc BL à 4 niveaux qui sera codé par deux pseudo-facteurs B 1 et B2 . On est donc
amené à chercher un plan 27−2 . Celui donné par la table ci dessous B 1 = ABCD
B2 = ABDE amène à confondre l’interaction CE avec B 1 B2 qui est un effet bloc.
On laisse à titre d’exercice, le soin au lecteur de vérifier que le plan suivant est
meilleur B1 = ABC B2 = CDE.
– Plus de deux niveaux. Les plans fractionnaires utilisent dans le cas général des tech-
nique de corps finis qui sont basés sur des nombres premiers. Pour ces raisons il
n’est possible que de travailler que sur un seul nombre premier : on peut donc faire
des plan pour des facteurs à 2, 4 ou 8 niveaux comme pour de facteurs à 3, 9 ou
27 niveaux, mais il est impossible de mélanger les deux. Pour des facteurs à 3 ni-
veaux il faut travailler non plus avec {−1, 1} mais avec les racines cubiques de l’unité
{1, j, j 2 }, ce qui amène à un présentation plus technique.

Exercice 5 Choisir un plan permettant d’expérimenter 7 facteurs A − G en estimant les

effets principaux en présence d’éventuelles interactions. Donner les relations complètes de
5. MÉTHODE TAGUSHI 27

définitions et l’alias du facteur A.

Exercice 6 On veut construire un plan pour 5 facteurs A − E, on veut pouvoir estimer

sans confusion les effets principaux et les interactions suivantes : A∗D, A∗E, B ∗D, B ∗E.
Proposez une solution de dimension minimale.

Exercice 7 Combien faut-il d’unités au moins pour expérimenter 12 facteurs en résolution

IV et V respectivement ?

Exercice 8 (repliement de plans) Un Expérimentateur réalise le plan 2 5−2 définie par

les relations D = ABC et E = BC. Plus tard il réalise un plan identique sauf que tous
les signes ont été inversés.
i. Quelle est la résolution du plan ainsi défini ?
ii. Aurait on choisi ce plan si, dès le départ, on avait su disposer de 16 unités ?

5 Méthode Tagushi
Certains facteurs d’un expérience que l’on peut controler en laboratoire peuvent ne
pas être contrôlés en utilisation normale. Un des apports de Tagushi est d’avoir proposé
des plans pour étudier l’influence de ces facteurs sur la variabilité du résultat.

Exemple 3 on veut régler un fraise pour une certaine performance donnée. Mais dans
l’utilisation future de cette fraise il y aura certain facteurs incontrôlés : la température de
l’atelier, le degré d’usure de la fraise, la température de l’huile de refroidissement.

On construit ainsi un plan fractionnaire pour les facteur contrôlés en utilisation normale
et un autre plan fractionnaire pour les facteurs non contrôlés en utilisation normale. On
“croise” ensuite les plans. Si le dispositif doit être réglé sur une spécification précise. On
va rechercher le réglage de facteurs contrôlés en utilisation normale qui minimise l’Erreur
Quadratique Moyenne (EQM), c’est à dire le carré du biais plus la variance où ce biais et
cette variance se calculent sous la loi de probabilité donnés par le plan sur les facteurs non
contrôlés en utilisation normale.
 28 CHAPITRE 3. PLANS FRACTIONNAIRES

Plans fractionnaire pour p facteurs à deux niveaux

Nombre Nombre de facteurs p
unités 3 4 5 6 7 8 9 10 11
III
23−1

4 6 6 6 6 6 6 6 6

C=AB
III III
IV III
25−2 26−3
24−1 27−4

Plan
8 6 6 6 6
Complet D=AB
E=AC
D=AB
F=BC
D=AB E= AC
D=ABC G=ABC
E=AC F=BC
V IV IV IV III III III
25−1 26−2 27−3 28−4 29−5 210−6 211−7
E=ABC
E=ABC F=BCD
Plan
16 E=ABC F=BCD G=ACD
Complet
E=BCD F=BCD G=ACD H=ABD
E=ABC F=ACD G=ACD H=ABD I=ABCD
E=ABC F = BCD G=ABC H=ABD I=ABCD J=AB
E=ABCD F=BCD G = ACD H=ABD I=ABCD J=AB K = AC
IV IV IV
VI IV IV
29−4 210−5 211−6
26−1 27−2 28−3

Plan F=ABC
32
Complet F=ABCD G=BCD
F=BCDE G=ACDE H=CDE
F=ABC
G=ACDE H=ABDE I=ACD
F=AB F=ABCD G=ABD
H=ABDE I=ACDE J=ADE
CDE G=ABDE H=BCDE
I=ABCE J=BCDE K =BDE
IV IV
VII V IV
210−4 211−5
27−1 28−2 29−3

Plan
64
Complet G=CDE
G=BCDF H=ABCD
G=ABCD
H=ACDE I=ABF
G=ABC G=ABCD H=ACEF
I=ABDE J=BDEF
DEF H=ABEF I=CDEF
J=ABCE K=ADEF
V
VIII VI V
211−4
28−1 29−2 210−3

Plan
128
Complet H=ABCG
I=BCDE
H=ACD H=ABCG
J=ACDF
H=ABC FG I=BCDE
K=ABC
DEFG I=BCEFG J=ACDF
DEFG

Tab. 3.1 – Tableau de plans fractionnaires de minimum d’aberration. En ligne, le nombre

d’unités 2p−q , en colonne le nombre de facteurs, une solution est donnée lorsqu’il existe
un plan de résolution III au moins. Chaque cellule indique de haut en bas : la résolution,
le nom du plan, les relations de définition.
Chapitre 4

Surfaces de Réponses, plans

isovariants

1 Cadre de l’étude
On considère une variable Y éventuellement influencée par m variables quantitatives
X1 , ..., Xm . On note x := (x1 , ..., xm ) un réalisation possible de l’ensemble des variables
explicatives. On suppose :
• E(Y ) au point x est une fonction polynomiale de degré q de x
• Toutes les observations de Y sont non-corrélées et de même variance (même si on
répète plusieurs fois le même x.
En ordonnant arbitrairement les unités du plan, on obtient un modèle linéaire classique
non-gaussien.
E(Y ) = Xβ ; Var(Y ) = σ 2 Id. (4.1)
On note p la dimension de β. Certaines colonnes de la matrice X sont liés fonctionnel-
lement : certaines sont des produits ou puissances d’autres. Nous supposerons qu’elle ne
sont pas liées linéairement de sorte que le modèle (4.1) est régulier. On note X(x) la ligne
de X qui correspondrait à l’observation au point x (Ce point ne fait pas forcément partie
des unités du plan correspondant au modèle (4.1) ).

Définition 8 On dit que le plan est isovariant (rotatable en anglais) si Var(X(x) β̂) ne
dépend de x qu’a travers |x| où β̂ est l’estimateur classique des moindres carrés.

2 Conditions d’isovariance
Pour examiner les conséquence de la définition nous introduisons les notations sui-
vantes.
-Si P est une transformation orthogonale de IR m il existe un endomorphisme (une
matrice si on préfère) unique QP tel que
 0
QP (X 0 (x)) = X(P (x)) . (4.2)

Comme P 0 = P −1 , il est facile d’en déduire que QP 0 = (QP )−1 .

29
30 CHAPITRE 4. SURFACES DE RÉPONSES, PLANS ISOVARIANTS

Lemme 1 Soit A une matrice symétrique définie positive de taille p, p. Alors il y a

équivalence entre
• T r(A) = T r(A−1 = p
et
• A = Id

Nous laissons la démonstration de ce lemme à titre d’exercice.

Nous notons M la matrice des moments, qui n’est rien d’autre que la matrice d’infor-
mation (qui sera elle notée N ) divisée par le nombre de données :
1 0 1
M= XX= N
n n
Les éléments de M sont du type
n
1X
[δ1 , ..., δm ] := X δ1 , ..., Xim
δm
,
n i=1 i1

où Xji est l’observation de la variable j sur l’unité i. La matrice M contient les éléments
ci-dessus pour δ1 + ... + δm ≤ 2q

Proposition 6 Le plan est isovariant si et seulement si M est isovariante c’est à dire

que pour toute transformation orthogonale P : M = Q P M Q0P

Démonstration : Condition suffisante : supposons M isovariante. Pour simplifier on

note X(x) écrit en vecteur ligne par : X et Q P : Q. On a alors

X(P (x)) = XQ0 ,

et également
1   1
2
Var (P (x)) β̂ ) = XQ0 2 Var(β̂)QX0 = XQ0 (X 0 X)−1 QX0 = X(X 0 X)−1 X0 .
σ σ

Passons maintenant à la condition suffisante : nous supposons donc maintenant que le

plan est isovariant. Remarquons d’abord que

T r(Var(X β̂)) = σ 2 T r(X(X 0 X)−1 X 0 ) = σ 2 T r((X 0 X)−1 X 0 X) = σ 2 p.

Soit maintenant le plan transformé par P . L’équation 4.2 montre que la matrice de ce plan
vaut X(QP )0 . Comme le plan est isovariant (vérifier ce dernier point en détail)

Var(X(QP )0 β̂) = Var(X β̂).

Donc

σ 2 p = T r(Var(X β̂)) = T r(Var(X(QP )0 β̂)) = σ 2 T r(X(QP )0 I −1 (QP )X 0 )

et
p = T r((QP )0 I −1 (QP )I) (4.3)
2. CONDITIONS D’ISOVARIANCE 31

On obtient exactement le même résultat en remplaçant I par M . On pose donc

A := (M 1/2 Q0P M −1 QP M 1/2 )
La relation (4.3) implique que T r(A) = p. Il est facile de vérifier par ailleurs que
A−1 = (M −1/2 (QP )−1 M (Q0P )−1 M −1/2 ).
Donc T r(A−1 ) = T r((Q0P )−1 M −1 (QP )−1 M −1 ) = p, car (QP )−1 = QP 0 et il suffit d’appli-
quer le lemme. 2.
Nous en venons maintenant au résultat principal de cette section.

Théorème 2 (Box et Hunter 1957) Avec nos notations, la matrice des moments est
isovariante si et seulement si tout moment d’ordre δ = m j=1 δj ≤ 2q
P

(i) est nul si l’un des δj est impair

(ii) est égal à
j δj !
Q
µδ
2, δ/2 j (δj /2)!
Q

sinon.
Les µδ sont certaines constantes.

L’interprétation de ce théorème peut se faire sur l’exemple incontournable des modèles

quadratiques qui correspondent à q = 2. Dans ce cas on doit considérer les moments
jusqu’à l’ordre 4. Compte tenu de (i), les moments non nuls ne peuvent être que ceux qui
correspondent à
- Un δj = 2, les autres δj 0 nuls. Il doivent être tous égaux et le théorème nous dit que
l’on note µ2 leur valeur commune.
- Deux δj égaux à 2, les autres nuls. Le théorème nous dit que l’on note maintenant
µ4 leur valeur commune .
- Un δj = 4, les autres δj 0 nuls. Il doivent être tous égaux et leur valeur est maintenant
3µ4 .
Nous reviendrons sur cet exemple ultérieurement. Démonstration du théorème :
On pose
1X
F (t) := (1 + xi1 t1 + ... + xim tm )2q
n i
Un petit peu de combinatoire montre que le monôme d’ordre (δ 1 , ..., δm ) a pour coefficient
(2q)!
[δ1 , ..., δm ]
(2q − δ)! j δj !
Q

Supposons la matrice des moments isovariante, alors les [δ 1 , ..., δm ] le sont également et
cela implique que pour toute transformation orthogonale P :
F (P (t)) = F (t)
F est en fait une fonction de la valeur absolue de t et comme c’est un polynôme
0
δ!
2 δ
t2δ 2δm
X X X
F (t) = aδ (ktk ) = aδ 1 ...tm
1

j=1,m δj !
Q
δ=0,q δ=0,q
P0
où la somme porte sur les δj ≤ 0 tels que δ1 + ... + δm = δ. En identifiant on obtient
le résultat. 2.
32 CHAPITRE 4. SURFACES DE RÉPONSES, PLANS ISOVARIANTS

3 Plans composites centrés de Box et Wilson

Nous allons nous concentrer notre attention sur les modèles quadratiques (q = 2) pour
m facteurs. Un plan composite est constitué de trois parties :
– une partie factorielle composée du plan factoriel complet {−1, +1} m (on peut mon-
trer qu’en fait une fraction régulière de résolution V suffit).
– d’une partie axiale composée de 2m points à distance α de l’origine et situés sur les
axes (α reste à déterminer).
– la partie centrale composée de n0 points au centre 0 (n0 reste à déterminer).
La distribution du plan est invariante
– par toute symétrie par rapport à une hyperplan d’équation x j = 0.
– Ainsi que par toute permutation des coordonnées.
Deux conséquences immédiates sont que [δ 1 , ...δm ] est nul dès qu’un δj est impair et que
[δ1 , ...δm ] = [δσ(1) , ...δσ(m) ]
Soient F le nombre de points de la partie factorielle du plan, il est immédiat que la
matrice d’information N vérifie (en notant dans l’ordre : la constante, les effets linéaires,
les carrés et enfin les produits
 
n 0. ... .0 b. ... .b 0. ... .0

 0 b 0... .0 0. ... .0 0. ... .0 

... ... ... ... ... ... ... ... ... ...
 
 
b 0. ... .0 c d... .d 0. ... .0
 
 
 
M =
 ... ... ... ... ... ... ... ... ... ... 
 (4.4)

 b 0 0... .0 d ...d c 0. ... .0 

0 0. ... .0 0. ... .0 a 0.. .0
 
 
... ... ... ... ... ... ... ... ... ...
 
 
0 0. ... .0 0. ... .0 0 ..0 .a
avec

x2i,1 = F + 2α2 ; d := a := x4i,1 = F + 2α4

X X X
b := xi,1 xi,2 = F ; c :=
1=1,n 1=1,n 1=1,n

le plan est donc isovariant si α = F 1/4

Remarquez que pour m = 2 et m = 4 les points de la partie axiale et factorielle
sont situés sur la même sphère. Pour les autres valeurs de m les points ont seulement
qu’approximativement la même norme.
Nombre de variables m 2 3 4 5 5 6 6
Taille p. fact. F 4 8 16 32 16 64 32
taille p. axiale 4 6 8 10 10 12 12
Point centraux
isovariance ≤1 ≤1 ≤1 ≤1 ≤1 ≤1 ≤1
précision uniforme 5 6 7 10 6 15 9
orthogonalité 8 12 12 17 10 24 15
α 1.41 1.68 2 2.38 2 2.83 2.38
Tableau 2 :Récapitulatif des paramètres et des tailles des plans composites centrés.
Chapitre 5

Bibliographie

Box, G., Hunter, W.G., Hunter J.S. (1978). Statistics for experimenters. Wiley, New-
York.
Coursol, J. (1980). Technique statistique des modèles linéaires. Cimpa, Nice.
Dacunha-Castelle, D. & Duflo, M. (1982, 1983). Probabilités et Statistiques , 1 Problèmes
à temps fixe , 2. Problèmes à temps mobiles. Masson, Paris.
Droesbeke, J-J., Fine J. & Saporta G. (1997) (Ed.). Plans d’expériences, Application à
l’entreprise. Technip.
Jiang, J. (1996). REML estimation : Asymptotic behaviour and related topics. Annals of
statistics, 24, 255-286.
Raghavarao , D. (1971). Construction and combinatorial problems in design of experiments.
Wiley, New-York.

33
34 CHAPITRE 5. BIBLIOGRAPHIE
Table des matières

1 Modèles mixtes 3
1 Analyse de la variance multivariable . . . . . . . . . . . . . . . . . . . . . . 3
2 Modèles à effets aléatoires et mixtes . . . . . . . . . . . . . . . . . . . . . . 5
2.1 Facteurs croisés et hiérarchisés . . . . . . . . . . . . . . . . . . . . . 5
2.2 Modèles mixtes équirépétés . . . . . . . . . . . . . . . . . . . . . . . 6
2.3 Modèles mixtes généraux . . . . . . . . . . . . . . . . . . . . . . . . 7
2.4 Estimation des effets fixes . . . . . . . . . . . . . . . . . . . . . . . . 8
2.5 Estimation par MIVQUE dans un modèle mixte . . . . . . . . . . . 8
2.6 Estimation par maximum de vraisemblance restreinte . . . . . . . . 9
2.7 Tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

2 Plans d’expériences randomisés 13

1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2 Nécessité de la randomisation . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3 Plans d’expériences classiques . . . . . . . . . . . . . . . . . . . . . . . . . . 15

3 Plans fractionnaires 21
1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2 Cadre général pour des facteurs à deux niveaux . . . . . . . . . . . . . . . . 22
3 Méthode des facteurs de base . . . . . . . . . . . . . . . . . . . . . . . . . . 24
4 Plan pour l’étude des effets principaux et des interaction doubles . . . . . . 25
5 Méthode Tagushi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

4 Surfaces de Réponses, plans isovariants 29

1 Cadre de l’étude . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2 Conditions d’isovariance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
3 Plans composites centrés de Box et Wilson . . . . . . . . . . . . . . . . . . 32

5 Bibliographie 33

35