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

## $$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:  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

## Replicated 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")

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

## Moving 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    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.

## 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 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 (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     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
##
## NULL

### Sequential 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 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

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

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)

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

## Mixture 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.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

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

## 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
## [13] parallel_3.3.1   rsconnect_0.7    htmltools_0.3.5  knitr_1.14