- \(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.

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 NA
```

The 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 object`

Analyzing this design is quite easy using the methods `effectPlot`

, `interactionPlot`

, `lm`

as well as `wirePlot`

and `contourPlot`

:

`effectPlot(fdo, classic = TRUE)`

`interactionPlot(fdo)`

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.003685
```

The 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
```

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)
```

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).

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: III
```

The 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 0
```

Fractional 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()`

```
##
## 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 NA
```

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)`

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")
```

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")
```

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.2518
```

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 193
```

`sao`

```
## 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 NA
```

Since 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)
```

At this point the step size was chosen quite small for illustration purposes.

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 NA
```

and 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 response
```

Looking 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`

(**r**esponse **s**urface **d**esign **o**bject) 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 NA
```

Using 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) = yield2
```

A 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)
```

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.

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()`

```
##
## Choose a response surface design by clicking into the appropriate field
## Waiting for your selection:
##
## NULL
```

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 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 NA
```

Switching 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 NA
```

Randomizing 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 NA
```

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 NA
```

**Blocking 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 NA
```

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] 37`

\subsection{Multiple Response Optimization - Desirabilites}
Many 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: 1
```

Besides 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)
```

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 values
```

The `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)=d4`

Fits 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.934
```

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.0
```

The 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)
```

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)
```

`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
```