\documentclass[fontsize=12pt, a4paper,bibliography=totoc]{scrartcl}
\usepackage[english]{babel} 
\usepackage[latin1]{inputenc}
\usepackage{natbib}
\usepackage{Sweavel}

\usepackage[pdftex]{graphicx}
\SweaveOpts{keep.source=TRUE}
%\SweaveOpts{prefix.string=graphics/plot, eps = FALSE, pdf = TRUE}
\def\Rcolor{\color{red}}
\def\Routcolor{\color[rgb]{0,0,0.5}}
\def\Rcommentcolor{\color[rgb]{0.259,0.631,0.259}}
\def\Rbackground{\color{white}}       %Rahmenfarbe Input
\def\Routbackground{\color{white}}    %Rahmenfarbe Output
%huebsche Schrift
\usepackage[T1]{fontenc}

\usepackage{lmodern}
\tolerance=2000
\usepackage{microtype}

\addtokomafont{sectioning}{\rmfamily}

\usepackage{url}
\usepackage{listings}
\usepackage{framed}
\usepackage{array,booktabs,dcolumn, multirow}  %schoene Tabelle
\usepackage{ragged2e}
\usepackage{amsmath,amsfonts,amssymb}
\usepackage[]{natbib}
\usepackage{marvosym}
\usepackage{keystroke}
\usepackage{blindtext}          % Erzeugung von Textpassagen
\usepackage{layout}
\usepackage{subfigure}
\usepackage{listings}

\usepackage[pdftitle={Working with the qualityTools package},pdfauthor={Thomas Roth}, pdfsubject={This document contains a brief illustration of the R package qualityTools for teaching Quality Science related statistical methods with regard to the Six Sigma methodology and Quality Management.},pdfkeywords={Process Capability Analysis for normal and non-normal distributions in R, Distribution Fitting and Identification in R, Anderson Darling Goodness of Fit Test for normal and non-normal distributions in R, Measurement Systems Analysis in R,  Gage R\&R in R, Repeatability and Reproducibility studies in R, Gage Capapbility in R, Design of Experiments in R, Desirability, factorial designs in R, fractional factorial designs in R, confounding and alias structure, response surface designs in R, sequential assembly of central composite designs in R, rotatibility and orthogonality and blocking, response surface methods, mixture designs in R, taguchi designs in R, simultaneous optimization of multiple response variables in R, QQPlot, Probability Plot, Six Sigma in R, OC-Curves in R, Operation Characteristic, Define, Measure, Analyze, Improve and Control of the Six Sigma approach, Methods for Measurement, analysis and improvement as part of the Continual Improvement of the quality management system}, colorlinks=true, linkcolor=red2]{hyperref}


\setkeys{Gin}{width=\textwidth}


\usepackage{fancyhdr,color,geometry,fancyvrb,ulem,xspace,longtable,colortbl,bibgerm} %scrextend wurde entfernt
\usepackage{ucs}
\usepackage{makeidx}
\usepackage[labelfont={bf}]{caption} % Bei Bild- und Tabellenbeschriftung wird nun Tabelle und Bild fett geschrieben
\usepackage[flushleft,neverdecrease]{paralist} % wird bentigt fr die rmische Liste
\usepackage{float}
\usepackage{rotating}

\addtokomafont{caption}{\small} %verkleinert die Bildunterschriften/Tabellenberschriften

\geometry{a4paper, top=30mm, left=25mm, right=25mm, bottom=20mm, headsep=12.5mm, footskip=12.5mm, marginparsep=6pt, marginparwidth=20mm, headheight=20pt}

%\renewcommand{\thefigure}{\arabic{section}.\arabic{figure}}
%\renewcommand{\thetable}{\arabic{section}.\arabic{table}}
\renewcommand{\em}{\it}
%\renewcommand{\emph}{\textit}
\renewcommand{\emph}{\textbf}
\renewcommand{\arraystretch}{1} % beachte, dass im Dokument Konfidenzintervalle auch arraystretch nach einer Tabelle gendert werden muss
\newcommand{\remph}[1]{\textcolor{navy}{\texttt{#1}}}
\newcommand{\Remph}[1]{\textcolor{red2}{\texttt{#1}}}
\newcommand{\R}{\textrm{R}\xspace}
\newcommand{\abe}{\\[1.0ex]} 

\definecolor{red2}{RGB}{255,0,0}
\definecolor{navy}{RGB}{0,0,128}
\definecolor{Gruen}{RGB}{66,161,66}


\newcommand{\rIn}[1]{\textcolor{red2}{\texttt{#1}}}			%Fuer R Output im Text
\newcommand{\rOut}[1]{\textcolor{navy}{\texttt{#1}}}		%Fuer R Output im Text
%
%%Umgebung fr R-Eingaben
%\lstnewenvironment{RI}{\footnotesize \color{red2} \lstset{language=R, basicstyle=\ttfamily, commentstyle=\color{Gruen}, breaklines=true}}{}
%%Umgebung fr R-Ausgaben
%\lstnewenvironment{RO}{\footnotesize \color{navy} \lstset{language=R, basicstyle=\ttfamily, commentstyle=\color{Gruen}, breaklines=true}}{}

\author{Thomas Roth \\Technical University Berlin}
%\Plainauthor{Thomas Roth}
\title{Working with the \pkg{qualityTools} package}
%\Plaintitle{Working with the qualityTools package}

%\subtitle{A short introduction}
%\date{\today}
%\thanks{DISCLAIMER: Used for teaching thus methods are optimized for illustrating special concepts and best practice.}

\fancypagestyle{standard}{\fancyhf{}
\rhead{\textbf{\textcolor{red2}{\textbar}} \copyright Roth}
\lhead{\nouppercase{\leftmark}}%\small \textcolor{navy} 
\cfoot{\thepage}
}

\pagestyle{standard}

\begin{document}
\SweaveOpts{concordance=FALSE}
%\maketitle
\begin{titlepage}

\begin{center}
\includegraphics[width=0.4\textwidth]{Bilder/qualityToolsGrafik.pdf}
\par
\vspace{1cm}
\huge{\textbf{Working with the qualityTools package}}
\\[1ex]
\large{\textbf{A short introduction\footnote{An updated version of this document can be found under \url{http://www.r-qualitytools.org}. A web application can be found under \url{http://webapps.r-qualitytools.org}}}}
\\[1ex]
\Large{Thomas Roth}
\\[1ex]
\large{\today}
\vspace{1cm}
\end{center}


%\VignetteIndexEntry{qualityTools: A Package for Teaching Statistics in Quality Science. R package version 1.31.}
%\VignetteDepends{qualityTools}
%\VignetteKeywords{Quality Science, Process Capability, Distribution Fitting, Gage R\&R, Measurement System Analysis, Design of Experiments, Desirability, factorial design, fractional factorial design, response surface design, response surface methods, simultaneous optimization of response variables, QQPlot, Probability Plot}
%\VignettePackage{qualityTools}

\begin{abstract}
This vignette is intended to give a short introduction into the methods of the qualityTools package. The qualityTools package contains methods associated with the \textbf{D}efine \textbf{M}easure \textbf{A}nalyze \textbf{I}mprove and \textbf{C}ontrol (i.e. \textbf{DMAIC}) problem solving cycle of the Six Sigma Quality Management methodology. Usage of these methods is illustrated with the help of artificially created datasets.

\begin{itemize}
\item \textbf{Define:} Pareto Chart
\item \textbf{Measure:} Probability and Quantile-Quantile Plots, Process Capability Ratios for various distributions and Gage R\&R
\item \textbf{Analyze:} Pareto Chart, Multi-Vari Chart, Dot Plot
\item \textbf{Improve:} Full and fractional factorial, response surface, mixture and taguchi designs as well as the desirability approach for simultaneous optimization of more than one response variable. Normal, Pareto and Lenth Plot of effects as well as Interaction Plots
\end{itemize}

\end{abstract}

\end{titlepage}
\pagebreak
\setcounter{page}{1}

\tableofcontents
\pagebreak





%%%%Introduction%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Working with the qualityTools package} \label{sec:intro}

Working with the qualityTools package is straightforward as you will see in the next few pages. The qualityTools package was implemented for teaching purposes in order to serve as a (Six-Sigma)-Toolbox and contains methods that are associated to a problem solving cycle. There are many problem solving cycles around with new ones emerging although most of these new ones take on special aspects. A very popular problem solving cycle is the \textbf{PDCA} cycle (i.e. plan $\curvearrowright$ do$\curvearrowright$ check$\curvearrowright$ act$\circlearrowleft$) which was made popular by Deming\footnote{William E. Deming} but goes back to Shewart\footnote{Walter A. Shewhart}. As part of the widely known and accepted Six-Sigma-Methodology some enhancements to this problem solving cycle were made and a problem solving cycle consisting of the five phases Define, Measure, Analyze, Improve and Control emerged.

\begin{description}
\item[Define]	Describe the problem and its (financial) consequences. Interdisciplinary workgroups contribute to the problem and its consequences which is the pivotal stage in narrowing down the problem. Process flow diagrams identify crucial process elements (i.e. activities), creativity techniques such as Brainwriting and Brainstorming as well as the SIPOC\footnote{Suppliers, inputs, process, outputs, customers} technique should lead, depending on the future size of the project, to possibly a project charter. Amongst other things, the project charter serves as a descripition of the process, customer\'s requirements in relation to corporate objectives.
\item[Measure] 	Come up with a reasonable plan for collecting the required data and make sure that the measurement systems are capable (i.e. no or known bias and as little system immanent variation contributing to the measurements as possible). \emph{Variation} and \emph{bias} are the enemy to finding effects. The bigger the background noise the less probable are the chances of success using limited resources for all kinds of experiments. Within the Measure phase a description of the situation is given with the help of process- or gage capability indices (MSA\footnote{Measurement Systems Analysis} Type I) or a Gage R\&R (MSA Type II)\cite{MSA.2010}.
\item[Analyze]	Try to find the root causes of the problem using various statistical methods such as histograms, regression, correlation, distribution identification, analysis of variance, multi-vari-charts.
\item[Improve]	Use designed experiments i.e. full and fractional factorials, response surface designs, mixture designs, taguchi designs and the desirability concept to find optimal settings or solutions for a problem.
\item[Control]	Once an improvement was achieved it needs to be secured, meaning arrangements need to be implemented in order to secure the level of improvement. Besides proper documentation, the use of statistical process control (i.e. quality-control-charts) can be used to monitor the behavior of a process. Although quite often referred to as \emph{Show Programm for Customers}, SPC is able to help to distinguish between \emph{common causes} and \emph{special causes} in the process behavior.
\end{description}
\pagebreak






%%%%DEFINE%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{qualityTools in DEFINE} \label{sec:define}
\marginpar{\includegraphics[width=1.5cm]{Bilder/D.pdf}}

Most techniques in the Define phase are not related to substantial use of statistical methods. The objective of the DEFINE phase is to bring together all parties concerned, grasp their knowledge and insights to the process involved, set a common objective and DEFINE how each party contributes(or the role each party takes) to the solving of the problem. In order not to get lost in subsequent meetings and ongoing discussion, this common objective, the contribution of each party, milestones and responsibilities need to be written down in what is known to be a Project Charter. Of course, problems with easy-to-identify causes are not subject of these kind of projects.

However, a classical visualization technique that is used in this phase and available in the qualityTools package is the pareto chart. Pareto charts are special forms of bar charts that help to separate the vital few from the trivial many causes for a given problem (e.g.the most frequent cause for a defective product). This way pareto charts visualize how much a cause contributes to a specific issue.

Suppose a company is investigating non compliant units (products). 120 units were investigated and 6 different types of defects (qualitative data) were found. The defects are named A to F. The defects data can be found in defects.

<<fig=FALSE, echo=FALSE>>=
library(qualityTools)
@

\begin{figure}[htbp]\begin{center}
<<fig=TRUE, echo=TRUE, results=hide, width=14, height=7>>=
#create artificial defect data set
defects = c(rep("E", 62), rep("B", 15), rep("F", 3), rep("A", 10),
rep("C",20), rep("D", 10))
paretoChart(defects)
@
\caption{Pareto Chart}\label{fig:paretochart}
\end{center}\end{figure}

This pareto chart might convey the message that in order to solve 68 percent of the problem 33 percent of the causes (\emph{vital few}\footnote{the vital few and the trivial many - 20 percent of the defects cause 80 percent of the problems}) need to be subject of an investigation.

Besides this use case, pareto charts are also used for visualizing the effect sizes of different factors for designed experiments (see \Remph{paretoPlot}).
\pagebreak






%%%MEASURE%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{qualityTools in MEASURE} \label{sec:measure}
\marginpar{\includegraphics[width=1.5cm]{Bilder/M.pdf}}

Collecting data involves the use of measurement systems often referred to as gages. In order to make a statement regarding the quality, i.e. the degree in which a set of inherent characteristics meets requirements, of a product\cite{ISO21747}, the capability of the measurement system used needs to be validated.

Gages can have two types of impairments:

\begin{itemize}
\item a bias (an assumed constant shift of values for measurements of equal magnitude)
\item variation
	\begin{itemize}
	\item introduced by other factors e.g. operators using these gages
	\item system immanent variation of the measurement system itself 
	\end{itemize}
\end{itemize}

These impairments lead to varying measurements for repeated measurements of the same unit (e.g. a product). The amount of tolerable variation of course depends on the number of distinctive categories you need to be able to identify in order to characterize the product. This tolerable amount of variation for a measurement system relates directly to the tolerance range of the characteristics of a product.

The capability of a measurement system is crucial for any conclusion based on data. Non-capable Measurement Systems due to a non adjusted bias, or a Measurement System immanent variation implicate two serious errors of judgement.

\begin{itemize}
\item Accepting items that are actually \emph{out of tolerance}
\item Declining items that are actually \emph{within tolerance}
\end{itemize}

\noindent Thus the capability of Measurement Systems is directly related to costs (see figure \ref{fig:gageCapability}).

\begin{figure}\begin{center}
	\includegraphics[width=\textwidth]{Bilder/gageCapability.pdf}
	\caption{Errors of judgement due to non-capable Measurement Systems}
	\label{fig:gageCapability}
\end{center}\end{figure}


\subsection{Gage Capability - MSA Type I}
Suppose an engineer wants to check the capability of an optical Measurement Device. An unit with known characteristic ($x_m = 10.033mm$) is repeatedly measured $n=25$ times. From the measurement values the mean $\overline{x_g}$ and standard deviation $s_g$\footnote{$\sigma_g$ denotes the standard deviation of the gage which is also referred to as repeatability} can be calculated.

Basically the calculation of an capability index comprises two steps. First a fraction of the tolerance width (i.e. $USL - LSL$)\footnote{Upper Specification Limit and Lower Specification Limit} is calculated. The fraction typically relates to 0.2. In a second step this fraction is set in relation with a measure of the process spread (i.e. the range in which 95.5\% or 99.73\% of the characteristics of a process are to be expected). For normal distributed measurement values this relates to $k=2\sigma_g$ and $k=3\sigma_g$ calculated from the measurement values.  For non-normal distributed data the corresponding quantiles can be taken. If there's no bias this calculation represents the capability index $c_g$ and reflects the true capability of the measurement device.

\begin{align}
c_g &= \frac{0.2\cdot(USL - LSL)}{6\cdot s_g}\\
	&= \frac{0.2\cdot(USL - LSL)}{X_{0.99865} - X_{0.00135}}
\end{align}

However, if there's a bias it is taken into account by substracting it from the numerator. In this case $c_g$ reflects only the potential capability (i.e. capability if bias is corrected) and $c_{gk}$ is an estimator of the actual capability. The bias is calculated as the difference between the known characteristic $x_m$ and the mean of the measurement values $x_g$

\begin{align}
c_{gk} = \frac{0.1\cdot(USL - LSL) - \left| x_m - x_g \right|}{3\cdot s_g}
\end{align}


Determining if the bias is due to chance or not can be done with the help of a t-test which has the general form:

\begin{align}
t = \frac{\text{difference in means}}{\text{standard error of the difference}} = \frac{Bias}{\frac{s_{Bias}}{\sqrt{n}}}
\end{align}

Besides bias and standard deviation it is important to check the run-chart of the measurement values. Using the qualityTools package, all this is easily achieved using the \Remph{cg} method. The output of the \Remph{cg} method is shown in figure \ref{fig:cgMethod}.


<<fig=FALSE, echo = TRUE>>=
x = c(9.991, 10.013, 10.001, 10.007, 10.010, 10.013, 10.008, 10.017, 10.005,
10.005, 10.002, 10.017, 10.005, 10.002,  9.996, 10.011, 10.009, 10.006,
10.008, 10.003, 10.002, 10.006, 10.010, 9.992, 10.013)
cg(x, target = 10.003, tolerance = c(9.903, 10.103))
@

\begin{figure}[htbp]
\begin{center}
\includegraphics[width=\textwidth]{Bilder/cg.pdf}
\caption{Potentially capable gage with bias}
\label{fig:cgMethod}
\end{center}\end{figure}


\subsection{Gage Repeatability\&Reproducibility - MSA Type II}
A common procedure applied in industry is to perform a Gage R\&R analysis to assess the repeatability and the reproducibility of a measurement system. R\&R stands for repeatability and reproducibility. Repeatability hereby refers to the precision of a measurement system (i.e. the standard deviation of subsequent measurements of the same unit). Reproducibility is the part of the overall variance that models the effect of different e.g. operators performing measurements on the same unit and a possible interaction between different operators and parts measured within this Gage R\&R. The overall model is given by
\begin{align}
	\sigma^2_{total} &= \sigma^2_{Parts} +\sigma^2_{Operator} +\sigma^2_{Parts\times Operator} +\sigma^2_{Error} \label{form:gageRR}
\end{align}
where $\sigma^2_{Parts}$ models the variation between different units of the same process. $\sigma^2_{Parts}$ is thus an estimate of the inherent process variability. Repeatability is modeled by $\sigma^2_{Error}$ and reproducibility by $\sigma^2_{Operator} +\sigma^2_{Parts\times Operator}$.

Suppose 10 randomly chosen units were measured by 3 randomly chosen operators. Each operator measured each unit two times in a randomly chosen order. The units were presented in a way they could not be distinguished by the operators.

The corresponding gage R\&R design can be created using the \Remph{gageRRDesign} method of the qualityTools package. The measurements are assigned to this design using the response method. Methods for analyzing this design are given by \Remph{gageRR} and \Remph{plot}.

<<echo=TRUE>>=
#create a gage RnR design
design = gageRRDesign(Operators=3, Parts=10, Measurements=2, randomize=FALSE)
#set the response	
response(design) = c(23,22,22,22,22,25,23,22,23,22,20,22,22,22,24,25,27,28,
23,24,23,24,24,22,22,22,24,23,22,24,20,20,25,24,22,24,21,20,21,22,21,22,21,
21,24,27,25,27,23,22,25,23,23,22,22,23,25,21,24,23)
#perform Gage RnR
gdo = gageRR(design)
@


<<GageRR, echo=TRUE>>=
#visualization of Gage RnR
plot(gdo)
@

\begin{figure}[htbp]\begin{center}
<<fig=TRUE, echo=FALSE, results=hide, width=7, height=7>>=
<<GageRR>>
@
\caption{Visualization of the Gage R\&R}
\label{fig:gageRR}
\end{center}\end{figure}

The standard graphical output of a Gage R\&R is given in figure \ref{fig:gageRR}.

The barplot gives a visual representation of the Variance Components. \rOut{totalRR} depicts the total \emph{R}epeatability and \emph{R}eproducibility. 48\% of the variance is due to 35\% \rOut{repeatability} (i.e. variation from the gage itself) and 13\% reproducibility (i.e. effect of operator and the interaction between operator and part). It can be seen from the AnOVa table that an interaction between parts and operators is not existing. The remaining 52\% (51.7 in column \rOut{VarCompContrib}) of variation stems from differences between parts taken from the process (i.e. process inherent variation) which can be seen also in the \emph{Measurement by Part} plot. The variaton for measurements taken by one operator is roughly equal for all three operators (\emph{Measurement by Operator}) although operator C seems to produce values that are most of the time larger than the values from the other operators (\emph{Interaction Operator: Part}).

Besides this interpretation of the results critical values (see table \ref{tab:gageCrit}) for totalRR also refered to as GRR\footnote{Gage Repeatability Reproducibility} are used within industry. However, a measurement system should never be judged by critical values alone.

\begin{table}[htbp]
\caption{critical values for judging the suitability of measurement system}
\label{tab:gageCrit}
\begin{center}
\begin{tabular}{cc}
\toprule
Contribution of \rOut{total RR} & Capabability\\
\midrule
$\leq 0.1$	& suitable\\
 $<0.1$ and $< 0.3$ & limited suitability depending upon circumstances\\
$\geq 0.3$  & not suitable\\
\bottomrule
\end{tabular}
\end{center}
\end{table}



\paragraph*{Checking for interaction}
The interaction plot provides a visual check of possible interactions between Operator and Part. For each Operator the average measurement value is shown as a function of the part number. Crossing lines indicate that operators are assigning different readings to identical depending on the combination of Operator and Part. Different readings means in the case of an interaction between Operator and Part that on average sometimes smaller or bigger values are assigned depending on the combination of Operator and Part. In this case, lines are practically not crossing but Operator C seems to systematically assign larger readings to the parts than his colleagues.

\paragraph*{Operators}
To check for an operator dependent effect, measurements are plotted grouped by operators in form of boxplots. Boxplots that differ in size or location might indicate e.g. possible different procedures within the measurement process, which then lead to a systematic difference in the readings. In this case one might discuss a possible effect for operator C which is also supported by the interaction plot.

\paragraph*{Inherent process variation}
Within this plot Measurements are grouped by operator. Due to the repeated measurements by different Operators per Part an insight into the process is given. A line connecting the mean of the measurements of each part provides an insight into the inherent process variation. Each part is measured number of operator times number of measurements per part.

\paragraph*{Components of variation}
In order to understand the output of a Gage R\&R study formula \ref{form:gageRR} should be referenced. The variance component \rOut{totalRR} (\rOut{VarComp} column) represents the total \emph{R}epeatability and \emph{R}eproducibility. Since variances are simply added 1.664 is the sum of 1.209 (\rOut{repeatability} given by $\sigma^2_{Error}$) and 0.455 (\rOut{reproducibility}). Reproducibility itself is the sum of \rOut{Operator} ($\sigma^2_{Operator}$) and \rOut{Operator:Part} ($\sigma^2_{Parts\times Operator}$). Since there's no interaction Reproducibility amounts to 0.455. \rOut{Part to Part} amounts to 1.781. Together with the total of repeatability and reproducibility this gives $\sigma^2_{Total} = 3.446$

\subsubsection{Relation to the Measurement Systems Terminology}
The Measurement Systems Analysis Manual \cite{MSA.2010}uses a specific Terminology for the terms \rOut{repeatability}, \rOut{reproducibility}, \rOut{Operator}, \rOut{Part to Part}, \rOut{totalRR} and the interaction \rOut{Operator:Part}. The objective of this paragraph is to give a short overview of these terms and how they relate to the terms used in the \rOut{gageRR} methods of the qualityTools package.
\begin{description}
\item[EV] stands for Equipment Variation which is the variation due to the \rOut{repeatability}
\item[AV] stands for Appraiser Variation which is the variation due to the \rOut{operator}s.
\item[INT] stands for the interaction Appraiser:Part which is the \rOut{Operator:Part} interaction
\item[GRR] stands for Gage Repeatability\&Reproducibility and refers to the variation introduced by the measurement system. The equivalent to this term is \rOut{totalRR} which is the sum of \rOut{repeatability} and \rOut{reproducibility}.
\item[PV] stands for Part Variation which relates to \rOut{Part to Part} 
\end{description}

\pagebreak






%%%%ANALYZE%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{qualityTools in ANALYZE} \label{sec:analyse}
\marginpar{\includegraphics[width=1.5cm]{Bilder/A.pdf}}

\subsection{Process Capability}

Besides the capability of a measurement system, often the capability of a process is of interest or needs to be assessed e.g. as part of a supplier customer relationship in industry. Process Capability Indices basically tells one how much of the tolerance range is being used by common cause variation of the considered process. Using these techniques one can state how many units (e. g. products) are expected to fall outside the tolerance range (i.e. defective regarding the requirements determined before) if for instance production continues without intervention. It also gives insights into where to center the process if shifting is possible and meaningful in terms of costs. There are three indices which are also defined in the corresponding ISO 21747:2006 document\cite{ISO21747} .

\begin{align}
c_p &= \frac{USL - LSL}{Q_{0.99865} - Q_{0.00135}}\\
\nonumber \\
c_{pkL} &= \frac{Q_{0.5} - LSL}{Q_{0.5} - Q_{0.00135}}\\
\nonumber\\
c_{pkU} &= \frac{USL - Q_{0.5}}{Q_{0.99865} - Q_{0.5}}
\end{align}

$c_p$ is the potential process capability giving one the process capability that could be achieved if the process can be centered within specification limits\footnote{USL - Upper Specification Limit\\LSL - Lower Specification Limit} and $c_{pk}$ is the actual process capability which incorporates the location of the distribution (i.e. the center) of the characteristic within the specification limits. For one sided specification limits $c_{pkL}$ and $c_{pkU}$ exist with $c_{pk}$ being equal to the smallest capability index. As one can imagine in addition the location of the distribution of the characteristic the shape of the distribution is relevant too. Assessing the fit of a specific distribution for given data can be done via probability plots (\Remph{ppPlot}) and quantile-quantile plots (\Remph{qqPlot}), as well as formal test methods like the Anderson Darling Test.

Process capabilities can be calculated with the \Remph{pcr} method of the qualityTools package. The \Remph{pcr} method plots a histogram of the data, the fitted distribution and returns the capability indices along with the estimated\footnote{Fitting the distribution itself is accomplished by the fitdistr method of the R-package MASS.} parameters of the distribution, an Anderson Darling Test for the specified distribution and the corresponding QQ-Plot.

<<NormPcr, echo=TRUE, results=hide>>=
set.seed(1234)
#generate some data
norm = rnorm(20, mean = 20)
#generate some data
weib = rweibull(20, shape = 2, scale = 8)
#process capability
pcr(norm, "normal", lsl = 17, usl = 23)
@
<<WeibPcr, echo=TRUE, results=hide>>=
#process cabapility
pcr(weib, "weibull", usl = 20)
@

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\setkeys{Gin}{width=0.45\textwidth}

\begin{figure}[htbp]\begin{center}
\subfigure[normal distribution]{
<<fig=TRUE, echo=FALSE, results=hide, height=7, width=7>>=
<<NormPcr>>
@
}
\hfill
\subfigure[weibull distribution]{
<<fig=TRUE, echo=FALSE, results=hide, height=7, width=7>>=
<<WeibPcr>>
@
}
\caption{Process Capability Ratios for weibull and normal distribution}
\label{fig:pcr}
\end{center}\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\setkeys{Gin}{width=\textwidth}

Along with the graphical representation an Anderson Darling Test for the corresponding distribution is returned.

<<echo=FALSE>>=
pcr(weib, "weibull", usl = 20)
@

Q-Q Plots can be calculated with the \Remph{qqPlot} function of the qualityTools package (figure \ref{fig:qqPlot}).

<<qqPlot, echo=TRUE, results=hide>>=
par(mfrow = c(1,2))
qqPlot(weib, "weibull"); qqPlot(weib, "normal")
@

\begin{figure}[htbp]\begin{center}
<<fig=TRUE, echo=FALSE, results=hide, height=7, width=14>>=
<<qqPlot>>
@
\caption{QQ-Plots for different distributions}
\label{fig:qqPlot}
\end{center}\end{figure}

Probability Plots can be calculated with the \Remph{ppPlot} function of the qualityTools package (figure \ref{fig:ppPlot}).

<<ppPlot, echo=TRUE, results=hide>>=
par(mfrow = c(1,2))
ppPlot(norm, "weibull"); ppPlot(norm, "normal")
@

\begin{figure}[htbp]\begin{center}
<<fig=TRUE, echo=FALSE, results=hide, height=7, width=14>>=
<<ppPlot>>
@
\caption{PP-Plots for different distributions}
\label{fig:ppPlot}
\end{center}\end{figure}
\pagebreak











%%%%IMPROVE%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{qualityTools in IMPROVE} \label{sec:improve}
\marginpar{\includegraphics[width=1.5cm]{Bilder/I.pdf}}

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 \textbf{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 (figure \ref{fig:blackbox}).

\begin{figure}[H]\begin{center}
	\includegraphics[width=0.65\textwidth]{Bilder/blackbox.pdf}
	\caption{\emph{Black Box} model of a process}
	\label{fig:blackbox}
\end{center}\end{figure}

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 \textbf{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 held 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 leads to variations in the response variables despite identical factor settings in an experiment.


\subsection{$2^k$ Factorial Designs}\label{subsec:facDesigns}
In order to find more about this black box model one can come up with a $2^k$ factorial design by using the method \Remph{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) is simulated by the method \Remph{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.

<<echo=TRUE>>=
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
@


The response of this fictional process is given by the \Remph{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 

<<echo=TRUE, results=hide>>=
#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.

% <<echo=TRUE, results=hide>>=
% yield = c(simProc(120,140, 2),simProc(80,120, 2),simProc(120,120, 2),
% simProc(80,140, 2),simProc(120,140, 1),simProc(80,120, 1),
% simProc(90,130, 1.5),simProc(120,120, 1), simProc(90,130, 1.5),
% simProc(90,130, 1.5),simProc(90,130, 1.5), simProc(80,140, 1))
% @

<<echo=TRUE, results=hide>>=
yield = c(simProc(120,140, 1),simProc(80,140, 1),simProc(120,140, 2),
simProc(120,120, 1),simProc(90,130, 1.5),simProc(90,130, 1.5),
simProc(80,120, 2),simProc(90,130, 1.5),simProc(90,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 \Remph{response} method.

<<echo=TRUE, results=hide>>=
response(fdo) = yield	#assign yield to the factorial design object
@

Analyzing this design is quite easy using the methods \Remph{effectPlot}, \Remph{interactionPlot}, \Remph{lm} as well as \Remph{wirePlot} and \Remph{contourPlot} (figure \ref{fig:effectplot})

<<effPlot, echo=TRUE, results=hide>>=
effectPlot(fdo, classic = TRUE)
@
<<iaPlot, echo=TRUE, results=hide>>=
interactionPlot(fdo)
@

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\setkeys{Gin}{width=0.45\textwidth}

\begin{figure}[htbp]\begin{center}
\subfigure{
<<fig=TRUE, echo=FALSE, results=hide, height=7, width=7>>=
<<effPlot>>
@
}
\hfill
\subfigure{
<<fig=TRUE, echo=FALSE, results=hide, height=7, width=7>>=
<<iaPlot>>
@
}
\caption{effect- and interaction plot for the factorial design}
\label{fig:effectplot}
\end{center}\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\setkeys{Gin}{width=\textwidth}

The factorial design in fdo can be handed without any further operations directly to the base \Remph{lm} method of R.

<<echo=TRUE>>=
lm.1 = lm(yield ~ A*B*C, data = fdo)
summary(lm.1)
@


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 \Remph{paretoPlot} method of the qualityTools package (figure \ref{fig:paretoNormalPlot}). Another visualization technique commonly found is a normal plot using the \Remph{normalPlot} method of the qualityTools package.

<<ParetPlot, echo=TRUE, results=hide>>=
par(mfrow = c(1,2))
paretoPlot(fdo)
normalPlot(fdo)
@

\begin{figure}[htbp]\begin{center}
<<fig=TRUE, echo=FALSE, results=hide, width=14, height=7>>=
<<ParetPlot>>
@
\caption{pareto plot of the standardized effects and normal plot of the coefficients}
\label{fig:paretoNormalPlot}
\end{center}\end{figure}

The relation between the factors A and B can be visualized as 3D representation in form of a wireframe or contour plot using the \Remph{wirePlot} and \Remph{contourPlot} method of the qualityTools package (figure \ref{fig:multipleResponsePlot1}). Again, no further transformation of the data is needed!

<<wirePlot, echo=TRUE, results=hide>>= %, width = 14, height = 7
par(mfrow = c(1,2))
wirePlot(A, B, yield, data = fdo)
contourPlot(A, B, yield, data = fdo)
@

\begin{figure}[htbp]\begin{center}
<<fig=TRUE, echo=FALSE, results=hide, width=14, height=7>>=
<<wirePlot>>
@
\caption{response surface and contour plot}
\label{fig:wirePlot}
\end{center}\end{figure}

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 the model is less wrong than others (i.e. we don't know of any better model).

\subsection{$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 \Remph{fracDesign} of the qualityTools 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.

<<echo=TRUE, results=hide>>=
fdo.frac = fracDesign(k = 3, gen = "C = AB", centerCube = 4)
@

In order to get more specific information about a design the \Remph{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. 

<<echo=TRUE>>=
summary(fdo.frac)
@

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 \Remph{aliasTable} or \Remph{confounds} of the qualityTools package.

<<echo=TRUE>>=
aliasTable(fdo.frac)
@

The latter gives a more human readable version of the first and adds the resolution and generator(s)of the design.

<<echo=TRUE>>=
confounds(fdo.frac)
@

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 \Remph{fracChoose} of the qualityTools package and simply clicking onto the desired design (figure \ref{fig:designTable}).

<<fracChoose, echo=TRUE, results=hide>>=
fracChoose()
@

\begin{figure}[htbp]\begin{center}
<<fig=TRUE, echo=FALSE, results=hide, width=14, height=7>>=
<<fracChoose>>
@
\caption{Choosing minimum aberration designs}
\label{fig:designTable}
\end{center}\end{figure}

\subsection{Replicated Designs and Center Points}

A replicated design with additional center points can be created by using the \Remph{replicates} and \Remph{centerCube} argument.

<<echo=TRUE, results=hide>>=
fdo1 = facDesign(k = 3, centerCube = 2, replicates = 2)
@

\subsection{Multiple Responses}\label{subsec:multipleResponses}

Once you have observed the response for the different factor combinations one can add one or more response vectors to the design with the \Remph{response} method of the qualityTools package. A second response to be named y2 is created, filled with random numbers and put together in a data.frame with \Remph{data.frame}. The method \Remph{response} is used again to add these values to the factorial design object \Remph{fdo}.

<<echo=TRUE, results=hide>>=
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 \Remph{wirePlot} and \Remph{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 $\sim$ 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 (figure \ref{fig:multipleResponsePlot1}).

<<multiresp, echo=TRUE, results=hide>>=
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")
@

\begin{figure}[htbp]\begin{center}
<<fig=TRUE, echo=FALSE, results=hide, width=14, height=7>>=
<<multiresp>>
@
\caption{wire plot with different formulas specified}
\label{fig:multipleResponsePlot1}
\end{center}\end{figure}

Using the \Remph{wirePlot} and \Remph{contourPlot} methods of the qualityTools package settings of the other n-2 factors can be set using the \Remph{factors} argument. A wireplot with the third factor C on -1 an C = 1 can be created as follows (figure \ref{fig:multipleResponsePlot2})

<<multiresp2, echo=TRUE, results=hide>>=
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")
@

\begin{figure}[htbp]\begin{center}
<<fig=TRUE, echo=FALSE, results=hide, width=14, height=7>>=
<<multiresp2>>
@
\caption{wire plot with formula and setting for factor C}
\label{fig:multipleResponsePlot2}
\end{center}\end{figure}

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 \Remph{fits} method of the qualityTools package and is especially useful when working with more than one response (see \ref{subsec:multipleResponses}). Of course \Remph{lm} can be used to analyze the fractional factorial designs.

<<echo=TRUE>>=
fits(fdo) = lm(yield ~ A+B, data = fdo)
fits(fdo) = lm(y2 ~ A*B*C, data = fdo)
fits(fdo)
@

\subsection{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 (figure \ref{fig:wirePlot}) can be figured out graphically or calculated using the \Remph{steepAscent} method of the qualityTools package.

<<echo=TRUE>>=
sao =steepAscent(factors=c("A","B"),response="yield",data=fdo,steps=20)
sao
@

Since we set the real values earlier using the \Remph{highs} and \Remph{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\footnote{steepest ascent object} can be set using the \Remph{response} method of the qualityTools package and then be plotted using the \Remph{plot} method. Of course one can easily use the base plot method itself. However for documentation purposes the \Remph{plot} method for a steepest ascent object might be more convenient (see figure \ref{fig:sao}).

<<steepAsc, echo=TRUE, results=hide>>=
predicted = simProc(sao[,5], sao[,6])
response(sao) = predicted 
plot(sao, type = "b", col = 2)
@

\begin{figure}[htbp]\begin{center}
<<fig=TRUE, echo=FALSE, results=hide, width=14, height=7>>=
<<steepAsc>>
@
\caption{predicted maximum at Delta = 11 (see sao)}\label{fig:sao}
\label{fig:steepestAscent}
\end{center}\end{figure}

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

\subsection{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, \dots, 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 in \ref{subsec:facDesigns} started off with a $2^k$ design containing center points. Sticking to a linear model we used the \Remph{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

<<echo=TRUE, results=hide>>=
#set the seed for randomization of the runs
set.seed(1234)
fdo2 = facDesign(k = 2, centerCube = 3)
names(fdo2) = c("Factor 1", "Factor 2") 
lows(fdo2) = c(134, 155)
highs(fdo2) = c(155, 175)
@

and the yield is calculated by using the \Remph{simProc} and assigned to the design with the help of the generic \Remph{response} method of the qualityTools package.

<<echo=TRUE, results=hide>>=
yield = c(simProc(134,175), simProc(144.5,165.5), simProc(155,155),
simProc(144.5,165.5), simProc(155,175), simProc(144.5,165.5),
simProc(134,155))
response(fdo2) = yield
@

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 \Remph{starDesign} method of the qualityTools package. By default the value of alpha is chosen so that both criteria, \emph{orthogonality} and \emph{rotatability} are approximately met. Simply call the \Remph{starDesign} method on the factorial design object fdo2. Calling rsdo\footnote{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.

<<echo=TRUE>>=
rsdo = starDesign(data = fdo2)
rsdo
@

Using the \Remph{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 \Remph{simProc} method of the qualityTools package.

<<echo=TRUE, results=hide>>=
yield2 = c(yield, simProc(130,165), simProc(155,165), simProc(144,155),
simProc(144,179),simProc(144,165),simProc(144,165),simProc(144,165))
response(rsdo) = yield2
@

A full quadratic model is fitted using the \Remph{lm} method

<<echo=TRUE, results=hide>>=
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 \Remph{wirePlot} and \Remph{contourPlot} method of the qualityTools package.

<<respSur, echo=TRUE, results=hide>>=
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)
@

\begin{figure}[htbp]\begin{center}
<<fig=TRUE, echo=FALSE, results=hide, width=14, height=7>>=
<<respSur>>
@
\caption{quadratic fit of the response surface design object rsdo}
\label{fig:responseSurfaceAndContour}
\end{center}\end{figure}

Figure \ref{fig:wholeProcessRegion} can be used to compare the outcomes of the factorial and response surface designs with the simulated process. The inactive Factor 3 was omitted.

\begin{figure}[htbp]\begin{center}
<<fig=TRUE, echo=FALSE, results=hide, width=14, height=7>>=
A = seq(40, 210, length = 100)
B = seq(90, 190, length = 100)
C = seq(90, 190, length = 100)
filled.contour(A, B,outer(A,B, simProc, noise = FALSE), xlab = "Factor 1", ylab = "Factor 2", color.palette = colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan","#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000"))) 
@
\caption{underlying black box process without noise}
\label{fig:wholeProcessRegion}
\end{center}\end{figure}

Besides this sequential strategy, response surface designs can be created using the method \Remph{rsmDesign} of the qualityTools package. A design with \Remph{alpha = 1.633}, 0 centerpoints in the cube portion and 6 center points in the star portion can be created with:



<<echo=TRUE, results=hide>>=
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 \Remph{so=TRUE} (i.e. standard order). \Remph{cc} stands for centerCube and \Remph{cs} for centerStar.

<<echo=TRUE, results=hide>>=
fdo = randomize(fdo, so = TRUE)
@

Response Surface Designs can also be chosen from a table by using the method \Remph{rsmChoose} of the qualityTools package (see figure \ref{fig:designTable2}).

<<desTab, echo=TRUE, results=hide>>=
rsdo = rsmChoose()
@

\begin{figure}[htbp]\begin{center}
<<fig=TRUE, echo=FALSE, results=hide, width=14, height=7>>=
<<desTab>>
@
\caption{choosing a predefined response surface design from a table}
\label{fig:designTable2}
\end{center}\end{figure}

\subsubsection{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 \Remph{starDesign} method of the qualityTools package. A star portion consists of axial runs and optional center points (\Remph{cs}) in the axial part as opposed to center points (\Remph{cc}) in the cube part.

<<echo=TRUE, results=hide>>=
fdo3 = facDesign(k = 6)
rsdo = starDesign(alpha = "orthogonal", data = fdo3)
@

In case no existing (fractional) factorial design is handed to the \Remph{starDesign} method a list with \Remph{data.frames} is returned which can be assigned to the existing (fractional) factorial design using the \Remph{star}, \Remph{centerStar} and \Remph{centerCube} methods of the qualityTools package.

\subsubsection{Randomization}

Randomization is achieved by using the \Remph{randomize} method of the qualityTools package. At this point randomization works for most of the designs types. A \Remph{random.seed} needs to be supplied which is helpful to have the same run order on any machine.

<<echo=TRUE, results=hide>>=
randomize(fdo, random.seed = 123)
@

The \Remph{randomize} method can also be used to obtain a design in standard order with the help of the \Remph{so} argument.
<<echo=TRUE, results=hide>>=
randomize(fdo, so = TRUE)
@

\subsubsection{Blocking}
Blocking is another relevant feature and can be achieved by the \Remph{blocking} method of the qualityTools package. At this point blocking a design afterwards is not always successful. However, it is unproblematic during the sequential assembly.


\subsection{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.

\begin{align}
\sqrt[n]{\prod_{i=1}^n {d_i}}
\end{align}

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$ and $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 \Remph{desires} method of the qualityTools package. The optimization direction for each response variable is defined via the \Remph{min}, \Remph{max} and \Remph{target} argument of the \Remph{desires} method. The \Remph{target} argument is set with \Remph{max} for maximization, \Remph{min} for minimization and a specific value for optimization towards a specific target. Three settings arise from this constellation

\begin{description}
\item[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.
\item[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.
\item[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]
\end{description}

Besides these settings the scale factor influences the shape of the \Remph{desirability} method.
Desirability methods can be created and plotted using the \Remph{desires} and \Remph{plot} method of the qualityTools package. Desirabilities are always attached to a response and thus should be assigned to factorial designs (figure \ref{fig:desirabilities}).

<<echo=TRUE>>=
d1 = desirability(y1, 120, 170, scale = c(1,1), target = "max")
d3 = desirability(y3, 400, 600, target = 500)
d1
@

Besides having a summary on the command line, the \Remph{desirability} method can be conveniently visualized using the \Remph{plot} method. With the desirabilities d1 and d3 one gets the following plots.

<<desirb, echo=TRUE, results=hide>>=
par(mfrow = c(1,2))
plot(d1, col = 2); plot(d3, col = 2)
@

\begin{figure}[htbp]\begin{center}
<<fig=TRUE, echo=FALSE, results=hide, width=14, height=7>>=
<<desirb>>
@
\caption{plotted desirabilities for y1 and y3}
\label{fig:desirabilities}
\end{center}\end{figure}

\subsection{Using desirabilities together with designed experiments}

The \Remph{desirability} methodology is supported by the factorial design objects. The output of the \Remph{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 \Remph{rsmDesign} of the qualityTools package. Then we use the \Remph{randomize} method to obtain the standard order of the design.

<<echo=TRUE, results=hide>>=
ddo = rsmDesign(k = 3, alpha = 1.633, cc = 0, cs = 6)
ddo = randomize(ddo, so = TRUE)
#optional
names(ddo) = c("silica", "silan", "sulfur")
#optional
highs(ddo) = c(1.7, 60, 2.8)
#optional
lows(ddo) = c(0.7, 40, 1.8)
@

The \Remph{summary} method gives an overview of the design. The values of the responses are incorporated with the \Remph{response} method of the qualityTools package.

<<echo=TRUE, results=hide>>=
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 \Remph{data.frame} of these 4 responses is assigned to the design object ddo.

<<echo=TRUE, results=hide>>=
response(ddo) = data.frame(y1, y2, y3, y4)[c(5,2,3,8,1,6,7,4,9:20),]
@

The desirabilities are incorporated with the \Remph{desires} method of the qualityTools package. y1 and y3 were already defined which leaves the desirabailities for y2 and y4 to be defined.

<<echo=TRUE, results=hide>>=
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 \Remph{desires} method is used as follows.

<<echo=TRUE, results=hide>>=
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.

<<echo=TRUE, results=hide>>=
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 \Remph{optimum} of the qualityTools package giving the same factor settings as stated in \cite{Derringer.1980} for an overall desirability of 0.58 and individual desirabilities of 0.187, 1, 0.664, 0.934 for y1, y2, y3 and y4.

<<echo=TRUE>>=
optimum(ddo, type = "optim")
@

\subsection{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 \Remph{mixDesign} can be used to e.g. create simplex lattice design and simplex centroid designs. The generic methods \Remph{response}, \Remph{names}, \Remph{highs}, \Remph{lows}, \Remph{units} and \Remph{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 \Remph{mixDesign} of the qualityTools package. mdo is an abbreviation of mix design object.

<<echo=TRUE, results=hide>>=
mdo = mixDesign(3,2, center = FALSE, axial = FALSE, randomize = FALSE,
replicates  = c(1,1,2,3))
names(mdo) = c("polyethylene", "polystyrene", "polypropylene")

#set response (i.e. yarn elongation)
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
@

Again the values of the response are associated with the method \Remph{response} of the qualityTools package. Calling mdo prints the design. The generic \Remph{summary} method can be used for a more detailed overview.

<<echo=TRUE>>=
mdo
@

The data can be visualized using the \Remph{wirePlot3} and \Remph{contourPlot3} methods (figure \ref{fig:response3AndContour3}). In addition to the \Remph{wirePlot} and \Remph{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 \Remph{form} argument with a formula.

<<respAcont, echo=TRUE, results=hide>>=
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)
@

\begin{figure}[htbp]\begin{center}
<<fig=TRUE, echo=FALSE, results=hide, width=14, height=7>>=
<<respAcont>>
@
\caption{ternary plots for the elongation example}
\label{fig:response3AndContour3}
\end{center}\end{figure}

\subsection{Taguchi Designs}

Taguchi Designs are available using the method \Remph{taguchiDesign} of the qualityTools package. There are two types of taguchi designs:

\begin{itemize}
\item Single level:	all factors have the same number of levels (e.g. two levels for a L4\_2)
\item Mixed level:	factors have different number of levels (e.g. two and three levels for L18\_2\_3)
\end{itemize}

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 \Remph{taguchiDesign} method of the qualityTools package. The generic method \Remph{names}, \Remph{units}, \Remph{values}, \Remph{summary}, \Remph{plot}, \Remph{lm} and other methods again are supported. This way the relevant information for each factor can be stored in the design object tdo\footnote{taguchi design object} itself.

<<echo=TRUE>>=
set.seed(1234)
tdo = taguchiDesign("L9_3")
values(tdo) = list(A  = c(20, 40, 60), B = c("mateial 1", "material 2",
"material 3"), C = c(1,2,3))
names(tdo) = c("Factor 1", "Factor 2", "Factor 3", "Factor 4")
summary(tdo)
@

The \Remph{response} method is used to assign the values of the response variables. \Remph{effectPlot} can be used once more to visualize the effect sizes of the factors (figure \ref{fig:taguchiEffectPlot}).

<<taguchi, echo=TRUE, results=hide>>=
response(tdo) = rnorm(9)
effectPlot(tdo, col = 2)
@

\begin{figure}[htbp]\begin{center}
<<fig=TRUE, echo=FALSE, results=hide, width=14, height=7>>=
<<taguchi>>
@
\caption{effect plot for the taguchi experiment}
\label{fig:taguchiEffectPlot}
\end{center}\end{figure}


\section{Web Application for the qualityTools package}
In some cases it is not convenient to install sosftware and sometimes it is simply not possible. In order to make use of the qualityTools package in these cases, methods that are routinely performed, subject to reports and especially associated to the different phases of Six Sigma Projects are also provided as a web application under \url{http://webapps.r-qualitytools.org}. Methods include:

\begin{itemize}
\item Statistical Process Control (SPC),
\begin{itemize}
\item \href{http://webapps.r-qualitytools.org/brew/PCR/pcr.html}{Process Capability Indices},
\item \href{http://webapps.r-qualitytools.org/brew/QCC/qcc.html}{Quality Control Charts}, 
\end{itemize}

\item Measurement Systems Analysis,
\begin{itemize}
\item \href{http://webapps.r-qualitytools.org/brew/Gage/gageRR.html}{Gage Reapeatability \& Reproducibility Studies (Gage R\&R)} according to \citet{MSA.2010},
\item \href{http://webapps.r-qualitytools.org/brew/CG/cg.html}{Gage Capability Indices} according to VDA 5,
\end{itemize}

\item Design of Experiments,
\begin{itemize}
\item \href{http://webapps.r-qualitytools.org/brew/facDesign/fracDesign.html}{Fractional Factorial Designs / Screening Designs},
\item \href{http://webapps.r-qualitytools.org/brew/facDesign/facDesign.html}{Full Factorial Designs},
\item \href{http://webapps.r-qualitytools.org/brew/facDesign/rsmDesign.html}{Response Surface Designs} including Desirabilities,
\end{itemize}
\end{itemize}

as well as other methods such as \href{http://webapps.r-qualitytools.org/brew/tTest/tTest.html}{t tests}, \href{http://webapps.r-qualitytools.org/brew/qqPlot/qqPlot.html}{quantile-quantile plots} and \href{http://webapps.r-qualitytools.org/brew/powerT/powerT.html}{power calculations} for one and two sample t tests. 

The web application builds upon the rApache project \citep{Horner.2010} which embeds the R interpreter inside the apache webserver and gives access to these methods within a web application.\footnote{For more information on rApache visit the  rApache website under \url{http://rapache.net/index.html}}. The generation of reports is accomplished by using \href{http://www.statistik.lmu.de/~leisch/Sweave/}{Sweave} \citep{Leisch.2003}. Data can be imported from .csv and .xls files, whole sessions can be saved and restored via \href{http://en.wikipedia.org/wiki/JSON}{JSON}.

\section{Session Information}
The version number of R and packages loaded for generating the vignette were:
<<echo=TRUE>>=
sessionInfo()
@

\section{R-Code in this Vignette}
All of the R-Code used in this vignette can be found in the RCode.R file.

<<echo=FALSE, results=hide,split=TRUE>>= 
Stangle(file="qualityTools.rnw", output="RCode.r", annotate=FALSE) 
@

% <<echo=TRUE, eval=FALSE>>=
% source("RCode.r")
% @


%\lstinputlisting[emptylines=0]{RCode.R} 



\clearpage

\appendix
%\input{anhang}

\clearpage


%einkommentieren fr Literaturverzeichnis
%\bibliographystyle{alphadin}
\bibliographystyle{plainnat}
\bibliography{vignette}

\nocite{*} %damit alle und nicht nur im Text zitierte Literaturangaben angegeben werden
\end{document}