1 Génération d’un Plan Factoriel Complet (PFC)

Au cours théorique vous avez appris ce qu’était un Plan Factoriel Complet. Faites tourner l’exemple ci-dessous et répondez aux questions qui suivent.

Design=FrF2(nruns = 16,nfactors = 4,replications =2,randomize=FALSE,repeat.only=TRUE)

creating full factorial with 16 runs ...

design.info(Design)

$type
[1] "full factorial"

$nruns
[1] 16

$nfactors
[1] 4

$nlevels
[1] 2 2 2 2

$factor.names
$factor.names$A
[1] "-1" "1" 

$factor.names$B
[1] "-1" "1" 

$factor.names$C
[1] "-1" "1" 

$factor.names$D
[1] "-1" "1" 


$replications
[1] 2

$repeat.only
[1] TRUE

$randomize
[1] FALSE

$seed
NULL

$creator
FrF2(nruns = 16, nfactors = 4, replications = 2, randomize = FALSE, 
    repeat.only = TRUE)

$quantitative
    A     B     C     D 
FALSE FALSE FALSE FALSE 

$FrF2.version
[1] "2.2-3"

DM=as.data.frame(desnum(Design))
DM  # DM = Design matrix

   A1 B1 C1 D1
1  -1 -1 -1 -1
2  -1 -1 -1 -1
3   1 -1 -1 -1
4   1 -1 -1 -1
5  -1  1 -1 -1
6  -1  1 -1 -1
7   1  1 -1 -1
8   1  1 -1 -1
9  -1 -1  1 -1
10 -1 -1  1 -1
11  1 -1  1 -1
12  1 -1  1 -1
13 -1  1  1 -1
14 -1  1  1 -1
15  1  1  1 -1
16  1  1  1 -1
17 -1 -1 -1  1
18 -1 -1 -1  1
19  1 -1 -1  1
20  1 -1 -1  1
21 -1  1 -1  1
22 -1  1 -1  1
23  1  1 -1  1
24  1  1 -1  1
25 -1 -1  1  1
26 -1 -1  1  1
27  1 -1  1  1
28  1 -1  1  1
29 -1  1  1  1
30 -1  1  1  1
31  1  1  1  1
32  1  1  1  1

1.1 Répondez aux questions suivantes concernant ce plan

Justifiez pourquoi il s’agit d’un Plan Factoriel Complet.
Combien de facteurs sont testés pour combien d’essais au total ?
Combien de niveaux a chacun des facteurs ?
Que représentent les valeurs -1 et 1 pour chacun des facteurs ?

Sur base du Design nous observons qu’il y a 4 colonnes et donc 4 facteurs ayant chacun 2 niveaux (-1 et 1). Pour un plan factoriel complet à 4 facteurs, il faut \(2^4 = 16\) essais pour avoir toutes les combinaisons de chacun des niveaux des facteurs or nous en avons 32. Cependant, notre Design est répliqué puique chacune des lignes est répétée une deuxième fois. Nous avons donc bien un Plan Factoriel Complet \(2^4\), avec 2 réplications.

Les valeurs -1 et 1 représentent les valeurs minimales et maximales du domaine de chacun des facteurs testés. Le domaine de tous les facteurs est donc compris entre -1 et 1, on parle de données standardisées (NB: Nous réfléchirons plus tard dans le TP à l’avantage de standardisé le domaine de chaque facteur).

1.2 Etude de l’orthogonalité du plan

La définition de l’orthogonalité n’est pas toujours claire dans la litterature. Nous allons dans ce cours toujours lier un plan à un modèle d’intérêt qu’on désire estimer. On dira alors qu’un plan est orthogonal pour ce modèle s’il permet d’estimer les paramètres de façon non corrélée.

Répondez sur cette base aux questions suivantes.

Quel est le plus grand modèle polymonial qui peut être estimé avec ce plan ?
On va se limiter ici au modèle comprenant tous les effets principaux et toutes les interactions 2 à 2 entre facteurs.
Comment vérifier l’orthogonalité du plan par rapport à ce dernier modèle ?
Etudiez le code ci-dessous et validez que le plan est bien orthogonal.

X=model.matrix(~A+B+C+D+A:B+A:C+A:D+B:C+B:D+C:D,Design)
# Matrice du modèle
X

   (Intercept) A1 B1 C1 D1 A1:B1 A1:C1 A1:D1 B1:C1 B1:D1 C1:D1
1            1 -1 -1 -1 -1     1     1     1     1     1     1
2            1 -1 -1 -1 -1     1     1     1     1     1     1
3            1  1 -1 -1 -1    -1    -1    -1     1     1     1
4            1  1 -1 -1 -1    -1    -1    -1     1     1     1
5            1 -1  1 -1 -1    -1     1     1    -1    -1     1
6            1 -1  1 -1 -1    -1     1     1    -1    -1     1
7            1  1  1 -1 -1     1    -1    -1    -1    -1     1
8            1  1  1 -1 -1     1    -1    -1    -1    -1     1
9            1 -1 -1  1 -1     1    -1     1    -1     1    -1
10           1 -1 -1  1 -1     1    -1     1    -1     1    -1
11           1  1 -1  1 -1    -1     1    -1    -1     1    -1
12           1  1 -1  1 -1    -1     1    -1    -1     1    -1
13           1 -1  1  1 -1    -1    -1     1     1    -1    -1
14           1 -1  1  1 -1    -1    -1     1     1    -1    -1
15           1  1  1  1 -1     1     1    -1     1    -1    -1
16           1  1  1  1 -1     1     1    -1     1    -1    -1
17           1 -1 -1 -1  1     1     1    -1     1    -1    -1
18           1 -1 -1 -1  1     1     1    -1     1    -1    -1
19           1  1 -1 -1  1    -1    -1     1     1    -1    -1
20           1  1 -1 -1  1    -1    -1     1     1    -1    -1
21           1 -1  1 -1  1    -1     1    -1    -1     1    -1
22           1 -1  1 -1  1    -1     1    -1    -1     1    -1
23           1  1  1 -1  1     1    -1     1    -1     1    -1
24           1  1  1 -1  1     1    -1     1    -1     1    -1
25           1 -1 -1  1  1     1    -1    -1    -1    -1     1
26           1 -1 -1  1  1     1    -1    -1    -1    -1     1
27           1  1 -1  1  1    -1     1     1    -1    -1     1
28           1  1 -1  1  1    -1     1     1    -1    -1     1
29           1 -1  1  1  1    -1    -1    -1     1     1     1
30           1 -1  1  1  1    -1    -1    -1     1     1     1
31           1  1  1  1  1     1     1     1     1     1     1
32           1  1  1  1  1     1     1     1     1     1     1
attr(,"assign")
 [1]  0  1  2  3  4  5  6  7  8  9 10
attr(,"contrasts")
attr(,"contrasts")$A
   [,1]
-1   -1
1     1

attr(,"contrasts")$B
   [,1]
-1   -1
1     1

attr(,"contrasts")$C
   [,1]
-1   -1
1     1

attr(,"contrasts")$D
   [,1]
-1   -1
1     1

# Matrice de X'X
t(X)%*%X

            (Intercept) A1 B1 C1 D1 A1:B1 A1:C1 A1:D1 B1:C1 B1:D1 C1:D1
(Intercept)          32  0  0  0  0     0     0     0     0     0     0
A1                    0 32  0  0  0     0     0     0     0     0     0
B1                    0  0 32  0  0     0     0     0     0     0     0
C1                    0  0  0 32  0     0     0     0     0     0     0
D1                    0  0  0  0 32     0     0     0     0     0     0
A1:B1                 0  0  0  0  0    32     0     0     0     0     0
A1:C1                 0  0  0  0  0     0    32     0     0     0     0
A1:D1                 0  0  0  0  0     0     0    32     0     0     0
B1:C1                 0  0  0  0  0     0     0     0    32     0     0
B1:D1                 0  0  0  0  0     0     0     0     0    32     0
C1:D1                 0  0  0  0  0     0     0     0     0     0    32

Le plus grand modèle estimable avec ce plan est un modèle à 16 termes : 1 (intercept) + 4 (effets principaux) + 6 (interactions des effets 2 à 2) + 4 (interactions des effets 3 à 3) + 1 (interaction avec tous les effets) = 16 termes.
Oui, car notre matrice X’X est diagonale. Cela implique que chaque facteur est estimé de manière non-corrélée (et indépendante en cas de normalité) des autres. Cela est mieux pour l’efficacité de nos estimations et une qualité appréciée pour les plans.

1.3 Présentation de l’exemple

Les vaccins comprenant un conjugué antigène de polysaccharide (PS) lié à une protéine ont de multiples applications. Dans le cadre de cette étude, on s’intéresse au Rendement en terme de PS couplé (%). Celui-ci dépend de :
C : Concentration initiale de PS [2, 4] mg/l
P : PH d’activation [7, 9]
R : Ratio PS/PR [0.8, 1.2]
T : Type de protéine (identiques mais 2 origines différents) (T1,T2)

1.4 Notion de valeurs standardisées

Le jeux de données avec lequel vous allez travailler est en valeurs standardisées. Travailler avec des données standardisées est intéressant lorsque nous travaillons avec des plans expérimentaux pour différentes raisons.

Cela permet de remettre le domaine des valeurs de tous les facteurs sur une même échelle entre -1 et 1 et donc d’interpréter les effets des facteurs directement par rapport au domaine d’intérêt sans trop s’occuper des unités de chacun.
Quand un facteur standardisé augmente d’une unité, cela veut dire que le facteur passe du milieu de son domaine à son maximum et la valeur des coefficients s’interprête alors en conséquence.
Nous pouvons alors directement comparer les coefficients des différents facteurs pour voir lequel à le plus d’impact sur la réponse.
D’autre part, cela modifie également l’interprétation du terme constant qui représente la valeur moyenne de Y quand tous les facteurs sont au milieu de leur domaine (en 0).
D’un point de vue statistique, en particulier dans les plans factoriels, la standardisation (et surtout le centrage) permet d’obtenir l’orthogonalité de la matrice du modèle et estimer les effets des facteurs indépendamment.

2 Analyse par régression des résultats d’un plan factoriel complet

2.1 Plan Factoriel Complet (PFC)

2.1.1 Lecture et visualisation des données

Téléchargez les données du fichier excel TP7_PFC.xlsx puis : * Visualiser le fichier de données et vérifier que le plan choisi est bien un plan factoriel complet. Est-il répété ?

vaccinPFC <- read_excel("TP7_PFC.xlsx") 
pander(vaccinPFC)

Cs	Ps	Rs	Ts	Rendement
-1	-1	-1	-1	59.6
-1	-1	-1	-1	58.6
1	-1	-1	-1	84.7
1	-1	-1	-1	88.2
-1	1	-1	-1	85.2
-1	1	-1	-1	84.6
1	1	-1	-1	92.7
1	1	-1	-1	93.3
-1	-1	1	-1	45.6
-1	-1	1	-1	43.4
1	-1	1	-1	74.1
1	-1	1	-1	71.7
-1	1	1	-1	70.8
-1	1	1	-1	70.1
1	1	1	-1	79.6
1	1	1	-1	80.4
-1	-1	-1	1	51.4
-1	-1	-1	1	51.3
1	-1	-1	1	78.1
1	-1	-1	1	77.9
-1	1	-1	1	77.7
-1	1	-1	1	76.6
1	1	-1	1	86.8
1	1	-1	1	84.8
-1	-1	1	1	37.5
-1	-1	1	1	35.2
1	-1	1	1	67.4
1	-1	1	1	66.3
-1	1	1	1	60.8
-1	1	1	1	61.5
1	1	1	1	72.8
1	1	1	1	72.4

2.1.2 Estimation du modèle complet de degré 2

Dites quel est le modèle le plus complet qu’il est possible d’estimer avec ce plan.
Estimez un modèle de régression lm d’ordre 2 (comprenant tous les effets principaux et toutes les interactions à 2) pour expliquer le Rendement en fonction des 4 facteurs.
Sortez le tableau résumé de la régression et une analyse des résidus.
Que pensez-vous de la qualité du modèle selon des différents critères que vous connaissez.
Quels termes du modèle ont un effet significatif sur la réponse ?
Sortez une table d’ANOVA et comparez là au tableau des coefficients. Est ce que cette table apporte quelque chose ?
Vérifier s’il y a des alias entre les termes du modèle. Pour cela, vous pouvez utiliser la fonction aliases ? Comment auriez-vous pu trouver cette réponse sans utiliser la fonction ?
Regardez aussi les écart-types des paramètres du modèle en calculant la matrice variance-covariance \(\sqrt{(X'X)^{-1}*S^2}\). Ils sont tous identiques. Pouvez-vous expliquer pourquoi ?

modPFC <- lm(Rendement ~ Cs+Ps+Rs+Ts+Cs:Ps+Cs:Rs+Cs:Ts+Ps:Rs+Ps:Ts+Rs:Ts, data = vaccinPFC)
summary(modPFC)


Call:
lm.default(formula = Rendement ~ Cs + Ps + Rs + Ts + Cs:Ps + 
    Cs:Rs + Cs:Ts + Ps:Rs + Ps:Ts + Rs:Ts, data = vaccinPFC)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.67813 -0.65000 -0.03437  0.45000  2.35938 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 70.03437    0.18668 375.155  < 2e-16 ***
Cs           9.41563    0.18668  50.437  < 2e-16 ***
Ps           8.09688    0.18668  43.373  < 2e-16 ***
Rs          -6.93437    0.18668 -37.146  < 2e-16 ***
Ts          -3.87812    0.18668 -20.774 1.75e-15 ***
Cs:Ps       -4.69688    0.18668 -25.160  < 2e-16 ***
Cs:Rs        0.57187    0.18668   3.063   0.0059 ** 
Cs:Ts        0.24062    0.18668   1.289   0.2114    
Ps:Rs       -0.14688    0.18668  -0.787   0.4402    
Ps:Ts       -0.07813    0.18668  -0.418   0.6798    
Rs:Ts        0.01562    0.18668   0.084   0.9341    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.056 on 21 degrees of freedom
Multiple R-squared:  0.997, Adjusted R-squared:  0.9955 
F-statistic: 688.1 on 10 and 21 DF,  p-value: < 2.2e-16

par(mfrow=c(2,2))
plot(modPFC)

les valeurs estimées ("leverages") sont tous = 0.34375
 et il n'y a aucune variable dépendante facteur ; aucun graphe no. 5

par(mfrow=c(1,1))
anova(modPFC)

Analysis of Variance Table

Response: Rendement
          Df  Sum Sq Mean Sq   F value    Pr(>F)    
Cs         1 2836.93 2836.93 2543.8885 < 2.2e-16 ***
Ps         1 2097.90 2097.90 1881.1986 < 2.2e-16 ***
Rs         1 1538.74 1538.74 1379.7945 < 2.2e-16 ***
Ts         1  481.28  481.28  431.5622 1.753e-15 ***
Cs:Ps      1  705.94  705.94  633.0205 < 2.2e-16 ***
Cs:Rs      1   10.47   10.47    9.3843    0.0059 ** 
Cs:Ts      1    1.85    1.85    1.6614    0.2114    
Ps:Rs      1    0.69    0.69    0.6190    0.4402    
Ps:Ts      1    0.20    0.20    0.1751    0.6798    
Rs:Ts      1    0.01    0.01    0.0070    0.9341    
Residuals 21   23.42    1.12                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

aliases(modPFC)

[1] no aliasing in the model

## Calcul des écarts types du modèle : 
Design=FrF2(nruns = 16,nfactors = 4,replications =2,randomize=FALSE,repeat.only=TRUE)

creating full factorial with 16 runs ...

X=model.matrix(modPFC,Design)

T1 <- solve(t(X)%*%X) #(X'X)^-1
T2 <- T1*1.1151 #(X'X)^-1 * S^2
T3 <- sqrt(T2) # sqrt((X'X)^-1 * S^2))
T3

            (Intercept)        Cs        Ps        Rs        Ts     Cs:Ps
(Intercept)   0.1866732 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
Cs            0.0000000 0.1866732 0.0000000 0.0000000 0.0000000 0.0000000
Ps            0.0000000 0.0000000 0.1866732 0.0000000 0.0000000 0.0000000
Rs            0.0000000 0.0000000 0.0000000 0.1866732 0.0000000 0.0000000
Ts            0.0000000 0.0000000 0.0000000 0.0000000 0.1866732 0.0000000
Cs:Ps         0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.1866732
Cs:Rs         0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
Cs:Ts         0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
Ps:Rs         0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
Ps:Ts         0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
Rs:Ts         0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
                Cs:Rs     Cs:Ts     Ps:Rs     Ps:Ts     Rs:Ts
(Intercept) 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
Cs          0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
Ps          0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
Rs          0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
Ts          0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
Cs:Ps       0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
Cs:Rs       0.1866732 0.0000000 0.0000000 0.0000000 0.0000000
Cs:Ts       0.0000000 0.1866732 0.0000000 0.0000000 0.0000000
Ps:Rs       0.0000000 0.0000000 0.1866732 0.0000000 0.0000000
Ps:Ts       0.0000000 0.0000000 0.0000000 0.1866732 0.0000000
Rs:Ts       0.0000000 0.0000000 0.0000000 0.0000000 0.1866732

Il est possible d’estimer un modèle d’ordre 1 avec toutes les interactions (d’ordre 2, 3 et 4).
Notre modèle semble intéressant pour décrire la variabilité de nos mesures (\(R^2 = 0.997\)). Les hypothèses sur les résidus sont respectées.
Tous les termes d’ordre 1 ont un effet significatif sur le Rendement. Par contre, uniquement les termes d’ordre 2 Cs:Ps et Cs:Rs sont significatifs.
L’ANOVA et le summary donnent des p-valeurs identiques. Effectivement, les deux tableaux testent si l’effet des différents termes/facteurs sur le rendement est significativement différent de zéro.
Aucun effet n’est confondu. C’est l’avantage des Plans Factoriels Complets. Nous avons suffisamment d’essais pour tester tous les niveaux de chaque facteur. Ainsi, nous avons toutes les combinaisons possibles qui nous permettent d’estimer tous les termes d’interactions.
Les écarts-types du modèles sont tous identiques puisqu’ils sont calculés à partir de \(S^2\) et de la matrice \((X'X)^{-1}\) qui sont les mêmes pour tous les paramètres puisque les valeurs sont standardisées.

2.1.3 Représentation graphique des résultats

Le package FrF2 propose deux manières, autres que l’anova, d’analyser l’effet des différents facteurs de votre modèle.

Le Daniel plot qui est un qq plot des coefficients du modèle. Si tous les paramètres théoriques valent 0, le graphique devrait ressembler à qqplot de données normales. Si pas, les coefficients plus importants ressortent comme des outliers. C’est une manière de visualiser les termes importants dans le modèle.
Le plot des interaction IAPlot permet de visualiser les interactions 2 à 2 entre facteurs.
La fonction MEPlot permet finalement de représenter les effets principaux des facteurs.

Questions :

Retrouvez vous bien les termes attendus sur le Daniel Plot ?
Retrouver vous bien les interactions attendues sur le graphe d’interaction ? Est ce que ces interactions sont fortes ou modérées.
Interprétez les effets principaux des facteurs et essayez de comprendre ce qu’ils représentent et quelle est leur relation avec les coefficients du modèle.

# DanielPlot
FrF2::DanielPlot(modPFC, half=TRUE, code = TRUE)

# InteractionPlot
FrF2::IAPlot(modPFC)

# MEPlot
FrF2::MEPlot(modPFC)

A l’aide de l’IAPlot, nous pouvons observer l’interaction Cs:Ps (significative à l’anova) car les droites ne sont pas parallèles. L’interaction Cs:Rs est plus difficile à distinguer (elle est également moins significative que l’autre dans l’anova). Le DanielPlot permet d’observer l’effet absolu des différents termes du modèle. On y retrouve les termes significatifs dans l’anova par ordre de significacité : Cs, Ps, Rs, Cs:Ps, Ts et Cs:Rs.

2.1.4 Simplification du modèle

Estimez un modèle comprenant uniquement les termes significatifs du modèle et comparez les tableaux résumés des deux modèles.

Les coefficients n’ont pas changé, pouvez-vous expliquer pourquoi ?
Essayez d’expliquer les autres changements qui sont assez minimes.

modPFC2 <- lm(Rendement ~ Cs+Ps+Rs+Ts+Cs:Ps+Cs:Rs, data = vaccinPFC)
summary(modPFC)


Call:
lm.default(formula = Rendement ~ Cs + Ps + Rs + Ts + Cs:Ps + 
    Cs:Rs + Cs:Ts + Ps:Rs + Ps:Ts + Rs:Ts, data = vaccinPFC)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.67813 -0.65000 -0.03437  0.45000  2.35938 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 70.03437    0.18668 375.155  < 2e-16 ***
Cs           9.41563    0.18668  50.437  < 2e-16 ***
Ps           8.09688    0.18668  43.373  < 2e-16 ***
Rs          -6.93437    0.18668 -37.146  < 2e-16 ***
Ts          -3.87812    0.18668 -20.774 1.75e-15 ***
Cs:Ps       -4.69688    0.18668 -25.160  < 2e-16 ***
Cs:Rs        0.57187    0.18668   3.063   0.0059 ** 
Cs:Ts        0.24062    0.18668   1.289   0.2114    
Ps:Rs       -0.14688    0.18668  -0.787   0.4402    
Ps:Ts       -0.07813    0.18668  -0.418   0.6798    
Rs:Ts        0.01562    0.18668   0.084   0.9341    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.056 on 21 degrees of freedom
Multiple R-squared:  0.997, Adjusted R-squared:  0.9955 
F-statistic: 688.1 on 10 and 21 DF,  p-value: < 2.2e-16

summary(modPFC2)


Call:
lm.default(formula = Rendement ~ Cs + Ps + Rs + Ts + Cs:Ps + 
    Cs:Rs, data = vaccinPFC)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.8656 -0.5531 -0.1031  0.5016  1.9094 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  70.0344     0.1808 387.252  < 2e-16 ***
Cs            9.4156     0.1808  52.063  < 2e-16 ***
Ps            8.0969     0.1808  44.771  < 2e-16 ***
Rs           -6.9344     0.1808 -38.343  < 2e-16 ***
Ts           -3.8781     0.1808 -21.444  < 2e-16 ***
Cs:Ps        -4.6969     0.1808 -25.971  < 2e-16 ***
Cs:Rs         0.5719     0.1808   3.162  0.00408 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.023 on 25 degrees of freedom
Multiple R-squared:  0.9966,    Adjusted R-squared:  0.9958 
F-statistic:  1222 on 6 and 25 DF,  p-value: < 2.2e-16

Les coefficients n’ont pas changés puisque les termes que nous avions enlevés n’étaient pas significatifs.
Le \(R^2\) a légèrement diminué du à la suppression des facteurs non significatifs, mais qui expliquaient quand même une partie de la variabilité. Par contre le \(R^2_{ajusté}\) a légèrement augmenté, ce qui montre bien que supprimer les termes du modèle était pertinent puisque ceux-ci n’expliquaient pas beaucoup de la variabilité du modèle.

#Check orthogonality
X <- model.matrix(modPFC)
t(X)%*%X

            (Intercept) Cs Ps Rs Ts Cs:Ps Cs:Rs Cs:Ts Ps:Rs Ps:Ts Rs:Ts
(Intercept)          32  0  0  0  0     0     0     0     0     0     0
Cs                    0 32  0  0  0     0     0     0     0     0     0
Ps                    0  0 32  0  0     0     0     0     0     0     0
Rs                    0  0  0 32  0     0     0     0     0     0     0
Ts                    0  0  0  0 32     0     0     0     0     0     0
Cs:Ps                 0  0  0  0  0    32     0     0     0     0     0
Cs:Rs                 0  0  0  0  0     0    32     0     0     0     0
Cs:Ts                 0  0  0  0  0     0     0    32     0     0     0
Ps:Rs                 0  0  0  0  0     0     0     0    32     0     0
Ps:Ts                 0  0  0  0  0     0     0     0     0    32     0
Rs:Ts                 0  0  0  0  0     0     0     0     0     0    32

Le plan est bien orthogonal.

3 Génération d’un Plan Factoriel Fractionaire (PFF)

4 Génération d’un Plan Factoriel Fractionnaire (PFF)

4.1 Qu’est-ce qu’un Plan Factoriel Fractionnaire ?

Un Plan Factoriel Fractionnaire (\(2^{k-r}\)) est un plan résultant de r fractionnements successifs d’un plan factoriel complet \(2^k\) à partir de r générateurs. Ce type de plan permet de réduire le nombre d’essais de manière importante tout en gardant un plan de qualité. Ce plan ne permettra par contre pas d’estimer des modèles d’ordre aussi élevé que les plans factoriels complets. Certains termes du modèle complet seront confondus entre eux. Le choix des générateurs est donc directement lié à la complexité du modèle qui sera estimable.

4.2 Plan factoriel fractionnaire - Version 1

Le code ci-dessous permet de générer un plan factoriel frationnaire. Testez le puis répondez aux questions qui suivent.

Design=FrF2(8,4,randomize=FALSE)
design.info(Design)

$type
[1] "FrF2"

$nruns
[1] 8

$nfactors
[1] 4

$factor.names
$factor.names$A
[1] -1  1

$factor.names$B
[1] -1  1

$factor.names$C
[1] -1  1

$factor.names$D
[1] -1  1


$catlg.name
[1] "catlg"

$catlg.entry
Design:  4-1.1 
   8  runs,  4  factors,  
   Resolution  IV 
   Generating columns:  7 
   WLP (3plus):  0 1 0 0 0 ,  0  clear 2fis

$aliased
$aliased$legend
[1] "A=A" "B=B" "C=C" "D=D"

$aliased$main
character(0)

$aliased$fi2
[1] "AB=CD" "AC=BD" "AD=BC"


$FrF2.version
[1] "2.2-3"

$replications
[1] 1

$repeat.only
[1] FALSE

$randomize
[1] FALSE

$seed
NULL

$creator
FrF2(8, 4, randomize = FALSE)

generators(Design)

$generators
[1] "D=ABC"

aliasprint(Design)

$legend
[1] A=A B=B C=C D=D

$main
character(0)

$fi2
[1] AB=CD AC=BD AD=BC

DM=as.data.frame(desnum(Design))
DM

   A  B  C  D
1 -1 -1 -1 -1
2  1 -1 -1  1
3 -1  1 -1  1
4  1  1 -1 -1
5 -1 -1  1  1
6  1 -1  1 -1
7 -1  1  1 -1
8  1  1  1  1

Combien d’essais comprend le plan ? combien de fois le plan complet a-t-il été fractionné ?
Sur base du script ci-dessous, retrouvez quels sont le/les générateur(s) utilisé(s) pour fractionner le PFC.
Vérifiez sur le plan généré que ce(s) générateur(s) sont bien vérifiés.
Quels effets pouvons-nous estimer avec ce plan et quels effets sont confondus ?

Il s’agit d’un Plan Factoriel Fractionaire avec 1 générateur et sans réplications :\(2^{4-1} = 8\). Le générateur utilisé est D=ABC. Ce plan permet d’estimer tous les effets d’ordre 1 sans qu’ils soient confondus avec d’autres effets d’ordre 1 ou 2. Par contre, les effets d’ordre 2 sont tous confondus entre eux. Les effets d’ordre 1 sont eux confondus avec des effets d’ordre 3 : D=ABC, A=BCD, B=ACD, D=ABC.

Exemple : Si nous souhaitons estimer A:B, l’estimation sera confondue avec celle de C:D. Cela n’est pas grave si l’effet de A:B est non nul et celui de C:D ne l’est pas. Par contre si les deux effets sont non nuls, l’estimation d’A:B sera faussée car elle comprendra également l’effet de C:D.

voici un code utile qui permet de trouver tous les alias :

alldata=data.frame(y=1,desnum(Design))
aliases(lm(y~(.)^3,data=alldata))

          
 A = B:C:D
 B = A:C:D
 C = A:B:D
 D = A:B:C
 A:B = C:D
 A:C = B:D
 A:D = B:C

#Representation des confusions entre les termes d'ordre 2 et ceux d'ordre 1
corrPlot(Design, scale = "corr")

kmax was reduced to number of factors

Le chargement a nécessité le package : RColorBrewer

4.2.1 Résolution du plan

Sur base de votre réponse à la question précédante, quelle est la résolution du plan ?

Il s’agit d’un plan de résolution IV car il permet d’estimer les effets d’ordre 1 et uniquement certains effets d’ordre 2 car ils sont confondus entre eux. La longueur du générateur indique aussi la résolution. On peut en effet réécrire le générateur comme suite ABCD=1.

4.2.2 Etude de l’orthogonalité du plan

Ce plan est-il orthogonal pour l’estimation du modèle d’ordre 1 ?
Et pour le modèle d’ordre 2 complet ?

# Check des matrices X'X
# Modèle d'ordre 1
X1=model.matrix(~A+B+C+D,DM)
t(X1)%*%X1

            (Intercept) A B C D
(Intercept)           8 0 0 0 0
A                     0 8 0 0 0
B                     0 0 8 0 0
C                     0 0 0 8 0
D                     0 0 0 0 8

# Modèle d'ordre 2
X2=model.matrix(~A+B+C+D+A:B+A:C+A:D+B:C+B:D+C:D,DM)
X2

  (Intercept)  A  B  C  D A:B A:C A:D B:C B:D C:D
1           1 -1 -1 -1 -1   1   1   1   1   1   1
2           1  1 -1 -1  1  -1  -1   1   1  -1  -1
3           1 -1  1 -1  1  -1   1  -1  -1   1  -1
4           1  1  1 -1 -1   1  -1  -1  -1  -1   1
5           1 -1 -1  1  1   1  -1  -1  -1  -1   1
6           1  1 -1  1 -1  -1   1  -1  -1   1  -1
7           1 -1  1  1 -1  -1  -1   1   1  -1  -1
8           1  1  1  1  1   1   1   1   1   1   1
attr(,"assign")
 [1]  0  1  2  3  4  5  6  7  8  9 10

t(X2)%*%X2

            (Intercept) A B C D A:B A:C A:D B:C B:D C:D
(Intercept)           8 0 0 0 0   0   0   0   0   0   0
A                     0 8 0 0 0   0   0   0   0   0   0
B                     0 0 8 0 0   0   0   0   0   0   0
C                     0 0 0 8 0   0   0   0   0   0   0
D                     0 0 0 0 8   0   0   0   0   0   0
A:B                   0 0 0 0 0   8   0   0   0   0   8
A:C                   0 0 0 0 0   0   8   0   0   8   0
A:D                   0 0 0 0 0   0   0   8   8   0   0
B:C                   0 0 0 0 0   0   0   8   8   0   0
B:D                   0 0 0 0 0   0   8   0   0   8   0
C:D                   0 0 0 0 0   8   0   0   0   0   8

Ce plan est bien orthogonal pour l’estimation d’un modèle d’ordre 1 car la matrice X1’X1 est diagonale. Ce n’est par contre pas le cas pour un modèle d’ordre 2. On voit que des effets sont confondus entre eux (les colonnes correspondantes sont identiques dans la matrice du modèle) et la matrice X2’X2 ne sera pas inversible et donc les paramètres pas estimables.

4.3 Plan factoriel fractionnaire - Version 2

Etudiez et exécutez le code ci-dessous, puis répondez aux questions qui suivent.

Design=FrF2(8,4,generators = c("AB")) #specification du generateur a utiliser
DM=as.data.frame(desnum(Design))
DM

   A  B  C  D
1  1  1  1  1
2  1  1 -1  1
3 -1 -1  1  1
4 -1  1 -1 -1
5 -1  1  1 -1
6 -1 -1 -1  1
7  1 -1 -1 -1
8  1 -1  1 -1

#reponse: 
design.info(Design)

$type
[1] "FrF2.generators"

$nruns
[1] 8

$nfactors
[1] 4

$factor.names
$factor.names$A
[1] -1  1

$factor.names$B
[1] -1  1

$factor.names$C
[1] -1  1

$factor.names$D
[1] -1  1


$generators
[1] "D=AB"

$aliased
$aliased$legend
[1] "A=A" "B=B" "C=C" "D=D"

$aliased$main
[1] "A=BD" "B=AD" "D=AB"

$aliased$fi2
character(0)


$FrF2.version
[1] "2.2-3"

$replications
[1] 1

$repeat.only
[1] FALSE

$randomize
[1] TRUE

$seed
NULL

$creator
FrF2(8, 4, generators = c("AB"))

aliasprint(Design)

$legend
[1] A=A B=B C=C D=D

$main
[1] A=BD B=AD D=AB

$fi2
character(0)

generators(Design)

$generators
[1] "D=AB"

X=model.matrix(~A+B+C+D+A:B+A:C+A:D+B:C+B:D+C:D,Design)
t(X)%*%X

            (Intercept) A1 B1 C1 D1 A1:B1 A1:C1 A1:D1 B1:C1 B1:D1 C1:D1
(Intercept)           8  0  0  0  0     0     0     0     0     0     0
A1                    0  8  0  0  0     0     0     0     0     8     0
B1                    0  0  8  0  0     0     0     8     0     0     0
C1                    0  0  0  8  0     0     0     0     0     0     0
D1                    0  0  0  0  8     8     0     0     0     0     0
A1:B1                 0  0  0  0  8     8     0     0     0     0     0
A1:C1                 0  0  0  0  0     0     8     0     0     0     0
A1:D1                 0  0  8  0  0     0     0     8     0     0     0
B1:C1                 0  0  0  0  0     0     0     0     8     0     0
B1:D1                 0  8  0  0  0     0     0     0     0     8     0
C1:D1                 0  0  0  0  0     0     0     0     0     0     8

#Representation des confusions entre les termes d'ordre 2 et ceux d'ordre 1
corrPlot(Design, scale = "corr")

kmax was reduced to number of factors

En regardant cet autre Plan Factoriel Fractionaire, qu’est-ce qui change par rapport au précédant ?
Quelle est sa résolution et pourquoi ?
Quel modèle peut être estimé avec ce plan ?

Ce Plan Factoriel Fractionaire n’utilise pas le même générateur. Le générateur utilisé (D=AB) est moins performant puisqu’avec celui-ci des effets d’ordre 1 sont confondus avec des effets d’ordre 2. Sa résolution est donc de III puiqu’il permet d’estimer uniquement un modèle d’ordre 1 et aucune interaction.

4.4 Plan factoriel fractionnaire - Version 3

Voici un code pour générer encore un plan factoriel frationnaire. Il y a ici 5 facteurs et on fractionne deux fois le plan \(2^5\). L’objectif que l’on poursuit est de trouver un plan qui permette d’estimer le modèle suivant : \[ Y = \beta_0 + \beta_1 A + \beta_2 B + \beta_3 C + \beta_4 D + \beta_5 E + \beta_6 A:B + \beta_7 A:C + \epsilon \] Etudiez le plan généré :

Quels sont les générateurs utilisés
Quelles sont les confusions liées
Dites si le plan répond à cet objectif.

Design=FrF2(8, nfactors=5)
MD=as.data.frame(desnum(Design))
MD

   A  B  C  D  E
1 -1 -1  1  1 -1
2  1  1 -1  1 -1
3 -1  1  1 -1 -1
4  1 -1  1 -1  1
5 -1 -1 -1  1  1
6  1  1  1  1  1
7  1 -1 -1 -1 -1
8 -1  1 -1 -1  1

# Reponse
generators(Design)

$generators
[1] "D=AB" "E=AC"

di=design.info(Design)
di$catlg.entry

Design:  5-2.1 
   8  runs,  5  factors,  
   Resolution  III 
   Generating columns:  3 5 
   WLP (3plus):  2 1 0 0 0 ,  0  clear 2fis

alldata=data.frame(y=1,desnum(Design))
aliases(lm(y~(.)^3,data=alldata))

                          
 A = C:E = B:D            
 B = C:D:E = A:D          
 C = B:D:E = A:E          
 D = B:C:E = A:B          
 E = B:C:D = A:C          
 B:C = D:E = A:B:E = A:C:D
 B:E = C:D = A:B:C = A:D:E

#Representation des confusions entre les termes d'ordre 2 et ceux d'ordre 1
corrPlot(Design, scale = "corr")

Non, ce modèle utilise les générateurs : D=AB et E=AC. L’effet de l’interaction A:B sera donc confondu avec l’effet de D.

4.4.1 Version 3 bis

Voici un code qui permet d’imposer à R de pouvoir estimer un modèle partiel donné.

Retrouvez les générateurs du plan généré
Retrouvez alias liés à ce plan et vérifiez que les termes du modèle que l’on désire estimer ne sont pas confondus entre eux.
Vérifiez aussi sur base de la matrice X’X que le plan généré permet bien d’estimer le modèle attendu de façon orthogonale.

Design=FrF2(8, nfactors=5, estimable = formula("~A+B+C+D+E+A:B+A:C"),clear=FALSE, res3=TRUE)
XM=as.data.frame(desnum(Design))
XM

   A  B  C  D  E
1 -1  1 -1 -1  1
2  1 -1 -1  1  1
3 -1 -1  1 -1  1
4 -1  1  1  1 -1
5 -1 -1 -1  1 -1
6  1 -1  1 -1 -1
7  1  1 -1 -1 -1
8  1  1  1  1  1

# Reponse
generators(Design)

$generators
[1] "C=DB" "E=DA"

alldata=data.frame(y=1,desnum(Design))
aliases(lm(y~(.)^3,data=alldata))

                          
 A = D:E = B:C:E          
 B = C:D = A:C:E          
 C = B:D = A:B:E          
 D = A:E = B:C            
 E = A:D = A:B:C          
 A:B = C:E = A:C:D = B:D:E
 A:C = B:E = A:B:D = C:D:E

X=model.matrix(~A+B+C+D+A:B+A:C,Design)
t(X)%*%X

            (Intercept) A1 B1 C1 D1 A1:B1 A1:C1
(Intercept)           8  0  0  0  0     0     0
A1                    0  8  0  0  0     0     0
B1                    0  0  8  0  0     0     0
C1                    0  0  0  8  0     0     0
D1                    0  0  0  0  8     0     0
A1:B1                 0  0  0  0  0     8     0
A1:C1                 0  0  0  0  0     0     8

#Representation des confusions entre les termes d'ordre 2 et ceux d'ordre 1
corrPlot(Design, scale = "corr")

Celui-ci le permet car il utilise deux autres générateurs : C=BD et E=DA. Les effets des interactions A:B et C:D ne sont donc pas confondus avec des effets d’ordre 1.

5 Analyse par régression des résultats d’un plan factoriel frationnaire

Ces analyses seront sur le même exemple que celui du PFC.

5.1 Plan Factoriel Fractionaire (PFF)

Pour le même problème que pour le PFC, on a réalisé un plan beaucoup moins coûteux :

Sans répétition
En fractionnant le plan à 16 essais en 8 essais uniquement à l’aide du générateur D=ABC (ou \(Ts=Cs*Ps*Rs\))

L’objectif ici est de vérifier si on parvient à obtenir autant d’information qu’avec le plan à 32 essais.

5.1.1 Téléchargement et visalisation des données

Téléchargez à présent les données du fichier Excel TP7_PFF.xlsx.

vaccinPFF <- read_excel("TP7_PFF.xlsx")
pander(vaccinPFF)

Cs	Ps	Rs	Ts	Rendement
-1	-1	-1	-1	58.5
-1	-1	1	1	36.65
-1	1	-1	1	75.25
-1	1	1	-1	71
1	-1	-1	1	78.75
1	-1	1	-1	73
1	1	-1	-1	93.5
1	1	1	1	73.5

5.1.2 Estimation du modèle d’ordre 2 complet

Estimez un modèle lm d’ordre 2 (toutes les interactions 2 à 2) et visualisez le résumé du modèle.

Que constatez-vous ? Vous attendiez-vous à ce problème
Etudiez les alias et également la matrice X’X pour encore mieux comprendre.

modPFF <- lm(Rendement ~ Cs+Ps+Rs+Ts+Cs:Ps+Cs:Rs+Cs:Ts+Ps:Rs+Ps:Ts+Rs:Ts, data = vaccinPFF)
summary(modPFF)


Call:
lm.default(formula = Rendement ~ Cs + Ps + Rs + Ts + Cs:Ps + 
    Cs:Rs + Cs:Ts + Ps:Rs + Ps:Ts + Rs:Ts, data = vaccinPFF)

Residuals:
ALL 8 residuals are 0: no residual degrees of freedom!

Coefficients: (3 not defined because of singularities)
            Estimate Std. Error t value Pr(>|t|)
(Intercept) 70.01875        NaN     NaN      NaN
Cs           9.66875        NaN     NaN      NaN
Ps           8.29375        NaN     NaN      NaN
Rs          -6.48125        NaN     NaN      NaN
Ts          -3.98125        NaN     NaN      NaN
Cs:Ps       -4.48125        NaN     NaN      NaN
Cs:Rs        0.04375        NaN     NaN      NaN
Cs:Ts        0.41875        NaN     NaN      NaN
Ps:Rs             NA         NA      NA       NA
Ps:Ts             NA         NA      NA       NA
Rs:Ts             NA         NA      NA       NA

Residual standard error: NaN on 0 degrees of freedom
Multiple R-squared:      1, Adjusted R-squared:    NaN 
F-statistic:   NaN on 7 and 0 DF,  p-value: NA

aliases(modPFF)

              
 Cs:Ps = Rs:Ts
 Cs:Rs = Ps:Ts
 Cs:Ts = Ps:Rs

X=model.matrix(modPFF)
t(X)%*%X

            (Intercept) Cs Ps Rs Ts Cs:Ps Cs:Rs Cs:Ts Ps:Rs Ps:Ts Rs:Ts
(Intercept)           8  0  0  0  0     0     0     0     0     0     0
Cs                    0  8  0  0  0     0     0     0     0     0     0
Ps                    0  0  8  0  0     0     0     0     0     0     0
Rs                    0  0  0  8  0     0     0     0     0     0     0
Ts                    0  0  0  0  8     0     0     0     0     0     0
Cs:Ps                 0  0  0  0  0     8     0     0     0     0     8
Cs:Rs                 0  0  0  0  0     0     8     0     0     8     0
Cs:Ts                 0  0  0  0  0     0     0     8     8     0     0
Ps:Rs                 0  0  0  0  0     0     0     8     8     0     0
Ps:Ts                 0  0  0  0  0     0     8     0     0     8     0
Rs:Ts                 0  0  0  0  0     8     0     0     0     0     8

#Representation des confusions entre les termes d'ordre 2 et ceux d'ordre 1
Design <- data2design(X[,2:5])
corrPlot(Design, scale = "corr")

kmax was reduced to number of factors

Car il n’est pas possible d’estimer tous les termes du modèle demandé. Effectivement, il s’agit d’un plan factoriel fractionaire avec seulement 8 essais. Cela ne permet pas d’estimer les termes d’interactions demandés.

5.1.3 Estimation d’un modèle simplifié

Essayez maintenant d’estimer le modèle simplifié, en gardant uniquement les termes principaux et l’interaction la plus significative identifiée pour le plan complet.

Comparez le tableau des coefficients de ce modèle et celui du modèle simplifié estimé pour le plan factoriel complet. Sont-ils comparables ?
Comparez également le graphe des effets principaux et d’interaction.
Comparez les écarts-types des modèles avec le PFC et le PFF. Pouvez-vous expliquer la différence entre les 2 ?
En terme d’interprétation, à quoi correspond la valeur donneé dans le summary pour Cs ?

modPFF <- lm(Rendement ~ Cs+Ps+Rs+Ts+Cs:Ps+Cs:Rs, data = vaccinPFF)
summary(modPFF)


Call:
lm.default(formula = Rendement ~ Cs + Ps + Rs + Ts + Cs:Ps + 
    Cs:Rs, data = vaccinPFF)

Residuals:
      1       2       3       4       5       6       7       8 
 0.4187 -0.4187 -0.4187  0.4187  0.4187 -0.4187 -0.4187  0.4187 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept) 70.01875    0.41875 167.209  0.00381 **
Cs           9.66875    0.41875  23.090  0.02755 * 
Ps           8.29375    0.41875  19.806  0.03212 * 
Rs          -6.48125    0.41875 -15.478  0.04107 * 
Ts          -3.98125    0.41875  -9.507  0.06671 . 
Cs:Ps       -4.48125    0.41875 -10.701  0.05932 . 
Cs:Rs        0.04375    0.41875   0.104  0.93373   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.184 on 1 degrees of freedom
Multiple R-squared:  0.9993,    Adjusted R-squared:  0.9949 
F-statistic: 228.3 on 6 and 1 DF,  p-value: 0.05062

summary(modPFC2)


Call:
lm.default(formula = Rendement ~ Cs + Ps + Rs + Ts + Cs:Ps + 
    Cs:Rs, data = vaccinPFC)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.8656 -0.5531 -0.1031  0.5016  1.9094 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  70.0344     0.1808 387.252  < 2e-16 ***
Cs            9.4156     0.1808  52.063  < 2e-16 ***
Ps            8.0969     0.1808  44.771  < 2e-16 ***
Rs           -6.9344     0.1808 -38.343  < 2e-16 ***
Ts           -3.8781     0.1808 -21.444  < 2e-16 ***
Cs:Ps        -4.6969     0.1808 -25.971  < 2e-16 ***
Cs:Rs         0.5719     0.1808   3.162  0.00408 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.023 on 25 degrees of freedom
Multiple R-squared:  0.9966,    Adjusted R-squared:  0.9958 
F-statistic:  1222 on 6 and 25 DF,  p-value: < 2.2e-16

MEPlot(modPFF)

MEPlot(modPFC2)

term2 <- summary(modPFF)$coefficients[2,1]

Le PFC donne un plus faible écart-type mais nous pouvons constater que l’estimation des paramètres avec le PFF est assez proche et pourtant beaucoup moins coûteuse en terme d’essais. Ceci illustre bien l’intérêt d’utiliser des PFF.
La valeur de 9.66875 signifie que lorsque tous les autres facteurs sont à la moitié de leur domaine, le fait d’augmenter le facteur Cs de la moitié jusqu’à sa valeur maximale, fait augmenter le Rendement de 9.66875%. C’est donc bien le facteur Cs qui à le plus d’impact sur le Rendement.

TP7 : Plans factoriels complets et fractionnaires à 2 niveaux

LBRAI2222

mai 09, 2023