qualityTools-package: Statistical Methods for Quality Sciences

Description Details Note Author(s) Examples

Description

This Package contains methods associated with the (D)efine (M)easure (A)nalyze (I)mprove and (C) ontrol (i.e. DMAIC) cycle of the Six Sigma Quality Management methodology.

  1. Define: Pareto Chart

  2. Measure: Probability and Quantile-Quantile Plots, Process Capability Ratios for various distributions and Gage R&R

  3. Analyze: Pareto Chart, Multi-Vari Chart, Dot Plot

  4. Improve: Full and fractional factorial, response surface and mixture 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 etc.

  5. Control: Quality Control Charts can be found in the qcc package

Details

Package: qualityTools
Type: Package
Version: 0.96.2
Date: 2012-07-01
URL: http://r-qualitytools.org
License: GPL-2
LazyLoad: yes

Note

This package is primarily used for teaching! The package vignette is available under http://www.r-qualitytools.org.

Author(s)

Thomas Roth [email protected]
Maintainer: Thomas Roth [email protected]

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19

Example output

Loading required package: Rsolnp
Loading required package: MASS

Attaching package: 'qualityTools'

The following object is masked from 'package:stats':

    sigma


prtChr> #artifical defects dataset
prtChr> defects = LETTERS[runif(80, 1, 10)]

prtChr> #paretoChart of the defects
prtChr> paretoChart(defects) 
                                                                      
Frequency          14    11    10    10     8     8     7     6      6
Cum. Frequency     14    25    35    45    53    61    68    74     80
Percentage      17.5% 13.8% 12.5% 12.5% 10.0% 10.0%  8.8%  7.5%   7.5%
Cum. Percentage 17.5% 31.2% 43.8% 56.2% 66.2% 76.2% 85.0% 92.5% 100.0%

                                                                   
Frequency       14.0 11.00 10.00 10.00  8.00  8.00  7.00  6.0   6.0
Cum. Frequency  14.0 25.00 35.00 45.00 53.00 61.00 68.00 74.0  80.0
Percentage      17.5 13.75 12.50 12.50 10.00 10.00  8.75  7.5   7.5
Cum. Percentage 17.5 31.25 43.75 56.25 66.25 76.25 85.00 92.5 100.0

mvPlot> #Example I
mvPlot> examp1 = expand.grid(c("Engine1","Engine2","Engine3"),c(10,20,30,40))                  

mvPlot> examp1 = as.data.frame(rbind(examp1, examp1, examp1))

mvPlot> examp1 = cbind(examp1, rnorm(36, 1, 0.02))

mvPlot> names(examp1) = c("factor1", "factor2", "response")

mvPlot> test1=mvPlot(response = examp1[,3], fac1 = examp1[,2],
mvPlot+              fac2 = examp1[,1],sort=FALSE,las=2,FUN=mean) 

mvPlot> #Example II
mvPlot> examp2=expand.grid(c("Op I","Op II","Op III"),c(1,2,3,4),
mvPlot+                    c("20.11.1987","21.11.1987"))

mvPlot> examp2=as.data.frame(rbind(examp2, examp2, examp2))

mvPlot> examp2=cbind(examp2, rnorm(72, 22, 2))

mvPlot> names(examp2) = c("factor1", "factor2", "factor3", "response")

mvPlot> test2=mvPlot(response = examp2[,4], fac1 = examp2[,1],
mvPlot+             fac2 = examp2[,2], fac3 = examp2[,3], sort=TRUE, FUN=mean, 
mvPlot+             labels=TRUE)

mvPlot> #Example III
mvPlot> examp3 = expand.grid(c("A","B","C"),c("I","II","III","IV"),c("H","I"),
mvPlot+                      c(1,2,3,4,5))

mvPlot> examp3 = as.data.frame(rbind(examp3, examp3, examp3))

mvPlot> examp3 = cbind(examp3, rnorm(360, 0, 2))

mvPlot> names(examp3) = c("factor1", "factor2", "factor3", "factor4", "response")

mvPlot> test3=mvPlot(response = examp3[,5], fac1 = examp3[,1],
mvPlot+              fac2 = examp3[,2], fac3 = examp3[,3], fac4 = examp3[,4], sort=TRUE, 
mvPlot+              quantile=TRUE, FUN=mean)

dotPlt> #create some data and grouping
dotPlt> x = rnorm(28)

dotPlt> g = rep(1:2, 14)

dotPlt> #dot plot with groups and no stacking
dotPlt> dotPlot(x, group = g, stacked = FALSE, pch = c(19, 20), 
dotPlt+         main = "Non stacked dot plot")

dotPlt> #dot plot with groups and stacking
dotPlt> x = rnorm(28)

dotPlt> dotPlot(x, group = g, stacked = TRUE, pch = c(19, 20), 
dotPlt+         main = "Stacked dot plot")

qqPlot> #set up the plotting window for 6 plots
qqPlot> par(mfrow = c(3,2))

qqPlot> #generate random data from weibull distribution
qqPlot> x = rweibull(20, 8, 2)

qqPlot> #Quantile-Quantile Plot for different distributions
qqPlot> qqPlot(x, "log-normal")

qqPlot> qqPlot(x, "normal")

qqPlot> qqPlot(x, "exponential", DB = TRUE)

qqPlot> qqPlot(x, "cauchy")

qqPlot> qqPlot(x, "weibull")

qqPlot> qqPlot(x, "logistic")        
Warning messages:
1: In densfun(x, parm[1], parm[2], ...) : NaNs produced
2: In densfun(x, parm[1], parm[2], ...) : NaNs produced
3: In densfun(x, parm[1], parm[2], ...) : NaNs produced
4: In densfun(x, parm[1], parm[2], ...) : NaNs produced
5: In densfun(x, parm[1], parm[2], ...) : NaNs produced
6: In densfun(x, parm[1], parm[2], ...) : NaNs produced

ppPlot> #set up the plotting window for 6 plots
ppPlot> par(mfrow = c(3,2))

ppPlot> #generate random data from weibull distribution
ppPlot> x = rweibull(20, 8, 2)

ppPlot> #Probability Plot for different distributions
ppPlot> ppPlot(x, "log-normal")
$meanlog
[1] 0.6193474

$sdlog
[1] 0.1378872

 [1] 0.2989610 0.6739502 0.9417645 1.1447345 1.2993570 1.4145314 1.4957853
 [8] 1.5468246 1.5702385 1.5678634 1.5409836 1.4904432 1.4167020 1.3198464
[15] 1.1995559 1.0550049 0.8846501 0.6857536 0.4531260 0.1741306

ppPlot> ppPlot(x, "normal")
$mean
[1] 1.874788

$sd
[1] 0.2449665

 [1] 0.2385840 0.5778563 0.8403335 1.0522880 1.2243003 1.3621126 1.4692402
 [8] 1.5479476 1.5996943 1.6253600 1.6253600 1.5996943 1.5479476 1.4692402
[15] 1.3621126 1.2243003 1.0522880 0.8403335 0.5778563 0.2385840

ppPlot> ppPlot(x, "exponential", DB = TRUE)
$rate
[1] 0.5333935

 [1] 0.52005869 0.49338902 0.46671934 0.44004966 0.41337999 0.38671031
 [7] 0.36004063 0.33337096 0.30670128 0.28003160 0.25336193 0.22669225
[13] 0.20002257 0.17335290 0.14668322 0.12001354 0.09334387 0.06667419
[19] 0.04000451 0.01333484

ppPlot> ppPlot(x, "cauchy")
$location
[1] 1.9406

$scale
[1] 0.1321547

 [1] 0.01482702 0.13126166 0.35273350 0.65756333 1.01591229 1.39270269
 [7] 1.75105165 2.05588148 2.27735332 2.39378795 2.39378795 2.27735332
[13] 2.05588148 1.75105165 1.39270269 1.01591229 0.65756333 0.35273350
[19] 0.13126166 0.01482702

ppPlot> ppPlot(x, "weibull")
$shape
[1] 9.022527

$scale
[1] 1.978806

 [1] 0.1691647 0.4362787 0.6659340 0.8686838 1.0480426 1.2054946 1.3415846
 [8] 1.4562917 1.5491788 1.6194442 1.6659197 1.6870235 1.6806657 1.6440838
[15] 1.5735652 1.4639542 1.3077046 1.0928044 0.7971100 0.3638550

ppPlot> ppPlot(x, "logistic")        
$location
[1] 1.890717

$scale
[1] 0.1374921

 [1] 0.1772829 0.5045743 0.7955001 1.0500601 1.2682544 1.4500830 1.5955459
 [8] 1.7046430 1.7773745 1.8137402 1.8137402 1.7773745 1.7046430 1.5955459
[15] 1.4500830 1.2682544 1.0500601 0.7955001 0.5045743 0.1772829
Warning messages:
1: In densfun(x, parm[1], parm[2], ...) : NaNs produced
2: In densfun(x, parm[1], parm[2], ...) : NaNs produced
3: In densfun(x, parm[1], parm[2], ...) : NaNs produced
4: In densfun(x, parm[1], parm[2], ...) : NaNs produced
5: In densfun(x, parm[1], parm[2], ...) : NaNs produced
6: In densfun(x, parm[1], parm[2], ...) : NaNs produced
7: In densfun(x, parm[1], parm[2], ...) : NaNs produced

pcr> x = rweibull(30, 2, 8) +100

pcr> #process capability for a weibull distribution
pcr> pcr(x, "weibull", lsl = 100, usl = 117)

	Anderson Darling Test for weibull distribution

data:  x 
A = 0.9503, shape = 30.604, scale = 108.499, p-value <= 0.01
alternative hypothesis: true distribution is not equal to weibull 


pcr> #box cox transformation and one sided
pcr> pcr(x, boxcox = TRUE, lsl = 1)

	Anderson Darling Test for normal distribution

data:  x 
A = 0.4198, mean = 0, sd = 0, p-value = 0.3058
alternative hypothesis: true distribution is not equal to normal 


pcr> #boxcox with lambda=2
pcr> pcr(x, boxcox = 2, lsl = 1)

	Anderson Darling Test for normal distribution

data:  x 
A = 0.4198, mean = 0, sd = 0, p-value = 0.3058
alternative hypothesis: true distribution is not equal to normal 


pcr> #process capability assuming a normal distribution
pcr> pcr(x, "normal", lsl = 0, usl = 17)

	Anderson Darling Test for normal distribution

data:  x 
A = 0.4117, mean = 106.859, sd = 3.320, p-value = 0.32
alternative hypothesis: true distribution is not equal to normal 


pcr> #process capability for a normal distribution and data in subgroups
pcr> #some artificial data with shifted means in subgroups
pcr> x = c(rnorm(5, mean = 1), rnorm(5, mean = 2), rnorm(5, mean = 0))

pcr> #grouping vector
pcr> group = c(rep(1,5), rep(2,5), rep(3,5))

pcr> #calculate process capability
pcr> pcr(x, grouping = group) #compare to sd(x)

	Anderson Darling Test for normal distribution

data:  x 
A = 0.279, mean = 0.756, sd = 1.320, p-value = 0.5949
alternative hypothesis: true distribution is not equal to normal 


gageRR> #create a crossed Gage R&R Design
gageRR> gdo = gageRRDesign(3,10, 2, randomize = FALSE)

gageRR> #set the response i.e. Measurements
gageRR> y = 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,
gageRR+       22,22,24,23,22,24,20,20,25,24,22,24,21,20,21,22,21,22,21,21,24,27,25,27,
gageRR+       23,22,25,23,23,22,22,23,25,21,24,23)

gageRR> response(gdo) = y

gageRR> #perform a Gage R&R
gageRR> gdo = gageRR(gdo, tolerance = 5)

AnOVa Table -  crossed Design
              Df Sum Sq Mean Sq F value   Pr(>F)    
Operator       2  20.63  10.317   8.597  0.00112 ** 
Part           9 107.07  11.896   9.914 7.31e-07 ***
Operator:Part 18  22.03   1.224   1.020  0.46732    
Residuals     30  36.00   1.200                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

----------
AnOVa Table Without Interaction -  crossed Design
            Df Sum Sq Mean Sq F value   Pr(>F)    
Operator     2  20.63  10.317   8.533 0.000675 ***
Part         9 107.07  11.896   9.840 2.39e-08 ***
Residuals   48  58.03   1.209                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

----------

Gage R&R
                 VarComp VarCompContrib Stdev StudyVar StudyVarContrib
totalRR            1.664          0.483 1.290     7.74           0.695
 repeatability     1.209          0.351 1.100     6.60           0.592
 reproducibility   0.455          0.132 0.675     4.05           0.364
   Operator        0.455          0.132 0.675     4.05           0.364
   Operator:Part   0.000          0.000 0.000     0.00           0.000
Part to Part       1.781          0.517 1.335     8.01           0.719
totalVar           3.446          1.000 1.856    11.14           1.000
                 P/T Ratio
totalRR               1.55
 repeatability        1.32
 reproducibility      0.81
   Operator           0.81
   Operator:Part      0.00
Part to Part          1.60
totalVar              2.23

---
 * Contrib equals Contribution in %
 **Number of Distinct Categories (truncated signal-to-noise-ratio) = 1 


gageRR> #summary
gageRR> summary(gdo)

Operators:	 3 	Parts:	 10
Measurements:	 2 	Total:	 60
----------

AnOVa Table -  crossed Design
              Df Sum Sq Mean Sq F value   Pr(>F)    
Operator       2  20.63  10.317   8.597  0.00112 ** 
Part           9 107.07  11.896   9.914 7.31e-07 ***
Operator:Part 18  22.03   1.224   1.020  0.46732    
Residuals     30  36.00   1.200                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

----------
AnOVa Table Without Interaction -  crossed Design
            Df Sum Sq Mean Sq F value   Pr(>F)    
Operator     2  20.63  10.317   8.533 0.000675 ***
Part         9 107.07  11.896   9.840 2.39e-08 ***
Residuals   48  58.03   1.209                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

----------

Gage R&R
                 VarComp VarCompContrib Stdev StudyVar StudyVarContrib
totalRR            1.664          0.483 1.290     7.74           0.695
 repeatability     1.209          0.351 1.100     6.60           0.592
 reproducibility   0.455          0.132 0.675     4.05           0.364
   Operator        0.455          0.132 0.675     4.05           0.364
   Operator:Part   0.000          0.000 0.000     0.00           0.000
Part to Part       1.781          0.517 1.335     8.01           0.719
totalVar           3.446          1.000 1.856    11.14           1.000
                 P/T Ratio
totalRR               1.55
 repeatability        1.32
 reproducibility      0.81
   Operator           0.81
   Operator:Part      0.00
Part to Part          1.60
totalVar              2.23

---
 * Contrib equals Contribution in %
 **Number of Distinct Categories (truncated signal-to-noise-ratio) = 1 


gageRR> #standard graphics for Gage R&R
gageRR> plot(gdo)

gageRR> ##create a crossed Gage R&R Design - 
gageRR> ##Vardeman, VanValkenburg 1999 - Two-Way Random-Effects Analyses and Gauge
gageRR> #gdo = gageRRDesign(Operators = 5, Parts = 2, Measurements = 3, randomize = FALSE)
gageRR> #
gageRR> ##Measurements
gageRR> #weight = c(3.481, 3.448, 3.485, 3.475, 3.472,
gageRR> #           3.258, 3.254, 3.256, 3.249, 3.241,
gageRR> #           3.477, 3.472, 3.464, 3.472, 3.470,
gageRR> #           3.254, 3.247, 3.257, 3.238, 3.250,
gageRR> #           3.470, 3.470, 3.477, 3.473, 3.474,
gageRR> #           3.258, 3.239, 3.245, 3.240, 3.254)
gageRR> #
gageRR> ##set the response i.e. Measurements
gageRR> #response(gdo) = weight
gageRR> #
gageRR> ##perform a Gage R&R
gageRR> #gdo = gageRR(gdo)
gageRR> #
gageRR> ##summary
gageRR> #summary(gdo)
gageRR> #
gageRR> ##standard graphics for Gage R&R
gageRR> #plot(gdo)
gageRR> #
gageRR> 
gageRR> 
gageRR> 
gageRR> 

fcDsgn> #returns a 2^3 full factorial design
fcDsgn> vp.full = facDesign(k = 3)        

fcDsgn> #generate some random response                        
fcDsgn> response(vp.full) = rnorm(2^3) 

fcDsgn> #summary of the full factorial design (especially no defining relation)                            
fcDsgn> summary(vp.full)                                           
Information about the factors:

           A       B       C
low       -1      -1      -1
high       1       1       1
name                        
unit                        
type numeric numeric numeric
-----------
  StandOrder RunOrder Block  A  B  C  rnorm.2.3.
4          4        1     1  1  1 -1 -0.91149314
3          3        2     1 -1  1 -1 -0.64604795
2          2        3     1  1 -1 -1 -1.21969321
5          5        4     1 -1 -1  1  0.29408749
6          6        5     1  1 -1  1 -0.53747162
7          7        6     1 -1  1  1 -0.04947574
8          8        7     1  1  1  1 -0.27433356
1          1        8     1 -1 -1 -1  0.24014825

fcDsgn> #------------
fcDsgn> 
fcDsgn> #returns a full factorial design with 3 replications per factor combination 
fcDsgn> #and 4 center points
fcDsgn> vp.rep = facDesign(k = 2, replicates = 3, centerCube = 4)  

fcDsgn> #set names
fcDsgn> names(vp.rep) = c("Name 1", "Name 2") 

fcDsgn> #set units                     
fcDsgn> units(vp.rep) = c("min", "F")         

fcDsgn> #set low and high factor values                     
fcDsgn> lows(vp.rep) = c(20, 40, 60)                               

fcDsgn> highs(vp.rep) = c(40, 60, 80)  

fcDsgn> #summary of the replicated full factorial Design        
fcDsgn> summary(vp.rep)                                                
Information about the factors:

           A       B
low       20      40
high      40      60
name  Name 1  Name 2
unit     min       F
type numeric numeric
-----------
   StandOrd RunOrder Block  A  B  y
4         4        1     1  1  1 NA
9         9        2     1 -1 -1 NA
3         3        3     1 -1  1 NA
14       14        4     1  0  0 NA
10       10        5     1  1 -1 NA
2         2        6     1  1 -1 NA
12       12        7     1  1  1 NA
11       11        8     1 -1  1 NA
15       15        9     1  0  0 NA
5         5       10     1 -1 -1 NA
16       16       11     1  0  0 NA
8         8       12     1  1  1 NA
13       13       13     1  0  0 NA
6         6       14     1  1 -1 NA
1         1       15     1 -1 -1 NA
7         7       16     1 -1  1 NA
Warning messages:
1: In `[<-`(`*tmp*`, i, value = <S4 object of class "doeFactor">) :
  implicit list embedding of S4 objects is deprecated
2: In `[<-`(`*tmp*`, i, value = <S4 object of class "doeFactor">) :
  implicit list embedding of S4 objects is deprecated
3: In `[<-`(`*tmp*`, i, value = <S4 object of class "doeFactor">) :
  implicit list embedding of S4 objects is deprecated
4: In `[<-`(`*tmp*`, i, value = new("doeFactor")) :
  implicit list embedding of S4 objects is deprecated
5: In `[<-`(`*tmp*`, i, value = new("doeFactor")) :
  implicit list embedding of S4 objects is deprecated

frcDsg> #returns a 2^3 full factorial design
frcDsg> vp.full = facDesign(k = 3)       

frcDsg> #design in 2 blocks                           
frcDsg> vp.full = blocking(vp.full, 2)   

frcDsg> #generate some random response                           
frcDsg> response(vp.full) = rnorm(2^3)   

frcDsg> #summary of the full factorial design (especially no defining relation)                           
frcDsg> summary(vp.full)                                           
Information about the factors:

           A       B       C
low       -1      -1      -1
high       1       1       1
name                        
unit                        
type numeric numeric numeric
-----------
  StandOrder RunOrder Block  A  B  C  rnorm.2.3.
6          6        1     1  1 -1  1 -0.91149314
4          4        2     1  1  1 -1 -0.64604795
7          7        3     1 -1  1  1 -1.21969321
1          1        4     1 -1 -1 -1  0.29408749
3          3        5     2 -1  1 -1 -0.53747162
5          5        6     2 -1 -1  1 -0.04947574
8          8        7     2  1  1  1 -0.27433356
2          2        8     2  1 -1 -1  0.24014825

frcDsg> #returns a 2^4-1 fractional factorial design. Factor D will be aliased with
frcDsg> vp.frac = fracDesign(k = 4, gen = "D=ABC") 

frcDsg> #the three-way-interaction ABC (i.e. I = ABCD)                 
frcDsg> response(vp.frac) = rnorm(2^(4-1))    

frcDsg> #summary of the fractional factorial design                      
frcDsg> summary(vp.frac)                                            
Information about the factors:

           A       B       C       D
low       -1      -1      -1      -1
high       1       1       1       1
name                                
unit                                
type numeric numeric numeric numeric
-----------
  StandOrder RunOrder Block  A  B  C  D rnorm.2..4...1..
4          4        1     1  1  1 -1 -1      -0.91149314
3          3        2     1 -1  1 -1  1      -0.64604795
2          2        3     1  1 -1 -1  1      -1.21969321
5          5        4     1 -1 -1  1  1       0.29408749
6          6        5     1  1 -1  1 -1      -0.53747162
7          7        6     1 -1  1  1 -1      -0.04947574
8          8        7     1  1  1  1  1      -0.27433356
1          1        8     1 -1 -1 -1 -1       0.24014825

---------

Defining relations:
I =  ABCD 		Columns: 1 2 3 4 

Resolution:  IV 


frcDsg> #returns a full factorial design with 3 replications per factor combination 
frcDsg> #and 4 center points
frcDsg> vp.rep = fracDesign(k = 3, replicates = 3, centerCube = 4)  

frcDsg> #summary of the replicated fractional factorial Design
frcDsg> summary(vp.rep)                                             
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
16       16        1     1  1  1  1 NA
4         4        2     1  1  1 -1 NA
14       14        3     1  1 -1  1 NA
9         9        4     1 -1 -1 -1 NA
3         3        5     1 -1  1 -1 NA
18       18        6     1  1 -1 -1 NA
24       24        7     1  1  1  1 NA
20       20        8     1  1  1 -1 NA
15       15        9     1 -1  1  1 NA
2         2       10     1  1 -1 -1 NA
28       28       11     1  0  0  0 NA
8         8       12     1  1  1  1 NA
7         7       13     1 -1  1  1 NA
23       23       14     1 -1  1  1 NA
22       22       15     1  1 -1  1 NA
12       12       16     1  1  1 -1 NA
25       25       17     1  0  0  0 NA
19       19       18     1 -1  1 -1 NA
5         5       19     1 -1 -1  1 NA
17       17       20     1 -1 -1 -1 NA
21       21       21     1 -1 -1  1 NA
10       10       22     1  1 -1 -1 NA
11       11       23     1 -1  1 -1 NA
27       27       24     1  0  0  0 NA
13       13       25     1 -1 -1  1 NA
6         6       26     1  1 -1  1 NA
1         1       27     1 -1 -1 -1 NA
26       26       28     1  0  0  0 NA
Warning messages:
1: In `[<-`(`*tmp*`, i, value = new("doeFactor")) :
  implicit list embedding of S4 objects is deprecated
2: In `[<-`(`*tmp*`, i, value = new("doeFactor")) :
  implicit list embedding of S4 objects is deprecated
3: In `[<-`(`*tmp*`, i, value = new("doeFactor")) :
  implicit list embedding of S4 objects is deprecated
4: In `[<-`(`*tmp*`, i, value = new("doeFactor")) :
  implicit list embedding of S4 objects is deprecated
5: In `[<-`(`*tmp*`, i, value = new("doeFactor")) :
  implicit list embedding of S4 objects is deprecated
6: In `[<-`(`*tmp*`, i, value = new("doeFactor")) :
  implicit list embedding of S4 objects is deprecated
7: In `[<-`(`*tmp*`, i, value = new("doeFactor")) :
  implicit list embedding of S4 objects is deprecated
8: In `[<-`(`*tmp*`, i, value = new("doeFactor")) :
  implicit list embedding of S4 objects is deprecated
9: In `[<-`(`*tmp*`, i, value = new("doeFactor")) :
  implicit list embedding of S4 objects is deprecated
10: In `[<-`(`*tmp*`, i, value = new("doeFactor")) :
  implicit list embedding of S4 objects is deprecated

effctP> #effectPlot for a 2^k factorial design
effctP> fdo = facDesign(k = 3)

effctP> #set response with generic response function
effctP> response(fdo) = rnorm(8)  

effctP> effectPlot(fdo)

effctP> #effectPlot for a taguchiDesign
effctP> tdo = taguchiDesign("L9_3")

effctP> response(tdo) = rnorm(9)

effctP> effectPlot(tdo, points = TRUE, col = 2, pch = 16, lty = 3)
Warning messages:
1: In `[<-`(`*tmp*`, i, value = new("doeFactor")) :
  implicit list embedding of S4 objects is deprecated
2: In `[<-`(`*tmp*`, i, value = new("doeFactor")) :
  implicit list embedding of S4 objects is deprecated
3: In `[<-`(`*tmp*`, i, value = new("doeFactor")) :
  implicit list embedding of S4 objects is deprecated
4: In `[<-`(`*tmp*`, i, value = <S4 object of class "taguchiFactor">) :
  implicit list embedding of S4 objects is deprecated
5: In `[<-`(`*tmp*`, i, value = <S4 object of class "taguchiFactor">) :
  implicit list embedding of S4 objects is deprecated
6: In `[<-`(`*tmp*`, i, value = <S4 object of class "taguchiFactor">) :
  implicit list embedding of S4 objects is deprecated
7: In `[<-`(`*tmp*`, i, value = <S4 object of class "taguchiFactor">) :
  implicit list embedding of S4 objects is deprecated

tgchDs> tdo = taguchiDesign("L9_3")

tgchDs> values(tdo) = list(A = c("material 1","material 2","material 3"), 
tgchDs+                    B = c(29, 30, 35))

tgchDs> names(tdo) = c("Factors", "Are", "Documented", "In The Design")

tgchDs> response(tdo) = rnorm(9)

tgchDs> summary(tdo)
Taguchi SINGLE Design
Information about the factors:

                 A       B          C             D
value 1 material 1      29          1             1
value 2 material 2      30          2             2
value 3 material 3      35          3             3
name       Factors     Are Documented In The Design
unit                                               
type       numeric numeric    numeric       numeric

-----------

  StandOrder RunOrder Replicate A B C D   rnorm(9)
1          9        1         1 3 3 2 1  0.4101212
2          7        2         1 3 1 3 2 -0.8929399
3          8        3         1 3 2 1 3 -1.7214333
4          2        4         1 1 2 2 2  1.3475997
5          6        5         1 2 3 1 2  0.5335563
6          5        6         1 2 2 3 1 -0.2017748
7          3        7         1 1 3 3 3  0.4606256
8          4        8         1 2 1 2 3 -0.2115619
9          1        9         1 1 1 1 1 -0.6531924

-----------


tgchDs> effectPlot(tdo)
Warning messages:
1: In `[<-`(`*tmp*`, i, value = new("taguchiFactor")) :
  implicit list embedding of S4 objects is deprecated
2: In `[<-`(`*tmp*`, i, value = new("taguchiFactor")) :
  implicit list embedding of S4 objects is deprecated
3: In `[<-`(`*tmp*`, i, value = new("taguchiFactor")) :
  implicit list embedding of S4 objects is deprecated
4: In `[<-`(`*tmp*`, i, value = new("taguchiFactor")) :
  implicit list embedding of S4 objects is deprecated

rsmDsg> #central composite design for 2 factors with 2 blocks, alpha = 1.41, 
rsmDsg> #5 centerpoints in the cube portion and 3 centerpoints in the star portion:
rsmDsg> rsmDesign(k = 2, blocks = 2, alpha = sqrt(2),cc = 5, cs = 3)
   StandOrd RunOrder Block      A      B  y
4         4        1     1  1.000  1.000 NA
3         3        2     1 -1.000  1.000 NA
2         2        3     1  1.000 -1.000 NA
5         5        4     1  0.000  0.000 NA
9         9        5     1  0.000  0.000 NA
6         6        6     1  0.000  0.000 NA
7         7        7     1  0.000  0.000 NA
8         8        8     1  0.000  0.000 NA
1         1        9     1 -1.000 -1.000 NA
13       13       10     2  0.000  1.414 NA
11       11       11     2  1.414  0.000 NA
16       16       12     2  0.000  0.000 NA
14       14       13     2  0.000  0.000 NA
10       10       14     2 -1.414  0.000 NA
15       15       15     2  0.000  0.000 NA
12       12       16     2  0.000 -1.414 NA

rsmDsg> #central composite design with both, orthogonality and near rotatability 
rsmDsg> rsmDesign(k = 2, blocks = 2, alpha = "both")
   StandOrd RunOrder Block      A      B  y
3         3        1     1 -1.000  1.000 NA
7         7        2     1  0.000  0.000 NA
2         2        3     1  1.000 -1.000 NA
6         6        4     1  0.000  0.000 NA
4         4        5     1  1.000  1.000 NA
5         5        6     1  0.000  0.000 NA
1         1        7     1 -1.000 -1.000 NA
13       13        8     2  0.000  0.000 NA
9         9        9     2  1.414  0.000 NA
12       12       10     2  0.000  0.000 NA
8         8       11     2 -1.414  0.000 NA
14       14       12     2  0.000  0.000 NA
11       11       13     2  0.000  1.414 NA
10       10       14     2  0.000 -1.414 NA

rsmDsg> #central composite design with
rsmDsg> #2 centerpoints in the factorial portion of the design i.e 2
rsmDsg> #1 centerpoint int the star portion of the design i.e. 1
rsmDsg> #2 replications per factorial point i.e. 2^3*2 = 16
rsmDsg> #3 replications per star points 3*2*3 = 18
rsmDsg> #makes a total of 37 factor combinations
rsmDsg> rsdo = rsmDesign(k = 3, blocks = 1, alpha = 2, cc = 2, cs = 1, fp = 2, sp = 3)

rsmDsg> nrow(rsdo) #37
[1] 37
Warning messages:
1: In `[<-`(`*tmp*`, i, value = new("doeFactor")) :
  implicit list embedding of S4 objects is deprecated
2: In `[<-`(`*tmp*`, i, value = new("doeFactor")) :
  implicit list embedding of S4 objects is deprecated
3: In `[<-`(`*tmp*`, i, value = new("doeFactor")) :
  implicit list embedding of S4 objects is deprecated
4: In `[<-`(`*tmp*`, i, value = new("doeFactor")) :
  implicit list embedding of S4 objects is deprecated
5: In `[<-`(`*tmp*`, i, value = new("doeFactor")) :
  implicit list embedding of S4 objects is deprecated
6: In `[<-`(`*tmp*`, i, value = new("doeFactor")) :
  implicit list embedding of S4 objects is deprecated
7: In `[<-`(`*tmp*`, i, value = new("doeFactor")) :
  implicit list embedding of S4 objects is deprecated

prtPlt> #factorial design with replications
prtPlt> #NA in response column and 2 replicates per factor combination
prtPlt> vp = fracDesign(k = 3, replicates = 2)  

prtPlt> #generate some data
prtPlt> y1 = 4*vp[,1] -7*vp[,2] + 2*vp[,2]*vp[,1] + 0.2*vp[,3] + rnorm(16)

prtPlt> y2 = 9*vp[,1] -2*vp[,2] + 5*vp[,2]*vp[,1] + 0.5*vp[,3] + rnorm(16)               

prtPlt> response(vp) = data.frame(y1,y2)

prtPlt> #show effects and interactions (nothing significant expected)
prtPlt> paretoPlot(vp)                         

prtPlt> #fractional factorial design --> Lenth Plot
prtPlt> vp = fracDesign(k = 4, gen = "D = ABC")

prtPlt> #generate some data
prtPlt> y = rnorm(8)                

prtPlt> response(vp) = y

prtPlt> #show effects and interactions (nothing significant expected)
prtPlt> paretoPlot(vp)                        
Warning messages:
1: In `[<-`(`*tmp*`, i, value = new("doeFactor")) :
  implicit list embedding of S4 objects is deprecated
2: In `[<-`(`*tmp*`, i, value = new("doeFactor")) :
  implicit list embedding of S4 objects is deprecated
3: In `[<-`(`*tmp*`, i, value = new("doeFactor")) :
  implicit list embedding of S4 objects is deprecated
4: In `[<-`(`*tmp*`, i, value = new("doeFactor")) :
  implicit list embedding of S4 objects is deprecated
5: In `[<-`(`*tmp*`, i, value = new("doeFactor")) :
  implicit list embedding of S4 objects is deprecated
6: In `[<-`(`*tmp*`, i, value = new("doeFactor")) :
  implicit list embedding of S4 objects is deprecated
7: In `[<-`(`*tmp*`, i, value = new("doeFactor")) :
  implicit list embedding of S4 objects is deprecated

wirPlt> #create a response surface design and assign random data to response y
wirPlt> fdo = rsmDesign(k = 3, blocks = 2)

wirPlt> response(fdo) = data.frame(y = rnorm(nrow(fdo)))

wirPlt> par(mfrow = c(3,2))

wirPlt> #I - display linear fit
wirPlt> wirePlot(A,B,y, data = fdo, form = "linear")

wirPlt> #II - display full fit (i.e. effect, interactions and quadratic effects
wirPlt> wirePlot(A,B,y, data = fdo, form = "full")

wirPlt> #III - display a fit specified before
wirPlt> fits(fdo) = lm(y ~ B + I(A^2) , data = fdo)

wirPlt> wirePlot(A,B,y, data = fdo, form = "fit")

wirPlt> #IV - display a fit given directly
wirPlt> wirePlot(A,B,y, data = fdo, form = "y ~ A*B + I(A^2)")

wirPlt> #V - display a fit using a different colorRamp
wirPlt> wirePlot(A,B,y, data = fdo, form = "full", col = 2)

wirPlt> #V - display a fit using a self defined colorRamp
wirPlt> myColour = colorRampPalette(c("green", "gray","blue"))

wirPlt> wirePlot(A,B,y, data = fdo, form = "full", col = myColour)
Warning messages:
1: In `[<-`(`*tmp*`, i, value = new("doeFactor")) :
  implicit list embedding of S4 objects is deprecated
2: In `[<-`(`*tmp*`, i, value = new("doeFactor")) :
  implicit list embedding of S4 objects is deprecated
3: In `[<-`(`*tmp*`, i, value = new("doeFactor")) :
  implicit list embedding of S4 objects is deprecated

cntrPl> #create a response surface design and assign random data to response y
cntrPl> fdo = rsmDesign(k = 3, blocks = 2)

cntrPl> response(fdo) = data.frame(y = rnorm(nrow(fdo)))

cntrPl> par(mfrow = c(2,3))

cntrPl> #I - display linear fit
cntrPl> contourPlot(A,B,y, data = fdo, form = "linear")

cntrPl> #II - display full fit (i.e. effect, interactions and quadratic effects
cntrPl> contourPlot(A,B,y, data = fdo, form = "full")

cntrPl> #III - display a fit specified before
cntrPl> fits(fdo) = lm(y ~ B + I(A^2) , data = fdo)

cntrPl> contourPlot(A,B,y, data = fdo, form = "fit")

cntrPl> #IV - display a fit given directly
cntrPl> contourPlot(A,B,y, data = fdo, form = "y ~ A*B + I(A^2)")

cntrPl> #V - display a fit using a different colorRamp
cntrPl> contourPlot(A,B,y, data = fdo, form = "full", col = 2)

cntrPl> #VI - display a fit using a self defined colorRamp
cntrPl> myColour = colorRampPalette(c("green", "gray","blue"))

cntrPl> contourPlot(A,B,y, data = fdo, form = "full", col = myColour)
Warning messages:
1: In `[<-`(`*tmp*`, i, value = new("doeFactor")) :
  implicit list embedding of S4 objects is deprecated
2: In `[<-`(`*tmp*`, i, value = new("doeFactor")) :
  implicit list embedding of S4 objects is deprecated
3: In `[<-`(`*tmp*`, i, value = new("doeFactor")) :
  implicit list embedding of S4 objects is deprecated

mxDsgn> #simplex lattice design with center (default response added with NA's)
mxDsgn> mixDesign(p = 3, n = 2, center = FALSE)
  StandOrder RunOrder    Type   A   B   C  y
1          1        6 2-blend 0.0 0.5 0.5 NA
2          2        3 1-blend 0.0 0.0 1.0 NA
3          3        5 2-blend 0.5 0.5 0.0 NA
4          4        4 1-blend 0.0 1.0 0.0 NA
5          5        1 2-blend 0.5 0.0 0.5 NA
6          6        2 1-blend 1.0 0.0 0.0 NA

mxDsgn> #simplex lattice design with a center (default response added with NA's)
mxDsgn> mixDesign(p = 3, n = 2, center = TRUE)
  StandOrder RunOrder    Type         A         B         C  y
1          1        4  center 0.3333333 0.3333333 0.3333333 NA
2          2        7 1-blend 0.0000000 1.0000000 0.0000000 NA
3          3        6 1-blend 0.0000000 0.0000000 1.0000000 NA
4          4        2 1-blend 1.0000000 0.0000000 0.0000000 NA
5          5        5 2-blend 0.0000000 0.5000000 0.5000000 NA
6          6        3 2-blend 0.5000000 0.0000000 0.5000000 NA
7          7        1 2-blend 0.5000000 0.5000000 0.0000000 NA

mxDsgn> #simplex lattice design with augmented points (default response added with NA's)
mxDsgn> mixDesign(p = 3, n = 2, center = FALSE, axial = TRUE)
  StandOrder RunOrder    Type         A         B         C  y
1          1        1 1-blend 1.0000000 0.0000000 0.0000000 NA
2          2        2 2-blend 0.5000000 0.5000000 0.0000000 NA
3          3        5 2-blend 0.0000000 0.5000000 0.5000000 NA
4          4        9 1-blend 0.0000000 0.0000000 1.0000000 NA
5          5        3 2-blend 0.5000000 0.0000000 0.5000000 NA
6          6        4   axial 0.6666667 0.1666667 0.1666667 NA
7          7        6   axial 0.1666667 0.1666667 0.6666667 NA
8          8        8   axial 0.1666667 0.6666667 0.1666667 NA
9          9        7 1-blend 0.0000000 1.0000000 0.0000000 NA

mxDsgn> #simplex lattice design with augmented points, center and 2 replications 
mxDsgn> #(default response added with NA's)
mxDsgn> mixDesign(p = 3, n = 2, center = TRUE, axial = TRUE, replicates = 2)
   StandOrder RunOrder    Type         A         B         C  y
1           1        9 2-blend 0.0000000 0.5000000 0.5000000 NA
2           2       14   axial 0.1666667 0.1666667 0.6666667 NA
3           3       13   axial 0.6666667 0.1666667 0.1666667 NA
4           4       18 1-blend 0.0000000 1.0000000 0.0000000 NA
5           5        7  center 0.3333333 0.3333333 0.3333333 NA
6           6       19 1-blend 0.0000000 1.0000000 0.0000000 NA
7           7        4 2-blend 0.5000000 0.0000000 0.5000000 NA
8           8        6   axial 0.1666667 0.6666667 0.1666667 NA
9           9       10 1-blend 1.0000000 0.0000000 0.0000000 NA
10         10        1 2-blend 0.0000000 0.5000000 0.5000000 NA
11         11       16   axial 0.1666667 0.1666667 0.6666667 NA
12         12       12 1-blend 0.0000000 0.0000000 1.0000000 NA
13         13        5 2-blend 0.5000000 0.5000000 0.0000000 NA
14         14       15 1-blend 1.0000000 0.0000000 0.0000000 NA
15         15       20  center 0.3333333 0.3333333 0.3333333 NA
16         16       17 1-blend 0.0000000 0.0000000 1.0000000 NA
17         17        2   axial 0.1666667 0.6666667 0.1666667 NA
18         18        3 2-blend 0.5000000 0.5000000 0.0000000 NA
19         19        8 2-blend 0.5000000 0.0000000 0.5000000 NA
20         20       11   axial 0.6666667 0.1666667 0.1666667 NA

mxDsgn> #simplex centroid design giving 2^(p-1) distinct design points 
mxDsgn> #(default response added with NA's)
mxDsgn> mixDesign(p = 3, n = 2, type = "centroid")
  StandOrder RunOrder    Type         A         B         C  y
1          1        7  center 0.3333333 0.3333333 0.3333333 NA
2          2        2 1-blend 0.0000000 1.0000000 0.0000000 NA
3          3        6 2-blend 0.5000000 0.0000000 0.5000000 NA
4          4        4 2-blend 0.0000000 0.5000000 0.5000000 NA
5          5        3 2-blend 0.5000000 0.5000000 0.0000000 NA
6          6        5 1-blend 1.0000000 0.0000000 0.0000000 NA
7          7        1 1-blend 0.0000000 0.0000000 1.0000000 NA

mxDsgn> #simplex centroid design with augmented points (default response added with NA's)
mxDsgn> mixDesign(p = 3, n = 2, type = "centroid", axial = TRUE)
   StandOrder RunOrder    Type         A         B         C  y
1           1        7   axial 0.6666667 0.1666667 0.1666667 NA
2           2        2 1-blend 0.0000000 1.0000000 0.0000000 NA
3           3        9 2-blend 0.0000000 0.5000000 0.5000000 NA
4           4        3   axial 0.1666667 0.6666667 0.1666667 NA
5           5       10 2-blend 0.5000000 0.5000000 0.0000000 NA
6           6        5  center 0.3333333 0.3333333 0.3333333 NA
7           7        6 1-blend 0.0000000 0.0000000 1.0000000 NA
8           8        1   axial 0.1666667 0.1666667 0.6666667 NA
9           9        4 1-blend 1.0000000 0.0000000 0.0000000 NA
10         10        8 2-blend 0.5000000 0.0000000 0.5000000 NA

mxDsgn> #yarn elongation example from Cornell (2002)
mxDsgn> mdo = mixDesign(3,2, center = FALSE, axial = FALSE, randomize = FALSE, 
mxDsgn+                 replicates  = c(1,1,2,3))

mxDsgn> #set names (optional)
mxDsgn> names(mdo) = c("polyethylene", "polystyrene", "polypropylene")

mxDsgn> #units(mdo) = "%" #set units (optional)
mxDsgn> #set response (i.e. yarn elongation)
mxDsgn> elongation = c(11.0, 12.4, 15.0, 14.8, 16.1, 17.7, 16.4, 16.6, 8.8, 10.0, 10.0,
mxDsgn+                9.7, 11.8, 16.8, 16.0)  

mxDsgn> response(mdo) = elongation
[1] "elongation"

mxDsgn> #print a summary of the design
mxDsgn> summary(mdo)
Simplex LATTICE Design
Information about the factors:

                A           B             C
low             0           0             0
high            1           1             1
name polyethylene polystyrene polypropylene
unit            %           %             %
type      numeric     numeric       numeric

-----------

Information about the Design Points:

           1-blend 2-blend
Unique           3       3
Replicates       2       3
Sub Total        6       9
Total           15        

-----------

Information about the constraints:

A >= 0 B >= 0 C >= 0

-----------

                              PseudoComponent _|_ Proportion _|_ Amount

   StandOrder RunOrder    Type |   A   B   C _ | _   A   B   C _ | _   A   B
1           1        1 1-blend | 1.0 0.0 0.0   |   1.0 0.0 0.0   |   1.0 0.0
2           2        2 1-blend | 1.0 0.0 0.0   |   1.0 0.0 0.0   |   1.0 0.0
3           3        3 2-blend | 0.5 0.5 0.0   |   0.5 0.5 0.0   |   0.5 0.5
4           4        4 2-blend | 0.5 0.5 0.0   |   0.5 0.5 0.0   |   0.5 0.5
5           5        5 2-blend | 0.5 0.5 0.0   |   0.5 0.5 0.0   |   0.5 0.5
6           6        6 2-blend | 0.5 0.0 0.5   |   0.5 0.0 0.5   |   0.5 0.0
7           7        7 2-blend | 0.5 0.0 0.5   |   0.5 0.0 0.5   |   0.5 0.0
8           8        8 2-blend | 0.5 0.0 0.5   |   0.5 0.0 0.5   |   0.5 0.0
9           9        9 1-blend | 0.0 1.0 0.0   |   0.0 1.0 0.0   |   0.0 1.0
10         10       10 1-blend | 0.0 1.0 0.0   |   0.0 1.0 0.0   |   0.0 1.0
11         11       11 2-blend | 0.0 0.5 0.5   |   0.0 0.5 0.5   |   0.0 0.5
12         12       12 2-blend | 0.0 0.5 0.5   |   0.0 0.5 0.5   |   0.0 0.5
13         13       13 2-blend | 0.0 0.5 0.5   |   0.0 0.5 0.5   |   0.0 0.5
14         14       14 1-blend | 0.0 0.0 1.0   |   0.0 0.0 1.0   |   0.0 0.0
15         15       15 1-blend | 0.0 0.0 1.0   |   0.0 0.0 1.0   |   0.0 0.0
     C | elongation
1  0.0 |       11.0
2  0.0 |       12.4
3  0.0 |       15.0
4  0.0 |       14.8
5  0.0 |       16.1
6  0.5 |       17.7
7  0.5 |       16.4
8  0.5 |       16.6
9  0.0 |        8.8
10 0.0 |       10.0
11 0.5 |       10.0
12 0.5 |        9.7
13 0.5 |       11.8
14 1.0 |       16.8
15 1.0 |       16.0

-----------

Mixture Total: 1 equals 1


mxDsgn> #show contourPlot and wirePlot
mxDsgn> dev.new(14, 14);par(mfrow = c(2,2))
dev.new(): using pdf(file="Rplots1.pdf")

mxDsgn> contourPlot3(A, B, C, elongation, data = mdo, form = "quadratic")

mxDsgn> wirePlot3(A, B, C, elongation, data = mdo, form = "quadratic", theta = 190)

mxDsgn> wirePlot3(A, B, C, elongation, data = mdo, form = "quadratic", phi = 390, 
mxDsgn+           theta = 0)

mxDsgn> wirePlot3(A, B, C, elongation, data = mdo, form = "quadratic", phi = 90)
There were 21 warnings (use warnings() to see them)

dsrblt> #Maximization of a response
dsrblt> #define a desirability for response y where higher values of y are better 
dsrblt> #as long as the response is smaller than high
dsrblt> d = desirability(y, low = 6, high = 18, target = "max")

dsrblt> #show and plot the desirability function
dsrblt> d; plot(d)
Target is to maximize y 
lower Bound:  6 
higher Bound:  18 
Scale factor is:  1 1 
importance:  1 


dsrblt> #Minimization of a response including a scaling factor
dsrblt> #define a desirability for response y where lower values of y are better 
dsrblt> #as long as the response is higher than low
dsrblt> d = desirability(y, low = 6, high = 18, scale = c(2),target = "min")

dsrblt> #show and plot the desirability function
dsrblt> d; plot(d)
Target is to minimize y 
lower Bound:  6 
higher Bound:  18 
Scale factor is:  2 
importance:  1 


dsrblt> #Specific target of a response is best including a scaling factor
dsrblt> #define a desirability for response y where desired value is at 8 
dsrblt> #and values lower than 6 as well as values higher than 18 are not acceptable
dsrblt> d = desirability(y, low = 6, high = 18, scale = c(0.5,2),target = 12)

dsrblt> #show and plot the desirability function
dsrblt> d; plot(d)
Target is  12  for y 
lower Bound:  6 
higher Bound:  18 
Scale factor is: low = 0.5 and high = 2 
importance:  1 


optimm> #BEWARE BIG EXAMPLE --Simultaneous Optimization of Several Response Variables--
optimm> #Source: DERRINGER, George; SUICH, Ronald: Simultaneous Optimization of Several 
optimm> #        Response Variables. Journal of Quality Technology Vol. 12, No. 4, 
optimm> #        p. 214-219
optimm> 
optimm> #Define the response suface design as given in the paper and sort via 
optimm> #Standard Order
optimm> fdo = rsmDesign(k = 3, alpha = 1.633, cc = 0, cs = 6)

optimm> fdo = randomize(fdo, so = TRUE)

optimm> #Attaching the 4 responses
optimm> y1 = c(102,120,117,198,103,132,132,139,102,154,96,163,116,153,133,133,140,142,
optimm+        145,142)

optimm> y2 = c(900,860,800,2294,490,1289,1270,1090,770,1690,700,1540,2184,1784,1300,
optimm+        1300,1145,1090,1260,1344)

optimm> y3 = c(470,410,570,240,640,270,410,380,590,260,520,380,520,290,380,380,430,
optimm+        430,390,390)

optimm> 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,
optimm+        70)

optimm> response(fdo) = data.frame(y1, y2, y3, y4)[c(5,2,3,8,1,6,7,4,9:20),]

optimm> #setting names and real values of the factors
optimm> names(fdo) = c("silica", "silan", "sulfur")

optimm> highs(fdo) = c(1.7, 60, 2.8)

optimm> lows(fdo) = c(0.7, 40, 1.8)

optimm> #summary of the response surface design
optimm> summary(fdo)
Information about the factors:

           A       B       C
low      0.7      40     1.8
high     1.7      60     2.8
name  silica   silan  sulfur
unit                        
type numeric numeric numeric
-----------
   StandOrd RunOrder Block      A      B      C  y1   y2  y3   y4
5         1        1     1 -1.000 -1.000 -1.000 103  490 640 62.5
2         2        2     1  1.000 -1.000 -1.000 120  860 410 65.0
3         3        3     1 -1.000  1.000 -1.000 117  800 570 77.5
8         4        4     1  1.000  1.000 -1.000 139 1090 380 70.0
1         5        5     1 -1.000 -1.000  1.000 102  900 470 67.5
6         6        6     1  1.000 -1.000  1.000 132 1289 270 67.0
7         7        7     1 -1.000  1.000  1.000 132 1270 410 78.0
4         8        8     1  1.000  1.000  1.000 198 2294 240 74.5
9         9        9     1 -1.633  0.000  0.000 102  770 590 76.0
10       10       10     1  1.633  0.000  0.000 154 1690 260 70.0
11       11       11     1  0.000 -1.633  0.000  96  700 520 63.0
12       12       12     1  0.000  1.633  0.000 163 1540 380 75.0
13       13       13     1  0.000  0.000 -1.633 116 2184 520 65.0
14       14       14     1  0.000  0.000  1.633 153 1784 290 71.0
15       15       15     1  0.000  0.000  0.000 133 1300 380 70.0
16       16       16     1  0.000  0.000  0.000 133 1300 380 68.5
17       17       17     1  0.000  0.000  0.000 140 1145 430 68.0
18       18       18     1  0.000  0.000  0.000 142 1090 430 68.0
19       19       19     1  0.000  0.000  0.000 145 1260 390 69.0
20       20       20     1  0.000  0.000  0.000 142 1344 390 70.0

optimm> #setting the desires
optimm> desires(fdo) = desirability(y1, 120, 170, scale = c(1,1), target = "max")

optimm> desires(fdo) = desirability(y2, 1000, 1300, target = "max")

optimm> desires(fdo) = desirability(y3, 400, 600, target = 500)

optimm> desires(fdo) = desirability(y4, 60, 75, target = 67.5)

optimm> desires(fdo)
$y1
Target is to maximize y1 
lower Bound:  120 
higher Bound:  170 
Scale factor is:  1 1 
importance:  1 


$y2
Target is to maximize y2 
lower Bound:  1000 
higher Bound:  1300 
Scale factor is:  1 1 
importance:  1 


$y3
Target is  500  for y3 
lower Bound:  400 
higher Bound:  600 
Scale factor is: low = 1 and high = 1 
importance:  1 


$y4
Target is  67.5  for y4 
lower Bound:  60 
higher Bound:  75 
Scale factor is: low = 1 and high = 1 
importance:  1 



optimm> #Have a look at some contourPlots
optimm> par(mfrow = c(2,2))

optimm> contourPlot(A, B, y1, data = fdo)

optimm> contourPlot(A, B, y2, data = fdo)

optimm> wirePlot(A, B, y1, data = fdo)

optimm> wirePlot(A, B, y2, data = fdo)

optimm> #setting the fits as in the paper
optimm> fits(fdo) = lm(y1 ~ A + B + C + A:B + A:C + B:C + I(A^2) + I(B^2) + I(C^2), 
optimm+                data = fdo)

optimm> fits(fdo) = lm(y2 ~ A + B + C + A:B + A:C + B:C + I(A^2) + I(B^2) + I(C^2), 
optimm+                data = fdo)

optimm> fits(fdo) = lm(y3 ~ A + B + C + A:B + A:C + B:C + I(A^2) + I(B^2) + I(C^2), 
optimm+                data = fdo)

optimm> fits(fdo) = lm(y4 ~ A + B + C + A:B + A:C + B:C + I(A^2) + I(B^2) + I(C^2),
optimm+                data = fdo)

optimm> #fits(fdo)
optimm> 
optimm> #calculate the same best factor settings as in the paper using type = "optim"
optimm> optimum(fdo, 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

optimm> #calculate (nearly) the same best factor settings as in the paper using type = "grid"
optimm> optimum(fdo, type = "grid")

composite (overall) desirability: 0.58

        A      B      C
coded 0.0  0.136 -0.817
real  1.2 51.361  1.892

                    y1       y2      y3     y4
Responses      130.660 1299.669 456.937 67.980
Desirabilities   0.213    0.999   0.569  0.936
Warning messages:
1: In `[<-`(`*tmp*`, i, value = new("doeFactor")) :
  implicit list embedding of S4 objects is deprecated
2: In `[<-`(`*tmp*`, i, value = new("doeFactor")) :
  implicit list embedding of S4 objects is deprecated
3: In `[<-`(`*tmp*`, i, value = new("doeFactor")) :
  implicit list embedding of S4 objects is deprecated

qualityTools documentation built on May 30, 2017, 1:43 a.m.