Quality Tools in IMPROVE
Thomas Roth - thomas [at] r-qualityTools [dot] org
2017-04-02
- \(2^k\) Factorial Designs
- \(2^{k-p}\) Fractional Factorial Designs
- Replicated Designs and Center Points
- Multiple Responses
- Moving to a process setting with an expected higher yield
- Response Surface Designs
- Using desirabilities together with designed experiments
- Mixture Designs
- Taguchi Designs
- Session Info
Each process has a purpose. The effectiveness of a process can be expressed with the help of (quality) characteristics. Those characteristics can be denoted as the responses of a process. In order to attain the desired values for the responses certain settings need to be arranged for the process. Those settings refer to the input variables of the process. Working with designed experiments it is helpful to refer to the (black box) process model.
In general input variables can be distinguished into controllable and disturbance variables. Input variables that can be controlled and have an assumed effect on the responses are denoted as factors. Input variables that are not factors are either hard to change (e.g. the hydraulic fluid in a machine) or varying them does not make good economic sense (e. g. the temperature or humidity in a factory building). These hard-to-change factors are also called uncontrollable input variables. It is attempted to hold those variables constant. Disturbance variables affect the outcomes of a process by introducing noise such as small variations in the controllable and uncontrollable input variables which in turn migth lead to variations in the response variables despite identical factor settings in an experiment.
\(2^k\) Factorial Designs
In order to find more about this black box model one can come up with a \(2^k\) factorial design by using the method facDesign of the qualityTools package. As used in textbooks k denotes the number of factors. A design with k factors and 2 combinations per factor gives you \(2^k\) different factor combinations and thus what is called runs.
Suppose a process has 5 factors A, B, C, D and E. The yield (i. e. response) of the process is measured in percent. Three of the five factors are assumed by the engineers to be relevant to the yield of the process. These three factors are to be named Factor 1, Factor 2 and Factor 3 (A, B and C). The (unknown relations of the factors of the) process are simulated by the method simProc of the qualityTools package. Factor 1 is to be varied from 80 to 120, factor B from 120 to 140 and factor C from 1 to 2 . Low factor settings are assigned a -1 and high values a +1.
set.seed(1234)
fdo = facDesign(k = 3, centerCube = 4) #fdo - factorial design object
names(fdo) = c("Factor 1", "Factor 2", "Factor 3") #optional
lows(fdo) = c(80, 120, 1) #optional
highs(fdo) = c(120, 140, 2) #optional
summary(fdo) #information about the factorial design## Information about the factors:
##
## A B C
## low 80 120 1
## high 120 140 2
## name Factor 1 Factor 2 Factor 3
## unit
## type numeric numeric numeric
## -----------
## StandOrd RunOrder Block A B C y
## 4 4 1 1 1 1 -1 NA
## 3 3 2 1 -1 1 -1 NA
## 8 8 3 1 1 1 1 NA
## 2 2 4 1 1 -1 -1 NA
## 12 12 5 1 0 0 0 NA
## 11 11 6 1 0 0 0 NA
## 5 5 7 1 -1 -1 1 NA
## 10 10 8 1 0 0 0 NA
## 9 9 9 1 0 0 0 NA
## 6 6 10 1 1 -1 1 NA
## 7 7 11 1 -1 1 1 NA
## 1 1 12 1 -1 -1 -1 NAThe response of this fictional process is given by the simProc method of the qualityTools package. The yield for Factor 1, Factor 2 and Factor 3 taking values of 80, 120 and 1 can be calculated using
#set first value
yield = simProc(x1 = 120, x2 = 140, x3 = 2)Setting all the yield of this artificial black box process gives a very long line of R-Code.
yield = c(simProc(120, 140, 1),simProc(80,140, 1),simProc(120,140, 2),
simProc(120,120, 1),simProc(100,130, 1.5),simProc(100,130, 1.5),
simProc(80,120, 2),simProc(100,130, 1.5), simProc(100,130, 1.5),
simProc(120,120, 2),simProc(80,140, 2), simProc(80,120, 1))Assigning the yield to the factorial design can be done using the response method.
response(fdo) = yield #assign yield to the factorial design objectAnalyzing this design is quite easy using the methods effectPlot, interactionPlot, lm as well as wirePlot and contourPlot:
effectPlot(fdo, classic = TRUE)
effect plot for the factorial design
interactionPlot(fdo)
interaction plot for the factorial design
The factorial design in fdo can be handed without any further operations directly to the base lm method of R.
lm.1 = lm(yield ~ A*B*C, data = fdo)
summary(lm.1)##
## Call:
## lm(formula = yield ~ A * B * C, data = fdo)
##
## Residuals:
## 1 2 3 4 5 6 7 8
## -0.01338 -0.01338 -0.01338 -0.01338 -0.01338 -0.01338 -0.01338 -0.01338
## 9 10 11 12
## 0.02893 0.03227 0.02076 0.02509
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.298e-01 9.543e-03 24.077 1.77e-05 ***
## A 7.242e-02 1.169e-02 6.197 0.003449 **
## B 1.134e-01 1.169e-02 9.702 0.000632 ***
## C -7.619e-05 1.169e-02 -0.007 0.995111
## A:B 7.834e-02 1.169e-02 6.703 0.002578 **
## A:C 1.823e-03 1.169e-02 0.156 0.883627
## B:C -2.139e-04 1.169e-02 -0.018 0.986277
## A:B:C -2.735e-03 1.169e-02 -0.234 0.826460
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03306 on 4 degrees of freedom
## Multiple R-squared: 0.978, Adjusted R-squared: 0.9394
## F-statistic: 25.36 on 7 and 4 DF, p-value: 0.003685The effects of A and B as well as the interaction A:B are identified to be significant. A Pareto plot of the standardized effects visualizes these findings and can be created with the paretoPlot method of the qualityTools package. Another visualization technique commonly found is a normal plot using the normalPlot method of the qualityTools package.
par(mfrow = c(1,2))
paretoPlot(fdo)
normalPlot(fdo)## A B C A:B A:C B:C
## 6.19653517 9.70212971 -0.00651905 6.70299438 0.15594853 -0.01829860
## A:B:C
## -0.23401870
pareto plot of the standardized effects and normal plot of the coefficients
The relation between the factors A and B can be visualized as 3D representation in form of a wireframe or contour plot using the wirePlot and contourPlot method of the qualityTools package. Again, no further transformation of the data is needed!
par(mfrow = c(1,2))
wirePlot(A, B, yield, data = fdo)
contourPlot(A, B, yield, data = fdo)
response surface and contour plot
One question that arises is whether this linear fit adequately describes the process. In order to find out, one can simply compare values predicted in the center of the design (i.e. A=0, B=0 and C=0) with the values observed in the center of the design. This difference could also be tested using a specialized t-Test. For now, let’s assume this model is less wrong than others (i. e. we don’t know of any better model).
\(2^{k-p}\) Fractional Factorial Designs
Imagine testing 5 different factors in a \(2^k\) design giving you \(2^5 = 32\) runs. This is likely to be quite expensive if run on any machine, process or setting within production, research or a similar environment. Before dismissing the design, it’s advisable to reflect what this design is capable of in terms of what types of interactions it can estimate. The highest interaction in a \(2^5\) design is the interaction between the five factors ABCDE. This interaction, even if significant, is really hard to interpret, and likely to be non-existent. The same applies for interactions between four factors and some of the interactions between 3 factors which is why most of the time fractional factorial designs are considered in the first stages of experimentation.
A fractional factorial design is denoted \(2^{k-p}\) meaning k factors are tested in \(2^{k-p}\) runs. In a \(2^{5-1}\) design five factors are tested in 24 runs (hence p=1 additional factor is tested without further runs). This works by confounding interactions with additional factors. This section will elaborate on this idea with the help of the methods of the qualityTools package.
For fractional factorial designs the method fracDesign of the package can be used. The generators can be given in the same notation that is used in textbooks on this matter. For a \(2^{3-1}\) design (i.e. 3 factors that are to be tested in a \(2^2\) by confounding the third factor with the interaction between the first two factors) this would be given by the argument gen = "C = AB" meaning the interaction between A and B is to be confounded with the effect of a third factor C. The effect estimated for C is then confounded with the interaction AB; they cannot be separately estimated, hence C = AB (alias) or the alias of C is AB.
fdo.frac = fracDesign(k = 3, gen = "C = AB", centerCube = 4)In order to get more specific information about a design the summary method can be used. For this example you will see on the last part the identity I = ABC of the design. The identity I of a design is the left part of the generator multiplied by the generator. The resolution is the (character-) length of the shortest identity.
summary(fdo.frac)## Information about the factors:
##
## A B C
## low -1 -1 -1
## high 1 1 1
## name
## unit
## type numeric numeric numeric
## -----------
## StandOrd RunOrder Block A B C y
## 4 4 1 1 1 1 1 NA
## 3 3 2 1 -1 1 -1 NA
## 2 2 3 1 1 -1 -1 NA
## 5 5 4 1 0 0 0 NA
## 6 6 5 1 0 0 0 NA
## 7 7 6 1 0 0 0 NA
## 8 8 7 1 0 0 0 NA
## 1 1 8 1 -1 -1 1 NA
##
## ---------
##
## Defining relations:
## I = ABC Columns: 1 2 3
##
## Resolution: IIIThe following rules apply
\begin{align} I\times A &= A\\ A\times A &= I\\ A\times B &= B\times A \end{align}By multiplying A, B and C you will find all confounded effects or aliases. A more convenient way to get an overview of the alias structure of a factorial design is to call the method aliasTable or confounds of the qualityTools package. The latter gives a more human readable version of the first.
aliasTable(fdo.frac)## C AC BC ABC
## Identity 0 0 0 1
## A 0 0 1 0
## B 0 1 0 0
## AB 1 0 0 0Fractional factorial designs can be generated by assigning the appropriate generators. However, most of the time standard fractional factorial designs known as minimum aberration designs (cite{Box05}) will be used. Such a design can be chosen from predefined tables by using the method fracChoose of the qualityTools package and simply clicking onto the desired design.
fracChoose()
Choosing minimum aberration designs
##
## Choose a fractional factorial design by clicking into the appropriate field
## Waiting for your selection:## StandOrder RunOrder Block A B C y
## 3 3 1 1 -1 1 -1 NA
## 2 2 2 1 1 -1 -1 NA
## 4 4 3 1 1 1 1 NA
## 1 1 4 1 -1 -1 1 NAReplicated Designs and Center Points
A replicated design with additional center points can be created by using the replicates and centerCube argument.
fdo1 = facDesign(k = 3, centerCube = 2, replicates = 2)Multiple Responses
Once you have observed the response for the different factor combinations one can add one or more response vectors to the design with the response method of the qualityTools package . A second response to be named y2 is created with the help of random numbers.
set.seed(1234)
y2 = rnorm(12, mean = 20)
response(fdo) = data.frame(yield, y2) A 3D visualization is done with the help of the methods wirePlot and contourPlot of the qualityTools package with no need to first create arrays of values or the like. Simply specify the formula you would like to fit with e. g. form = "yield ~ A+B". Specifying this fit for response yield one can see that there’s actually no practical difference to the fit that included an interaction term.
par(mfrow = c(1,2))
wirePlot(A, B, yield, data = fdo, form = "yield~A+B+C+A*B")
contourPlot(A, B, y2, data = fdo, form = "y2~A+B+C+A*B")
wire plot with different formulas specified
Using the wirePlot and contourPlot methods of the qualityTools package settings of the other n-2 factors can be set using the factors argument. A wireplot with the third factor C on -1 an C = 1 can be created as follows:
par(mfrow = c(1,2))
wirePlot(A,B,y2, data = fdo, factors = list(C=-1), form = "y2~A*B*C")
wirePlot(A,B,y2, data = fdo, factors = list(C=1), form = "y2~A*B*C")
wire plot with formula and setting for factor C
If no formula is explicitly given the methods default to the full fit or the fit stored in the factorial design object fdo. Storing a fit can be done using the fits method of the qualityTools package and is especially useful when working with more than one response. Again lm can be used to analyze the fractional factorial designs.
fits(fdo) = lm(yield ~ A+B, data = fdo)
fits(fdo) = lm(y2 ~ A*B*C, data = fdo)
fits(fdo)## $yield
##
## Call:
## lm(formula = yield ~ A + B, data = fdo)
##
## Coefficients:
## (Intercept) A B
## 0.22975 0.07242 0.11339
##
##
## $y2
##
## Call:
## lm(formula = y2 ~ A * B * C, data = fdo)
##
## Coefficients:
## (Intercept) A B C A:B
## 19.5577 -0.1982 0.5608 0.4270 0.2175
## A:C B:C A:B:C
## 0.5098 -0.0428 0.2518Moving to a process setting with an expected higher yield
Since our process can be adequately modeled by a linear relationship the direction in which to go for an expected higher yield is easy to determine. A contour plot of factor A and B illustrate that we simply need to “step up the stairs”. The shortest way to get up these stairs can be figured out graphically or calculated using the steepAscent method of the qualityTools package.
sao =steepAscent(factors=c("A","B"),response="yield",data=fdo,steps=20)##
## Steepest Ascent for fdo
##
## Run Delta A.coded B.coded A.real B.real
## 1 1 0 0.0 0.000 100 130
## 2 2 1 0.2 0.313 104 133
## 3 3 2 0.4 0.626 108 136
## 4 4 3 0.6 0.939 112 139
## 5 5 4 0.8 1.253 116 143
## 6 6 5 1.0 1.566 120 146
## 7 7 6 1.2 1.879 124 149
## 8 8 7 1.4 2.192 128 152
## 9 9 8 1.6 2.505 132 155
## 10 10 9 1.8 2.818 136 158
## 11 11 10 2.0 3.131 140 161
## 12 12 11 2.2 3.445 144 164
## 13 13 12 2.4 3.758 148 168
## 14 14 13 2.6 4.071 152 171
## 15 15 14 2.8 4.384 156 174
## 16 16 15 3.0 4.697 160 177
## 17 17 16 3.2 5.010 164 180
## 18 18 17 3.4 5.323 168 183
## 19 19 18 3.6 5.637 172 186
## 20 20 19 3.8 5.950 176 189
## 21 21 20 4.0 6.263 180 193sao## Run Delta A.coded B.coded A.real B.real "yield"
## 1 1 0 0.0 0.0000000 100 130.0000 NA
## 2 2 1 0.2 0.3131469 104 133.1315 NA
## 3 3 2 0.4 0.6262939 108 136.2629 NA
## 4 4 3 0.6 0.9394408 112 139.3944 NA
## 5 5 4 0.8 1.2525877 116 142.5259 NA
## 6 6 5 1.0 1.5657346 120 145.6573 NA
## 7 7 6 1.2 1.8788816 124 148.7888 NA
## 8 8 7 1.4 2.1920285 128 151.9203 NA
## 9 9 8 1.6 2.5051754 132 155.0518 NA
## 10 10 9 1.8 2.8183223 136 158.1832 NA
## 11 11 10 2.0 3.1314693 140 161.3147 NA
## 12 12 11 2.2 3.4446162 144 164.4462 NA
## 13 13 12 2.4 3.7577631 148 167.5776 NA
## 14 14 13 2.6 4.0709100 152 170.7091 NA
## 15 15 14 2.8 4.3840570 156 173.8406 NA
## 16 16 15 3.0 4.6972039 160 176.9720 NA
## 17 17 16 3.2 5.0103508 164 180.1035 NA
## 18 18 17 3.4 5.3234977 168 183.2350 NA
## 19 19 18 3.6 5.6366447 172 186.3664 NA
## 20 20 19 3.8 5.9497916 176 189.4979 NA
## 21 21 20 4.0 6.2629385 180 192.6294 NASince we set the real values earlier using the highs and lows methods of the qualityTools package factors settings are displayed in coded as well as real values. Again the values of the response of sao (__s__teepest __a__scent __o__bject) can be set using the response method of the qualityTools package and then be plotted using the plot method. Of course one can easily use the base plot method itself. However for documentation purposes the plot method for a steepest ascent object might be more convenient.
predicted = simProc(sao[,5], sao[,6])
response(sao) = predicted
plot(sao, type = "b", col = 2)
predicted maximum at Delta = 11 (see sao)
At this point the step size was chosen quite small for illustration purposes.
Response Surface Designs
Not all relations are linear and thus in order to detect and model non-linear relationships sometimes more than two combinations per factor are needed. At the beginning all a black box might need is a \(2^k\) or \(2^{k-p}\) design. In order to find out whether a response surface design (i.e. a design with more than two combination per factors) is needed one can compare the expected value of one’s response variable(s) with the observed one(s) using centerpoints (i.e. the 0, 0, , 0 setting). The bigger the difference between observed and expected values, the more unlikely this difference is the result of random noise.
For now, let’s return to the initial simulated process. The project started off with a \(2^k\) design containing center points. Sticking to a linear model we used the steepAscent method of the qualityTools package to move to a better process region. The center of the new process region is defined by 144 and 165 in real values. This region is the start of a new design. Again one starts by using a factorial design
fdo2 = facDesign(k = 2, centerCube = 3) #create a factorial design with 2 factors and 3 centerpoints
fdo2 = randomize(fdo2, random.seed = 1234) #randomize the factoriald design
names(fdo2) = c("Factor 1", "Factor 2")
lows(fdo2) = c(134, 155)
highs(fdo2) = c(154, 175)
summary(fdo2)## Information about the factors:
##
## A B
## low 134 155
## high 154 175
## name Factor 1 Factor 2
## unit
## type numeric numeric
## -----------
## StandOrd RunOrder Block A B y
## 1 1 1 1 -1 -1 NA
## 6 6 2 1 0 0 NA
## 4 4 3 1 1 1 NA
## 2 2 4 1 1 -1 NA
## 5 5 5 1 0 0 NA
## 3 3 6 1 -1 1 NA
## 7 7 7 1 0 0 NAand the yield is calculated by using the simProc and assigned to the design with the help of the generic response method of the qualityTools package.
yield = c(simProc(134,155), simProc(144,165), simProc(154,175), simProc(154,155), simProc(144,165), simProc(134,175), simProc(144,165)) #simulate yield
response(fdo2) = yield #set yield as responseLooking at the residual graphics one will notice a substantial difference between expected and observed values (a test for lack of fit could of course be performed and will be significant). To come up with a model that describes the relationship one needs to add further points which are referred to as the star portion of the response surface design.
Adding the star portion is easily done using the starDesign method of the qualityTools package. By default the value of alpha is chosen so that both criteria, orthogonality and rotatability are approximately met. Simply call the starDesign method on the factorial design object fdo2. Calling rsdo (response surface design object) will show you the resulting response surface design. It should have a cube portion consisting of 4 runs, 3 center points in the cube portion, 4 axial and 3 center points in the star portion.
rsdo = starDesign(data = fdo2) #turn factorial design into response surface design adding a star portion
rsdo## StandOrd RunOrder Block A B yield
## 1 1 1 1 -1.000 -1.000 0.7554
## 6 6 2 1 0.000 0.000 0.7954
## 4 4 3 1 1.000 1.000 0.6727
## 2 2 4 1 1.000 -1.000 0.8005
## 5 5 5 1 0.000 0.000 NA
## 3 3 6 1 -1.000 1.000 NA
## 7 7 7 1 0.000 0.000 NA
## 8 8 8 2 -1.414 0.000 NA
## 9 9 9 2 1.414 0.000 NA
## 10 10 10 2 0.000 -1.414 NA
## 11 11 11 2 0.000 1.414 NA
## 12 12 12 2 0.000 0.000 NA
## 13 13 13 2 0.000 0.000 NA
## 14 14 14 2 0.000 0.000 NAUsing the star method of the qualityTools package one can easily assemble designs sequentially. This sequential strategy saves resources since compared to starting off with a response surface design from the very beginning, the star portion is only run if really needed. The yields for the process are still given by the simProc method of the qualityTools package.
yield2 = c(yield, simProc(130,165), simProc(158,165), simProc(144,151), simProc(144,179), simProc(144,165), simProc(144,165), simProc(144,165))
response(rsdo) = yield2A full quadratic model is fitted using the lm method
lm.3 = lm(yield2 ~ A*B + I(A^2) + I(B^2), data = rsdo)and one sees that there are significant quadratic components. The response surface can be visualized using the wirePlot and contourPlot method of the qualityTools package.
par(mfrow=c(1,2))
wirePlot(A,B,yield2,form="yield2~A*B+I(A^2)+I(B^2)",data=rsdo,theta=-70)
contourPlot(A,B,yield2,form="yield2~A*B+I(A^2)+I(B^2)",data=rsdo)
quadratic fit of the response surface design object rsdo
The following image can be used to compare the outcomes of the factorial and response surface designs with the simulated process. The inactive Factor 3 was omitted.
underlying black box process without noise
Besides this sequential strategy, response surface designs can be created using the method rsmDesign of the qualityTools package. A design with alpha = 1.633, 0 centerpoints in the cube portion and 6 center points in the star portion can be created with:
fdo = rsmDesign(k = 3, alpha = 1.633, cc = 0, cs = 6)and the design can be put in standard order using the randomize method with argument so=TRUE (i. e. standard order). cc stands for centerCube and cs for centerStar.
fdo = randomize(fdo, so = TRUE)Response Surface Designs can also be chosen from a table by using the method rsmChoose of the qualityTools package.
rsdo = rsmChoose()
choosing a predefined response surface design from a table
##
## Choose a response surface design by clicking into the appropriate field
## Waiting for your selection:
##
## NULLSequential Assembly of Response Surface Designs
Sequential assembly is a very important feature of Response Surface Designs. Depending on the features of the (fractional) factorial design a star portion can be augmented using the starDesign method of the qualityTools package. A star portion consists of axial runs and optional center points (cs) in the axial part as opposed to center points (cc) in the cube part.
fdo3 = facDesign(k = 6)
rsdo = starDesign(alpha = "orthogonal", data = fdo3)In case no existing (fractional) factorial design is handed to the starDesign method a list with data.frames is returned which can be assigned to the existing (fractional) factorial design using the star, centerStar and centerCube methods of the qualityTools package.
Randomization
Randomization is achieved by using the randomize method. At this point randomization works for most of the designs types. A random.seed needs to be supplied which is helpful to have the same run order and thus making code and findings reproducible.
Generating a factorial design object with randomized order:
fdo = facDesign(k = 3)
fdo## StandOrder RunOrder Block A B C y
## 4 4 1 1 1 1 -1 NA
## 3 3 2 1 -1 1 -1 NA
## 2 2 3 1 1 -1 -1 NA
## 5 5 4 1 -1 -1 1 NA
## 6 6 5 1 1 -1 1 NA
## 7 7 6 1 -1 1 1 NA
## 8 8 7 1 1 1 1 NA
## 1 1 8 1 -1 -1 -1 NASwitching to standard order using randomize:
randomize(fdo, so = TRUE)## StandOrder RunOrder Block A B C y
## 1 1 1 1 -1 -1 -1 NA
## 2 2 2 1 1 -1 -1 NA
## 3 3 3 1 -1 1 -1 NA
## 4 4 4 1 1 1 -1 NA
## 5 5 5 1 -1 -1 1 NA
## 6 6 6 1 1 -1 1 NA
## 7 7 7 1 -1 1 1 NA
## 8 8 8 1 1 1 1 NARandomizing using random.seed:
randomize(fdo, random.seed = 1234)## StandOrder RunOrder Block A B C y
## 1 1 1 1 -1 -1 -1 NA
## 6 6 2 1 1 -1 1 NA
## 8 8 3 1 1 1 1 NA
## 3 3 4 1 -1 1 -1 NA
## 2 2 5 1 1 -1 -1 NA
## 4 4 6 1 1 1 -1 NA
## 5 5 7 1 -1 -1 1 NA
## 7 7 8 1 -1 1 1 NABlocking
Blocking is another relevant feature and can be achieved by the blocking method of the qualityTools package. In designed experiments, uncontrollable input variables (i. e. they cannot be held constant throughout the design) are identified (see Black Box process modell). The idea of blocking is to run the design in subsets called blocks. Changes in the level of the uncontrollable input variables are assigned to subsets of the designs i. e. blocks were this variable can be held constant (e. g. different lots). Then the blocking variable can be included in the model and is isolated. At this point blocking a design afterwards is not always successful. However, it is unproblematic during the sequential assembly.
Blocking a \(2^k\) Full Factorial Design: In \(2^k\) full factorial designs the parameter blocks can be used to specify the number of blocks. Blocking is supported for \(k \geq 3\) factors.
bfdo = facDesign(k = 3, blocks = 2)
bfdo## StandOrder RunOrder Block A B C y
## 6 6 1 1 1 -1 1 NA
## 4 4 2 1 1 1 -1 NA
## 7 7 3 1 -1 1 1 NA
## 1 1 4 1 -1 -1 -1 NA
## 3 3 5 2 -1 1 -1 NA
## 5 5 6 2 -1 -1 1 NA
## 8 8 7 2 1 1 1 NA
## 2 2 8 2 1 -1 -1 NABlocking a Response Surface Design: In response surface designs the parameter can be used to specify the number of blocks. Blocking is supported for \(k \geq 2\) factors.
brsdo = rsmDesign(k = 2, blocks = 2)
brsdo## StandOrd RunOrder Block A B y
## 3 3 1 1 -1.000 1.000 NA
## 2 2 2 1 1.000 -1.000 NA
## 4 4 3 1 1.000 1.000 NA
## 5 5 4 1 0.000 0.000 NA
## 1 1 5 1 -1.000 -1.000 NA
## 9 9 6 2 0.000 1.414 NA
## 6 6 7 2 -1.414 0.000 NA
## 7 7 8 2 1.414 0.000 NA
## 10 10 9 2 0.000 0.000 NA
## 8 8 10 2 0.000 -1.414 NAReplications
Replications in response surface designs are assigned with the replication parameters cc (centerCube), cs (centerStar), fp (factorialPoints) and sp (starPoints). For example, a #central composite design with:
- 2 centerpoints in the factorial portion of the design i. e. \(2\)
- 1 centerpoint int the star portion of the design i. e. \(1\)
- 2 replications per factorial point i. e. \(2^3 \times 2 = 16\)
- 3 replications per star points \(3\times 2\times 3 = 18\)
makes a total of 37 factor combinations
rsdo = rsmDesign(k=3, blocks=1, alpha=2, cc=2, cs=1, fp=2, sp=3)
nrow(rsdo) #37## [1] 37Many problems involve the simultaneous optimization of more than one response variable. Optimization can be achieved by either maximizing or minimizing the value of the response or by trying to set the response on a specific target. Optimization using the Desirabilities approach (cite{Derringer.1980}), the (predicted) values of the response variables are transformed into values within the interval [0,1] using three different desirability methods for the three different optimization criterias (i. e. minimize, maximize, target). Each value of a response variable can be assigned a specific desirability, optimizing more than one response variable. The geometric mean of the specific desirabilities characterizes the overall desirability.
\[\sqrt[n]{\prod_{i=1}^n {d_i}}\]
This means, for the predicted values of the responses, each factor combination has a corresponding specific desirability and an overall desirability can be calculated. Suppose we have three responses. For a specific setting of the factors the responses have desirabilities such as
- \(d_1=0.7\) for \(y_1\),
- \(d_2=0.8\) for \(y_2\)
- \(d_3=0.2\) for \(y_3\).
The overall desirability \(d_{all}\) is then given by the geometric mean:
\begin{align} d_{all} &= \sqrt[n]{d_1\cdot d_2, \ldots, d_n}\\ \nonumber\\ &= \sqrt[3]{d_1\cdot d_2\cdot d_3}\\ \nonumber\\ & = \sqrt[3]{0.7 \cdot 0.8 \cdot 0.2} \end{align}Desirability methods can be defined using the desires method of the qualityTools package. The optimization direction for each response variable is defined via the min, max and target argument of the desires method. The target argument is set with max for maximization, min for minimization and a specific value for optimization towards a specific target. Three settings arise from this constellation
target = max: min is the lowest acceptable value. If the response variable takes values below min the corresponding desirability will be zero. For values equal or greater than min the desirability will be greater zero.
target = min: max is the highest acceptable value. If the response variable takes values above max the corresponding desirability will be zero. For values equal or less than max the desirability will be greater zero.
target = value: a response variable with a value of value relates to the highest achievable desirability of 1. Values outside min or max lead to a desirability of zero, inside min and max to values within (0,1]
Besides these settings the scale factor influences the shape of the desirability method. Desirability methods can be created and plotted using the desires and plot method of the qualityTools package. Desirabilities are always attached to a response and thus should be assigned to factorial designs.
d1 = desirability(y1, 120, 170, scale = c(1,1), target = "max")
d3 = desirability(y3, 400, 600, target = 500)
d1## Target is to maximize y1
## lower Bound: 120
## higher Bound: 170
## Scale factor is: 1 1
## importance: 1Besides having a summary on the command line, the desirability method can be conveniently visualized using the plot method. With the desirabilities d1 and d3 one gets the following plots.
par(mfrow = c(1,2))
plot(d1, col = 2); plot(d3, col = 2)
plotted desirabilities for y1 and y3
Using desirabilities together with designed experiments
The desirability methodology is supported by the factorial design objects. The output of the desirability method can be stored in the design object, so that information that belongs to each other is stored in the same place (i. e. the design itself). In the following few R lines a designed experiment that uses desirabilities will be shown. The data used comes from (cite{Derringer.1980}). Four responses y1, y2, y3, and y4 were defined. Factors used in this experiment were silica, silan, and sulfur with high factor settings of 1.7, 60, 2.8 and low factor settings of 0.7, 40, 1.8. It was desired to have y1 and y2 maximized and y3 and y4 set on a specific target (see below).
First of all the corresponding design that was used in the paper is created using the method rsmDesign of the qualityTools package. Then we use the randomize method to obtain the standard order of the design.
ddo = rsmDesign(k = 3, alpha = 1.633, cc = 0, cs = 6)
ddo = randomize(ddo, so = TRUE)
names(ddo) = c("silica", "silan", "sulfur") #optional setting of names
highs(ddo) = c(1.7, 60, 2.8) #optional setting of high values
lows(ddo) = c(0.7, 40, 1.8) #optional low valuesThe summary method gives an overview of the design. The values of the responses are attached with the response method of the qualityTools package.
y1 = c(102, 120, 117, 198, 103, 132, 132, 139, 102, 154, 96, 163, 116, 153, 133, 133, 140, 142, 145, 142)
y2 = c(900, 860, 800, 2294, 490, 1289, 1270, 1090, 770, 1690, 700, 1540, 2184, 1784, 1300, 1300, 1145, 1090, 1260, 1344)
y3 = c(470, 410, 570, 240, 640, 270, 410, 380, 590, 260, 520, 380, 520, 290, 380, 380, 430, 430, 390, 390)
y4 = c(67.5, 65, 77.5, 74.5, 62.5, 67, 78, 70, 76, 70, 63, 75, 65, 71, 70, 68.5, 68, 68, 69, 70)The sorted data.frame of these 4 responses is assigned to the design object ddo.
response(ddo) = data.frame(y1, y2, y3, y4)[c(5,2,3,8,1,6,7,4,9:20),]The desirabilities are attached with the desires method of the qualityTools package. y1 and y3 were already defined which leaves the desirabailities for y2 and y4 to be defined.
d2 = desirability(y2, 1000, 1300, target = "max")
d4 = desirability(y4, 60, 75, target = 67.5)The desirabilities need to be defined with the names of the response variables in order to use them with the responses of the design object. The desires method is used as follows.
desires(ddo)=d1; desires(ddo)=d2; desires(ddo)=d3; desires(ddo)=d4Fits are set as in (cite{Derringer.1980}) using the fits methods of the qualityTools package.
fits(ddo) = lm(y1 ~ A+B+C+A:B+A:C+B:C+I(A^2)+I(B^2)+I(C^2), data = ddo)
fits(ddo) = lm(y2 ~ A+B+C+A:B+A:C+B:C+I(A^2)+I(B^2)+I(C^2), data = ddo)
fits(ddo) = lm(y3 ~ A+B+C+A:B+A:C+B:C+I(A^2)+I(B^2)+I(C^2), data = ddo)
fits(ddo) = lm(y4 ~ A+B+C+A:B+A:C+B:C+I(A^2)+I(B^2)+I(C^2), data = ddo)The overall optimum can now be calculated using the method optimum giving the same factor settings as stated in (cite{Derringer.1980}) for an overall desirability of 0.58 and individual desirabilities of \(y1= 0.187\), \(y2 = 1\), \(y3 = 0.664\) and \(y4 = 0.934\).
optimum(ddo, type = "optim")##
## composite (overall) desirability: 0.583
##
## A B C
## coded -0.0533 0.144 -0.872
## real 1.1733 51.442 1.864
##
## y1 y2 y3 y4
## Responses 129.333 1300 466.397 67.997
## Desirabilities 0.187 1 0.664 0.934Mixture Designs
At this time the generation of the different kinds of mixture designs is fully supported including a ternary contour and 3D plot. Analyzing these designs however needs to be done without any specific support by a method of the qualityTools package.
Following the introduced name convention of the qualityTools package the method mixDesign can be used to e. g. create simplex lattice design and simplex centroid designs. The generic methods response, names, highs, lows, units and types are again supported. A famous data set (cite{Cor02}) is given by the elongation of yarn for various mixtures of three factors. This example can be reconstructed using the method mixDesign of the qualityTools package. mdo is an abbreviation of mix design object.
mdo = mixDesign(3, 2, center = FALSE, axial = FALSE, randomize = FALSE, replicates = c(1, 1, 2, 3))
names(mdo) = c("polyethylene", "polystyrene", "polypropylene")
elongation = c(11.0, 12.4, 15.0, 14.8, 16.1, 17.7, 16.4, 16.6, 8.8, 10.0, 10.0, 9.7, 11.8, 16.8, 16.0)
response(mdo) = elongation #set response (i.e. yarn elongation)## [1] "elongation"Again the values of the response are attached with the method response. Calling mdo prints the design. The generic summary method can be used for a more detailed overview.
mdo## StandOrder RunOrder Type A B C elongation
## 1 1 1 1-blend 1.0 0.0 0.0 11.0
## 2 2 2 1-blend 1.0 0.0 0.0 12.4
## 3 3 3 2-blend 0.5 0.5 0.0 15.0
## 4 4 4 2-blend 0.5 0.5 0.0 14.8
## 5 5 5 2-blend 0.5 0.5 0.0 16.1
## 6 6 6 2-blend 0.5 0.0 0.5 17.7
## 7 7 7 2-blend 0.5 0.0 0.5 16.4
## 8 8 8 2-blend 0.5 0.0 0.5 16.6
## 9 9 9 1-blend 0.0 1.0 0.0 8.8
## 10 10 10 1-blend 0.0 1.0 0.0 10.0
## 11 11 11 2-blend 0.0 0.5 0.5 10.0
## 12 12 12 2-blend 0.0 0.5 0.5 9.7
## 13 13 13 2-blend 0.0 0.5 0.5 11.8
## 14 14 14 1-blend 0.0 0.0 1.0 16.8
## 15 15 15 1-blend 0.0 0.0 1.0 16.0The data can be visualized using the wirePlot3 and contourPlot3 methods . In addition to the wirePlot and contourPlot methods the name of the third factor (i. e. C) and the type of standard fit must be given. Of course it is possible to specify a fit manually using the form argument with a formula.
par(mfrow=c(1,2))
contourPlot3(A, B, C, elongation, data = mdo, form = "quadratic")
wirePlot3(A, B, C, elongation, data=mdo, form="quadratic", theta=-170)
ternary plots for the elongation example
Taguchi Designs
Taguchi Designs are available using the method taguchiDesign of the qualityTools package. There are two types of taguchi designs:
Single level: all factors have the same number of levels (e.g. two levels for a L4_2)
Mixed level: factors have different number of levels (e.g. two and three levels for L18_2_3)
Most of the designs that became popular as taguchi designs however are simple \(2^k\) fractional factorial designs with a very low resolution of III (i. e. main effects are confounded with two factor interactions) or other mixed level designs and are originally due to contributions by other e. g. Plackett and Burman, Fisher, Finney and Rao (cite{Box88}).
A design can be created using the taguchiDesign method of the qualityTools package. The generic method names, units, values, summary, plot, lm and other methods again are supported. This way the relevant information for each factor can be stored in the design object tdo (taguchi design object) itself.
set.seed(1234)
tdo = taguchiDesign("L9_3")
values(tdo) = list(A = c(20, 40, 60), B = c("material 1", "material 2", "material 3"), C = c(1,2,3))
names(tdo) = c("Factor 1", "Factor 2", "Factor 3", "Factor 4")
summary(tdo)## Taguchi SINGLE Design
## Information about the factors:
##
## A B C D
## value 1 20 material 1 1 1
## value 2 40 material 2 2 2
## value 3 60 material 3 3 3
## name Factor 1 Factor 2 Factor 3 Factor 4
## unit
## type numeric numeric numeric numeric
##
## -----------
##
## StandOrder RunOrder Replicate A B C D y
## 1 7 1 1 3 1 3 2 NA
## 2 1 2 1 1 1 1 1 NA
## 3 6 3 1 2 3 1 2 NA
## 4 4 4 1 2 1 2 3 NA
## 5 2 5 1 1 2 2 2 NA
## 6 8 6 1 3 2 1 3 NA
## 7 5 7 1 2 2 3 1 NA
## 8 3 8 1 1 3 3 3 NA
## 9 9 9 1 3 3 2 1 NA
##
## -----------The response method is used to assign the values of the response variables. effectPlot can be used once more to visualize the effect sizes of the factors:
response(tdo) = rnorm(9)
effectPlot(tdo, col = 2)
effect plot for the taguchi experiment
Session Info
sessionInfo()## R version 3.3.1 (2016-06-21)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.1 LTS
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] qualityTools_1.55 MASS_7.3-45 Rsolnp_1.16
##
## loaded via a namespace (and not attached):
## [1] Rcpp_0.12.6 codetools_0.2-14 digest_0.6.10 truncnorm_1.0-7
## [5] formatR_1.4 magrittr_1.5 evaluate_0.9 stringi_1.1.1
## [9] rmarkdown_1.0 tools_3.3.1 stringr_1.1.0 yaml_2.1.13
## [13] parallel_3.3.1 rsconnect_0.7 htmltools_0.3.5 knitr_1.14