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

source('_format.R')

In Chapter 4 the text introduces designs that can accomodate investigations where complete randomization isn't feasible. These designs are useful when you want to estimate the effect of a factor and avoid bias due to other nuissance factors.

Example 4.1

Example 4.1 describes an experiment where batches of resin are extruded at different pressures and the proportion of defective grafts are recorded. The goal of the experiment is to maximize yield by minimizing the number of defects (Flicks). The structure of the data is presented below:

str(Table4.3)

\noindent In generating the ANOVA table we can either include BatchOfResin and ignore the $p$-value or we can mark it with Error to indicate each level is a block:

model <- aov(Flicks ~ Error(BatchOfResin) + ExtrusionPressure, data=Table4.3)
print(summary(model))

\noindent What if we hadn't recognized that this was a completely randomized block design? That would yield the following ANOVA table:

summary(aov(Flicks ~ ExtrusionPressure, data=Table4.3))

\noindent To estimate the means we turn to the emmeans package [@emmeans]:

library(emmeans)
model <- aov(Flicks ~ ExtrusionPressure + BatchOfResin, data=Table4.3)
exmodel <- emmeans(model, ~ ExtrusionPressure)
print(exmodel)
blmodel <- emmeans(model, ~ BatchOfResin)
print(blmodel)

\noindent We follow this with a boxplot for the defects:

boxplot(Flicks ~ ExtrusionPressure, data=Table4.3)

\noindent Consider the diagnostic plots for Example 4.1:

op <- par(mfrow=c(1,2))
qqnorm(resid(model), main='Normal Probability Plot')
qqline(resid(model), col='orange')
plot(x=predict(model), y=resid(model), main='Fitted vs. Residual')
par(op)

\noindent These plots do not show strong signs of any problem. The plot of residual versus each of the factors would be informative if some model misspecificaiton had occurred, but as in Figure 4.6, we find no strong evidence of problems from such a plot:

op <- par(mfrow=c(1,2))
plot(y=resid(model), x=as.numeric(Table4.3$ExtrusionPressure), main='Residual vs. Extrusion Pressure')
plot(y=resid(model), x=as.numeric(Table4.3$BatchOfResin), main='Residual vs. Catch')
par(op)

Let's assume the blocks are random effects, then we arrive at the following model:

library(lme4)
library(lmerTest)
model <- lmer(
    Flicks ~ -1 + ExtrusionPressure + (1|BatchOfResin),
    data=Table4.3,
    REML=TRUE
)
print(summary(model))
print(rand(model))

\noindent Confidence intervals on the variance components and the fixed effects can also be evaluated:

confint(model, oldNames=FALSE)

---

Example 4.2

In Example 4.2 an experimenter is studying the effects five different formulation of rocket propellant on the burn rate of that propellant (see Table 4.9 in the textbook). Formulations are constructed from batches that can only supply five formulations in all. Also, the formulations are prepared by different operators. The two levels of blocking --- the batch and the operator --- require the use of a Latin square design. First, let's look at the structure of the data:

str(Table4.9)

\noindent Next, let's conduct the ANOVA test to see if Formulation has an effect. Again, wrapping Batch and Operator in Error is optional and just prevents rendering $p$-values for them.

y <- Table4.9$BurnRate - 25
model <- aov(y ~ Error(Batch + Operator) + Formulation, data=Table4.9)
summary(model)

In this situation the design was provided, but you can use R to find Latin squares and Graeco-Latin squares. The package agricolae [@agricolae] contains functions for finding such designs. Constructing a Latin square can be done with design.lsd:

library(agricolae)

trt <- c("A", "B", "C", "D")
outdesign <- design.lsd(trt)
print(outdesign$sketch)

\noindent The agricolae package can also be used to construct Graeco-Latin squares with design.graeco and a host of other designs (balanced incomplete block designs, split plots, strip plots, etc.).

---

Example 4.3

In Example 4.3 the rocket propellant data is analyzed with an additional blocking factor: TestAssembly. The ANOVA table for this model is presented below:

y <- Table4.20$BurnRate - 25
model <- aov(y ~ Error(Batch + Operator + TestAssembly) + Formulation, data=Table4.20)
summary(model)

---

Section 4.4 --- Balanced Incomplete Block Design

The data in Table 4.22 in the textbook is from an experiment when an engineer tests four types of catalyst to optimize reaction time. Batches of material are only enough for three runs, so we must treat the batches as blocks. This is a balanced incomplete block design and is not orthogonal. Because of this, ReactionTime ~ Batch + Catalyst and ReactionTime ~ Catalyst + Batch give different results. Taking Batch before Catalyst adjusts Catalyst for the effect of Batch. This is what we want:

model <- aov(ReactionTime ~ Batch + Catalyst, data=Table4.22)
summary(model)

\noindent Compare this to the other order and we get the unadjusted treatment and adjusted block in Table 4.25:

model <- aov(ReactionTime ~ Catalyst + Batch, data=Table4.22)
summary(model)

\noindent Interestingly, if you wrap Batch with Error then R evaluates it first and the order does not matter.

Let's fit a model for parameter estimation. I'm going to use emmeans so we can't use the Error function:

library(emmeans)
model <- aov(ReactionTime ~ Batch + Catalyst, Table4.22)
print(summary(model))
camodel <- emmeans(model, ~Catalyst)
print(camodel)
bamodel <- emmeans(model, ~Batch)
print(bamodel)
print(pairs(camodel, adjust='tukey'))

---

References

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