Cours plan pour surface de réponse (plan continu)

Bonjour, donc on va étudier aujourd'hui les plans pour surface de réponse ou plans pour variables quantitatives. Donc l'objectif, ça va être toujours d'étudier une variable quantitative Y qui est une variable réponse et cette fois de choisir des expériences pour voir les effets de variables quantitatives X. Et on va chercher à modéliser la surface de réponse, la réponse de Y en fonction de cette variable quantitative. Et pour modéliser cette surface de réponse, eh bien on va avoir besoin de faire quelques expériences, donc pas trop d'expériences puisqu'on est dans un cadre de planification expérimentale, mais des essais bien choisis pour qu'on puisse modéliser la réponse en fonction des variables quantitatives de façon la plus précise possible. Donc petit rappel d'abord de régression linéaire simple, donc le cas où on s'intéresse à 2 variables quantitatives, une variable quantitative Y qu'on cherche à expliquer en fonction d'une variable quantitative X. Donc on écrit le modèle sous la forme Yi= bêta 0+ bêta 1 Xi+ epsilon i. Donc on modélise la relation linéaire entre X et Y et puis il y a un résidu epsilon i et avec des hypothèses sur les résidus, ils sont indépendants, de même variance et d'espérance nulle. Et quand on va estimer les paramètres de ce modèle, eh bien on va les estimer par grâce au critère des moindres carrés, donc en essayant de minimiser l'écart entre les observations Yi et la modélisation par le modèle, c'est-à-dire les bêta 0 chapeau plus bêta 1 chapeau Xi. Bêta 0 chapeau et bêta 1 chapeau, ce sont les estimations des paramètres du modèle, des estimations qui justement vont chercher à minimiser l'écart entre observations et prévisions. Plus exactement, vont minimiser le carré des écarts, la somme des carrés des écarts entre les observations et les prévisions. Si on veut estimer ces paramètres, eh bien on peut dériver cette somme et on va avoir que bêta 0 chapeau est égal à Y barre moins bêta 1 chapeau X barre. Et bêta 1 chapeau, c'est une formule qui nous donne en gros la covariance entre X et Y divisée par la variance de X. Alors je dis en gros parce que normalement on devrait mettre 1 sur n- 1 devant les 2 signes somme. Mais ces quantités 1 sur n- 1 étant au numérateur et au dénominateur, elles s'éliminent et il reste juste la formule qu'on a écrite ici. La variance de bêta chapeau, elle peut s'écrire comme sigma carré divisé par la somme des Xi moins X barre au carré et ça c'est là encore la variance de X multipliée par n- 1. Donc qu'est-ce qu'on peut voir grâce à cette formule? Eh bien que la précision sur l'estimation du coefficient, c'est-à-dire la variance coefficient est d'autant plus petite que dénominateur est grand. Donc n, c'est quoi? C'est le nombre d'expériences. Donc on se doute bien que plus on fera d'expériences, plus on sera précis dans l'estimation des paramètres. Donc ça, c'est assez logique. Et on voit aussi que quand la variance de X est grande, alors l'estimation des paramètres est plus précise. Et donc là, ça veut dire quoi variance de X grande? Ça veut dire qu'on va prendre des X qui sont très différents. Alors je parle pas du numérateur sigma carré parce que ce sigma carré c'est la variance résiduelle et on peut pas vraiment jouer dessus. Mais là on voit que on peut jouer sur le choix des expérimentations, sur le choix des X pour faire en sorte que la variance des bêtas soit la plus petite possible et on peut choisir des X qui soient suffisamment différents. Si on se place dans le cadre de la régression multiple, bien le modèle s'écrit quasiment de la même façon. On va écrire réponse Y pour l'individu i, c'est bêta 0 plus bêta 1 fois la valeur prise par la variable 1 pour l'individu i, donc qu'on va noter Xi1, plus bêta 2 Xi2, plus etc., jusqu'à plus bêta p Xip. Donc ça, c'est pour les variables, la p-ième variable. Et puis on a toujours un résidu epsilon i, et on fait toujours les mêmes hypothèses sur ces résidus. Donc des hypothèses de, ils sont indépendants, de espérance nulle et de même variance sigma carré. D'un point de vue matriciel, on peut écrire ce modèle sous la forme Y= X bêta+ E avec X une matrice, donc une grosse matrice avec une colonne de 1, colonne de 1 qui va servir à estimer la constante, donc le bêta 0, et puis après toutes les valeurs de X. Donc à partir de cette formule matricielle, on va chercher à estimer toujours les paramètres, donc les bêtas, et les bêtas c'est toujours estimé par moindre carré ici, on va chercher à minimiser la somme des carrés des écarts entre les observations Yi et les prévisions par le modèle. Et en dérivant donc cette formule-là matriciellement, on peut trouver que c'est égal à X' que multiplie X'. Donc ça c'est l'estimation du paramètre. Et puis on a les propriétés d'estimation de ce paramètre, c'est-à-dire tout d'abord l'estimateur est sans biais, l'espérance de bêta chapeau est égale à bêta et la variance de bêta chapeau est égale à X prime X moins 1 sigma carré. Donc X prime X moins 1 sigma carré, eh bien là encore on voit que la variance du paramètre, la variance de bêta dépend des X et donc des expériences qu'on va réaliser. Si on choisit bien les X, on va pouvoir avoir une estimation précise des paramètres de notre modèle. Donc tout l'enjeu des plans est là, c'est trouver les X, les valeurs de X, les expériences à réaliser pour pouvoir ensuite estimer les paramètres de façon très précise. Alors ensuite, on peut, à partir des estimations des coefficients, prédire des valeurs de Y, soit pour les valeurs passées, soit pour de nouvelles valeurs de X et donc il suffit de remplacer dans la formule les coefficients par leur estimation donc bêta 0 chapeau, bêta 1 chapeau, etc. On a la variance de prédiction qui est égale à sigma carré que multiplie 1 plus X prime 0, X prime X moins 1, X0. Donc X0 c'est le vecteur des X là où on fait la prédiction et on multiplie par X prime X moins 1, X prime X moins 1 c'est toujours là où on a fait nos expérimentations. Et donc on voit que la précision sur la prédiction de Y dépend là encore des X qui ont servi à estimer les paramètres du modèle. Alors je vais prendre un petit exemple où je m'intéresse à 2 facteurs: un facteur température de cuisson, je sais que cette température de cuisson varie entre 120 et 140 degrés, et la durée de la cuisson qui varie entre 40 et 60 minutes. Et ce qui m'intéresse, c'est d'étudier une réponse y qui est par exemple le moelleux de pain de mie. Donc on va mesurer une quantité y qui mesure le moelleux de pain de mie. Et l'objectif de cette problématique, ça va être de trouver finalement le couple x1, x2 qui va permettre d'optimiser le moelleux du pain de mie. Donc les questions derrière, c'est toujours on a notre variable y réponse, on a nos variables x1, x2 et quels sont les effets de ces 2 facteurs, les effets de X1 et X2, quel est le rôle de chacune de ces variables dans la variabilité de Y, est-ce qu'il y a un facteur multiplicatif de X1 et X2 sur la réponse Y. Et puis ensuite, eh bien ça va être d'essayer de trouver un optimum, donc trouver le couple X1, X2 qui permet d'optimiser la variable Y. Et puis bien évidemment, quand on va construire ce genre de modèle, on veut être le plus précis possible, donc avoir les estimations les plus précises possibles des coefficients et par suite les estimations les plus précises possibles des prédictions pour de nouvelles valeurs de X. Alors le problème souvent c'est que nos variables X1, X2 ne jouent pas forcément de façon linéaire seulement sur Y. Par exemple, on peut imaginer que le moelleux dépend de la température de cuisson. La variable Y prend des valeurs faibles quand la température est très faible puis La valeur de Y augmente mais elle atteint un maximum et ensuite elle rediminue. Donc en pratique ici ce qu'on voit c'est qu'on a un effet qui est quadratique ou qui ressemble en tout cas à un effet quadratique. Donc on n'a pas juste un effet linéaire qui serait, eh bien, la variable Y augmente quand le X augmente. Et donc on va vouloir prendre en compte cet effet quadratique dans le modèle qui permet d'exprimer Y en fonction des X1, X2. Donc pour ce faire, bien on va construire un modèle avec à la fois des effets linéaires, donc des effets de X1 et de X2 seuls, puis des effets quadratiques, c'est-à-dire un effet de X1 carré, un effet de X2 carré, et puis un effet d'interaction entre X1 et X2. Donc le modèle, je vais l'écrire bien YI, donc la réponse pour l'expérience I, c'est bêta 0 plus bêta 1 fois X1 Xi1, la valeur prise par la première variable pour l'individu i, plus bêta 2 fois Xi2. Donc ça c'est le modèle linéaire pour l'instant. Ensuite je vais rajouter les termes quadratiques, donc je vais mettre bêta 11 par exemple, donc c'est un coefficient que je note bêta 11 Xi1 carré, bêta 22 Xi2 carré. Donc ça ce sont les, la prise en compte des termes, des 2 termes quadratiques. Et puis enfin un terme d'interaction, donc un coefficient bêta 1,2 que multiplie Xi1 Xi2. Ça, ça va être l'effet du produit de Xi1 par Xi2. Et puis on a toujours un résidu avec toujours les mêmes hypothèses sur les résidus. Alors l'interaction ici entre 2 variables quantitatives, ça signifie que l'effet d'une variable X1 sur Y dépend des autres valeurs de la variable X2. Par exemple, ça signifie que si je fixe X2 et que je cherche l'optimum de Y en fonction de X1 pour un X2 fixé, je vais avoir une valeur pour X1. Si je prends une autre valeur de X2, l'optimum de Y ne sera pas pour la même valeur de X1. Ça, ça voudra dire qu'il y a une interaction entre X1 et X2. Donc ici on va visualiser la modélisation de Y selon la nature de la relation entre X1, X2 et Y. Donc le graphe en haut à gauche, vous avez uniquement des effets linéaires de X1 et X2 et donc vous avez un plan qui est incliné. Le schéma suivant, vous avez un effet linéaire de X1, donc on voit un plan qui monte et puis vous avez un effet linéaire de X2 mais aussi un effet quadratique de X2. Ici, on a un peu comme une feuille de papier qui serait incurvée. Est sur le côté. Le graphe suivant, bien vous avez un effet linéaire de x1, donc les valeurs ont tendance à monter par rapport à x1, un effet linéaire de x2, mais également une interaction entre x1 et x2. Donc vous voyez une forme qui est incurvée et donc un optimum de x1 qui va dépendre des valeurs de x2. Même chose, un optimum de x2 qui dépend des valeurs de x1. Vous avez ensuite un un plan donc en creux ici avec des effets quadratiques de x1 et quadratiques de x2. Donc il y a un minimum ici puisque les coefficients devant x1² et x2² sont positifs. Ensuite vous avez une surface de réponse qui est en forme de selle de cheval où on a un effet quadratique qui a un effet positif et un effet quadratique, celui de x2, qui a un effet négatif. Et puis enfin, vous avez le dernier graphe qui est l'opposé de celui qu'il y a en bas à gauche où vous avez un optimum avec cette fois des effets des termes quadratiques qui sont négatifs. Donc la modélisation qu'on propose permet de prendre en compte tous ces types de surfaces de réponse. Donc par rapport à l'effet de 2 variables, on a un peu tous les types de surfaces de réponse possibles ici. Alors vous avez un petit exercice qu'il est possible de faire où on va vous demander d'optimiser là encore une recette de galette et d'optimiser 2 paramètres. Donc d'optimiser la quantité de farine sachant que cette quantité varie entre 45% et 55% et d'optimiser la température de cuisson qui varie entre 180 et 220 degrés. Donc vous avez ces 2 facteurs qui peuvent varier simultanément et le but c'est d'essayer de trouver 10 essais qui vont permettre de modéliser au mieux la réponse de Y en fonction de ces 2 facteurs. Donc là encore, vous avez une feuille Excel que vous pouvez télécharger et puis comparer le plan que vous allez construire par rapport à un plan qui est optimal. Et si vous n'arrivez pas à construire le plan, eh bien vous allez suivre le cours et après vous pourrez construire un plan assez facilement. Alors quelques idées d'abord sur la construction du plan. Donc on a dit que la qualité de la modélisation va dépendre de bêta chapeau et de la précision de bêta chapeau plus exactement, donc de la variance de bêta chapeau. Et la variance de bêta chapeau c'est X prime X moins un sigma carré. Donc sigma carré on peut pas y toucher, mais X prime X moins un on va pouvoir jouer dessus pour trouver les bonnes valeurs de X qui permettent de minimiser la variance de bêta chapeau. Donc on va pouvoir connaître la qualité du plan avant même de faire les expériences parce qu'une fois qu'on a défini les expériences, on connaît X et donc on peut calculer X prime X moins un et donc avoir une idée de la qualité de notre plan. Alors première idée qu'on a déjà évoquée tout à l'heure, c'est il faut prendre des X qui soient les plus dispersés possibles. Alors on le voyait sur la formule, mais comment on peut le dire autrement? Bien, ici j'ai fait un schéma à gauche avec des points bleus et puis à droite avec des points rouges. Il y a une même variabilité dans la position des points, des points bleus qui sont à gauche par rapport aux points rouges qui sont à gauche et des points bleus qui sont à droite par rapport aux points rouges qui sont à droite. Donc c'est la même variabilité mais dans un cas j'ai séparé beaucoup plus les points bleus que les points rouges. Et si je trace des droites entre 2 points pris au hasard, entre 2 points bleus ou entre 2 points rouges, eh bien je vois que les pentes de mes droites sont beaucoup plus homogènes à partir des points bleus qu'avec les points rouges. C'est très logique, les points sont plus éloignés et donc on va avoir une précision qui est plus importante. C'est exactement ça quand on va dire qu'on veut une variance des X la plus grande possible, on veut plus de précision dans l'estimation de la pente de la droite. C'est exactement ça qu'on est en train de faire là. Maintenant, si je prends des points uniquement à gauche et à droite dans la figure bleue ici, eh bien je peux avoir une relation entre Y et X qui est linéaire, mais je peux avoir n'importe quel autre type de relation. Donc ici, j'ai mis une relation quadratique et on voit qu'on peut avoir différentes formes de modèles. Et donc pour tester la linéarité, j'ai envie d'avoir des essais au centre et pas uniquement des essais aux extrémités. Donc là on voit bien l'idée d'avoir aussi des essais au centre et pas seulement aux extrémités de X. Et puis bien si on veut s'intéresser à 2 facteurs, on va vouloir ce qu'on appelle l'orthogonalité entre 2 facteurs. Alors si j'ai 2 facteurs, la variance de bêta 1 chapeau, la variance de l'estimateur de bêta 1, c'est sigma carré divisé par n que multiplie 1 moins r12 que multiplie la variance de X1. R12, c'est la valeur absolue ici de la corrélation entre X1 et X2. Et donc ce qu'on voit bien ici, c'est que si la valeur absolue de la corrélation est proche de 1, 1 moins r12 va être proche de 0 et donc la variance va augmenter sensiblement. Et donc la variance de bêta 1 chapeau sera petite si la corrélation entre X1 et X2 est proche de 0. Si elle est égale à 0, la variance de bêta 1 chapeau est égale en fait à la formule qu'on avait quand on avait la régression simple, c'est-à-dire à la variance de bêta 1 chapeau si j'avais fait la régression uniquement en fonction de la variable X1. Ça, ça veut dire quoi? Ça veut dire que je vais pouvoir estimer 2 coefficients, le coefficient bêta 1 et le coefficient bêta 2, avec la même précision que si j'avais fait mes expérimentations pour n'estimer que bêta 1. Donc ça veut dire que je vais avoir une estimation qui sera précise. Et s'il y a une corrélation entre X1 et X2, eh bien la variance de bêta 1 chapeau augmente et donc je suis moins précis. Alors, pour construire les plans d'expériences avec des variables quantitatives, eh bien, on va d'abord s'affranchir des unités. Donc tout à l'heure, on avait des températures, on avait des pourcentages d'humidité. On va en fait réduire chacune des variables pour construire des plans qui soient sans unité et donc proposer des plans qui soient toujours les mêmes. Alors, comment est-ce qu'on fait pour transformer les variables et les ramener sans unité, eh bien en fait on va s'arranger pour que le minimum soit égal à -1 et le maximum soit égal à 1. Donc on va construire des plans où x1 va varier entre -1 et 1 et x2 va varier entre -1 et 1 également. Donc comment est-ce qu'on va faire ce changement de variable? On va prendre x et puis soustraire la moyenne, le centre du domaine, xmax+ xmin/ 2, ça c'est la moyenne du domaine et puis xmax moins xmin sur 2, c'est l'étendue divisée par 2. Donc en fait, on peut facilement calculer les nouvelles valeurs de x une fois qu'on s'affranchit des unités avec cette formule. Une fois qu'on aura fait ça, on va construire toujours en fait le même plan, toujours la même formule et il faudra juste après recoder pour retrouver les unités. Et ce qui est intéressant aussi avec ce codage, c'est qu'on a une interprétation des coefficients du modèle qui est assez intéressante. Donc les interprétations sont les suivantes: si mon modèle s'écrit Y= bêta 0+ bêta 1 X+ bêta 11 X², j'ai la même chose si j'ai plusieurs variables, mais là je vais m'intéresser à une variable, c'est exactement pareil pour les autres variables. Donc bêta 0, c'est quoi? Si j'ai fait le changement de variable proposé, bêta 0 c'est égal à y au centre du domaine, donc c'est égal à y au point x= 0. Bêta 1 ce sera égal à quoi? Pour calculer bêta 1, en fait je vais juste avoir besoin de calculer y pour le point 1 et y pour le point -1. Quand x= 1, y vaut bêta 0+ bêta 1+ bêta 1 1. Quand x égale -1, y en -1 c'est égal à bêta 0 moins bêta 1 plus bêta 11. Et donc bêta 0 c'est la valeur de y au centre du domaine, bêta 1, il suffit que je fasse y en +1 moins y en -1, les bêta 0 et les bêta 11 vont s'éliminer, je vais me retrouver avec 2 bêta 1, j'ai plus qu'à diviser par 2 et j'ai bêta 1 qui sera égal à la la moyenne entre ce qui se passe en Y quand X vaut 1 moins ce qui se passe en Y quand X vaut -1. Et puis bêta 1, 1, en fait je vais sommer Y1 et Y-1. En faisant ça, je vais avoir 2 bêta 0 plus 2 bêta 1. Et donc bêta 1, 1, c'est la moyenne de ce qui se passe en Y quand X vaut 1 et quand X vaut -1, la moyenne de cette quantité moins Beta zéro. Alors, un petit schéma pour reprendre tout ça. Donc j'ai dit que beta zéro c'était la valeur de y quand x vaut zéro. Donc là, ici, si j'ai le point, si j'ai mes trois expériences pour x égale moins un, x égale zéro et x égale un, beta zéro c'est juste la valeur de y quand x vaut zéro. Beta 1, c'est quoi? C'est l'écart qu'il y a entre y quand x vaut 1 et y quand x vaut -1, mais c'est cet écart-là divisé par 2. Donc j'ai pris l'écart entre ce qui se passe en 1 et ce qui se passe en -1 et je divise par 2 pour voir l'effet de beta 1. Et puis là, a11, c'est quoi? C'est l'écart entre la moyenne de ce qui se passe en 1 et en -1 moins bêta zéro, donc moins le y zéro ici. Donc la moyenne de ce qui se passe entre un et moins un, c'est ce point-là, moins y zéro, ça va me donner une idée de bêta un un. Donc on visualise en fait directement grâce aux coefficients la relation qu'il y a entre y et x. Si bêta un un est petit, ça veut dire que le y zéro va être quasiment sur la droite entre y en 1 et y en -1. Et si le bêta 1,1 est très grand, ça veut dire qu'au contraire, on a une courbe qui est très en cloche. Et donc voilà quelques exemples d'interprétation. Donc avec bêta 0 dans le graphe de gauche, bêta 1 qui est donné et bêta 1,1 qui vaut 0. Donc pour le deuxième schéma, bêta 0 vaut -0.8 puisque l'Y0 vaut -0.8. Bêta 1, c'est l'écart entre 1.2 et -0.8, donc ça fait un écart de 2 divisé par 2, ça fait 1. Et puis bêta 1,1, c'est l'écart entre la moyenne entre 1.2 et -0.8, c'est-à-dire 0,2 et puis je fais ça moins moins 0,8 qui est le bêta 0 ou le y0. Voilà et on peut retrouver comme ça tous les coefficients. Donc pour le troisième schéma vous avez y0 qui vaut 0 donc bêta 0 vaut 0, y1 égal y moins 1 donc l'écart entre les deux c'est 0 divisé par 2 c'est toujours 0 donc bêta 1 vaut 0 et puis ensuite bêta 11 je vais avoir 0,5, la moyenne entre Y1 et Y-1, c'est 0,5, moins 0, moins bêta 0 ou moins Y0. Donc bêta 1, 1 qui vaut 0,5. Et puis enfin, bien, pour le dernier schéma, vous avez bêta 0 qui vaut 0, bêta 1 c'est -2 moins -2, ça fait 0, et puis bêta 1, 1 c'est la moyenne de -2 et -2, Ça fait -2- 0, donc c'est -2. -2X². Voilà, donc on retrouve tout de suite les coefficients et la forme du modèle quand on voit les coefficients. Alors comment construire maintenant un plan d'expériences? Eh bien, il y a un plan qui s'appelle les plans composites centrés qui se construisent de la façon suivante. Donc quand on étudie K facteurs, on va construire un plan factoriel complet éventuellement un plan fractionnaire 2^k ou 2^. Donc 2^k, c'est quand on prend un plan factoriel complet et 2^ quand il est possible de prendre un plan fractionnaire. Donc le plan factoriel complet va nous permettre de construire un carré en 2 dimensions, un cube en 3 dimensions, etc. Donc ici, ça correspond aux points bleus qui sont sur le carré. Donc on retrouve les 4 points du plan factoriel 2^2. Ensuite, on met ce qu'on appelle des points en étoile. Point en étoile, je vais prendre alpha qui est la racine quatrième du nombre d'essais du plan que j'ai construit précédemment. Donc dans mon plan complet, j'avais 4 points. Donc là, je vais prendre la racine quatrième de 4. La racine quatrième de 4, c'est la racine de racine de 4, c'est-à-dire racine de 2. Donc les points en étoile, eh bien je mets tous les facteurs à 0 sauf 1. Que je mets à racine de 2 et moins racine de 2. Donc ça, ça correspond au point rouge. Donc en fait, on va construire en mettant des 0 pour tous les facteurs sauf un qui prend la racine quatrième du nombre d'essais du plan factoriel. Et puis on rajoute des points au centre, c'est-à-dire des points en 0. Et donc là, on rajoute un certain nombre de points au centre selon la qualité du plan qu'on souhaite avoir. Donc après, le nombre d'expériences totales, eh bien c'est dans mon exemple 4 pour le plan factoriel plus 4 pour les points en étoile plus un certain nombre de points au centre. Donc c'est de cette façon qu'on va construire les plans composites centrés. Alors si on veut construire des plans composites centrés avec R, on va utiliser le package RSM et il suffit d'utiliser la fonction C. C, D et de préciser le nombre de facteurs qu'on souhaite étudier. Donc CCD de 2, ça me construit le plan composite centré, Central Composite Design, à 2 facteurs. Et puis bien la ligne d'en dessous, vous voyez CCD de 2 et je précise en fait les variables, donc en précisant que la première variable c'est une variable X1 qui est la température avec je donne le centre du domaine, donc le centre c'est 130 et divisé par 10, ça c'est la moitié de l'étendue puisque je voulais que les températures varient entre 120 et 140. Et puis même chose, la deuxième variable c'est la durée qui est centrée autour de 50 avec une variabilité de 10, donc ça peut aller de 40 à 60. Et donc quand je regarde le plan d'expérience qui a été construit, eh bien vous voyez que j'ai les essais avec les points au centre. Donc j'ai un certain nombre de points avec 2 points au centre avec 130 en température et 50 pour le temps. Et puis ensuite, j'ai les points du plan factoriel 2-2. Donc par exemple, ce sont les valeurs 120, 40, 120, 60, 140, 60 et 140, 40. Ça, ce sont les points du plan factoriel. Et puis vous avez les points en étoile. Avec 50, 115, 144, 50, 130, 64,14 et 130, 35,8. Alors quelles sont les propriétés de ce plan composite centré? Eh bien une première propriété qui est importante qui est l'isovariance par rotation. Cette propriété on l'a si on choisit pour les points en étoile de les mettre à la racine quatrième du nombre de d'essais qu'on a réalisés précédemment dans le plan factoriel. Cette isovariance par rotation, ça veut dire que la précision du plan dépend de la distance au centre du domaine. J'ai le centre du domaine qui est ici et la précision dépend uniquement de la distance au centre du domaine. Ça ne dépend pas du tout de la direction, que j'aille à droite, en haut, à gauche, en bas, la qualité du plan, la qualité de la précision est la même quelle que soit la distance au centre du domaine. Ça, c'est ce qu'on appelle l'isovariance par rotation. On a aussi une autre propriété qui est la précision uniforme, c'est-à-dire que la précision est identique à une distance de 1 du centre du domaine, quelle que soit la direction et quelle que soit la distance au centre, à une distance inférieure à 1. Autre propriété qui est très importante, c'est que les différents effets sont orthogonaux entre eux, c'est-à-dire qu'il n'y a pas de corrélation entre les effets principaux. Entre eux. Mais il y a une corrélation entre les effets quadratiques en fonction du nombre de points au centre qu'on va choisir. Donc si on choisit bien le nombre de points au centre N0, on va avoir une corrélation qui sera très faible entre les effets quadratiques. Donc en pratique, comment on construit le plan d'expériences? Eh bien, on définit les expériences qu'on doit réaliser et puis après, les points au centre, eh bien, on les met un petit peu partout On va pas faire tous les points au centre au début ou à la fin par exemple. Il faut qu'on en mette un petit peu partout. Maintenant, il faut aussi éventuellement s'adapter à la réalité du terrain. Par exemple, on va faire toutes les expériences à 140 degrés à peu près au même moment pour éviter d'avoir une fois à chauffer le four, une fois à le faire refroidir, ensuite à le faire réchauffer, etc. Donc il faut aussi prendre en compte ces contraintes pratiques quand on construit des plans d'expériences. Idéalement, statistiquement, il faudrait choisir au hasard l'ordre dans lequel on réalise nos essais. Mais après, il faut aussi prendre en compte les contraintes techniques. Alors voici un petit tableau qui permet de voir les plans composites centrés en fonction du nombre de facteurs qu'on étudie. Donc là, vous avez un tableau avec de 2 à 6 facteurs qu'on va étudier. Et donc je vous donne le plan factoriel complet ou le plan fractionnaire qu'on va utiliser pour construire le plan composite centré. Donc vous voyez que pour 2, 3, ou 4 facteurs, on prend le plan complet, donc le plan 2 puissance 2, 2 puissance 3 ou 2 puissance 4. Avec 5 facteurs, on va prendre le plan 2 puissance 5 moins 1. Et avec 6 facteurs, on va prendre le plan 2 puissance 6 moins 1. Ça nous donne un nombre d'essais pour ce plan fractionnaire ou ce plan complet. Donc avec 3 facteurs par exemple, on aura 8 essais ici. Ensuite, on va déterminer la valeur du alpha pour les points en étoile. Donc je vais prendre la racine quatrième de 4, la racine quatrième de 8, la racine quatrième de 16. Donc racine quatrième de 16, c'est la racine de racine de 16. Racine de 16, c'est 4. Racine de 4, c'est 2. Donc la racine de racine de 16, c'est 2. Ensuite, vous avez le nombre de points axiaux, donc avec les points, ce qu'on appelle les points en étoile. Eh bien, c'est le nombre de facteurs multiplié par 2 puisqu'il y a une valeur minimum, une valeur maximum. Donc pour 3 facteurs par exemple, eh bien, on aura 6 points en étoile. Et puis ensuite, on rajoute des points au centre. Et puis bien là, selon qu'on veut l'orthogonalité ou qu'on veut une précision uniforme dans notre plan, eh bien, on va choisir un nombre de points au centre qui sera un petit peu différent. Alors, une fois qu'on a construit un plan d'expériences, idéalement, avant de se lancer dans la réalisation des expériences, le mieux, c'est d'évaluer la qualité de ce plan. Pour évaluer la qualité de ce plan, eh bien on va calculer X prime X moins 1. Donc on va calculer la valeur qui, multipliée par sigma carré, donnera la variance des paramètres. Et si ce X prime X moins 1 est petit, eh bien ça veut dire qu'on aura des plans qui seront de bonne qualité et donc des estimations pour les modèles qui seront précises avec des variances qui seront petites sur les estimateurs des paramètres. Ici j'utilise le package RSM pour construire un plan et je construis un plan composite centré avec 2 facteurs. La fonction model.matrix me permet de calculer la matrice X pour des essais réalisés et pour un modèle donné. Donc là le modèle je l'écris X1+ X2+ alors on écrit I entre parenthèses X1 carré. Ça, ça signifie qu'on va calculer x1 carré et qu'on va définir un nouveau facteur x1 carré. Parce que sinon, dans R, le x1 carré n'est pas compris. Il est compris comme l'interaction entre x1 et x1 et c'est confondu complètement avec x1. Donc là, le grand I majuscule qui est une fonction va faire en sorte qu'on calcule x1 au carré et qu'on le considère comme un facteur. Même chose, je calcule x2 carré, je le considère comme un facteur et je calcule le produit de x1 et x2 et je le considère comme un facteur. Et puis je dis que je veux utiliser ce modèle-là sur les données du plan que je viens de construire. Et donc en calculant transposé de X multiplié par X, donc c'est-à-dire X prime X, je vois que j'ai une matrice qui est presque diagonale. Et ce qui m'intéresse surtout, c'est l'inverse de cette matrice. Donc la fonction qui permet d'inverser une matrice dans R, c'est la fonction solve. Donc si je fais solve de transposé de X pourcentage étoile pourcentage Donc ça, ça me calcule X prime X moins un. Eh bien, je vois que j'ai une matrice X prime X moins un qui est presque parfaite parce que bien la constante est indépendante des effets principaux parfaitement. Elle est légèrement liée aux effets quadratiques, mais très très peu. On voit que les effets principaux X1, X2 sont indépendants de tous les autres effets, de la constante et de tous les autres effets. Les effets quadratiques sont légèrement liés à la constante et puis le produit de X1, X2 est indépendant des autres effets. Donc on a construit une matrice X, on a choisi des essais qui vont permettre d'estimer tous les coefficients sans qu'il y ait de corrélation, de liaison entre les coefficients. Donc on ne peut pas faire mieux ici que cette matrice. Ou alors on fait plus d'essais et dans ce cas-là, en faisant plus d'essais, on sera plus précis et on aura des variances, donc des valeurs sur la diagonale qui seront plus petites. Mais pour un nombre d'essais relativement limité, eh bien ici on a une matrice de qualité. Et donc après, eh bien ce qu'on va faire, c'est qu'on va réaliser les expériences et puis on va étudier les données à partir d'un modèle de régression Donc le modèle de régression qui s'écrit avec la constante plus la somme des termes linéaires plus la somme des termes quadratiques plus la somme des produits 2 à 2 des effets plus évidemment un résidu. Donc en fait ici il y a des termes somme mais c'est vraiment uniquement les effets linéaires, les effets quadratiques et les effets des interactions seuls. La résiduelle peut également se décomposer en 2 termes. Pourquoi elle peut se décomposer en 2 termes? Parce qu'en fait on va pouvoir décomposer la résiduelle quand on répète exactement une même expérience. Pourquoi? Parce qu'ici on a des points au centre. Donc les points au centre mesurent la répétabilité de l'expérimentation. Et c'est ce qu'on va appeler l'erreur pure. Donc là, on réalise exactement la même expérience et on voit quelle est la variabilité des données pour ces expériences. Donc cette variabilité, c'est l'erreur pure. L'erreur pure, j'ai juste à calculer la variance des Y pour tous les points au centre. Donc cette variance, elle a N0- 1 degrés de liberté si j'ai N0 points au centre. Et puis dans la résiduelle du modèle, j'ai aussi une erreur d'ajustement, le fait que mon modèle ne s'ajuste pas parfaitement au plan. Cette erreur d'ajustement, en fait, on va la calculer par différence. On va dire que c'est l'erreur résiduelle du modèle qu'on calcule quand on estime les coefficients du modèle, moins l'erreur pure qu'on vient de calculer. Et donc l'erreur d'ajustement c'est l'erreur résiduelle moins l'erreur pure et les degrés de liberté associés à cette erreur d'ajustement, eh bien c'est les degrés de liberté de la résiduelle moins les degrés de liberté de l'erreur pure. Grâce à cette décomposition de l'erreur résiduelle, on va pouvoir savoir si l'erreur d'ajustement est du même ordre de grandeur ou est beaucoup plus importante que l'erreur pure. Si elle est du même ordre de grandeur que l'erreur pure, ça va vouloir dire que notre modèle s'ajuste bien. Il fait la même erreur quasiment que quand on répète exactement 2 fois la même expérience. Et au contraire, si l'erreur d'ajustement est beaucoup plus grande que l'erreur pure, ça voudra dire que notre modèle n'est pas très bon, il s'ajuste pas très bien aux données. Et donc ensuite, eh bien, ce qu'on peut faire, c'est construire des tests, tester est-ce qu'il y a des effets linéaires, est-ce qu'il y a des effets quadratiques, est-ce qu'il y a des effets de l'interaction. Donc c'est des tests très globaux avec, bien globalement, est-ce qu'il y a l'effet de l'ensemble des interactions? Est-ce qu'il y a au moins une interaction qui est significative ou pas? Et puis éventuellement, on peut aussi tester est-ce que j'ai cette interaction qui est significative ou pas? Donc on peut s'intéresser soit globalement à l'ensemble des interactions, ou bien globalement à l'ensemble des effets quadratiques, ou encore globalement à l'ensemble des effets linéaires, ou bien on peut s'intéresser aussi uniquement à un effet linéaire ou un effet quadratique. Ça, ce sont les tests sur les effets du modèle. On a ensuite le test d'ajustement du modèle, c'est ce que j'ai dit tout à l'heure, on va tester si le modèle est bien ajusté ou au contraire s'il est mal ajusté. Et donc on va comparer la variabilité liée à l'ajustement par rapport à la variabilité de l'erreur pure. Si l'erreur d'ajustement est du même ordre de grandeur que l'erreur pure, alors on va dire que le modèle s'ajuste bien. Si au contraire la variabilité d'ajustement est beaucoup plus grande que la variabilité de l'erreur pure, c'est que le modèle n'est pas très bon et il faut changer de modèle. Donc une erreur d'ajustement significative va inciter à par exemple rajouter des effets quadratiques si on ne les avait pas mis dans le modèle ou rajouter des interactions si on ne les avait pas mises dans le modèle. Alors voici les lignes de code R qui permettent de construire le modèle. Donc tout d'abord ici j'ai construit le plan, ensuite je donne le vecteur de valeurs Y, donc là ça veut dire que j'ai réalisé mes expériences Et donc je récupère les valeurs de Y et puis je vais construire le modèle. Donc pour construire le modèle, on va utiliser la fonction RSM, RSM qui est une fonction qui permet d'étudier les surfaces de réponse, donc des modèles avec des effets quadratiques directement. Donc là on dit qu'on explique la variable Y en fonction de SO, ça signifie second order, de X1, X2. Ça veut dire qu'on va construire un modèle de second ordre, c'est-à-dire en mettant X1, X2 X1 carré, X2 carré et l'interaction X1 X2. Si j'avais mis FO de X1 X2, ça voulait dire que je prenais les first order, c'est-à-dire les effets linéaires de X1 X2. TWI, ça donne l'interaction entre X1 X2 et PQ, c'est les effets quadratiques de X1 et de X2. Et donc en mettant directement second order, vous avez l'effet global des effets linéaires. Donc vous avez FO de X1, X2, ça signifie que je suis en train d'étudier l'ensemble des effets linéaires. Et donc vous avez ici une probabilité critique de 10⁻⁶ qui signifie que les effets linéaires sont significativement différents de 0. Il y a au moins un des facteurs qui a un effet linéaire sur Y. Ensuite vous avez TWI qui est l'interaction, donc il y a un degré de liberté associé Et vous avez une probabilité critique qui est inférieure à 5%, ça veut dire qu'il y a un effet de l'interaction X1 X2 et que cet effet est significatif. Et puis ensuite PQ de X1 X2, c'est les effets quadratiques, donc c'est, il y a 2 effets quadratiques, donc c'est 2 degrés de liberté ici. Et quand on le teste, eh bien la probabilité critique est inférieure à 5%, et là encore ça signifie qu'il y a au moins un effet quadratique qui a un effet significatif sur Y. Et puis, puisque j'ai utilisé cette fonction RSM, vous avez aussi un test sur l'erreur d'ajustement, la covfit, l'erreur d'ajustement qui ici a une probabilité critique supérieure à 0,05. Ça veut dire que l'erreur d'ajustement n'est pas significative et donc elle est du même ordre de grandeur que l'erreur pure. Elle est testée ici par rapport à l'erreur pure, on calcule le carré moyen de l'erreur d'ajustement sur le carré moyen de l'erreur pure. Et ici, la probabilité critique étant supérieure à 5%, on accepte l'hypothèse qu'il n'y a pas d'erreur d'ajustement. Et donc on a un modèle qui est plutôt bien choisi. Alors ensuite, puisque j'ai vu qu'il y avait des effets linéaires, des effets quadratiques et des effets de l'interaction, eh bien on va pouvoir tester les effets linéaires, quadratiques, un par un. Donc là on voit que pour chaque test, les coefficients sont significativement différents de 0 et donc X1, X2 ont un effet sur Y, X1², X2² aussi, et l'interaction entre X1 et X2 aussi. Une fois que j'ai construit ce modèle-là, eh bien ce qu'on peut faire c'est rechercher l'optimum. Alors là le modèle s'écrit 4,62+ 2,23X1 1,10X2- 1,5X1X2+ 0,5X1²+ 0,75X2². Et donc je peux dériver par rapport à x1 et dériver par rapport à x2. Donc là j'ai dérivé par rapport à x1 et par rapport à x2. Je me retrouve avec 2 équations, 2 inconnues, donc je vais trouver la solution. Il suffit que je résolve ce système d'équations et en résolvant ce système d'équations, eh bien je vais trouver un optimum qui annule les dérivées par rapport à x1 et par rapport à x2. Et donc quand j'ai trouvé l'optimum, eh bien c'est ce qu'on appelle un point stationnaire et dans R, le point stationnaire est donné, c'est 6,68 pour x1 et 5,94 pour x2. Alors comment savoir maintenant si c'est un maximum ou un minimum ou éventuellement un point selle? D'accord que sur une selle de cheval, on va avoir ce qu'on appelle un point selle. Là où les dérivées s'annulent, c'est ni un minimum ni un maximum. Pour savoir si c'est un maximum, il faut que les valeurs propres qui sont données dans R soient toutes négatives. Donc si les valeurs propres données dans R sont toutes négatives, on va avoir un point stationnaire qui est un maximum. Si toutes les valeurs propres sont positives, le point stationnaire sera un minimum. Et si les valeurs propres sont de signes différents, alors on aura un point selle. Donc ça, c'est pour vous dire comment on va retrouver avec R si on a un maximum, un minimum ou un point selle. L'autre solution, évidemment, c'est de construire des graphes, de construire la surface de réponse et de voir si on est en train d'avoir un maximum ou un minimum. Sauf qu'on peut construire des graphes si on a 2 variables explicatives. Si on en a plus que 2, ça va être impossible de construire des graphes. Pour construire les graphes, il y a 2 fonctions: la fonction contour ou la fonction persp. Donc on va appliquer la fonction contour sur les résultats de notre modélisation. Donc les résultats de la modélisation ici c'est cr.rsm, je l'avais appelé comme ça, et je dis que je veux dessiner en fonction de X1 et X2. Donc ici je donne en fait le ~X1+ X2, c'est juste pour me dire qu'est-ce que je mets en abscisse, qu'est-ce que je mets en ordonnée. Donc là ça veut dire que j'ai X1 en abscisse et X2 en ordonnée. Et puis pour Persp, on a une surface en 3 dimensions avec toujours X1 et X2 et puis on donne des niveaux de couleurs. Alors problème, c'est que si on a plus que 2 variables, on ne va pas pouvoir visualiser facilement, on va être obligé de faire les graphes en considérant 2 variables uniquement à la fois. Donc par exemple construire X1 X2 puis X1 X3 puis X2 X3. Le problème est que si je veux construire X1, X2, il faudra que je fixe X3. Et à quelle valeur est-ce que je fixe X3? Donc en général, ce qu'on fait, c'est qu'on va fixer X3 à la valeur de l'optimum pour voir comment jouent X1, X2. Parce qu'évidemment, en bougeant X3, l'optimum sur X1 et X2 sera pas le même. Donc le mieux, c'est de fixer X3 à l'optimum et de tracer la surface de réponse ou l'image avec les lignes de niveau avec X1, X2. Puis après on fixe X1 à l'optimum et on trace X2, X3. En pratique, si les essais coûtent très cher, on peut construire le plan de façon séquentielle. Donc on peut commencer par construire le plan factoriel, c'est-à-dire soit un plan complet soit un plan fractionnaire, et rajouter des points au centre. Et à partir de ces essais, on va construire un modèle linéaire sans effet quadratique Et grâce au point au centre, on va pouvoir tester la qualité d'ajustement du modèle en comparant l'erreur d'ajustement du modèle avec l'erreur pure. Si l'erreur d'ajustement du modèle est faible, si elle est du même ordre de grandeur que l'erreur pure, alors le modèle linéaire convient et on va s'arrêter là et considérer qu'on n'a pas besoin des effets quadratiques. Si maintenant l'erreur d'ajustement est forte, ça voudra dire que le modèle ne convient pas et on va rajouter les points en étoile pour pouvoir construire par la suite un modèle plus complet avec les effets linéaires mais aussi les effets quadratiques et les effets d'interaction. L'intérêt est de pouvoir construire peu d'essais dans un premier temps et vérifier que le modèle convient ou ne convient pas. Et si le modèle convient pas, on fait plus d'essais et on construit un modèle plus complet. Alors pour terminer, je vais proposer un dernier type de plan autre que les plans composites centrés. Ce sont des plans qu'on appelle plans de Box-Benken. Qui permettent aussi de s'intéresser à des variables quantitatives et qui permettent éventuellement d'avoir un petit peu moins d'essais et des plans qui sont de qualité assez bonne. Donc on fera des TD pour comparer la qualité des plans composite centré et des plans de Box-Benken. Alors comment on construit ces plans? Eh bien, ce qu'on fait, c'est qu'on va construire le plan complet pour un couple de 2 facteurs. Donc par exemple, en bleu ici, j'ai construit le plan complet pour le premier facteur et le deuxième facteur. Donc j'ai un plan en 4 essais et ce que je fais c'est que je rajoute que des zéros pour le troisième facteur. Ensuite je refais la même chose, je construis le plan complet pour la variable 1 et la variable 3 et donc là ce sont les points rouges et je mets x2 à la valeur 0. Et puis ensuite je construis à nouveau le plan complet pour les facteurs 2 et 3 et et je rajoute des points au centre pour la première variable. Donc je mets la première variable à 0. Et puis je rajoute des points au centre du domaine, donc avec pour X1, X2 et X3 que des 0. Alors quel est l'intérêt de ce plan d'expérience? Eh bien c'est un plan d'expérience où en fait on voit qu'on va tester pour chaque facteur uniquement 3 niveaux: le niveau 1, le niveau -1 et le niveau 0. Ce qui est assez intéressant parce qu'on va avoir à tester que ces 3 niveaux alors qu'avec le plan composite centré, on va avoir besoin de tester 5 niveaux puisqu'on a les points en étoile en plus des points sur le carré et du point au centre. Donc le fait d'avoir moins à changer de niveau permet d'un point de vue pratique d'avoir des plans qui sont plus faciles en général à réaliser. La qualité du plan est un petit peu moins bonne mais par rapport aux contraintes pratiques, dans certaines situations, ça peut être intéressant. Et l'autre intérêt de ce plan de Box-Benken, c'est que je peux éventuellement rajouter un facteur si je veux. Donc là par exemple, j'ai étudié 3 facteurs, je veux rajouter un 4e facteur, eh bien il suffit que je mette une 4e colonne avec des 0 pour les points en bleu, des 0 pour les points rouges, des 0 pour les points verts et puis après je vais rajouter des essais avec X1 X4, X2 X4, X3 X4 et les autres facteurs à 0. Et donc je peux construire mon plan de façon séquentielle. Alors qu'avec le plan composite centré, je ne peux pas rajouter un facteur une fois que j'ai construit un plan composite centré avec 3 facteurs par exemple. Donc sur R, pour construire ces plans de Box-Benken, c'est toujours le package RSM et c'est la fonction BBD, Box-Benken Design, et il suffit de préciser le nombre de facteurs qu'on veut utiliser.