In ehassler/MontgomeryDAE: MontgomeryDAE - Data and examples from Montgomery's "Design and Analysis of Experiments"

source('_format.R')

Chapter 9 in the text introduces designs with three-level factors and some non-regular fractional factorial (two-level) designs. The chapter also introduces optimal designs which (as of this writing) are not very well supported in R.

Example 9.1

The syrup loss experiment of Table 9.1 investigated how nozzle design, the filling speed, and the operating pressure affect the amount of syrup lost while filling containers. Each factor is investigated at three levels using a $3^3$ design. Note that we want to treat Speed and Pressure as factors for the ANOVA model:

model <- aov(SyrupLoss ~ NozzleType * factor(Speed) * factor(Pressure), data=Table9.1)
print(summary(model))

\noindent Next we'll fit regression models for each nozzle type to fit a contour plot for each nozzle type:

code.speed <- function(x){
    (x - 120) / (140 - 120)
}

code.pressure <- function(x){
    (x - 15) / (20 - 15)
}

df.reg <- data.frame(
    'NozzleType'=Table9.1$NozzleType,
    'Speed'=code.speed(Table9.1$Speed),
    'Pressure'=code.pressure(Table9.1$Pressure),
    'SyrupLoss'=Table9.1$SyrupLoss
)

model1 <- lm(SyrupLoss ~ Speed * Pressure + I(Speed^2) + I(Pressure^2), data=df.reg[df.reg$NozzleType == 1,])
model2 <- lm(SyrupLoss ~ Speed * Pressure + I(Speed^2) + I(Pressure^2), data=df.reg[df.reg$NozzleType == 2,])
model3 <- lm(SyrupLoss ~ Speed * Pressure + I(Speed^2) + I(Pressure^2), data=df.reg[df.reg$NozzleType == 3,])

print(coef(model1))
print(coef(model2))
print(coef(model3))

\noindent Next, we generate contour plots with the following code:

x <- seq(from=100, to=140, length.out=100)
y <- seq(from=10, to=20, length.out=100)
df0 <- expand.grid('Speed'=code.speed(x), 'Pressure'=code.pressure(y))

z1 <- matrix(predict(model1, df0), ncol=100)
z2 <- matrix(predict(model2, df0), ncol=100)
z3 <- matrix(predict(model3, df0), ncol=100)

par(mfrow=c(2,2))
contour(x, y, z1, main='Nozzle Type 1')
contour(x, y, z2, main='Nozzle Type 2')
contour(x, y, z3, main='Nozzle Type 3')
plot.new()  # This just uses up the last cell.

---

Example 9.2

This example illustrates the design of a $3^2$ in three blocks (aliasing the $AB^2$ effect with blocks):

model <- aov(Response ~ Error(Block) + A + B + A:B, data=Example9.2)
print(summary(model))

---

Figure 9.10

Constructing correlation plots is easiest to do using the corrplot package [@corrplot]. Here, I construct the correlation plot for the $2^{6-2}$ Resolution IV design of Figure 9.10:

library(FrF2)
library(corrplot)

X.d <- FrF2(nfactors=6, nruns=16)
X.m <- model.matrix(~.*., data=X.d)
X.corr <- t(X.m) %*% X.m / 16
corrplot(X.corr, method='color', tl.col="black", tl.srt=45, tl.cex=0.7)

---

Example 9.3 --- The Spin Coating Experiment

Recall the spin coating experiment from Section 8.7.2. Again we use stepwise regression.

model <- p.stepwise(
    lm(Thickness ~ 1, Table9.24),
    ~ (A + B + C + D + E + F)^2,
    alpha.to.enter=0.20,
    alpha.to.leave=0.20
)
print(summary(model))
print(anova(model))

---

Example 9.4

At present there's no easy way to construct $I$-optimal designs with blocks using R. We will construct an I-optimal design for six factors at levels $(-1, 0, 1)$. The design from the text has balanced blocks with three runs in each block, but the AlgDesign package does not support this constraint for $I$-optimal designs. We construct the $I$-optimal design (without blocking) as follows:

library(AlgDesign)

candidate.list <- expand.grid(
    'A'=c(-1, 0, 1),
    'B'=c(-1, 0, 1),
    'C'=c(-1, 0, 1),
    'D'=c(-1, 0, 1),
    'E'=c(-1, 0, 1),
    'F'=c(-1, 0, 1)
)

opt.i <- optFederov(~A+B+C+D+E+F, data=candidate.list, nTrials=12, criterion='I', approximate=FALSE, nRepeats=1000)

\noindent The optimal design discovered has I-criterion r opt.i$I and the following layout:

kable(opt.i$design) %>% kable_styling(bootstrap_options = "striped", full_width = F)

\noindent Compare this to the design in Table 9.36. The value of the $I$-criterion for this design is as follows:

eval.design(~ A + B + C + D + E + F, design=Table9.36, X=candidate.list)$I

---

References

ehassler/MontgomeryDAE documentation built on March 20, 2021, 7:08 p.m.