TP8 : Plans pour surfaces de réponses
LBIR2222
avril 22, 2024
1 Importation de la fonction plotDesign
La ligne de code ci-dessous vous permet d’importer la fonction plotDesign dans votre Rmarkdown. Cette fonction est utile pour pouvoir représenter graphiquement vos plans d’expérience.
Attention, pour que cela marche, vous devez impérativement avoir téléchargé la fonction plotDesign sur Moodle et l’avoir chargée dans le même dossier que celui du Rmarkdown.
2 Plan Factoriel Complet
Nous souhaitons estimer les paramètres d’un modèle polynomial quadratique complet à l’aide d’un Plan Factoriel Complet.
- Combien de niveaux différents devrions-nous tester au minimum pour chaque facteur ?
- Combien d’essais comprendra le plan à 3 facteurs ? et à k facteurs ?
- Combien de paramètres comprend le modèle à estimer pour 3 facteurs et pour k facteurs ?
- Le code ci-dessous permet de générez ce plan à 3 facteurs en définissant les niveaux comme -1, 0 et 1.
Réponse :
Pour pouvoir estimer des effets quadratiques, il faut au minimum 3 valeurs par facteur.
Dans un plan factoriel complet, toutes les combinaisons de chacun des niveaux des facteurs sont testées : \(3^3 = 27\) essais. \(3^k\) essais pour \(k\) facteurs.
Le modèle a 3 facteurs a 10 paramètres : terme constant, 3 termes principaux \(X_i\), 3 termes quadratiques \(X_i^2\), 3 termes d’interaction \(X_i*X_j\). \(1+k+k+k*(k-1)/2\) pour \(k\) facteurs.
#Generation du design du PFC
PFC <- DoE.base::fac.design(nfactors = 3, factor.names=list(x1=c(-1,0,1),x2=c(-1,0,1),x3=c(-1,0,1)),randomize=FALSE)creating full factorial with 27 runs ...adfPFC <- as.data.frame(PFC) #transformation en dataframe
#Representation du design en 2D
plotDesign(adfPFC, x="x1", y="x2", cols = "x3", title = "Plan factoriel complet à 3 facteurs")Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
ℹ Please use tidy evaluation idioms with `aes()`.
ℹ See also `vignette("ggplot2-in-packages")` for more information.
This warning is displayed once every 8 hours.
Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
generated.
3 Génération d’un Plan Box & Behnken
Le plan de Box & Behnken permet aussi d’estimer un modèle quadratique complet à un coût moindre que le plan factoriel complet.
- A l’aide du script ci-dessous, générez un plan Box-Behnken pour estimer un modèle quadratique à 3 facteurs.
- Justifiez graphiqment pourquoi il s’agit d’un Plan Box & Behnken.
- Combien d’essais y a-t-il au total ?
Fraction d’un plan factoriel complet \(3^k\) qui permet d’estimer un modèle quadratique. Il se construit en combinant de manière particulière un plan factoriel \(2^{k-r}\) avec un plan de type “bloc balancé incomplet” dans le but d’obtenir un plan à 3 niveaux le plus rotatable possible qui évite les sommets du cube. Pour 3 facteurs et 1 point au centre, cela fait 13 essais.
#Generation du design du PBB
PBB <- rsm::bbd(k = 3, n0 = 1, randomize = F) #fonction specifique du package rsm pour les PBB
adfPBB <- as.data.frame(PBB) #transformation en dataframe
#Representation du design en 2D
plotDesign(adfPBB, x="x1", y="x2", cols = "x3", title = "Plan Box & Behnken à 3 facteurs")
y numeric
4 Génération d’un Plan Composite Centré
Voici comment générer des plans composites centrés. On va en générer un avec les points sur les faces du cube (alpha=1) et un rotatables (isovariance par rotation) pour 3 facteurs et un essai au centre.
4.1 Plan composite avec points sur face du cube :
#Generation du design du PCC1
PCC1 <- rsm::ccd(basis = 3, n0 = c(1,0), randomize = F, alpha = "faces") #fonction specifique du package rsm pour les PCC
adfPCC1 <- as.data.frame(PCC1) #transformation en dataframe
pander(adfPCC1)| run.order | std.order | x1 | x2 | x3 | Block |
|---|---|---|---|---|---|
| 1 | 1 | -1 | -1 | -1 | 1 |
| 2 | 2 | 1 | -1 | -1 | 1 |
| 3 | 3 | -1 | 1 | -1 | 1 |
| 4 | 4 | 1 | 1 | -1 | 1 |
| 5 | 5 | -1 | -1 | 1 | 1 |
| 6 | 6 | 1 | -1 | 1 | 1 |
| 7 | 7 | -1 | 1 | 1 | 1 |
| 8 | 8 | 1 | 1 | 1 | 1 |
| 9 | 9 | 0 | 0 | 0 | 1 |
| 1 | 1 | -1 | 0 | 0 | 2 |
| 2 | 2 | 1 | 0 | 0 | 2 |
| 3 | 3 | 0 | -1 | 0 | 2 |
| 4 | 4 | 0 | 1 | 0 | 2 |
| 5 | 5 | 0 | 0 | -1 | 2 |
| 6 | 6 | 0 | 0 | 1 | 2 |
#Representation du design en 2D
plotDesign(adfPCC1, x="x1", y="x2", cols = "x3", title = "Plan Composite Centré à k=3 et alpha = 1")
y numeric
4.2 Plan composite rotatable :
#Generation du design du PCC2 -> alpha = sqrt^4(2^3)
PCC2 <- rsm::ccd(basis = 3, n0 = c(1,0), randomize = F, alpha = "rotatable") #fonction specifique du package rsm pour les PCC
adfPCC2 <- as.data.frame(PCC2) #transformation en dataframe
pander(adfPCC2)| run.order | std.order | x1 | x2 | x3 | Block |
|---|---|---|---|---|---|
| 1 | 1 | -1 | -1 | -1 | 1 |
| 2 | 2 | 1 | -1 | -1 | 1 |
| 3 | 3 | -1 | 1 | -1 | 1 |
| 4 | 4 | 1 | 1 | -1 | 1 |
| 5 | 5 | -1 | -1 | 1 | 1 |
| 6 | 6 | 1 | -1 | 1 | 1 |
| 7 | 7 | -1 | 1 | 1 | 1 |
| 8 | 8 | 1 | 1 | 1 | 1 |
| 9 | 9 | 0 | 0 | 0 | 1 |
| 1 | 1 | -1.682 | 0 | 0 | 2 |
| 2 | 2 | 1.682 | 0 | 0 | 2 |
| 3 | 3 | 0 | -1.682 | 0 | 2 |
| 4 | 4 | 0 | 1.682 | 0 | 2 |
| 5 | 5 | 0 | 0 | -1.682 | 2 |
| 6 | 6 | 0 | 0 | 1.682 | 2 |
#Representation du design en 2D
plotDesign(adfPCC2, x="x1", y="x2", cols = "x3", title = "Plan Composite Centré à 3 facteurs et rotatable")
y numeric
5 Comparaison des propriétés des 3 classes de plan
5.1 Nombre d’essais
Quel plan comprend le moins d’essais ?
C’est le plan de Box & Behnken : 13 essais contre 27 pour le plan factoriel complet et 15 pour les plans composites.
5.2 Précision et indépendance des paramètres
Calculez la matrice \((X'X)^{-1}\) pour chacun des plans générés et comparez les en termes des écart-types des paramètres et des corrélations entre paramètres.
- Vous devez pour cela d’abord créer les matrices du modèle X pour chaque plan (standardisé pour que ce soit comparable)
- Ensuite vous calculez \((X'X)^{-1}\)
- Puis vous comparez les écart-types et corrélations
# Génération des matrices du modèle
ModQuad=formula(~x1+x2+x3+I(x1^2)+I(x2^2)+I(x3^2)+x1:x2+x1:x3+x2:x3)
adfPFC=apply(PFC,2,as.numeric) # necessaire de convertir la dataframe en numeric pour le PFC. Les autres le sont deja
adfPFC <- as.data.frame(adfPFC)
XPFC=model.matrix(ModQuad,data=adfPFC)
XPBB=model.matrix(ModQuad,data=adfPBB)
XPCC1=model.matrix(ModQuad,data=adfPCC1)
XPCC2=model.matrix(ModQuad,data=adfPCC2)
# Calcul des matrices de VC (X'X)^-1
MVCPFC=solve(t(XPFC)%*%XPFC)
MVCPBB=solve(t(XPBB)%*%XPBB)
MVCPCC1=solve(t(XPCC1)%*%XPCC1)
MVCPCC2=solve(t(XPCC2)%*%XPCC2)
# Comparaison des écart-types des paramètres
dfse=cbind(PFC=sqrt(diag(MVCPFC)),PBB=sqrt(diag(MVCPBB)),PCC1=sqrt(diag(MVCPCC1)),PCC2=sqrt(diag(MVCPCC2)))
pander(dfse)| PFC | PBB | PCC1 | PCC2 | |
|---|---|---|---|---|
| (Intercept) | 0.5092 | 1 | 0.5375 | 0.9942 |
| x1 | 0.2357 | 0.3536 | 0.3162 | 0.2706 |
| x2 | 0.2357 | 0.3536 | 0.3162 | 0.2706 |
| x3 | 0.2357 | 0.3536 | 0.3162 | 0.2706 |
| I(x1^2) | 0.4082 | 0.6614 | 0.6236 | 0.4065 |
| I(x2^2) | 0.4082 | 0.6614 | 0.6236 | 0.4065 |
| I(x3^2) | 0.4082 | 0.6614 | 0.6236 | 0.4065 |
| x1:x2 | 0.2887 | 0.5 | 0.3536 | 0.3536 |
| x1:x3 | 0.2887 | 0.5 | 0.3536 | 0.3536 |
| x2:x3 | 0.2887 | 0.5 | 0.3536 | 0.3536 |
# Matrices des corrélations (attention on utilise la fonction corrPlot du package corrPlot et pas du package DoE.base)
corrplot::corrplot(cov2cor(MVCPFC),main="Plan factoriel complet")
Le plan Box & Behnken est celui qui utilise le moins d’essais mais il est celui pour lequel les écarts-types des paramètres sont les plus grands. Le premier PCC est plus intéressant en terme d’orthogonalité (estimation des paramètres) mais moins intéressant que le PCC2 qui est rotatable et pour qui les écarts-types sont plus faibles de manière générale (meilleures qualités de prédiction).
5.3 Etude des variances de prédiction liées aux plans
Sans considérer le PFC, quel plan a la meilleure variance de prédiction ?
5.3.1 Variances de prédiction avec contour plots
Attention ces graphes ne regardent ce qui se passe que pour 2 variables en supposant que la troisième =0.
PBB.df <- as.data.frame((PBB[,3:5]))
###
par(mfrow=c(2,2))
varfcn(adfPFC, ~SO(x1,x2,x3), contour=TRUE) #NB : SO(x1,x2,x3) = modele quadratique complet
varfcn(PBB, ~SO(x1,x2,x3), contour=TRUE)
varfcn(PCC1, ~SO(x1,x2,x3), contour=TRUE)
varfcn(PCC2, ~SO(x1,x2,x3), contour=TRUE)
On peut observer que le seul plan rotatable est le PCC2. Effectivement, lorsqu’un plan est rotatable, il produit des informations qui prédisent ŷ avec la même précision à tous les points équidistants de l’origine codée du plan.
5.3.2 Variances de prédiction directionnelles
La fonction varfcn() génère des directions par défaut le
long d’un axe, et sur une diagonale passant par un coin dans chaque
dimension. Ainsi, le graphique produit montre comment la variance évolue
le long de chacun de ces vecteurs, pour les distances fournies. Dans une
conception rotative, ces courbes seront toutes les mêmes (car peu
importe la direction dans laquelle on va, l’évolution de la variance
devrait être la même).
Le seul plan pour lequel c’est le cas est le PCC2. Il s’agit bien d’un plan rotatable comme mentionné précédamment. Nous pouvons également remarquer que le PBB semble avoir de bonnes propriétés proches de la rotatabilité également.
par(mfrow=c(2,2))
varfcn(adfPFC, ~SO(x1,x2,x3), contour=FALSE)
varfcn(PBB, ~SO(x1,x2,x3), contour=FALSE)
varfcn(PCC1, ~SO(x1,x2,x3), contour=FALSE)
varfcn(PCC2, ~SO(x1,x2,x3), contour=FALSE)
5.3.3 Fraction of Design Space Plot
La fonction suivante prend des points au hasard dans le domaine et calcule leur variance de prédiction. Cela permet de réaliser un graphe de la fréquence cumulée de la variance de prédiction dans le domaine. On peut observer que le PBB et le PCC2 sont tous les deux meilleurs en terme de variance de prédiction.
PBB.df <- as.data.frame((PBB[,3:5]))
form <- formula(~x1 + x2 + x3)
out <- vdg::spv(n = 1000, design = PBB.df, type = "spherical", formula = form)
out2 <- vdg::spv(n = 1000, design = PCC1, type = "spherical", formula = form)
out3 <- vdg::spv(n = 1000, design = PCC2, type = "spherical", formula = form)
plot(out)$fds / plot(out2)$fds / plot(out3)$fdsPlot 5 = 'boxplots' cannot be produced: 'at' is FALSE
Plot 5 = 'boxplots' cannot be produced: 'at' is FALSE
Plot 5 = 'boxplots' cannot be produced: 'at' is FALSE
6 Surface de réponse - Blanchiment du papier
6.1 Contexte
Le blanchiment du papier est le traitement chimique de la pâte de bois ou du papier recyclé dans le but d’éclaircir sa couleur et de blanchir la pâte, produisant ainsi un papier à la blancheur adaptée à l’usage final du papier. De nos jours, ce procédé est toujours en évolution suite aux efforts pour minimiser son impact environnemental et pour réduire sa dépendance à la chlorine.
6.2 Etude du plan expérimental
Importez le fichier
TP8_test.csvReprésentez le plan graphiquement et/ou sous forme de tableau
Quels sont les facteurs testés dans l’expérience ? Sont-ils quantitatifs ou qualitatifs ?
Quel est le type de plan d’expérience ?
#Importation des donnees
mydata <- read.csv("TP8_test.csv")
mydata$TypAlka <- as.factor(mydata$TypAlka) #specifie que TypAlka est un facteur qualitatif
#Representation du design en 2D
plotDesign(as.data.frame(mydata[,1:3]),x="QtH2O2",y="QtAlka",cols="TypAlka")
y numeric
#Representation du design en 3D
plotly::plot_ly(mydata, x=~QtH2O2,y=~QtAlka,z=~TypAlka, type= "scatter3d", mode = "markers") %>% layout(title = "Blanchiment de papier")%>% add_markers(color=~TypAlka)Les facteurs testés sont les suivants :
- QtH202 : facteur quantitatif - la quantité de H2O2
- QtAlka : facteur quantitatif - la quantité d’alcali
- TypAlka : facteur qualitatif - le type d’alcali
Le plan est un plan factoriel complet 3x3x3 avec le point central qui a été répété trois fois pour chaque type d’Alcali.
6.3 Représentation graphique des réponses
- Sur base du script suivant, que pouvez-vous dire de l’effet des différents facteurs sur la Brightness (blancheur) ?
ggplot(data = mydata, aes(x = QtH2O2, y = Brightness, colour=factor(QtAlka))) +
geom_point() +
theme_bw()+ labs(color = "QtAlka") +
geom_smooth(se = FALSE, method = lm, formula = y ~ poly(x, 2),linewidth=0.5) +
facet_wrap( ~ TypAlka) 
Le NaOH est l’alcali qui donne la meilleure blancheur quelque soient les valeurs des quantités d’Alcali et de H2O2
Au niveau des effets des quantités d’Alcali et de H2O2, nous pouvons constater qu’il semble y avoir une interaction entre ces deux facteurs. Celle-ci semble avoir un effet similaire pour les trois TypAlka.
6.4 Modélisation
6.4.1 Standardisation des données
- Nous avons vu au TP précédant qu’il était intéressant de travailler avec des données standardisées. Voici un script qui permet de standardiser les données.
Ajout des colonnes:
#Standardisation de QtH202
milieuQtH2O2=(max(mydata$QtH2O2)+min(mydata$QtH2O2))/2
etendueQtH2O2=(max(mydata$QtH2O2)-min(mydata$QtH2O2))/2
mydata$QtH2O2s <- (mydata$QtH2O2-milieuQtH2O2)/etendueQtH2O2
#Standardisation de QtAlka
milieuQtAlka=(max(mydata$QtAlka)+min(mydata$QtAlka))/2
etendueQtAlka=(max(mydata$QtAlka)-min(mydata$QtAlka))/2
mydata$QtAlkas <- (mydata$QtAlka-milieuQtAlka)/etendueQtAlka
pander(mydata)| QtH2O2 | QtAlka | TypAlka | Brightness | QtH2O2s | QtAlkas |
|---|---|---|---|---|---|
| 0.4 | 0.6 | CaO | 68.33 | -1 | -1 |
| 1 | 0.6 | CaO | 73.21 | 0 | -1 |
| 1.6 | 0.6 | CaO | 76.77 | 1 | -1 |
| 0.4 | 1 | CaO | 73.52 | -1 | 0 |
| 1 | 1 | CaO | 76.05 | 0 | 0 |
| 1 | 1 | CaO | 75.59 | 0 | 0 |
| 1 | 1 | CaO | 76.57 | 0 | 0 |
| 1.6 | 1 | CaO | 77.08 | 1 | 0 |
| 0.4 | 1.4 | CaO | 73.95 | -1 | 1 |
| 1 | 1.4 | CaO | 73.67 | 0 | 1 |
| 1.6 | 1.4 | CaO | 72.45 | 1 | 1 |
| 0.4 | 0.6 | MgO | 68.23 | -1 | -1 |
| 1 | 0.6 | MgO | 73.05 | 0 | -1 |
| 1.6 | 0.6 | MgO | 76.49 | 1 | -1 |
| 0.4 | 1 | MgO | 74.25 | -1 | 0 |
| 1 | 1 | MgO | 76.37 | 0 | 0 |
| 1 | 1 | MgO | 75.92 | 0 | 0 |
| 1 | 1 | MgO | 76.08 | 0 | 0 |
| 1.6 | 1 | MgO | 76.78 | 1 | 0 |
| 0.4 | 1.4 | MgO | 74.82 | -1 | 1 |
| 1 | 1.4 | MgO | 74.62 | 0 | 1 |
| 1.6 | 1.4 | MgO | 71.94 | 1 | 1 |
| 0.4 | 0.6 | NaOH | 70.55 | -1 | -1 |
| 1 | 0.6 | NaOH | 75.79 | 0 | -1 |
| 1.6 | 0.6 | NaOH | 79.04 | 1 | -1 |
| 0.4 | 1 | NaOH | 76.82 | -1 | 0 |
| 1 | 1 | NaOH | 78.94 | 0 | 0 |
| 1 | 1 | NaOH | 78.51 | 0 | 0 |
| 1 | 1 | NaOH | 79.13 | 0 | 0 |
| 1.6 | 1 | NaOH | 79.21 | 1 | 0 |
| 0.4 | 1.4 | NaOH | 76.77 | -1 | 1 |
| 1 | 1.4 | NaOH | 76.24 | 0 | 1 |
| 1.6 | 1.4 | NaOH | 74.07 | 1 | 1 |
6.4.2 Estimation du modèle
Estimez maintenant votre modèle quadratique complet (interactions + termes quadratiques). Faudra-t-il utiliser
lmoulmer?Notez aussi que dans les plans multifacteurs, il est d’usage de coder les facteurs catégoriels autrement qu’en choisissant le premier niveau comme niveau de référence. Il est conseillé d’utiliser l’option ‘contrasts = list(TypAlka = “contr.sum”)’. Avec ce codage, le niveau de référence doit être considéré comme un Alcali moyen fictif.
Utilisez la fonction
model.matrixpour voir comment ce codage est réaliséRéalisez une Anova pour voir quels effets sont significatifs
Réponse
- Nous utilisons la fonction
lmcar celle-ci est adaptée aux modèles linéaires à facteurs fixes, sans effets aléatoires. Effectivement, le seul facteur qui aurait pu être traité en aléatoire est le facteur TypAlka (les autres étant quantitatifs). Comme nous nous intéressons à 3 types d’alcali bien définis (et ne souhaitons pas étudier l’effet des autres facteurs sur n’importe quel type d’alcali), nous traiterons ce facteur comme fixe.
Estimation du modèle:
# Estimation du modele
mod_sd <- lm(Brightness ~ QtH2O2s * QtAlkas + TypAlka + QtH2O2s:TypAlka + QtAlkas:TypAlka + I(QtH2O2s^2) + I(QtAlkas^2), data = mydata, contrasts = list(TypAlka = "contr.sum"))
summary(mod_sd)
Call:
lm.default(formula = Brightness ~ QtH2O2s * QtAlkas + TypAlka +
QtH2O2s:TypAlka + QtAlkas:TypAlka + I(QtH2O2s^2) + I(QtAlkas^2),
data = mydata, contrasts = list(TypAlka = "contr.sum"))
Residuals:
Min 1Q Median 3Q Max
-0.55675 -0.11854 -0.02718 0.09962 0.44871
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 77.03281 0.08536 902.415 < 2e-16 ***
QtH2O2s 1.47722 0.06793 21.745 7.00e-16 ***
QtAlkas 0.39278 0.06793 5.782 9.71e-06 ***
TypAlka1 -0.88606 0.07095 -12.488 3.48e-11 ***
TypAlka2 -0.76242 0.07095 -10.745 5.42e-10 ***
I(QtH2O2s^2) -0.77868 0.10455 -7.448 2.54e-07 ***
I(QtAlkas^2) -2.62535 0.10455 -25.111 < 2e-16 ***
QtH2O2s:QtAlkas -2.68917 0.08320 -32.321 < 2e-16 ***
QtH2O2s:TypAlka1 0.27278 0.09607 2.839 0.00982 **
QtH2O2s:TypAlka2 -0.15889 0.09607 -1.654 0.11303
QtAlkas:TypAlka1 -0.09944 0.09607 -1.035 0.31240
QtAlkas:TypAlka2 0.20889 0.09607 2.174 0.04126 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.2882 on 21 degrees of freedom
Multiple R-squared: 0.9929, Adjusted R-squared: 0.9892
F-statistic: 268.2 on 11 and 21 DF, p-value: < 2.2e-16Anova Table (Type III tests)
Response: Brightness
Sum Sq Df F value Pr(>F)
(Intercept) 67648 1 8.1435e+05 < 2.2e-16 ***
QtH2O2s 39 1 4.7285e+02 6.996e-16 ***
QtAlkas 3 1 3.3429e+01 9.714e-06 ***
TypAlka 45 2 2.7039e+02 1.030e-15 ***
I(QtH2O2s^2) 5 1 5.5474e+01 2.540e-07 ***
I(QtAlkas^2) 52 1 6.3059e+02 < 2.2e-16 ***
QtH2O2s:QtAlkas 87 1 1.0447e+03 < 2.2e-16 ***
QtH2O2s:TypAlka 1 2 4.0673e+00 0.03214 *
QtAlkas:TypAlka 0 2 2.3655e+00 0.11844
Residuals 2 21
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 16.4.3 Simplification du modèle
Sur base des résultats de l’ANOVA, pouvez-vous supprimer des termes/effets dans votre modèle ? Si oui, mettez à jour votre modèle (pour un \(\alpha = 0.01\)).
Sur base des résultats de l’ANOVA, nous constatons que les interactions entre QtH2O2s:TypAlka et QtAlkas:TypAlka ne sont pas significatives. Nous pouvons les enlever du modèle :
# Estimation du modele
mod_sd <- lm(Brightness ~ QtH2O2s * QtAlkas + TypAlka + I(QtH2O2s^2) + I(QtAlkas^2), data = mydata, contrasts = list(TypAlka = "contr.sum"))
summary(mod_sd)
Call:
lm.default(formula = Brightness ~ QtH2O2s * QtAlkas + TypAlka +
I(QtH2O2s^2) + I(QtAlkas^2), data = mydata, contrasts = list(TypAlka = "contr.sum"))
Residuals:
Min 1Q Median 3Q Max
-0.55675 -0.19038 -0.09675 0.23552 0.58219
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 77.03281 0.09935 775.347 < 2e-16 ***
QtH2O2s 1.47722 0.07907 18.683 3.35e-16 ***
QtAlkas 0.39278 0.07907 4.968 4.05e-05 ***
TypAlka1 -0.88606 0.08258 -10.729 7.60e-11 ***
TypAlka2 -0.76242 0.08258 -9.232 1.57e-09 ***
I(QtH2O2s^2) -0.77868 0.12168 -6.399 1.06e-06 ***
I(QtAlkas^2) -2.62535 0.12168 -21.576 < 2e-16 ***
QtH2O2s:QtAlkas -2.68917 0.09684 -27.770 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.3355 on 25 degrees of freedom
Multiple R-squared: 0.9886, Adjusted R-squared: 0.9854
F-statistic: 309.8 on 7 and 25 DF, p-value: < 2.2e-16Anova Table (Type III tests)
Response: Brightness
Sum Sq Df F value Pr(>F)
(Intercept) 67648 1 601162.770 < 2.2e-16 ***
QtH2O2s 39 1 349.060 3.350e-16 ***
QtAlkas 3 1 24.678 4.053e-05 ***
TypAlka 45 2 199.606 4.261e-16 ***
I(QtH2O2s^2) 5 1 40.952 1.063e-06 ***
I(QtAlkas^2) 52 1 465.505 < 2.2e-16 ***
QtH2O2s:QtAlkas 87 1 771.174 < 2.2e-16 ***
Residuals 3 25
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 16.4.4 Analyse des résidus
Analysez les résidus de votre modèle. Respectent-ils les hypothèses du modèle ?
Les hypothèses semblent bien respectées.
6.4.5 Comparaison des alcalis
Sur base de votre modèle, quel type d’alcali semble le plus intéressant à utiliser pour avoir une Brightness maximale ? Les différences de Brightness par rapport aux autres alkalis sont-elles significativement différentes ?
L’alcali NaOH obtient une Brightness significativement supérieure aux deux autres alcali.
contrast estimate SE df t.ratio p.value
CaO - MgO -0.124 0.143 25 -0.864 0.6672
CaO - NaOH -2.535 0.143 25 -17.719 <.0001
MgO - NaOH -2.411 0.143 25 -16.855 <.0001
P value adjustment: tukey method for comparing a family of 3 estimates 6.5 Analyses avec le
package rsm
Pour continuer la suite des analyses nous utiliserons un autre
package, il s’agit du package rsm. Sur base du script
réalisé ci-dessous, quelles sont les analyses supplémentaires réalisées
avec la fonction rsm? Observez-vous des différences avec
les résultats de votre modèle précédant ?
Estimation du modèle:
### Methode rsm
# SO : Second Order model
res.model2 <- rsm(Brightness ~ SO(QtH2O2s, QtAlkas) + TypAlka, data = mydata, contrast=list(TypAlka= "contr.sum"))
summary(res.model2)
Call:
rsm(formula = Brightness ~ SO(QtH2O2s, QtAlkas) + TypAlka, data = mydata,
contrasts = list(TypAlka = "contr.sum"))
Estimate Std. Error t value Pr(>|t|)
(Intercept) 77.032807 0.099353 775.3469 < 2.2e-16 ***
QtH2O2s 1.477222 0.079067 18.6831 3.350e-16 ***
QtAlkas 0.392778 0.079067 4.9677 4.053e-05 ***
QtH2O2s:QtAlkas -2.689167 0.096837 -27.7700 < 2.2e-16 ***
QtH2O2s^2 -0.778684 0.121682 -6.3994 1.063e-06 ***
QtAlkas^2 -2.625351 0.121682 -21.5756 < 2.2e-16 ***
TypAlka1 -0.886061 0.082583 -10.7293 7.602e-11 ***
TypAlka2 -0.762424 0.082583 -9.2322 1.567e-09 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Multiple R-squared: 0.9886, Adjusted R-squared: 0.9854
F-statistic: 309.8 on 7 and 25 DF, p-value: < 2.2e-16
Analysis of Variance Table
Response: Brightness
Df Sum Sq Mean Sq F value Pr(>F)
FO(QtH2O2s, QtAlkas) 2 42.056 21.028 186.8687 9.240e-16
TWI(QtH2O2s, QtAlkas) 1 86.779 86.779 771.1741 < 2.2e-16
PQ(QtH2O2s, QtAlkas) 2 70.275 35.137 312.2513 < 2.2e-16
TypAlka 2 44.923 22.461 199.6058 4.261e-16
Residuals 25 2.813 0.113
Lack of fit 3 0.050 0.017 0.1332 0.9392
Pure error 22 2.763 0.126
Stationary point of response surface:
QtH2O2s QtAlkas
7.085221 -3.553918
Eigenanalysis:
eigen() decomposition
$values
[1] -0.07092932 -3.33310576
$vectors
[,1] [,2]
QtH2O2s -0.8848967 0.4657874
QtAlkas 0.4657874 0.8848967Anova Table (Type III tests)
Response: Brightness
Sum Sq Df F value Pr(>F)
(Intercept) 67648 1 601162.77 < 2.2e-16 ***
FO(QtH2O2s, QtAlkas) 42 2 186.87 9.240e-16 ***
TWI(QtH2O2s, QtAlkas) 87 1 771.17 < 2.2e-16 ***
PQ(QtH2O2s, QtAlkas) 70 2 312.25 < 2.2e-16 ***
TypAlka 45 2 199.61 4.261e-16 ***
Residuals 3 25
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Les résultats du summary sont équivalents (les effets quadratiques sont regroupés sous “PQ()”; les effets d’ordre 1 sous”FO()” et l’effet d’interaction sous “TWI()”). Un test supplémentaire est réalisé, il s’agit du test de Lack of Fit. Celui-ci test les hypothèses suivantes :
-\(H_0 :\) il n’y a pas de problème
d’ajustement du modèle par rapport aux mesures réelles
-\(H_1 :\) il y a un problème
d’ajustement du modèle par rapport aux mesures réelles
NB : Le test de Lack of Fit nécessite d’avoir des répétitions au centre de notre plan d’expérience. Cela permet de calculer la ‘Pure error’ et de la comparer à l’erreur résiduelle.
6.5.1 Voici une manière de représenter les courbes d’iso-réponses par type d’alcali
Sur base de ces figures, identifiez (graphiquement) des valeurs de type d’alcali, QtH2O2s et QtAlkas qui permettraient d’obtenir une Brigntness proche de 77,5 (plusieurs réponses possibles). Calculez aussi l’intervalle de prédictions.
Retransformez ces valeurs en valeurs réelles.
par (mfrow = c(1,3))
for (i in c("CaO","MgO","NaOH")){
contour.lm(res.model2, ~ QtH2O2s+QtAlkas , image = T, at = list(TypAlka=i), main = i)
}
# Réponse
## valeurs trouvées graphiquement
H202_opts = -0.85
alka_opts = 0.5
predictxy <- data.frame("QtH2O2s" = c(H202_opts),
"QtAlkas" = c(alka_opts),
"TypAlka" = "NaOH")
## modèle linéaire (possibilité de réutilisé celui plus haut directement)
mod_sd <- lm(Brightness ~ QtH2O2s * QtAlkas + I(QtH2O2s^2) + I(QtAlkas^2) + TypAlka, data = mydata, contrasts = list(TypAlka = "contr.sum"))
## prédiction des valeurs et de l'intervale de confiance sur les predictions
res_predict <- predict(mod_sd,newdata = predictxy, interval = 'prediction')
print(res_predict) fit lwr upr
1 77.546 76.79575 78.29626## transformation en valeurs réelles
H202_opt <- H202_opts*(max(mydata$QtH2O2)-min(mydata$QtH2O2))/2 + mean(mydata$QtH2O2)
alka_opt <- alka_opts*(max(mydata$QtAlka)-min(mydata$QtAlka))/2 + mean(mydata$QtAlka)
print(H202_opt)[1] 0.49[1] 1.2# Remontrer où se trouve le point choisi
contour.lm(res.model2, ~ QtH2O2s+QtAlkas, image = T, at = list(TypAlka="NaOH"), main = "NaOH")
points(predictxy,col="red", pch=16)