Andiilla_Persiapan_Rancob_UAS
Andi Illa Erviani Nensi
2025-06-10
#Fraksional Faktorial 2k-p ## Menyusun Rancangan Fraksional Faktorial 2^k-p Cara paling mudah untuk membuat rancangan faktorial 2^k−1 di R adalah menggunakan fungsi FrF2. Pada contoh berikut ini disusun rancangan 2^5−1 dengan generator 𝐸= ABCD
library(FrF2) ## 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':
##
## lengthsdesign<-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.generatorsFungsi alias dapat digunakan untuk menginformasikan struktur alias pada rancangan yang telah disusun sebelumnya. Namun fungsi ini memerlukan vektor peubah respon,𝑦, yang dapat pula 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:COutput 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 Percobaan dilakukan di suatu pabrik yang membuat produk sirkuit terintegrasi. Tujuan percobaan: mengetahui pengaruh 5 buah faktor (A, B, C, D, E) dalam meningkatkan hasil produksi
library(FrF2)
Data <- FrF2(16, 5, generators = "ABCD",
factor.names =list(A=c(-1,1), B=c(-1,1), C=c(-1,1), D=c(-1,1),E=c(-1,1)),
randomize = FALSE)
Data$y <- c(8, 9, 34, 52, 16, 22, 45, 60, 6, 10, 30, 50, 15, 21, 44, 63)
Data## A B C D E y
## 1 -1 -1 -1 -1 1 8
## 2 1 -1 -1 -1 -1 9
## 3 -1 1 -1 -1 -1 34
## 4 1 1 -1 -1 1 52
## 5 -1 -1 1 -1 -1 16
## 6 1 -1 1 -1 1 22
## 7 -1 1 1 -1 1 45
## 8 1 1 1 -1 -1 60
## 9 -1 -1 -1 1 -1 6
## 10 1 -1 -1 1 1 10
## 11 -1 1 -1 1 1 30
## 12 1 1 -1 1 -1 50
## 13 -1 -1 1 1 1 15
## 14 1 -1 1 1 -1 21
## 15 -1 1 1 1 -1 44
## 16 1 1 1 1 1 63
## class=design, type= FrF2.generatorsaliases(lm( y~ (.)^4, data = Data)) ##
## 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:Clibrary(DoE.base)
model <- lm(y ~ (.)^2, data = Data)
effect <- 2 * coef(model)[-1]
result <- data.frame(
Faktor = names(effect),
Effect = as.vector(effect)
)
print(result)## Faktor Effect
## 1 A1 11.125
## 2 B1 33.875
## 3 C1 10.875
## 4 D1 -0.875
## 5 E1 0.625
## 6 A1:B1 6.875
## 7 A1:C1 0.375
## 8 A1:D1 1.125
## 9 A1:E1 1.125
## 10 B1:C1 0.625
## 11 B1:D1 -0.125
## 12 B1:E1 -0.125
## 13 C1:D1 0.875
## 14 C1:E1 0.375
## 15 D1:E1 -1.375library(DoE.base)
mod1 <- lm( y ~ (.)^2, data = Data)
summary(mod1)##
## Call:
## lm.default(formula = y ~ (.)^2, data = Data)
##
## Residuals:
## ALL 16 residuals are 0: no residual degrees of freedom!
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 30.3125 NaN NaN NaN
## A1 5.5625 NaN NaN NaN
## B1 16.9375 NaN NaN NaN
## C1 5.4375 NaN NaN NaN
## D1 -0.4375 NaN NaN NaN
## E1 0.3125 NaN NaN NaN
## A1:B1 3.4375 NaN NaN NaN
## A1:C1 0.1875 NaN NaN NaN
## A1:D1 0.5625 NaN NaN NaN
## A1:E1 0.5625 NaN NaN NaN
## B1:C1 0.3125 NaN NaN NaN
## B1:D1 -0.0625 NaN NaN NaN
## B1:E1 -0.0625 NaN NaN NaN
## C1:D1 0.4375 NaN NaN NaN
## C1:E1 0.1875 NaN NaN NaN
## D1:E1 -0.6875 NaN NaN NaN
##
## 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: NAFormula 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.
library(daewr) LGB(coef(mod1)[-1], rpt = TRUE)
## Effect Report
##
## Label Half Effect Sig(.05)
## A1 5.5625 yes
## B1 16.9375 yes
## C1 5.4375 yes
## D1 -0.4375 no
## E1 0.3125 no
## A1:B1 3.4375 yes
## A1:C1 0.1875 no
## A1:D1 0.5625 no
## A1:E1 0.5625 no
## B1:C1 0.3125 no
## B1:D1 -0.0625 no
## B1:E1 -0.0625 no
## C1:D1 0.4375 no
## C1:E1 0.1875 no
## D1:E1 -0.6875 no
##
## Lawson, Grimshaw & Burt Rn Statistic = 6.098508
## 95th percentile of Rn = 1.201Anova
model <- aov(y~A+B+C+A:B, data = Data)
summary(model)## Df Sum Sq Mean Sq F value Pr(>F)
## A 1 495 495 193.19 2.53e-08 ***
## B 1 4590 4590 1791.24 1.56e-13 ***
## C 1 473 473 184.61 3.21e-08 ***
## A:B 1 189 189 73.78 3.30e-06 ***
## Residuals 11 28 3
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Faktorial 3k
Faktorial design 3^2
library(readxl)
Data <- read_excel("D:/Andi Illa Erviani Nensi/data_rancangan_percobaan.xlsx",sheet="data_faktorial3p2")
Data## # A tibble: 18 × 3
## Angle Speed Y
## <dbl> <dbl> <dbl>
## 1 0 0 -2
## 2 0 0 -1
## 3 1 0 0
## 4 1 0 2
## 5 2 0 -1
## 6 2 0 0
## 7 0 1 -3
## 8 0 1 0
## 9 1 1 1
## 10 1 1 3
## 11 2 1 5
## 12 2 1 6
## 13 0 2 2
## 14 0 2 3
## 15 1 2 4
## 16 1 2 6
## 17 2 2 0
## 18 2 2 -1Data$Angle = as.factor(Data$Angle)
Data$Speed = as.factor(Data$Speed)
str(Data)## tibble [18 × 3] (S3: tbl_df/tbl/data.frame)
## $ Angle: Factor w/ 3 levels "0","1","2": 1 1 2 2 3 3 1 1 2 2 ...
## $ Speed: Factor w/ 3 levels "0","1","2": 1 1 1 1 1 1 2 2 2 2 ...
## $ Y : num [1:18] -2 -1 0 2 -1 0 -3 0 1 3 ...anova = aov(Y~Angle*Speed, data=Data)
summary(anova)## Df Sum Sq Mean Sq F value Pr(>F)
## Angle 2 24.33 12.167 8.423 0.00868 **
## Speed 2 25.33 12.667 8.769 0.00770 **
## Angle:Speed 4 61.33 15.333 10.615 0.00184 **
## Residuals 9 13.00 1.444
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Menggunakan Kontras
# Fungsi kontras linier dan kuadratik
contrast_L <- function(x) {
ifelse(x == 0, -1, ifelse(x == 1, 0, 1))
}
contrast_Q <- function(x) {
ifelse(x == 1, -2, 1)
}
# kolom kontras
Data$A_L <- contrast_L(Data$Angle)
Data$A_Q <- contrast_Q(Data$Angle)
Data$B_L <- contrast_L(Data$Speed)
Data$B_Q <- contrast_Q(Data$Speed)
Data## # A tibble: 18 × 7
## Angle Speed Y A_L A_Q B_L B_Q
## <fct> <fct> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0 0 -2 -1 1 -1 1
## 2 0 0 -1 -1 1 -1 1
## 3 1 0 0 0 -2 -1 1
## 4 1 0 2 0 -2 -1 1
## 5 2 0 -1 1 1 -1 1
## 6 2 0 0 1 1 -1 1
## 7 0 1 -3 -1 1 0 -2
## 8 0 1 0 -1 1 0 -2
## 9 1 1 1 0 -2 0 -2
## 10 1 1 3 0 -2 0 -2
## 11 2 1 5 1 1 0 -2
## 12 2 1 6 1 1 0 -2
## 13 0 2 2 -1 1 1 1
## 14 0 2 3 -1 1 1 1
## 15 1 2 4 0 -2 1 1
## 16 1 2 6 0 -2 1 1
## 17 2 2 0 1 1 1 1
## 18 2 2 -1 1 1 1 1model <- lm(Y ~ A_L + A_Q + B_L + B_Q + A_L*B_L + A_L*B_Q + A_Q*B_L + A_Q*B_Q, data = Data)
anova(model)## Analysis of Variance Table
##
## Response: Y
## Df Sum Sq Mean Sq F value Pr(>F)
## A_L 1 8.333 8.333 5.7692 0.0397723 *
## A_Q 1 16.000 16.000 11.0769 0.0088243 **
## B_L 1 21.333 21.333 14.7692 0.0039479 **
## B_Q 1 4.000 4.000 2.7692 0.1304507
## A_L:B_L 1 8.000 8.000 5.5385 0.0430650 *
## A_L:B_Q 1 42.667 42.667 29.5385 0.0004137 ***
## A_Q:B_L 1 2.667 2.667 1.8462 0.2073056
## A_Q:B_Q 1 8.000 8.000 5.5385 0.0430650 *
## Residuals 9 13.000 1.444
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1OPTIMAL DESIGN
library(AlgDesign)
library(ggplot2)
library(hrbrthemes)
library(mixexp)## Loading required package: latticeOptimal Design untuk kasus response surface
dat<-gen.factorial(levels=3,nVars=3,varNames=c("A","B","C"))Kriteria D-Optimal mengguakan algoritma fedorov
set.seed(123)
desD<-optFederov(~.,dat,nTrials=14,eval=TRUE,crit="D")Rancangan D-Optimal:
desD$design## A B C
## 1 -1 -1 -1
## 3 1 -1 -1
## 6 1 0 -1
## 7 -1 1 -1
## 8 0 1 -1
## 9 1 1 -1
## 10 -1 -1 0
## 12 1 -1 0
## 19 -1 -1 1
## 21 1 -1 1
## 24 1 0 1
## 25 -1 1 1
## 26 0 1 1
## 27 1 1 1cat("Nilai Determinan:", desD$D, "\n")## Nilai Determinan: 0.8854694plot rancangan D-Optimal untuk faktor A dan faktor B
datadesignD<-desD$design
ggplot(datadesignD, aes(x=A, y=B)) +
geom_point()+
theme_ipsum()## Warning in grid.Call(C_stringMetric, as.graphicsAnnot(x$label)): font family
## not found in Windows font database
## Warning in grid.Call(C_stringMetric, as.graphicsAnnot(x$label)): font family
## not found in Windows font database## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
plot rancangan D-Optimal untuk faktor A dan faktor C
datadesignD<-desD$design
ggplot(datadesignD, aes(x=A, y=C)) +
geom_point()+
theme_ipsum()## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
### Kriteria I-Optimal mengguakan algoritma fedorov
set.seed(123)
desI<-optFederov(~.,dat,nTrials=14,eval=TRUE,crit="I")Rancangan D-Optimal:
desI$design## A B C
## 1 -1 -1 -1
## 2 0 -1 -1
## 3 1 -1 -1
## 7 -1 1 -1
## 8 0 1 -1
## 9 1 1 -1
## 16 -1 1 0
## 18 1 1 0
## 19 -1 -1 1
## 21 1 -1 1
## 22 -1 0 1
## 24 1 0 1
## 25 -1 1 1
## 27 1 1 1cat("Average Variances:", desI$I, "\n")## Average Variances: 3.376694plot rancangan I-Optimal untuk faktor A dan faktor B
datadesignI<-desI$design
ggplot(datadesignI, aes(x=A, y=B)) +
geom_point()+
theme_ipsum()## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
plot rancangan I-Optimal untuk faktor A dan faktor B
ggplot(datadesignI, aes(x=A, y=C)) +
geom_point()+
theme_ipsum()## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Optimal Design untuk kasus mixture design
dat<-gen.mixture(4,3)Kriteria D-Optimal mengguakan algoritma fedorov
set.seed(123)
desM_D<-optFederov(~X1+X2+X3-1,dat,6,nullify=TRUE,crit = "D")Rancangan D-Optimal:
desM_D$design## X1 X2 X3
## 1 1.0000000 0.0000000 0.0000000
## 3 0.3333333 0.6666667 0.0000000
## 4 0.0000000 1.0000000 0.0000000
## 5 0.6666667 0.0000000 0.3333333
## 9 0.0000000 0.3333333 0.6666667
## 10 0.0000000 0.0000000 1.0000000cat("Nilai Determinan:", desM_D$D, "\n")## Nilai Determinan: 0.2543809plot rancangan D-Optimal untuk mixture design
DesignPoints(desM_D$design)
Kriteria I-Optimal mengguakan algoritma fedorov
set.seed(123)
desM_I<-optFederov(~X1+X2+X3-1,dat,7,nullify=TRUE,crit = "I")Rancangan I-Optimal:
desM_I$design## X1 X2 X3
## 1 1.0000000 0.0000000 0.0000000
## 2 0.6666667 0.3333333 0.0000000
## 3 0.3333333 0.6666667 0.0000000
## 4 0.0000000 1.0000000 0.0000000
## 8 0.3333333 0.0000000 0.6666667
## 9 0.0000000 0.3333333 0.6666667
## 10 0.0000000 0.0000000 1.0000000cat("Average Variances:", desM_I$I, "\n")## Average Variances: 2.654545plot rancangan I-Optimal untuk mixture design
DesignPoints(desM_I$design)
RESPON SURFACE METHODS
Central Composite Design (CCD)
library(readxl)
code <- read_xlsx("D:/Andi Illa Erviani Nensi/data_rancangan_percobaan.xlsx",sheet="CCD")
data_nat <- read_xlsx("D:/Andi Illa Erviani Nensi/data_rancangan_percobaan.xlsx",sheet="respon_surface")
code## # A tibble: 13 × 3
## x1 x2 y
## <dbl> <dbl> <dbl>
## 1 -1 -1 0.0141
## 2 1 -1 0.0149
## 3 -1 1 0.0137
## 4 1 1 0.0132
## 5 -1.41 0 0.0117
## 6 1.41 0 0.0108
## 7 0 -1.41 0.014
## 8 0 1.41 0.0132
## 9 0 0 0.026
## 10 0 0 0.0266
## 11 0 0 0.0274
## 12 0 0 0.0268
## 13 0 0 0.0276data_nat## # A tibble: 13 × 3
## x1 x2 y
## <dbl> <dbl> <dbl>
## 1 30 6 0.0141
## 2 44 6 0.0149
## 3 30 8 0.0137
## 4 44 8 0.0132
## 5 27.1 7 0.0117
## 6 46.9 7 0.0108
## 7 37 5.59 0.014
## 8 37 8.41 0.0132
## 9 37 7 0.026
## 10 37 7 0.0266
## 11 37 7 0.0274
## 12 37 7 0.0268
## 13 37 7 0.0276library(rsm)ccd <- rsm(y ~ SO(x1,x2),data=code)
summary(ccd)##
## Call:
## rsm(formula = y ~ SO(x1, x2), data = code)
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.02687611 0.00045976 58.4573 1.125e-10 ***
## x1 -0.00012149 0.00036401 -0.3337 0.7483
## x2 -0.00040470 0.00036401 -1.1118 0.3030
## x1:x2 -0.00032500 0.00051403 -0.6323 0.5473
## x1^2 -0.00744825 0.00039145 -19.0273 2.756e-07 ***
## x2^2 -0.00626622 0.00039145 -16.0077 9.017e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Multiple R-squared: 0.9875, Adjusted R-squared: 0.9785
## F-statistic: 110.3 on 5 and 7 DF, p-value: 1.688e-06
##
## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## FO(x1, x2) 2 0.00000142 7.1200e-07 0.6737 0.54002
## TWI(x1, x2) 1 0.00000042 4.2200e-07 0.3998 0.54730
## PQ(x1, x2) 2 0.00058083 2.9041e-04 274.7795 2.231e-07
## Residuals 7 0.00000740 1.0570e-06
## Lack of fit 3 0.00000575 1.9170e-06 4.6523 0.08581
## Pure error 4 0.00000165 4.1200e-07
##
## Stationary point of response surface:
## x1 x2
## -0.00745505 -0.03209919
##
## Eigenanalysis:
## eigen() decomposition
## $values
## [1] -0.006244286 -0.007470185
##
## $vectors
## [,1] [,2]
## x1 0.1337578 -0.9910141
## x2 -0.9910141 -0.1337578# stationary point
sp_ccd <- canonical(ccd)$xs
print(sp_ccd)## x1 x2
## -0.00745505 -0.03209919calculate_optim <- function(sp_ccd, data_code, var_names) {
optim_vals <- sapply(1:length(var_names), function(i) {
var <- var_names[i]
sp_ccd[i] * ((max(data_nat[[var]]) - min(data_nat[[var]]))/2) + median(data_nat[[var]])
})
names(optim_vals) <- var_names
return(round(optim_vals, 2))
}
var_names <- c("x1", "x2")
# Nilai optimal
optim_vals <- calculate_optim(sp_ccd, data_nat, var_names)
print(optim_vals)## x1 x2
## 36.93 6.95ccd1 <- rsm(y ~ SO(x1,x2),data=data_nat)
contour(ccd1, ~ x1 + x2,
image = TRUE, main="second-order model",
axes=T)
persp(ccd, ~ x1 + x2, col=terrain.colors(50), contours="color", zlab = "y", main="second-order model")
Box-Behnken Design (BBD)
library(readxl)
df <- read_xlsx("D:/Andi Illa Erviani Nensi/data_rancangan_percobaan.xlsx",sheet="rsm_bbd1")
df.ori <- read_xlsx("D:/Andi Illa Erviani Nensi/data_rancangan_percobaan.xlsx",sheet="rsm_bbd2")
df## # A tibble: 15 × 4
## x1 x2 x3 y
## <dbl> <dbl> <dbl> <dbl>
## 1 -1 0 1 19.4
## 2 -1 0 -1 19.5
## 3 0 0 0 27.8
## 4 1 0 1 19
## 5 0 0 0 27.1
## 6 0 -1 -1 23.2
## 7 0 0 0 27.3
## 8 1 1 0 18.2
## 9 0 1 1 15.2
## 10 0 -1 1 19.3
## 11 1 0 -1 22.3
## 12 1 -1 0 18.5
## 13 0 1 -1 20.6
## 14 -1 -1 0 21.4
## 15 -1 1 0 15.5df.ori## # A tibble: 15 × 4
## x1 x2 x3 y
## <dbl> <dbl> <dbl> <dbl>
## 1 12 170 15 19.4
## 2 12 170 10 19.5
## 3 15 170 0 27.8
## 4 18 170 1 19
## 5 15 170 0 27.1
## 6 15 150 -1 23.2
## 7 15 170 0 27.3
## 8 18 190 0 18.2
## 9 15 190 1 15.2
## 10 15 150 1 19.3
## 11 18 170 -1 22.3
## 12 18 150 0 18.5
## 13 15 190 -1 20.6
## 14 12 150 0 21.4
## 15 12 190 0 15.5bbd <- rsm(y ~ SO(x1,x2,x3),data=df)
summary(bbd)##
## Call:
## rsm(formula = y ~ SO(x1, x2, x3), data = df)
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 27.40000 0.60484 45.3011 9.897e-08 ***
## x1 0.27500 0.37039 0.7425 0.4911922
## x2 -1.61250 0.37039 -4.3535 0.0073347 **
## x3 -1.58750 0.37039 -4.2860 0.0078184 **
## x1:x2 1.40000 0.52381 2.6727 0.0442033 *
## x1:x3 -0.80000 0.52381 -1.5273 0.1872248
## x2:x3 -0.37500 0.52381 -0.7159 0.5060792
## x1^2 -4.26250 0.54520 -7.8183 0.0005489 ***
## x2^2 -4.73750 0.54520 -8.6895 0.0003339 ***
## x3^2 -3.08750 0.54520 -5.6631 0.0023872 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Multiple R-squared: 0.975, Adjusted R-squared: 0.93
## F-statistic: 21.68 on 9 and 5 DF, p-value: 0.001726
##
## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## FO(x1, x2, x3) 3 41.568 13.856 12.6249 0.0090579
## TWI(x1, x2, x3) 3 10.962 3.654 3.3295 0.1141539
## PQ(x1, x2, x3) 3 161.620 53.873 49.0873 0.0003934
## Residuals 5 5.487 1.097
## Lack of fit 3 5.227 1.742 13.4038 0.0702220
## Pure error 2 0.260 0.130
##
## Stationary point of response surface:
## x1 x2 x3
## 0.03028555 -0.15575398 -0.25154991
##
## Eigenanalysis:
## eigen() decomposition
## $values
## [1] -2.872400 -3.972702 -5.242398
##
## $vectors
## [,1] [,2] [,3]
## x1 -0.3746510 0.7139732 0.59150568
## x2 -0.2308845 0.5460306 -0.80532162
## x3 0.8979582 0.4382840 0.03972553# stationary point
sp <- canonical(bbd)$xs
print(sp)## x1 x2 x3
## 0.03028555 -0.15575398 -0.25154991calculate_optim <- function(sp, df, var_names) {
optim_vals <- sapply(1:length(var_names), function(i) {
var <- var_names[i]
sp[i] * ((max(df.ori[[var]]) - min(df.ori[[var]]))/2) + median(df.ori[[var]])
})
names(optim_vals) <- var_names
return(round(optim_vals, 2))
}
var_names <- c("x1", "x2", "x3")
# Nilai optimal
optim_vals <- calculate_optim(sp, df.ori, var_names)
print(optim_vals)## x1 x2 x3
## 15.09 166.88 -2.01contour(bbd, ~ x1 + x2 + x3,
image = TRUE, main="second-order model",
axes=T)
persp(bbd, ~ x1 + x2 + x3, col=terrain.colors(50), contours="color", zlab = "y", main="second-order model")
