5 Mixture experiments
All designs discussed in the previous chapter 4 belong to the class of orthogonal designs in which all factors are varied independently. However, in many applications and especially in the chemical sciences there are cases where the factors cannot and should not be varied independently. If, say, there are three compounds x1,x2,x3, all in the range of 0-1 with 1 denoting the pure component, the design from blending these components must necessarily obey the mixture rule \(\sum_{i}x_{i}=1\). Strictly speaking, all problems in chemistry are mixture problem, but in most cases the solvent is treated as an inert “filler” not contributing to the mixture constraint and thereby allowing the remaining concentration factors to vary independently.
Due to the mixture constraints, mixture problems require special parametric models and corresponding to these models special designs. Mixture models and mixture designs will be introduced in the following sections and their use in R will be demonstrated with some examples borrowed from the literature.
A comprehensive overview on mixture experiments and a good source for further reading is the book (John A. Cornell 1990). (John Lawson, Cameron Willden 2016) gives an overview of the package mixexp, a collection of functions for creating and analyzing mixture designs in R.
5.1 Mixture models
Trying to identify the parametric OLS model \(\hat y = a_{0} + a_{1}\cdot x_{1} + a_{2}\cdot x_{2} + a_{3}\cdot x_{3}\) given the mixture constraint \(\sum_{i=1}^3x_{i}=1\) will fail because the mixture constraint renders one variable redundant. One way of dealing with the singularity is to drop one variable, say x3, and proceed with the reduced model \(\hat y =a_{0}^{'} + a_{1}\cdot x_{1}^{'} + a_{2}\cdot x_{2}^{'}\) and the auxiliary condition x3 = 1 - x1 - x2. This is called the slack variable approach with the slack variable, here x3, usually choosen to be the most inert component. Note that any other variable could have been excluded without affecting model statistics, as all slack variable models are statistically degenerated and will give the same predictions.
A more natural way of dealing with the linear constraints is to integrate the mixture contraints \(\sum_{i}x_{i}=1\) into the regression equation, leading to the so called Scheffe models, named after the scientist having first derived these equations. From the parametric model \[\hat y = a_{0} + a_{1}\cdot x_{1} + a_{2}\cdot x_{2} + a_{3}\cdot x_{3}\] and the mixture contraint \[x_{1}+x_{2}+x_{3}=1\] follows by substitution \[\hat y = a_{0}(x_{1}+x_{2}+x_{3}) + a_{1}\cdot x_{1} + a_{2}\cdot x_{2} + a_{3}\cdot x_{3} \Leftrightarrow \] \[\hat y = (a_{0} + a_{1})\cdot x_{1} + (a_{0} + a_{2})\cdot x_{2} + (a_{0} + a_{3})\cdot x_{3}\] the Scheffe main effect model \[\hat y = a_{1}^{'}\cdot x_{1} + a_{2}^{'}\cdot x_{2} + a_{3}^{'}\cdot x_{3}\] which extends for K mixture components to the Scheffe main effect model, equation (5.1)
From equation (5.1) and the mixture constraint follows the identity \(\hat {y}_{x_{i}=1} = a_{i}^{'}\), so the regression coefficients are the expected values \(\hat y_{i}\) of the pure components \(x_{i}\) with intermediate mixtures simply becoming weighted sums of the pure components. In the linear Scheffe model, the components do not interact and this might be an assumtion too tight for many applications thus motivating the quadratic Scheffe models.
Starting with the quadratic model \(\hat y= a_{0} + a_{1}\cdot x_{1} + a_{2}\cdot x_{2} + a_{3}\cdot x_{3} + a_{4}\cdot x_{1}^2 + a_{5}\cdot x_{2}^2 + a_{6}\cdot x_{3}^2 + a_{7} \cdot x_{1} \cdot x_{2} + a_{8} \cdot x_{1} \cdot x_{3}+ a_{9} \cdot x_{2} \cdot x_{3}\) by integrating the mixture constraint the quadratic Scheffe model is obtained \[\hat y= a_{1}^{'}\cdot x_{1} + a_{2}^{'}\cdot x_{2} + a_{3}^{'}\cdot x_{3} + a_{4}^{'} \cdot x_{1} \cdot x_{2} + a_{5}^{'} \cdot x_{1} \cdot x_{3} + a_{6}^{'} \cdot x_{2} \cdot x_{3}\] with the general quadratic Scheffe model for K components given by equation (5.2)
For any binary mixture xi=0.5, xj=0.5 follows from equation (5.2) the equality \[\hat y_{x_{i}=0.5;x_{j}=0.5} = \frac {a_{i} + a_{j}} {2} + \frac {a_{ij}} {4} \] which is depicted in figure 5.1 graphically for two binary mixtures with aij>0 (synergism) and aij<0 (antagonism). The straight line indicates the linear Scheffe model with aij=0.
Figure 5.1: Response plot of a quadratic Scheffe model for a binary mixture x1,x2 with synergistic and antagonistic effects labelled red and blue, respectively
Accounting for higher order effects such as ternary interactions can be achieved by augmenting the quadratic Scheffe model with higher order interaction terms as in the full cubic Scheffe model, see equation (5.3)
\[\begin{equation} y = \sum_{i=1}^{K}a_{i}\cdot x_{i} + \sum_{j>i}^{K} a_{ij} \cdot x_{i} \cdot x_{j} + \sum_{j>i}^{K} \delta_{ij} \cdot x_{i} \cdot x_{j} \cdot (x_{i} - x_{j}) + \sum_{k>j>i}^{K} a_{ijk} \cdot x_{i} \cdot x_{j} \cdot x_{k} + \epsilon \tag{5.3} \end{equation}\]Figure 5.2 shows a binary mixture plot of a full cubic model, which now can account for synergistic and antagonistic effects in mixture space. The cubic model underlying figure 5.2 is of the parametric form \(\hat y=a_{1} \cdot x_{1} + a_{2} \cdot x_{2} + a_{12} \cdot x_{1} \cdot x_{2} + \delta_{12} \cdot x_{1} \cdot x_{2} \cdot (x_{1}-x_{2})\) with parameters set to a1=50, a2=100, a12=0 and \(\delta_{12}=-200\).
Figure 5.2: Response plot of a cubic Scheffe model for a binary mixture x1,x2 with synergistic and antagonistic effects in mixture space
The full cubic model, formula (5.3), is very flexible due to the large number of regression parameters, hence it tends to overfit the data easily. Often the reduced (or special) cubic model, equation (5.4), is sufficiently complex to adequately describe the experimental situation.
\[\begin{equation} y = \sum_{i=1}^{K}a_{i}\cdot x_{i} + \sum_{j>i}^{K} a_{ij} \cdot x_{i} \cdot x_{j} + \sum_{k>j>i}^{K} a_{ijk} \cdot x_{i} \cdot x_{j} \cdot x_{k} + \epsilon \tag{5.4} \end{equation}\]The mixture models described above - linear, quadratic, full and the reduced (special) cubic model - are sufficient for modeling most mixture sytems appropriately and can be used as building blocks for modeling mixture-process designs, in which K mixture components x1, x2,…xK are combined with P process factors z1, z2,…zP. The first mixture-process model results from fully crossing the quadratic Scheffe model with a factorial (bilinear) process model, i.e. \[y =\Bigg ( \sum_{i=1}^{K}a_{i}\cdot x_{i} + \sum_{j>i}^{K} a_{ij} \cdot x_{i} \cdot x_{j} \Bigg ) \cdot \Bigg (\alpha_{0} + \sum_{l=1}^P \alpha_{l} \cdot z_{l} + \sum_{m>l}^R \alpha_{lm} \cdot z_{l}\cdot z_{m}\Bigg) + \epsilon\] and from product expansion follows directly the parametric model, equation (5.5)
\[\begin{align} y = & \sum_{i=1}^{K}b_{i}\cdot x_{i} + \sum_{j>i}^{K} b_{ij} \cdot x_{i} \cdot x_{j} + \sum_{i=1}^K \sum_{l=1}^P b_{il}^{'} \cdot x_{i} \cdot z_{l} + \sum_{j>i}^K \sum_{l=1}^P b_{ijl}^{'} \cdot x_{i} \cdot x_{j} \cdot z_{l} \ + \notag \\ & \sum_{i}^K \sum_{m>l}^P b_{ilm}^{'} \cdot x_{i} \cdot z_{l} \cdot z_{m} + \sum_{j>i}^K \sum_{m>l}^P b_{ijlm}^{'} \cdot x_{i} \cdot x_{j} \cdot z_{l} \cdot z_{m} + \epsilon \tag{5.5} \end{align}\]The number of regression coefficient in formula (5.5) as a function of K and P grows large rapidly making crossed designs of limited use for most applications. A more parsimonious model is obtained by specifically linking a quadratic Scheffe model with an RSM model leading to equation (5.6)
\[\begin{equation} y = \sum_{i=1}^{K}b_{i}\cdot x_{i} + \sum_{j>i}^{K} b_{ij} \cdot x_{i} \cdot x_{j} + \sum_{i=1}^K \sum_{l=1}^P b_{il}^{'} \cdot x_{i} \cdot z_{l} + \sum_{m>l}^P b_{lm}^{'} \cdot z_{l} \cdot z_{m}+ \sum_{l=1}^P b_{ll} \cdot z_{l}^2 + \epsilon \tag{5.6} \end{equation}\]The four mixture models, equations #1:(5.1), #2:(5.2), #3:(5.3), #4:(5.4) and the two mixture-process designs, fully and partially crossed mixture-process designs, equations #5:(5.5) and #6:(5.6), can be conveniently referenced by the model argument in the R-function mixexp::MixModel(…, model=#) with # denoting the model number in the above list.
When analyzing mixture experiments all principles of OLS model validation apply such as statistical testing and residual analysis as will be demonstrated with some example data provided with the R package mixexp from (R.H. Myers, D.C. Montgomery 2002). The aim of the experiment was to find an optimal blend of the components x1,x2,x3 maximizing the etching rate of silicon wafers (response erate). The design of this experiment was a partially replicated Simplex-Centroid Design (explained later) augmented by some interior points. The design maximally supports a reduced cubic model and will serve here as an example to demonstrate some functions from mixexp. Different from orthogonal models, it is usually not possible in mixture models to individually remove non-significant model terms as this will violate the mixture constraints and render the model singular. It is usually sufficient to choose between the models #1-#4 to appropriately describing the data within the model domain41.
The example data is shown in figure 5.3 as a ternary plots and supports either model #1, #2 or #4. Replicates are marked by circles and well balanced over the design space so as not to disturb the orthogonality of the design.
Figure 5.3: Ternary plot of the etching data with replicates marked by circles
Following the principle of parsimony, the data is first analysed with a quadratic Scheffe model and next with a reduced cubic model, and residuals are being compared in figure 5.4 (note that R2 in the lm-output of the no-intercept models are biased and should be ignored. R2 values reported by MixModel are correct and should be used instead).
library(mixexp)
data("etch")
summary(lm(erate ~ -1 + x1 + x2 + x3 +
x1:x2 + x1:x3 + x2:x3,
data = etch)) ##
## Call:
## lm(formula = erate ~ -1 + x1 + x2 + x3 + x1:x2 + x1:x3 + x2:x3,
## data = etch)
##
## Residuals:
## Min 1Q Median 3Q Max
## -157.68 -30.13 11.62 38.04 177.90
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## x1 534.64 83.35 6.414 0.000206 ***
## x2 329.16 83.35 3.949 0.004241 **
## x3 252.73 83.35 3.032 0.016254 *
## x1:x2 1343.11 469.58 2.860 0.021145 *
## x1:x3 644.53 469.58 1.373 0.207133
## x2:x3 711.68 469.58 1.516 0.168101
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 120.3 on 8 degrees of freedom
## Multiple R-squared: 0.9721, Adjusted R-squared: 0.9511
## F-statistic: 46.38 on 6 and 8 DF, p-value: 8.741e-06# note: R2 is wrong because when there is no intercept
# then SS.model is not adjusted by the mean
# Normally, the intercept a0 ensures mean adjustment
(res1 <- MixModel(etch, "erate",
mixcomps = c("x1", "x2", "x3"), model = 2))##
## coefficients Std.err t.value Prob
## x1 534.6383 83.34866 6.414481 0.000205933
## x2 329.1622 83.34866 3.949220 0.004240658
## x3 252.7336 83.34866 3.032245 0.016253810
## x2:x1 1343.1061 469.58404 2.860204 0.021144784
## x3:x1 644.5346 469.58404 1.372565 0.207132664
## x2:x3 711.6775 469.58404 1.515549 0.168101200
##
## Residual standard error: 120.261 on 8 degrees of freedom
## Corrected Multiple R-squared: 0.7583617##
## Call:
## lm(formula = mixmodnI, data = frame)
##
## Coefficients:
## x1 x2 x3 x1:x2 x1:x3 x2:x3
## 534.6 329.2 252.7 1343.1 644.5 711.7# quadratic Scheffe R2 correct
(res2 <- MixModel(etch, "erate",
mixcomps = c("x1", "x2", "x3"), model = 4)) ##
## coefficients Std.err t.value Prob
## x1 550.199515 23.22446 23.69051468 6.067419e-08
## x2 344.723325 23.22446 14.84311192 1.509476e-06
## x3 268.294753 23.22446 11.55224716 8.203511e-06
## x2:x1 689.537037 146.51489 4.70625916 2.192427e-03
## x3:x1 -9.034392 146.51489 -0.06166193 9.525557e-01
## x2:x3 58.108466 146.51489 0.39660451 7.034720e-01
## x2:x3:x1 9243.333333 940.85336 9.82441444 2.404146e-05
##
## Residual standard error: 33.43177 on 7 degrees of freedom
## Corrected Multiple R-squared: 0.9836603##
## Call:
## lm(formula = mixmodnI, data = frame)
##
## Coefficients:
## x1 x2 x3 x1:x2 x1:x3 x2:x3 x1:x2:x3
## 550.200 344.723 268.295 689.537 -9.034 58.108 9243.333# reduced cubic
kappa(res1)## [1] 9.353397kappa(res2)## [1] 57.53584par(mfrow=c(1,2))
colcol <- rep("black")
colcol[duplicated(etch[,1:3])]<- "red" #color label replicates
par(mfrow=c(1,2))
plot(predict(res1), rstudent(res1),ylim=c(-4,3),
xlab=expression( paste("predicted response ", hat(y)) ),
ylab="studentized residuals",type="n",
main=expression(paste("(",frac(paste("y - ",
hat(y)),sigma), ") ~ ", hat(y) )) )
text(predict(res1), rstudent(res1),label=1:nrow(etch), col=colcol)
abline(h=c(-2,0,2),lty=c(2,1,2))
grid()
# check residuals of reduced cubic model
plot(predict(res2), rstudent(res2),ylim=c(-4,3),
xlab=expression( paste("predicted response ", hat(y)) ),
ylab="studentized residuals",type="n",
main=expression(paste("(",frac(paste("y - ",
hat(y)),sigma), ") ~ ", hat(y) )) )
text(predict(res2), rstudent(res2),label=1:nrow(etch),col=colcol)
abline(h=c(-2,0,2),lty=c(2,1,2))
grid()
Figure 5.4: Studentized residuals versus model predictions from quadratic (left panel) and reduced cubic Scheffe models of the etching data. The replicate pairs (3,10), (2,9), (8,1) and (7,11) are better described within relication error by the cubic model (right panel) compared with the quadratic model (left panel).
The corrected R2 value of the quadratic Scheffe model, R2=0.76 is small compared with the cubic model, R2=0.98, indicating that the term x1:x2:x3 contributes substantially to the variance of the response. Similar conclusions can be drawn from figure 5.4 which reveals the quadratic Scheffe model to be biased by, e.g., not appropriately describing the replicate pair #7,#11. In the cubic Scheffe model the ternary term x1:x2:x3 is found significant and the residuals of the cubic model are now in better accordance with the normality assumption of the error term \(\epsilon\).
The cubic model can be visualized with a ternary contour plot shown in figure 5.5. The maximum lies close to the centroid (x1=0.33,x2=0.33,x3=0.33)42.
library(mixexp)
data("etch")
invisible(capture.output(res <- MixModel(etch, "erate",
mixcomps = c("x1", "x2", "x3"), model = 4)) )
# plot cubic mixture model as ternary contour plot
ModelPlot(model = res,
dimensions = list(x1 = "x1", x2 = "x2", x3 = "x3"),
contour = TRUE, cuts = 6, fill = TRUE)
Figure 5.5: Contour plot of the reduced cubic model of the etching data
Contour plots can also be used for dimensions higher than K=3 by conditioning (slicing) the plots on variables set constant. However, this can become confusing when there are many slices to compare, and therefore another popular, and less confusing way of visualizing high dimensional mixture models are the Cox effect plots. First, the Cox direction is defined as the line that connects the centroid point with the vertices (see figure 5.6) of the mixture space. When moving in the Cox direction of, say, x1 the ratio of the remaining mixture components stays constant, here \(\frac {x_{2}}{x_{3}}=1\), and this is next to equivalent to an “independent” effect of x1 in mixture space. The Cox effects are the traces obtained from sclicing the response surface along the Cox directions. Figure 5.6 shows the Cox direction in three dimensional mixture space along with a Cox effect plot of the reduced cubic model of the etching data.
rm(list=ls())
library(mixexp)
library(Ternary)
data(etch)
invisible(capture.output(res <- MixModel(etch, "erate",
mixcomps = c("x1", "x2", "x3"), model = 4)) )
# plot cubic mixture model as ternary contour plot
par(mfrow=c(1,2))
TernaryPlot(atip="x1", btip="x2", ctip="x3", axis.cex=1.2,lab.cex = 1.3,
axis.labels=seq(0,1,0.1))
TernaryArrows(c(0,0.5,0.5),c(1,0,0),col="red",length=0.15, lwd=2)
TernaryArrows(c(0.5,0,0.5),c(0,1,0),col="red",length=0.15, lwd=2)
TernaryArrows(c(0.5,0.5,0),c(0,0,1),col="red",length=0.15, lwd=2)
AddToTernary(points, list(x1=rep(0.33,3),x2=rep(0.33,3),
x3=rep(0.33,3)),pch=16,
col="black",cex=1.5)
title(main="Cox direction",cex.main=1.5)
ModelEff(nfac = 3, mod = 4, dir = 2, ufunc = res,
dimensions = list("x1", "x2", "x3"))
Figure 5.6: Cox direction in 3D mixture space (left panel) and effect plot along the Cox direction of the reduced cubic etching model (right panel)
5.2 Mixture designs
The previous section introduced six different parametric models for modeling mixture and mixture process systems. Associated with each model are mixture and mixture process designs optimally supporting the chosen model. Models and designs are two sides of the same coin: The parametric model determines the design and the design determines and limits the model which can be estimated from the data.
Mixture designs can be broadly divided into regular designs in which the design covers the whole mixture space (\(0 \leq x_{i} \leq 1\)) and irregular (or constrained) designs in which only a subspace, \(LB_{i}(>0) \leq x_{i} \leq UB_{i}(<1)\), of the whole mixture space is populated. Figure 5.7 shows examples of a regular and irregular mixture design in three dimensions.
Figure 5.7: Examples of regular (left panel) and irregular (constrained) designs (right panel)
Linear and quadratic Scheffe models (model #1, #2) can be estimated with (regular) Simplex Lattice Designs (SLD) which are shown for K=3 in figure 5.8 and created with the mixexp R-function SLD(fac,lev). The agument fac denotes the number of factors and lev specifies the number of levels with xi>0.
Figure 5.8: Simplex Lattice Designs for estimating linear (SLD(3,1), left panel) and quadratic Scheffe models (SLD(3,2), right panel)
The pure components (the vertex points in figure 5.8) are sufficient for estimating the linear Scheffe model. The quadratic Scheffe model is supported by augmenting the former design with all binary mixtures at the center of the three edges. Both designs are saturated43 and need therefore additional runs for estimating the error variance and to check for the presence of higher order effects.
The reduced cubic model (model #4) can be estimated with Simplex Centroid Designs (SCD). Again, the SCD is saturated and need additional runs for estimating error variance and Lack of Fit.
Finally, three level Simplex Lattice Designs can be used for estimating the full cubic model (model #3). All designs can be created conveniently with the functions mixexp::SLD() and mixexp::SCD() and examples for K=3 of the former and the latter are depicted in the figures 5.8 and 5.9.
Figure 5.9: Simplex Centroid Design SCD(3) (left panel) for reduced cubic model and three level Simplex Lattice Design SLD(3,3) (right panel) for the full cubic Scheffe model
Classical mixture models and designs are rigid in terms of model structure and design size similar to what has been said about classical orthogonal designs. A second order Simplex-Lattice Design for instance will render all two-way interaction estimable, irrespective of whether they are of interest or physically possible. Model structure and design size can be specified more flexibly with optimal designs as shown with the following R-code. Here the aim was to parsimoniously identify a non-standard Scheffe model, in R-notation ~ -1 + x1 + x2 + x3 + x1:x2 + x1:x2:x3, with six runs from a candidate set of 60 regular mixture candidates. Figure 5.10 depicts the D-optimal design points in mixture space comprising pure components (referring to main effects x1,x2,x3), the binary mixture x1,x2 (for the interaction term x1:x2) and two points close to the centroid for identifying the three way interaction x1:x2:x3.
rm(list=ls())
library(AlgDesign)
library(Ternary)
set.seed(1234)
xx <- expand.grid(x1=seq(0,1,0.1),
x2=seq(0,1,0.1),
x3=seq(0,1,0.1) ) # grid design
x <- xx[apply(xx,1,sum)==1,] # sum_x.i=1 only
mix <- optFederov(~ -1 + x1 + x2 + x3 +
x1:x2 + x1:x2:x3, data=x,nTrials=6,
nRepeats=1000)$design
TernaryPlot(atip="x1", btip="x2", ctip="x3",
axis.labels=seq(0,1,0.1), lab.cex=1)
AddToTernary(points, mix,pch=16,col="red",cex=1)
Figure 5.10: D-optimal mixture design of the model ~ -1 + x1 + x2 + x3 + x1:x2 + x1:x2:x3 written in R-formula notation
Designs for the fully crossed mixture-process model (#5), equation (5.5), results from crossing mixture and process designs. Let fac be a full or fractional factorial design with N1 rows and P columns and mix a mixture design of dimension (N2)x(K), the Cartesian product of dimension (N1 * N2)x(P+K) supports the fully crossed mixture-process design as the following R-code shows for K=3 and P=3, respectively
rm(list=ls())
set.seed(1234)
fac <- FrF2::FrF2(8,3,randomize=F) # 2^3 full factorial design; N=8
mix <- mixexp::SLD(3,3) # SLD design N=10
# 3vertex+3x2edges+centroid
fac <- data.frame(sapply(fac,function(x)
as.numeric(as.character(x))))
# convert factor to numerical
colnames(fac) <- paste("z",1:3,sep="")
x <- merge(mix,fac) # Cartesian product
dim(x)## [1] 80 6x$y <- rnorm(nrow(x))
(res <- mixexp::MixModel(x, "y",
mixcomps = c("x1", "x2", "x3"),
procvars = c("z1", "z2", "z3"),
model = 5)) # Fully crossed mixture-process model #5##
## coefficients Std.err t.value Prob
## x1 0.05990870 0.3082853 0.19432872 0.84695368
## x2 -0.57818736 0.3082853 -1.87549413 0.06842376
## x3 -0.20611270 0.3082853 -0.66857769 0.50780621
## x2:x1 0.05555502 1.3647251 0.04070785 0.96774194
## x3:x1 -1.00098144 1.3647251 -0.73346746 0.46777458
## x2:x3 0.85963536 1.3647251 0.62989636 0.53253190
## x1:z1 -0.17240678 0.3082853 -0.55924415 0.57927649
## x1:z2 0.48903461 0.3082853 1.58630507 0.12095767
## x1:z3 0.12188272 0.3082853 0.39535684 0.69479036
## x2:z1 0.22176294 0.3082853 0.71934310 0.47633042
## x2:z2 0.26971744 0.3082853 0.87489542 0.38712724
## x2:z3 -0.03665631 0.3082853 -0.11890384 0.90597790
## x3:z1 0.18778266 0.3082853 0.60911964 0.54606984
## x3:z2 0.04506181 0.3082853 0.14616916 0.88456052
## x3:z3 -0.12555989 0.3082853 -0.40728467 0.68608424
## x2:x1:z1 -1.82763803 1.3647251 -1.33919867 0.18846181
## x2:x1:z2 -0.75822632 1.3647251 -0.55558905 0.58174831
## x2:x1:z3 0.06996290 1.3647251 0.05126520 0.95938274
## x3:x1:z1 0.13423785 1.3647251 0.09836256 0.92216144
## x3:x1:z2 -0.90892644 1.3647251 -0.66601431 0.50942517
## x3:x1:z3 1.30367180 1.3647251 0.95526330 0.34548243
## x2:x3:z1 -0.72218119 1.3647251 -0.52917705 0.59976052
## x2:x3:z2 -1.33291497 1.3647251 -0.97669119 0.33489917
## x2:x3:z3 1.88583297 1.3647251 1.38184092 0.17509079
## x1:z1:z2 0.26702234 0.3082853 0.86615321 0.39184221
## x1:z1:z3 -0.58110294 0.3082853 -1.88495153 0.06710090
## x1:z2:z3 -0.14089678 0.3082853 -0.45703367 0.65024895
## x2:z1:z2 -0.13425104 0.3082853 -0.43547654 0.66568022
## x2:z1:z3 0.01472591 0.3082853 0.04776715 0.96215199
## x2:z2:z3 -0.17024018 0.3082853 -0.55221625 0.58403377
## x3:z1:z2 -0.30992865 0.3082853 -1.00533047 0.32109681
## x3:z1:z3 -0.69175568 0.3082853 -2.24388117 0.03074588
## x3:z2:z3 0.45376014 0.3082853 1.47188357 0.14928650
## x2:x1:z1:z2 -0.45968552 1.3647251 -0.33683378 0.73809548
## x2:x1:z1:z3 0.93480706 1.3647251 0.68497829 0.49751482
## x2:x1:z2:z3 2.83909884 1.3647251 2.08034487 0.04428803
## x3:x1:z1:z2 -0.64270562 1.3647251 -0.47094146 0.64037490
## x3:x1:z1:z3 3.11414546 1.3647251 2.28188481 0.02818482
## x3:x1:z2:z3 2.41847375 1.3647251 1.77213255 0.08439084
## x2:x3:z1:z2 -1.50904885 1.3647251 -1.10575299 0.27578515
## x2:x3:z1:z3 3.01676321 1.3647251 2.21052813 0.03316260
## x2:x3:z2:z3 -1.35284331 1.3647251 -0.99129365 0.32781272
##
## Residual standard error: 0.926512 on 38 degrees of freedom
## Corrected Multiple R-squared: 0.5948409##
## Call:
## lm(formula = mixmod, data = frame)
##
## Coefficients:
## x1 x2 x3 x2:x1 x3:x1 x2:x3 x1:z1 x1:z2
## 0.05991 -0.57819 -0.20611 0.05556 -1.00098 0.85964 -0.17241 0.48903
## x1:z3 x2:z1 x2:z2 x2:z3 x3:z1 x3:z2 x3:z3 x2:x1:z1
## 0.12188 0.22176 0.26972 -0.03666 0.18778 0.04506 -0.12556 -1.82764
## x2:x1:z2 x2:x1:z3 x3:x1:z1 x3:x1:z2 x3:x1:z3 x2:x3:z1 x2:x3:z2 x2:x3:z3
## -0.75823 0.06996 0.13424 -0.90893 1.30367 -0.72218 -1.33291 1.88583
## x1:z1:z2 x1:z1:z3 x1:z2:z3 x2:z1:z2 x2:z1:z3 x2:z2:z3 x3:z1:z2 x3:z1:z3
## 0.26702 -0.58110 -0.14090 -0.13425 0.01473 -0.17024 -0.30993 -0.69176
## x3:z2:z3 x2:x1:z1:z2 x2:x1:z1:z3 x2:x1:z2:z3 x3:x1:z1:z2 x3:x1:z1:z3 x3:x1:z2:z3 x2:x3:z1:z2
## 0.45376 -0.45969 0.93481 2.83910 -0.64271 3.11415 2.41847 -1.50905
## x2:x3:z1:z3 x2:x3:z2:z3
## 3.01676 -1.35284kappa(res)## [1] 15.35277Fully crossed designs and models are very bulky and hardly feasible in ordinary laboratory setups, and the partially crossed designs, model #6 (equation (5.6))), might be a more realistic alternative for most projects. Partially crossed and parsimonious mixture-process designs can be generated by D-optimally selecting points from a fully crossed candidate set with AlgDesign::optFederov() as shown by the following code
rm(list=ls())
set.seed(1234)
rsm <- rsm::bbd(3, n0=1,randomize=F)[,-(1:2)] # Box-Behnken design;N=13
sld <- mixexp::SLD(3,3) # cubic Simplex-lattice design;N=10
colnames(rsm) <- paste("z",1:3,sep="")
candidate <- merge(sld,rsm) # Cartesian product dim=(130 x 6)
mix <- AlgDesign::optFederov(~ -1 + (x1+x2+x3)^3 +
x1:z1 + x1:z2 + x1:z3 +
x2:z1 + x2:z2 + x2:z3 +
x3:z1 + x3:z2 + x3:z3 +
z1:z2 + z1:z3 + z2:z3 +
I(z1^2) + I(z2^2) + I(z3^2) ,
data=candidate,nTrials=30,
criterion="D", nRepeats=1000)$design
mix$y <- rnorm(nrow(mix))
(res <- mixexp::MixModel(mix, "y",
mixcomps = c("x1", "x2", "x3"),
procvars = c("z1", "z2", "z3"),
model = 6))## [1] "Warning, when using Model 6 the design in the process variables allow fitting the full quadratic model"
##
## coefficients Std.err t.value Prob
## x1 1.1568064 1.1661825 0.9919600 0.3471356
## x2 1.4896200 1.0769369 1.3832007 0.1999491
## x3 0.5302152 1.2014390 0.4413168 0.6693956
## I(z1^2) -0.4855809 0.5771805 -0.8412983 0.4219706
## I(z2^2) -0.2406845 0.5947270 -0.4046974 0.6951494
## I(z3^2) -0.2628178 0.6226630 -0.4220868 0.6828656
## x1:x2 -3.0211103 3.3953217 -0.8897862 0.3967412
## x1:x3 0.3875330 3.5732353 0.1084544 0.9160147
## x2:x3 -2.6004735 3.4484189 -0.7541060 0.4700549
## x1:z1 -0.5310433 0.5583275 -0.9511321 0.3663794
## x1:z2 -0.4529325 0.5666571 -0.7993062 0.4446950
## x1:z3 -0.2559378 0.4999943 -0.5118815 0.6210569
## x2:z1 -0.3321763 0.5764581 -0.5762367 0.5785786
## x2:z2 -0.3683906 0.5352710 -0.6882320 0.5086556
## x2:z3 -0.2032059 0.5336694 -0.3807712 0.7122050
## x3:z1 0.4478906 0.6118973 0.7319702 0.4828108
## x3:z2 0.4189375 0.6242892 0.6710632 0.5190298
## x3:z3 -0.8799695 0.5389427 -1.6327700 0.1369515
## z1:z2 0.1701386 0.4439679 0.3832228 0.7104494
## z1:z3 0.4408563 0.4035429 1.0924644 0.3030079
## z2:z3 0.1727396 0.4084885 0.4228750 0.6823111
##
## Residual standard error: 1.188226 on 9 degrees of freedom
## Corrected Multiple R-squared: 0.5440434##
## Call:
## lm(formula = mixmod, data = frame)
##
## Coefficients:
## x1 x2 x3 I(z1^2) I(z2^2) I(z3^2) x1:x2 x1:x3 x2:x3 x1:z1 x1:z2 x1:z3
## 1.1568 1.4896 0.5302 -0.4856 -0.2407 -0.2628 -3.0211 0.3875 -2.6005 -0.5310 -0.4529 -0.2559
## x2:z1 x2:z2 x2:z3 x3:z1 x3:z2 x3:z3 z1:z2 z1:z3 z2:z3
## -0.3322 -0.3684 -0.2032 0.4479 0.4189 -0.8800 0.1701 0.4409 0.1727kappa(res)## [1] 24.00369Constrained or irregular designs can be constructed with the extreme vertices algorithm implemented in the function Xvert() from the package mixexp by specifying dimension and constrained region \((LB_{i} \ge 0) \leq x_{i} \leq UB_{i} \le 1\) with LBi and UBi denoting lower and upper bounds of xi, respectively. Consistency requires the sum of the upper bounds to exceed 1, that is \(\sum_{i} UB_{i}> 1\).
Regular and irregular designs are often sparse in supporting higher order terms. Additional interior points can be created with the function Fillv() by averaging and thereby creating midpoints between all possible pairs.
Conversely, if a constrained design has too many runs given project objectives, it can be reduced by D-optimally selecting a subset of design points with the R-function AlgDesign::optFederov().
Augmentation of a constrained design with interior points is demonstrated with the following R-code, and figure 5.11 shows the design prior and after augmentation.
rm(list=ls())
x1 <- mixexp::Xvert(3,uc=c(.3,.3,1),lc=c(0,0,0),ndm=1,plot=F)[,-4]
x2 <- mixexp::Fillv(3,x1)
par(mfrow=c(1,2))
TernaryPlot(atip="x1", btip="x2", ctip="x3",axis.cex=1.2,lab.cex = 1.5,
axis.labels=seq(0,1,0.1))
AddToTernary(points, x1,pch=16,col="red",cex=1.5)
TernaryPlot(atip="x1", btip="x2", ctip="x3",axis.cex=1.2,lab.cex = 1.5,
axis.labels=seq(0,1,0.1))
AddToTernary(points, x2,pch=16,col="red",cex=1.5)
Figure 5.11: Constrained design \(x_{1} < 0.3; x_{2} < 0.3; \sum_{i} x_{i} = 1\) (left panel) and corresponding design augmented by interior points (right panel);
Table 5.1 gives an overview of the so far discussed mixture models and designs combined with short hints on how to create these designs in R.
| # | Model | Design | R-function |
|---|---|---|---|
| 1 | linear Scheffe | one-level Simplex-lattice | SLD(K,1) |
| 2 | quadratic Scheffe | two-level Simplex-lattice | SLD(K,2) |
| 3 | full cubic Scheffe | three-level Simplex-lattice | SLD(K,3) |
| 4 | special cubic Scheffe | Simplex-centroid | SCD(K) |
| 5 | fully crossed mix-process | Cartesian(mix*fac) | merge(mix,fac) |
| 6 | partially crossed mix-process | D-optimal from (mix*rsm) | optFederov(merge(mix,rsm)) |
| 7 | constrained mix design | extreme vertices | Xvert(K,LB,UB) |
Finally and similar to orthogonal designs, mixture designs can be blocked if there are nuisance parameters such as different raw material batches, technicians or lab equipment. Although there are a couple of optimally blocked mixture designs derived from first principles (for details see (John A. Cornell 1990)), algorithmic blocking with the R-function AlgDesign::optBlock() will suffice for most practical purposes. The following example code will block a third order Simplex-lattic design with ten runs into two blocks and depict both blocks in figure 5.12. In addition the example code shows a few work-arounds to circumvent a couple of problems not covered by the standard functionality. First, define the variable xij=xi-xj for specifying the cubic term xi:xj:(xi-xj) in the formula expression. Second, define a new variable according to block2=as.integer(block==“B2”) as a unique block indicator. This is required because lm() by default drops the first level of a factor (here B1 in block={B1,B2}) when the model has an intercept. However, when there is no intercept both blockB1 and blockB2 are kept in the model matrix thus rendering the model singular because of blockB1+blockB2=1 and x1+x2+x3=1. This problem is circumvented with the new block variable block2=as.integer(block==“B2”). Then the location parameters of the full cubic lm-model agree well with the true model despite the fact that no term is found significant. However, this is not surpsising as there is only 1 degree of freedom rendering the test statistics very conservative.
rm(list=ls())
set.seed(12345)
x <- mixexp::SLD(3,3)
x$x12 <- x$x1-x$x2 # first workaround
x$x13 <- x$x1-x$x3
x$x23 <- x$x2-x$x3
res <- AlgDesign::optBlock(~ -1 +(x1+x2+x3)^3 +
x1:x2:x12 + x1:x3:x13 + x2:x3:x23 , x,
c(6,6),criterion="Dpc", nRepeats=1000)
x <- res$design
x$block <- factor(c(rep("B1",6),rep("B2",6) ))
x$y <- 20*(x$block=="B2") + 3*x$x1 + 2*x$x2 +
5*x$x3 + 10*x$x1*x$x2*x$x12 + rnorm(nrow(x),0,0.1)
x$block2 <- as.integer(x$block=="B2") # second workaround
# that won't work as blockB1+blockB2=1
# => confounded with x1+x2+x3=1
summary(res <- (lm(y ~ -1 + block +(x1+x2+x3)^3 +
x1:x2:x12 + x1:x3:x13 + x2:x3:x23 ,data=x)))##
## Call:
## lm(formula = y ~ -1 + block + (x1 + x2 + x3)^3 + x1:x2:x12 +
## x1:x3:x13 + x2:x3:x23, data = x)
##
## Residuals:
## 1 2 4 5 6 7 3 51 61 8
## 1.251e-17 -1.431e-17 2.339e-17 1.079e-01 -1.079e-01 -1.002e-17 -7.683e-18 -1.079e-01 1.079e-01 -2.503e-17
## 9 10
## 6.674e-18 -1.067e-17
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## blockB1 5.0703 0.3052 16.613 0.0383 *
## blockB2 24.9924 0.2158 115.811 0.0055 **
## x1 -2.0946 0.3738 -5.604 0.1124
## x2 -2.9911 0.3738 -8.002 0.0791 .
## x3 NA NA NA NA
## x1:x2 -0.3166 1.0857 -0.292 0.8194
## x1:x3 0.3885 0.9403 0.413 0.7505
## x2:x3 -0.5591 0.9711 -0.576 0.6674
## x1:x2:x3 1.5672 5.7401 0.273 0.8303
## x1:x2:x12 10.8875 2.6148 4.164 0.1501
## x1:x3:x13 0.3372 1.9270 0.175 0.8897
## x2:x3:x23 -1.3537 2.3787 -0.569 0.6706
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2158 on 1 degrees of freedom
## Multiple R-squared: 1, Adjusted R-squared: 0.9998
## F-statistic: 6666 on 11 and 1 DF, p-value: 0.009553kappa(res) # 1.90344e+17 singular design## [1] 1.90344e+17## but this trick will help
x$block2 <- as.integer(x$block=="B2")
summary(res <- (lm(y ~ -1 + block2 +(x1+x2+x3)^3 +
x1:x2:x12 + x1:x3:x13 + x2:x3:x23 ,data=x)))##
## Call:
## lm(formula = y ~ -1 + block2 + (x1 + x2 + x3)^3 + x1:x2:x12 +
## x1:x3:x13 + x2:x3:x23, data = x)
##
## Residuals:
## 1 2 4 5 6 7 3 51 61 8
## -1.739e-17 2.315e-17 -1.655e-17 1.079e-01 -1.079e-01 9.307e-18 -7.520e-18 -1.079e-01 1.079e-01 -1.215e-18
## 9 10
## 6.300e-18 -3.478e-17
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## block2 19.9221 0.2158 92.316 0.0069 **
## x1 2.9757 0.2158 13.789 0.0461 *
## x2 2.0792 0.2158 9.635 0.0658 .
## x3 5.0703 0.3052 16.613 0.0383 *
## x1:x2 -0.3166 1.0857 -0.292 0.8194
## x1:x3 0.3885 0.9403 0.413 0.7505
## x2:x3 -0.5591 0.9711 -0.576 0.6674
## x1:x2:x3 1.5672 5.7401 0.273 0.8303
## x1:x2:x12 10.8875 2.6148 4.164 0.1501
## x1:x3:x13 0.3372 1.9270 0.175 0.8897
## x2:x3:x23 -1.3537 2.3787 -0.569 0.6706
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2158 on 1 degrees of freedom
## Multiple R-squared: 1, Adjusted R-squared: 0.9998
## F-statistic: 6666 on 11 and 1 DF, p-value: 0.009553kappa(res) # 47.5957## [1] 47.5957
Figure 5.12: Full cubic mixture design in two blocks B1 (left panel) and B2 (right panel)
Reference
John A. Cornell. 1990. Experiments with Mixtures. 2nd ed. Wiley, New York.
John Lawson, Cameron Willden. 2016. “Mixture Experiments in R Using Mixexp.” Journal of Statistical Software 72. https://www.jstatsoft.org/article/view/v072c02.
R.H. Myers, D.C. Montgomery. 2002. Response Surface Methodology: Process and Product Optimization Using Designed Experiments. 2nd ed. John Wiley & Sons.
More precisely it is sufficient to reduce the cubic model #4 to the quadratic model #2 if third order effects are not significant, or to reduce #2 to #1 if quadratic effects turn out not significant.↩
The precise location of the optimum remains vague in contour plots. With the aid of optimization, the subject of chapter 6, the maximum can be located to be \(x_{1}^*=0.41,x_{2}^*=0.34,x_{3}^*=0.25\) with an expected etching rate at this point of \(\hat y=832\).↩
The number of runs equal the number of regression. This will give a perfect fit with R2=1.↩