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

source('_format.R')
set.seed(98477112)

A key feature of many response surface designs is splitting the residual error into a pure error term and a lack of fit term. The pure error term is invariant to the specified model as it compares the mean of runs at a specific location with the observed values. This allows testing model misspecification, but isn't supported by native R, so we will use the rsm package [@rsm] in this chapter. The rsm package also has facilities for design construction --- the vignette for the package is a good resource on the many features.

Example 11.1

An engineer investigating the effect of reaction time and reaction temperature on yield collected the data of Table 11.1. We'll use rsm's coded.data function to encode our data on the $[-1, 1]$ scale, and then fit a first order model using the FO helper in the formula. Note that we can use code2val to convert coded values to natural units.

library(rsm)

df <- Table11.1
cf <- coded.data(df, x1 ~ (Time - 35) / 5, x2 ~ (Temperature - 155) / 5)
model <- rsm(Yield ~ FO(x1, x2), data=cf)
summary(model)

\noindent The lack of fit is not significant so there's no need to add terms, but let's add the two-factor interaction term into the model using the TWI helper anyway.

model <- rsm(Yield ~ FO(x1, x2) + TWI(x1, x2), data=cf)
summary(model)

\noindent We could fit a second order model using the SO helper, but there is no present need.

---

Example 11.2

We conduct the test for lack of fit over the linear, linear with interaction, and the full quadratic model, for each of the responses. First, we begin with Yield. For Yield, the first order model shows lack of fit:

cf <- coded.data(Table11.6, x1 ~ (Time - 85) / 5, x2 ~ (Temperature - 175) / 5)
model <- rsm(Yield ~ FO(x1, x2), data=cf)
summary(model)$lof

\noindent The first order model with interactions shows lack of fit as well:

model <- rsm(Yield ~ FO(x1, x2) + TWI(x1, x2), data=cf)
summary(model)$lof

\noindent The second order model does not show significant lack of fit:

model <- rsm(Yield ~ SO(x1, x2), data=cf)
summary(model)$lof

\noindent Next, we examine Viscosity. Lack of fit is no longer significant at the second order model:

model.fo <- rsm(Viscosity ~ FO(x1, x2), data=cf)
model.in <- rsm(Viscosity ~ FO(x1, x2) + TWI(x1, x2), data=cf)
model.so <- rsm(Viscosity ~ SO(x1, x2), data=cf)

rbind(
    'First Order LOF'=summary(model.fo)$lof[3,],
    'First Order Interaction LOF'=summary(model.in)$lof[4,],
    'Second Order LOF'=summary(model.so)$lof[5,]
)

\noindent Finally, we examine MolecularWeight. Lack of fit is not significant for the first order model, so this is the model we shall use:

model.fo <- rsm(MolecularWeight ~ FO(x1, x2), data=cf)
model.in <- rsm(MolecularWeight ~ FO(x1, x2) + TWI(x1, x2), data=cf)
model.so <- rsm(MolecularWeight ~ SO(x1, x2), data=cf)

rbind(
    'First Order LOF'=summary(model.fo)$lof[3,],
    'First Order Interaction LOF'=summary(model.in)$lof[4,],
    'Second Order LOF'=summary(model.so)$lof[5,]
)

Constructing the response surface plots is pretty easy with rsm. Figure 11.11 gives the response surface for the second order model of Yield, which we present below:

 contour(model, ~ x1 + x2, image = TRUE)

---

Example 11.3

Note that the AlgDesign package uses the Federov exchange algorithm and chooses from a candidate set of points, whereas software like JMP and Design Expert use an efficient search algorithm that does not require enumeration of points. As of this time I'm unaware of an R package for constructing such designs, however the existing optimization routines in R are sufficient for such construction.

For both the $D$-optimal design in Table 11.13 and the $I$-optimal design in Table 11.14, a fraction of design space plot was created. The package Vdgraph [@Vdgraph] provides the ability to construct such plots. The following is the fraction of design space plot for the $D$-optimal design of Table 11.13:

library(Vdgraph)

FDSPlot(Table11.13[,c('X1','X2','X3','X4')], mod=2)

\noindent Next, we present the fraction of design space plot for the $I$-optimal design of Table 11.14:

FDSPlot(Table11.14[,c('X1','X2','X3','X4')], mod=2)

---

Example 11.4

Computational fluid dynamics (CFD) models of exhaust from a jet turbine was studied to determine the effect of location in the exhaust plume on temperature. The R packages DiceDesign [@DiceDesign] and DiceEval [@DiceEval] are explicitly for the design and analysis of computer experiments (https://www.jstatsoft.org/article/view/v065i11/v65i11.pdf). We can construct many space filling designs with DiceDesign, though the specific sphere-packing design from JMP is not supported in this package. Let's construct a Latin Hypercube design for the space $x \in (0.05, 0.095)$ and $R \in (0, 0.062)$:

library(DiceDesign)
df <- lhsDesign(10, 2)$design
colnames(df) <- c('x', 'R')
# Translate the points which are on [0,1] scale to the required scale
df[,'x'] <- df[,'x'] * 0.9 + 0.05
df[,'R'] <- df[,'R'] * 0.062
plot(df, main='LHS Design')

\noindent Let us now move to the design in Table 11.16. After fitting the model we will create a plot of the estimated mean response surface:

library(DiceEval)
model <- modelFit(Table11.16[,c('x','R')], Table11.16[,'Temperature'], type='Kriging')

\noindent The contour plot of the response surface is given below. Compare this to Table 11.17's Contour Profiler result:

x.range <- range(Table11.16[,'x'])
R.range <- range(Table11.16[,'R'])
x <- seq(from=x.range[1], to=x.range[2], length.out=100)
R <- seq(from=R.range[1], to=R.range[2], length.out=100)
z <- matrix(
    modelPredict(model, expand.grid('x'=x, 'R'=R)),
    nrow=length(x)
)
contour(x, R, z, main='Gaussian Process Mean Surface')

---

Example 11.5 --- A Three-Component Mixture

The yarn elongation experiment investigates how Polyethylene, Polystyrene, and Polypropylene affect yarn elongation in kilograms of force applied. A ${3,2}$ simplex lattice design was used to investigate this mixture experiment. The package mixexp [@mixexp] provides a way to construct mixture experiments with R (https://www.jstatsoft.org/index.php/jss/article/view/v072c02/v72c02.pdf), and does the package qualityTools [@qualityTools] ( https://CRAN.R-project.org/package=qualityTools). Analysis of mixture experiments can use the regular lm function:

Note: qualityTools was recently removed from CRAN. You can install the last hosted version by downloading it manually from https://cran.r-project.org/src/contrib/Archive/qualityTools/qualityTools_1.55.tar.gz.

model <- lm(Elongation ~ -1 + (Polyethylene + Polystyrene + Polypropylene)^2, data=Table11.19)
print(summary(model))

\noindent Here, I use the qualityTools package to create a visualization of the mixture model:

library(qualityTools)

Y <- Table11.19$Elongation
mdo <- mixDesign(3, 2, center=FALSE, axial=FALSE, randomize=FALSE, replicates=c(1,1,2,3))
response(mdo) <- Y
contourPlot3(A, B, C, Y, data=mdo, form='quadratic')

---

Example 11.6

Construction of a $D$-optimal constrained mixture experiment can be achieved using the Xvert function in the mixexp package to generate a host of candidate points, then using the AlgDesign package to select a $D$-optimal subset of these points. Analysis of the optimal experiment in Example 11.6 proceeds as follows.

We cannot use the rsm package because it has trouble with dropping the intercept. So instead we'll use a helper function in this package pe.summary. First, we need to code the units:

df <- Table11.20
for(s in c('Monomer', 'Crosslinker', 'Resin')){
    df[,s] <- df[,s] / 100
}

\noindent Next, we examine the hardness model:

library(mixexp)

model <- MixModel(df, "Hardness", c("Monomer", "Crosslinker", "Resin"), 2)

\noindent Fitting the ANOVA table we find no lack of fit:

model <- aov(Hardness ~ -1 + (Monomer + Crosslinker + Resin)^2, data=df)
pe.summary(model)

\noindent Next, we examine the Solids model:

model <- MixModel(df, "Solids", c("Monomer", "Crosslinker", "Resin"), 2)

\noindent Fitting the ANOVA table we find evidence of lack of fit:

model <- aov(Solids ~ -1 + (Monomer + Crosslinker + Resin)^2, data=df)
pe.summary(model)

\noindent Since the quadratic model leaves no pure error we can't estimate lack of fit, so we'll stop here.

---

References

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