library(BsMD)## LenthPlot function from this package
library (FrF2)## FrF2 function is provided by this package
## Loading required package: DoE.base
## Loading required package: grid
## Loading required package: conf.design
## Registered S3 method overwritten by 'DoE.base':
## method from
## factorize.factor conf.design
##
## Attaching package: 'DoE.base'
## The following objects are masked from 'package:stats':
##
## aov, lm
## The following object is masked from 'package:graphics':
##
## plot.design
## The following object is masked from 'package:base':
##
## lengths
##
## Attaching package: 'FrF2'
## The following object is masked from 'package:BsMD':
##
## DanielPlot
qqnorm.lenth <- function (effects)
{
abseffects<- abs(effects)
m <- length (abseffects)
# Psuedo sd based all abs(effect)
s0<-1.5*median(abseffects)
non.sig <- abseffects<2.5*s0
# Refined estimate of sd
pse<-1.5*median(abseffects[non.sig])
sme<-qt(1- (1-0.95^(1/m))/2,m/3)*pse;
sig <- abseffects>sme
hqq<-qqnorm(effects, type = "n")
text(hqq$x, hqq$y, labels = names (effects), col = sig + 1)
qqline(effects)
cat ("Significant Factors Selected by Lenth's Method:\n",
names(abseffects)[sig], "\n")
names(effects)[sig]
}
## creating full factorial with 32 runs ...
## A B C D E
## 1 1 1 1 -1 1
## 2 1 -1 -1 -1 1
## 3 1 1 1 -1 -1
## 4 1 1 1 1 -1
## 5 1 1 1 1 1
## 6 -1 -1 -1 -1 -1
## 7 -1 -1 -1 -1 1
## 8 -1 1 1 -1 1
## 9 -1 -1 1 1 -1
## 10 -1 -1 -1 1 1
## 11 -1 1 1 1 -1
## 12 -1 -1 1 -1 1
## 13 1 -1 -1 1 1
## 14 -1 -1 -1 1 -1
## 15 1 -1 -1 -1 -1
## 16 -1 -1 1 -1 -1
## 17 1 -1 1 -1 1
## 18 -1 1 -1 1 -1
## 19 -1 -1 1 1 1
## 20 -1 1 -1 1 1
## 21 1 1 -1 -1 1
## 22 1 -1 -1 1 -1
## 23 -1 1 -1 -1 -1
## 24 1 -1 1 -1 -1
## 25 1 1 -1 1 1
## 26 1 1 -1 -1 -1
## 27 -1 1 -1 -1 1
## 28 -1 1 1 -1 -1
## 29 1 -1 1 1 1
## 30 -1 1 1 1 1
## 31 1 1 -1 1 -1
## 32 1 -1 1 1 -1
## class=design, type= full factorial
## look at the full design matrix
model.des1 <- data.frame( model.matrix(~A*B*C*D*E, des1))
model.des1 [with(model.des1, order (-A1.B1.C1.D1.E1, D1, C1, B1, A1)),
c("A1","B1", "C1", "D1", "E1", "A1.B1.C1.D1", "A1.C1.D1.E1",
"A1.B1.C1", "D1.E1", "A1.B1.C1.D1.E1" ) ]
## A1 B1 C1 D1 E1 A1.B1.C1.D1 A1.C1.D1.E1 A1.B1.C1 D1.E1 A1.B1.C1.D1.E1
## 7 -1 -1 -1 -1 1 1 -1 -1 -1 1
## 15 1 -1 -1 -1 -1 -1 -1 1 1 1
## 23 -1 1 -1 -1 -1 -1 1 1 1 1
## 21 1 1 -1 -1 1 1 1 -1 -1 1
## 16 -1 -1 1 -1 -1 -1 -1 1 1 1
## 17 1 -1 1 -1 1 1 -1 -1 -1 1
## 8 -1 1 1 -1 1 1 1 -1 -1 1
## 3 1 1 1 -1 -1 -1 1 1 1 1
## 14 -1 -1 -1 1 -1 -1 -1 -1 -1 1
## 13 1 -1 -1 1 1 1 -1 1 1 1
## 20 -1 1 -1 1 1 1 1 1 1 1
## 31 1 1 -1 1 -1 -1 1 -1 -1 1
## 19 -1 -1 1 1 1 1 -1 1 1 1
## 32 1 -1 1 1 -1 -1 -1 -1 -1 1
## 11 -1 1 1 1 -1 -1 1 -1 -1 1
## 5 1 1 1 1 1 1 1 1 1 1
## 6 -1 -1 -1 -1 -1 1 1 -1 1 -1
## 2 1 -1 -1 -1 1 -1 1 1 -1 -1
## 27 -1 1 -1 -1 1 -1 -1 1 -1 -1
## 26 1 1 -1 -1 -1 1 -1 -1 1 -1
## 12 -1 -1 1 -1 1 -1 1 1 -1 -1
## 24 1 -1 1 -1 -1 1 1 -1 1 -1
## 28 -1 1 1 -1 -1 1 -1 -1 1 -1
## 1 1 1 1 -1 1 -1 -1 1 -1 -1
## 10 -1 -1 -1 1 1 -1 1 -1 1 -1
## 22 1 -1 -1 1 -1 1 1 1 -1 -1
## 18 -1 1 -1 1 -1 1 -1 1 -1 -1
## 25 1 1 -1 1 1 -1 -1 -1 1 -1
## 9 -1 -1 1 1 -1 1 1 1 -1 -1
## 29 1 -1 1 1 1 -1 1 -1 1 -1
## 30 -1 1 1 1 1 -1 -1 -1 1 -1
## 4 1 1 1 1 -1 1 -1 1 -1 -1
des2 <- FrF2 (nfactors = 5, nruns=16, generators = c("ABCD"), alias.info = 3)
des2 [with (des2, order (D,C,B,A)), ]
## A B C D E
## 13 -1 -1 -1 -1 1
## 11 1 -1 -1 -1 -1
## 14 -1 1 -1 -1 -1
## 10 1 1 -1 -1 1
## 15 -1 -1 1 -1 -1
## 4 1 -1 1 -1 1
## 5 -1 1 1 -1 1
## 2 1 1 1 -1 -1
## 9 -1 -1 -1 1 -1
## 3 1 -1 -1 1 1
## 1 -1 1 -1 1 1
## 7 1 1 -1 1 -1
## 8 -1 -1 1 1 1
## 6 1 -1 1 1 -1
## 16 -1 1 1 1 -1
## 12 1 1 1 1 1
## class=design, type= FrF2.generators
## $legend
## [1] "A=A" "B=B" "C=C" "D=D" "E=E"
##
## $main
## character(0)
##
## $fi2
## [1] "AB=CDE" "AC=BDE" "AD=BCE" "AE=BCD" "BC=ADE" "BD=ACE" "BE=ACD" "CD=ABE"
## [9] "CE=ABD" "DE=ABC"
##
## $fi3
## character(0)
Aliase Structure: E = ABCD, A = BCDE, B=ACDE, C = ABDE, D = ABCE,
AE = BCD, BE = ACD, CE = ABD, AB = CDE, AC = BDE, AD = BCE, BC = ADE, BD = ACE, CD = ABE, that is, all third order interactions are confounded with 2nd order interactions. Fourth-order interactions are also confounded with 1st-order main effect.
des3 <- FrF2 (nfactors = 5, nruns=16, generators = c("-ABCD"), alias.info = 3)
des3[with (des3, order (D,C,B,A)), ]
## A B C D E
## 5 -1 -1 -1 -1 -1
## 7 1 -1 -1 -1 1
## 1 -1 1 -1 -1 1
## 14 1 1 -1 -1 -1
## 2 -1 -1 1 -1 1
## 11 1 -1 1 -1 -1
## 3 -1 1 1 -1 -1
## 15 1 1 1 -1 1
## 10 -1 -1 -1 1 1
## 4 1 -1 -1 1 -1
## 6 -1 1 -1 1 -1
## 9 1 1 -1 1 1
## 13 -1 -1 1 1 -1
## 16 1 -1 1 1 1
## 12 -1 1 1 1 1
## 8 1 1 1 1 -1
## class=design, type= FrF2.generators
## $legend
## [1] "A=A" "B=B" "C=C" "D=D" "E=E"
##
## $main
## character(0)
##
## $fi2
## [1] "AB=-CDE" "AC=-BDE" "AD=-BCE" "AE=-BCD" "BC=-ADE" "BD=-ACE" "BE=-ACD"
## [8] "CD=-ABE" "CE=-ABD" "DE=-ABC"
##
## $fi3
## character(0)
Aliase Structure: E = -ABCD, A = -BCDE, B = -ACDE, C = -ABDE, D = -ABCE,
AE = -BCD, BE = -ACD, CE = -ABD, AB = -CDE, AC = -BDE, AD = -BCE, BC = - ADE, BD = - ACE, CD = - ABE, that is, all third order interactions are confounded with 2nd order interactions. Fourth-order interactions are also confounded with 1st-order main effect.
A \(2^{5-1}\) design used for process improvement Five factors in a manufacturing process for an integrated circuit were investigated in a \(2^{5−1}\) design with the objective of improving the process yield. The five factors.
## y
## 1 8
## 2 9
## 3 34
## 4 52
## 5 16
## 6 22
## 7 45
## 8 60
## 9 6
## 10 10
## 11 30
## 12 50
## 13 15
## 14 21
## 15 44
## 16 63
# define the design matrix
A<-rep(c(-1,1),8)
B<-rep(c(-1,-1,1,1),4)
C<-rep(c(rep(-1,4),rep(1,4)),2)
D<-c(rep(-1,8),rep(1,8))
E<-A*B*C*D
# alternatively, we can generate these design with FrF2:
des2 <- FrF2 (nfactors = 5, nruns=16, generators = c("ABCD"), alias.info = 3)
des2 [with (des2, order (D,C,B,A)), ]
## A B C D E
## 3 -1 -1 -1 -1 1
## 11 1 -1 -1 -1 -1
## 7 -1 1 -1 -1 -1
## 5 1 1 -1 -1 1
## 15 -1 -1 1 -1 -1
## 14 1 -1 1 -1 1
## 1 -1 1 1 -1 1
## 8 1 1 1 -1 -1
## 12 -1 -1 -1 1 -1
## 4 1 -1 -1 1 1
## 9 -1 1 -1 1 1
## 6 1 1 -1 1 -1
## 2 -1 -1 1 1 1
## 10 1 -1 1 1 -1
## 13 -1 1 1 1 -1
## 16 1 1 1 1 1
## class=design, type= FrF2.generators
## $legend
## [1] "A=A" "B=B" "C=C" "D=D" "E=E"
##
## $main
## character(0)
##
## $fi2
## [1] "AB=CDE" "AC=BDE" "AD=BCE" "AE=BCD" "BC=ADE" "BD=ACE" "BE=ACD" "CD=ABE"
## [9] "CE=ABD" "DE=ABC"
##
## $fi3
## character(0)
##
## Call:
## lm.default(formula = y ~ (A + B + C + D + E)^2, data = process)
##
## Residuals:
## ALL 16 residuals are 0: no residual degrees of freedom!
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 30.3125 NA NA NA
## A 5.5625 NA NA NA
## B 16.9375 NA NA NA
## C 5.4375 NA NA NA
## D -0.4375 NA NA NA
## E 0.3125 NA NA NA
## A:B 3.4375 NA NA NA
## A:C 0.1875 NA NA NA
## A:D 0.5625 NA NA NA
## A:E 0.5625 NA NA NA
## B:C 0.3125 NA NA NA
## B:D -0.0625 NA NA NA
## B:E -0.0625 NA NA NA
## C:D 0.4375 NA NA NA
## C:E 0.1875 NA NA NA
## D:E -0.6875 NA NA NA
##
## Residual standard error: NaN on 0 degrees of freedom
## Multiple R-squared: 1, Adjusted R-squared: NaN
## F-statistic: NaN on 15 and 0 DF, p-value: NA
## Warning in anova.lm(g): ANOVA F-tests on an essentially perfect fit are
## unreliable
## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## A 1 495.1 495.1
## B 1 4590.1 4590.1
## C 1 473.1 473.1
## D 1 3.1 3.1
## E 1 1.6 1.6
## A:B 1 189.1 189.1
## A:C 1 0.6 0.6
## A:D 1 5.1 5.1
## A:E 1 5.1 5.1
## B:C 1 1.6 1.6
## B:D 1 0.1 0.1
## B:E 1 0.1 0.1
## C:D 1 3.1 3.1
## C:E 1 0.6 0.6
## D:E 1 7.6 7.6
## Residuals 0 0.0
## Significant Factors Selected by Lenth's Method:
## A B C A:B
## [1] "A" "B" "C" "A:B"
## alpha PSE ME SME
## 0.050000 0.468750 1.204960 2.446243
## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## A 1 495.1 495.1 193.195 2.535e-08 ***
## B 1 4590.1 4590.1 1791.244 1.560e-13 ***
## C 1 473.1 473.1 184.610 3.214e-08 ***
## A:B 1 189.1 189.1 73.781 3.302e-06 ***
## Residuals 11 28.2 2.6
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Call:
## lm.default(formula = y ~ A + B + C + A * B, data = process)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.812 -0.875 0.125 1.188 2.188
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 30.3125 0.4002 75.74 2.64e-16 ***
## A 5.5625 0.4002 13.90 2.53e-08 ***
## B 16.9375 0.4002 42.32 1.56e-13 ***
## C 5.4375 0.4002 13.59 3.21e-08 ***
## A:B 3.4375 0.4002 8.59 3.30e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.601 on 11 degrees of freedom
## Multiple R-squared: 0.9951, Adjusted R-squared: 0.9933
## F-statistic: 560.7 on 4 and 11 DF, p-value: 1.252e-12
des5<- FrF2 (nfactors = 6, nruns=16, generators = c("ABC", "BCD"),
alias.info = 3);
des5[with (des5,order(D,C,B,A)),]
## A B C D E F
## 4 -1 -1 -1 -1 -1 -1
## 16 1 -1 -1 -1 1 -1
## 1 -1 1 -1 -1 1 1
## 12 1 1 -1 -1 -1 1
## 7 -1 -1 1 -1 1 1
## 3 1 -1 1 -1 -1 1
## 9 -1 1 1 -1 -1 -1
## 2 1 1 1 -1 1 -1
## 15 -1 -1 -1 1 -1 1
## 11 1 -1 -1 1 1 1
## 14 -1 1 -1 1 1 -1
## 13 1 1 -1 1 -1 -1
## 8 -1 -1 1 1 1 -1
## 5 1 -1 1 1 -1 -1
## 6 -1 1 1 1 -1 1
## 10 1 1 1 1 1 1
## class=design, type= FrF2.generators
## $legend
## [1] "A=A" "B=B" "C=C" "D=D" "E=E" "F=F"
##
## $main
## [1] "A=BCE=DEF" "B=ACE=CDF" "C=ABE=BDF" "D=AEF=BCF" "E=ABC=ADF" "F=ADE=BCD"
##
## $fi2
## [1] "AB=CE" "AC=BE" "AD=EF" "AE=BC=DF" "AF=DE" "BD=CF" "BF=CD"
##
## $fi3
## [1] "ABD=ACF=BEF=CDE" "ABF=ACD=BDE=CEF"
## creating full factorial with 64 runs ...
model_des6 <- data.frame(model.matrix (~A*B*C*D*E*F, des6))
model_des6[ with(model_des6, order (-B1.C1.D1.F1, -A1.B1.C1.E1, D1, C1, B1, A1)),
c("A1","B1", "C1","D1", "A1.B1.C1.E1", "B1.C1.D1.F1")]
## A1 B1 C1 D1 A1.B1.C1.E1 B1.C1.D1.F1
## 45 -1 -1 -1 -1 1 1
## 28 1 -1 -1 -1 1 1
## 58 -1 1 -1 -1 1 1
## 3 1 1 -1 -1 1 1
## 60 -1 -1 1 -1 1 1
## 57 1 -1 1 -1 1 1
## 24 -1 1 1 -1 1 1
## 39 1 1 1 -1 1 1
## 19 -1 -1 -1 1 1 1
## 12 1 -1 -1 1 1 1
## 33 -1 1 -1 1 1 1
## 7 1 1 -1 1 1 1
## 14 -1 -1 1 1 1 1
## 9 1 -1 1 1 1 1
## 55 -1 1 1 1 1 1
## 35 1 1 1 1 1 1
## 1 -1 -1 -1 -1 -1 1
## 44 1 -1 -1 -1 -1 1
## 2 -1 1 -1 -1 -1 1
## 36 1 1 -1 -1 -1 1
## 43 -1 -1 1 -1 -1 1
## 40 1 -1 1 -1 -1 1
## 11 -1 1 1 -1 -1 1
## 20 1 1 1 -1 -1 1
## 5 -1 -1 -1 1 -1 1
## 8 1 -1 -1 1 -1 1
## 31 -1 1 -1 1 -1 1
## 51 1 1 -1 1 -1 1
## 38 -1 -1 1 1 -1 1
## 18 1 -1 1 1 -1 1
## 29 -1 1 1 1 -1 1
## 56 1 1 1 1 -1 1
## 4 -1 -1 -1 -1 1 -1
## 52 1 -1 -1 -1 1 -1
## 37 -1 1 -1 -1 1 -1
## 23 1 1 -1 -1 1 -1
## 63 -1 -1 1 -1 1 -1
## 27 1 -1 1 -1 1 -1
## 50 -1 1 1 -1 1 -1
## 26 1 1 1 -1 1 -1
## 15 -1 -1 -1 1 1 -1
## 49 1 -1 -1 1 1 -1
## 22 -1 1 -1 1 1 -1
## 32 1 1 -1 1 1 -1
## 64 -1 -1 1 1 1 -1
## 17 1 -1 1 1 1 -1
## 21 -1 1 1 1 1 -1
## 54 1 1 1 1 1 -1
## 59 -1 -1 -1 -1 -1 -1
## 62 1 -1 -1 -1 -1 -1
## 30 -1 1 -1 -1 -1 -1
## 16 1 1 -1 -1 -1 -1
## 25 -1 -1 1 -1 -1 -1
## 48 1 -1 1 -1 -1 -1
## 6 -1 1 1 -1 -1 -1
## 47 1 1 1 -1 -1 -1
## 53 -1 -1 -1 1 -1 -1
## 41 1 -1 -1 1 -1 -1
## 13 -1 1 -1 1 -1 -1
## 61 1 1 -1 1 -1 -1
## 10 -1 -1 1 1 -1 -1
## 46 1 -1 1 1 -1 -1
## 42 -1 1 1 1 -1 -1
## 34 1 1 1 1 -1 -1
# read in the data
shrinkage<- read.table("data/shrinkage.txt", header=T)
knitr::kable(shrinkage, align="c")
y |
---|
6 |
10 |
32 |
60 |
4 |
15 |
26 |
60 |
8 |
12 |
34 |
60 |
16 |
5 |
37 |
52 |
attach(shrinkage)
# define the design matrix
A<-rep(c(-1,1),8)
B<-rep(c(-1,-1,1,1),4)
C<-rep(c(rep(-1,4),rep(1,4)),2)
D<-c(rep(-1,8),rep(1,8))
E<-A*B*C
F<-B*C*D
g <- lm(y~(A+B+C+D+E+F)^2+A*B*D+A*B*F,data=shrinkage)
anova(g)
## Warning in anova.lm(g): ANOVA F-tests on an essentially perfect fit are
## unreliable
## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## A 1 770.1 770.1
## B 1 5076.6 5076.6
## C 1 3.1 3.1
## D 1 7.6 7.6
## E 1 0.6 0.6
## F 1 0.6 0.6
## A:B 1 564.1 564.1
## A:C 1 10.6 10.6
## A:D 1 115.6 115.6
## A:E 1 14.1 14.1
## A:F 1 1.6 1.6
## B:D 1 0.1 0.1
## B:F 1 0.1 0.1
## A:B:D 1 0.1 0.1
## A:B:F 1 95.1 95.1
## Residuals 0 0.0
##
## Call:
## lm.default(formula = y ~ (A + B + C + D + E + F)^2 + A * B *
## D + A * B * F, data = shrinkage)
##
## Residuals:
## ALL 16 residuals are 0: no residual degrees of freedom!
##
## Coefficients: (8 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 27.3125 NA NA NA
## A 6.9375 NA NA NA
## B 17.8125 NA NA NA
## C -0.4375 NA NA NA
## D 0.6875 NA NA NA
## E 0.1875 NA NA NA
## F 0.1875 NA NA NA
## A:B 5.9375 NA NA NA
## A:C -0.8125 NA NA NA
## A:D -2.6875 NA NA NA
## A:E -0.9375 NA NA NA
## A:F 0.3125 NA NA NA
## B:C NA NA NA NA
## B:D -0.0625 NA NA NA
## B:E NA NA NA NA
## B:F -0.0625 NA NA NA
## C:D NA NA NA NA
## C:E NA NA NA NA
## C:F NA NA NA NA
## D:E NA NA NA NA
## D:F NA NA NA NA
## E:F NA NA NA NA
## A:B:D 0.0625 NA NA NA
## A:B:F -2.4375 NA NA NA
##
## Residual standard error: NaN on 0 degrees of freedom
## Multiple R-squared: 1, Adjusted R-squared: NaN
## F-statistic: NaN on 15 and 0 DF, p-value: NA
## Significant Factors Selected by Lenth's Method:
## A B A:B A:D
## [1] "A" "B" "A:B" "A:D"
## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## A 1 770.1 770.1 55.1961 3.978e-05 ***
## B 1 5076.6 5076.6 363.8751 1.378e-08 ***
## D 1 7.6 7.6 0.5421 0.4803288
## A:B 1 564.1 564.1 40.4306 0.0001315 ***
## A:D 1 115.6 115.6 8.2832 0.0182356 *
## B:D 1 0.1 0.1 0.0045 0.9480994
## Residuals 9 125.6 14.0
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## look at the fitting results
par (mfrow = c(2,2))
plot (y~factor(A))
plot (y~factor(B))
interaction.plot(factor (A), factor (B), y, ylim = range (y))
points (factor(A), y, pch = B + 2)
interaction.plot(factor (A), factor (D), y, ylim = range (y))
points (factor(A), y, pch = D + 2)