Menyusun Rancangan Faktorial 2^k-1

Misal akan disusun rancangan 2^5-1 dengan generator 𝐸 =𝐴𝐡𝐢𝐷

#Menggunakan fungsi FrF2
library(FrF2)
## Warning: package 'FrF2' was built under R version 4.0.5
## Loading required package: DoE.base
## Warning: package 'DoE.base' was built under R version 4.0.3
## Loading required package: grid
## Loading required package: conf.design
## Warning: package 'conf.design' was built under R version 4.0.3
## 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
design<-FrF2(16, 5, generators="ABCD", randomize=F)
design
##     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
## class=design, type= FrF2.generators

Dalam prakteknya, kombinasi perlakuan sebaiknya diacak. Hal ini dapat diatur dengan menggunakan argumen randomize=TRUE.

library(FrF2)
set.seed(1)
design<-FrF2(16, 5, generators="ABCD", randomize=T)
design
##     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
## class=design, type= FrF2.generators

Untuk menentukan alias, diperlukan vektor peubah respon y, ini dapat dimisalkan sebagai suatu vektor peubah dari sebaran uniform.

library(FrF2)
y <- runif(16, 0, 1)
aliases( lm( y~ (.)^4, data = design))
##             
##  A = B:C:D:E
##  B = A:C:D:E
##  C = A:B:D:E
##  D = A:B:C:E
##  E = A:B:C:D
##  A:B = C:D:E
##  A:C = B:D:E
##  A:D = B:C:E
##  A:E = B:C:D
##  B:C = A:D:E
##  B:D = A:C:E
##  B:E = A:C:D
##  C:D = A:B:E
##  C:E = A:B:D
##  D:E = A:B:C

Output di atas menunjukkan bahwa pada pola alias pada rancangan 2^5βˆ’1, pengaruh utama berbaur dengan interaksi 4-faktor, dan interaksi 2-faktor berbaur dengan interaksi 3-faktor. Oleh karenanya, jika interaksi 3-faktor dan 4-faktor diasumsikan dapat diabaikan, maka pendugaan untuk semua pengaruh utama dan interaksi 2-faktor dapat dilakukan.

Ilustrasi

Misalkan suatu percobaan untuk menemukan cara membuat soup dengan rincian faktor dan taraf sebagai berikut.

library(readxl)
## Warning: package 'readxl' was built under R version 4.0.4
tabel1 <- read_excel("D:/1. Kuliah S2 (STATISTIKA DAN SAINS DATA)/Materi Kuliah/SEMESTER 2/Materi kuliah/3. Analisis dan Perancangan Percobaan/STA521 (S2 NOW)/Pertemuan 12/tabel1.xlsx")
knitr::kable(tabel1)
Factor Label Name Low Level High Level
A Number of Ports 1 3
B Temperature Cooling Water Ambient
C Mixing Time 60 sec. 80 sec.
D Batch Weight 1500 lb 2000 lb
E Delay Days 7 1

Misalkan pada percobaan tersebut diperoleh data seperti pada tabel berikut ini.

library(readxl)
tabel2 <- read_excel("D:/1. Kuliah S2 (STATISTIKA DAN SAINS DATA)/Materi Kuliah/SEMESTER 2/Materi kuliah/3. Analisis dan Perancangan Percobaan/STA521 (S2 NOW)/Pertemuan 12/tabel2.xlsx")
knitr::kable(tabel2)
Random Run Order (A) Number of Ports (B) Temperture (C) Mixing Time (sec) (D) Batch Weight (lb) (E)Delay (days) Response
12 1 Cool Water 60 1500 1 1.13
13 3 Cool Water 60 1500 7 1.25
5 1 Ambient 60 1500 7 0.97
3 3 Ambient 60 1500 1 1.70
6 1 Cool Water 80 1500 7 1.47
4 3 Cool Water 80 1500 1 1.28
16 1 Ambient 80 1500 1 1.18
14 3 Ambient 80 1500 7 0.98
1 1 Cool Water 60 2000 7 0.78
15 3 Cool Water 60 2000 1 1.36
7 1 Ambient 60 2000 1 1.85
10 3 Ambient 60 2000 7 0.62
11 1 Cool Water 80 2000 1 1.09
2 3 Cool Water 80 2000 7 1.10
9 1 Ambient 80 2000 7 0.76
8 3 Ambient 80 2000 1 2.10

R code berikut ini menunjukkan penggunaan fungsi add.response pada package DoE.base untuk menambahkan peubah respon. Setelah itu, dilakukan pemodelan (mod1) menggunakan fungsi lm.

library(FrF2)
soup <- FrF2(16, 5, generators = "ABCD",
 factor.names =list(Ports=c(1,3),
 Temp=c("Cool","Ambient"),
 MixTime=c(60,80),
 BatchWt=c(1500,2000),
 delay=c(7,1)), randomize = FALSE)
y <- c(1.13, 1.25, .97, 1.70, 1.47, 1.28, 1.18, .98, .78, 1.36, 1.85, .62, 1.09, 1.10, .76, 2.10 )
library(DoE.base)
soup <- add.response( soup , y )
mod1 <- lm( y ~ (.)^2, data = soup)
summary(mod1)
## 
## Call:
## lm.default(formula = y ~ (.)^2, data = soup)
## 
## Residuals:
## ALL 16 residuals are 0: no residual degrees of freedom!
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)
## (Intercept)        1.22625         NA      NA       NA
## Ports1             0.07250         NA      NA       NA
## Temp1              0.04375         NA      NA       NA
## MixTime1           0.01875         NA      NA       NA
## BatchWt1          -0.01875         NA      NA       NA
## delay1             0.23500         NA      NA       NA
## Ports1:Temp1       0.00750         NA      NA       NA
## Ports1:MixTime1    0.04750         NA      NA       NA
## Ports1:BatchWt1    0.01500         NA      NA       NA
## Ports1:delay1      0.07625         NA      NA       NA
## Temp1:MixTime1    -0.03375         NA      NA       NA
## Temp1:BatchWt1     0.08125         NA      NA       NA
## Temp1:delay1       0.20250         NA      NA       NA
## MixTime1:BatchWt1  0.03625         NA      NA       NA
## MixTime1:delay1   -0.06750         NA      NA       NA
## BatchWt1:delay1    0.15750         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

Formula yang digunakan pada fungsi di atas, formula=y(.)^2, menghasilkan saturated model pada pengaruh utama dan interaksi 2-faktor. Oleh karena tidak ada ulangan pada percobaan tersebut, normal probability plot dari koefisien regresi dapat digunakan sebagai alat bantu untuk menilai pengaruh yang signifikan. Perhatikan pula bahwa setiap pengaruh pada model ini berbaur dengan satu interaksi lain, sehingga hanya 15 pengaruh saja yang dapat diduga.

soupc<-FrF2(16,5,generators="ABCD",randomize=FALSE)
soupc<-add.response(soupc, y)
modc<-lm(y~(.)^2, data=soupc)
library(daewr)
## Warning: package 'daewr' was built under R version 4.0.3
## 
## Attaching package: 'daewr'
## The following object is masked _by_ '.GlobalEnv':
## 
##     soup
LGB(coef(modc)[-1], rpt = FALSE)

Plot tersebut menunjukkan bahwa pengaruh utama E (Delay Time), BE (interaksi antara Temperature dan Delay Time), dan DE (interaksi antara Batch Weight dan Delay Time) terlihat signifikan. Plot interaksi antar faktor juga dapat dihasilkan dengan menggunakan fungsi IAPlot().

IAPlot(soupc)

interaction.plot(soupc$E, soupc$B, soupc$y, cex=0.5)

interaction.plot(soupc$E, soupc$D, soupc$y, cex=0.5)

Dari kedua output di atas terlihat bahwa cenderung terdapat interaksi antara Delay Time dengan Temperature, dan dengan Batch Weight.

Uji Anova

fit<-aov(y ~ A+B+C+D+E, data=soupc)
summary(fit)
##             Df Sum Sq Mean Sq F value Pr(>F)  
## A            1 0.0841  0.0841   0.599  0.457  
## B            1 0.0306  0.0306   0.218  0.651  
## C            1 0.0056  0.0056   0.040  0.845  
## D            1 0.0056  0.0056   0.040  0.845  
## E            1 0.8836  0.8836   6.292  0.031 *
## Residuals   10 1.4044  0.1404                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
dat<-soupc
dat$B<-as.factor(dat$B); dat$D<-as.factor(dat$D);dat$E<-as.factor(dat$E)
fit1<-aov(y~B+D+E+B:E+D:E, data=dat)
summary(fit1)
##             Df Sum Sq Mean Sq F value  Pr(>F)   
## B            1 0.0306  0.0306   0.694 0.42418   
## D            1 0.0056  0.0056   0.128 0.72844   
## E            1 0.8836  0.8836  20.031 0.00119 **
## B:E          1 0.6561  0.6561  14.873 0.00318 **
## D:E          1 0.3969  0.3969   8.997 0.01335 * 
## Residuals   10 0.4411  0.0441                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
library("car")
## Warning: package 'car' was built under R version 4.0.4
## Loading required package: carData
qqPlot(residuals(fit1))

## [1]  7 13
library(nortest)
## Warning: package 'nortest' was built under R version 4.0.3
ad.test(residuals(fit1))
## 
##  Anderson-Darling normality test
## 
## data:  residuals(fit1)
## A = 0.34888, p-value = 0.4294
plot(fit1$fitted.values, fit1$residuals)

Fraksi 1/4 atau lebih pada Rancangan 2^k

Misalkan ingin disusun rancangan faktorial 2^6βˆ’3 yaitu dengan generatornya adalah AB, AC dan BC

library(FrF2)
frac <- FrF2( 8, 6, generators = c("AB", "AC", "BC"))
frac
##    A  B  C  D  E  F
## 1  1  1 -1  1 -1 -1
## 2 -1  1 -1 -1  1 -1
## 3  1 -1  1 -1  1 -1
## 4  1 -1 -1 -1 -1  1
## 5 -1 -1  1  1 -1 -1
## 6  1  1  1  1  1  1
## 7 -1 -1 -1  1  1  1
## 8 -1  1  1 -1 -1  1
## class=design, type= FrF2.generators
library(FrF2)
y <- runif(8, 0, 1)
aliases( lm( y~ (.)^4, data = frac))
##                                                   
##  A = B:D = C:E = B:E:F = C:D:F = A:B:C:F = A:D:E:F
##  B = C:F = A:E:F = C:D:E = A:B:C:E = B:D:E:F = A:D
##  C = B:F = A:D:F = B:D:E = A:B:C:D = C:D:E:F = A:E
##  D = E:F = A:C:F = B:C:E = A:C:D:E = B:C:D:F = A:B
##  E = D:F = A:B:F = B:C:D = A:B:D:E = B:C:E:F = A:C
##  F = B:C = D:E = A:B:E = A:C:D = A:B:D:F = A:C:E:F
##  A:F = B:E = C:D = A:B:C = A:D:E = B:D:F = C:E:F

Output di atas dapat membantu kita menentukan definining relation berikut: 𝐼 = 𝐴𝐡𝐷 = 𝐴𝐢𝐸 = 𝐴𝐡𝐸𝐹 = 𝐴𝐢𝐷𝐹 = 𝐡𝐢𝐹 = 𝐷𝐸𝐹.