Week 13 Optimization: Mixture Experiments
Imelda Sianturi
1 Mixture Experiments (Eksperimen Campuran)
Dalam banyak aplikasi, terutama dalam ilmu kimia, ada kasus dimana faktor tidak dapat dan tidak boleh di variasikan secara independen. Contoh ada tiga senyawa x1, x2, x3 dan semuanya dalam kisaran 0-1, dengan 1 menunjukkan komponen murni, maka desain dari pencampuran komponen ini harus mematuhi aturan campuran \(\sum_{i}x_{i}=1\).
Sebenarnya, semua masalah dalam kimia adalah masalah campuran, tetapi dalam banyak kasus, pelarut di perlakukan sebagai “pengisi” inert yang tidak berkontribusi pada batasan campuran, dan dengan demikian memungkinkan faktor konsentrasi yang tersisa untuk bervariansi secara independen.
Karena kendala campuran, masalah campuran memerlukan model parametrik khusus dan sesuai dengan model ini desain khusus. Model campuran dan desain campuran akan diperkenalkan pada bagian berikut dan penggunaannya di R akan ditunjukkan dengan beberapa contoh yang dipinjam dari literatur.
Sebuah tinjauan komprehensif tentang eksperimen campuran dan sumber yang baik untuk bacaan lebih lanjut adalah buku (John A. Cornell 1990). (John Lawson, Cameron Willden 2016) memberikan gambaran umum tentang packages mixexp, kumpulan fungsi untuk membuat dan menganalisis desain campuran di R.
1.1 Mixture Models (Model Campuran)
Mari kita coba mengidentifikasi model OLS parametrik, \(\hat y = a_{0} + a_{1}\cdot x_{1} + a_{2}\cdot x_{2} + a_{3}\cdot x_{3}\) dengan kendala campurannya adalah \(\sum_{i=1}^3x_{i}=1\) yang akan gagal karena batasan campuran membuat satu variabel menjadi berlebihan.
Salah satu cara untuk menangani singularitas adalah dengan menghilangkan satu variabel, katakanlah x3, dan lanjutkan dengan model tereduksi \(\hat y =a_{0}^{'} + a_{1}\cdot x_{1}^{'} + a_{2}\cdot x_{2}^{'}\) dengan kondisi tambahan \(x_3=1-x_1-x_2\). Ini disebut pendekatan variabel slack dengan variabel slack (dalam kasus adalah x3), biasanya dipilih sebagai komponen yang paling inert.
Perhatikan bahwa variabel lain dapat dikecualikan tanpa memengaruhi statistik model, karena semua model variabel slack secara statistik merosot dan akan memberikan prediksi yang sama. Cara yang lebih alami untuk menangani batasan linier adalah dengan mengintegrasikan batasan campuran \(\sum_{i}x_{i}=1\) ke dalam persamaan regresi, yang disebut model Scheffe (nama ilmuwan yang pertama kali memperoleh persamaan ini).
Dari model parametrik \[\hat y = a_{0} + a_{1}\cdot x_{1} + a_{2}\cdot x_{2} + a_{3}\cdot x_{3}\] dan batasan campurannya \[x_{1}+x_{2}+x_{3}=1\] diikuti dengan substitusi \[\hat y = a_{0}(x_{1}+x_{2}+x_{3}) + a_{1}\cdot x_{1} + a_{2}\cdot x_{2} + a_{3}\cdot x_{3} \Leftrightarrow\] \[\hat y = (a_{0} + a_{1})\cdot x_{1} + (a_{0} + a_{2})\cdot x_{2} + (a_{0} + a_{3})\cdot x_{3}\] model efek utama Scheffe \[\hat y = a_{1}^{'}\cdot x_{1} + a_{2}^{'}\cdot x_{2} + a_{3}^{'}\cdot x_{3}\] yang meluas untuk komponen campuran K ke model efek utama Scheffe, persamaan (5.1): \[\begin{equation} y = \sum_{i=1}^{K}a_{i}^{'}\cdot x_{i} + \epsilon \tag{5.1} \end{equation}\]
Dari persamaan 5.1 dan batasan campuran mengikuti identitas \(\hat {y}_{x_{i}=1} = a_{i}^{'}\), jadi koefisien regresi adalah nilai \(\hat y_{i}\) yang diharapkan dari komponen murni \(x_{i}\), dengan campuran menengah hanya menjadi jumlah tertimbang dari komponen murni. Dalam model Scheffe linier, komponen tidak berinteraksi dan ini mungkin merupakan asumsi yang terlalu ketat untuk banyak aplikasi sehingga memotivasi model kuadrat Scheffe.
Dimulai dengan model kuadrat \[\hat y= a_{0} + a_{1}\cdot x_{1} + a_{2}\cdot x_{2} + a_{3}\cdot x_{3} + a_{4}\cdot x_{1}^2 + a_{5}\cdot x_{2}^2 + a_{6}\cdot x_{3}^2 + a_{7} \cdot x_{1} \cdot x_{2} + a_{8} \cdot x_{1} \cdot x_{3}+ a_{9} \cdot x_{2} \cdot x_{3}\] dengan mengintegrasikan batasan campuran maka diperoleh model kuadrat Scheffe, \[\hat y= a_{1}^{'}\cdot x_{1} + a_{2}^{'}\cdot x_{2} + a_{3}^{'}\cdot x_{3} + a_{4}^{'} \cdot x_{1} \cdot x_{2} + a_{5}^{'} \cdot x_{1} \cdot x_{3} + a_{6}^{'} \cdot x_{2} \cdot x_{3}\] dengan model kuadrat Scheffe umum untuk komponen K yang diberikan oleh persamaan (5.2): \[\begin{equation} y = \sum_{i=1}^{K}a_{i}^{'}\cdot x_{i} + \sum_{j>i}^{K} a_{ij}^{'} \cdot x_{i} \cdot x_{j} + \epsilon \tag{5.2} \end{equation}\] Untuk campuran biner apa pun \(x_{i}=0.5\), \(x_{j}=0.5\) mengikuti dari persamaan (5.2), maka persamaannya adalah; \[\hat y_{x_{i}=0.5;x_{j}=0.5} = \frac {a_{i} + a_{j}} {2} + \frac {a_{ij}} {4}\] yang digambarkan dalam gambar 5.1 secara grafis untuk dua campuran biner dengan \(x_{ij}>0\) (sinergisme) dan \(x_{ij}<0\) (antagonisme). Garis lurus menunjukkan model Scheffe linier dengan \(x_{ij}=0\)
Gambar 5.1 diatas: Plot respons model kuadrat Scheffe untuk campuran biner x1, x2 dengan efek sinergis dan antagonis berlabel merah dan biru.
Perhitungan untuk efek orde tinggi seperti interaksi terner dapat dicapai dengan menambah model kuadrat Scheffe dengan istilah interaksi orde tinggi seperti dalam model kubik penuh Scheffe, lihat persamaan (5.3) dibawah ini; \[\begin{equation} y = \sum_{i=1}^{K}a_{i}\cdot x_{i} + \sum_{j>i}^{K} a_{ij} \cdot x_{i} \cdot x_{j} + \sum_{j>i}^{K} \delta_{ij} \cdot x_{i} \cdot x_{j} \cdot (x_{i} - x_{j}) + \sum_{k>j>i}^{K} a_{ijk} \cdot x_{i} \cdot x_{j} \cdot x_{k} + \epsilon \tag{5.3} \end{equation}\]
Gambar 5.2 dibawah ini, menunjukkan plot campuran biner dari model kubik penuh, yang sekarang dapat menjelaskan efek sinergis dan antagonistik dalam ruang campuran. Model kubik yang mendasari gambar 5.2 adalah bentuk parametrik dari, \[\hat y=a_{1} \cdot x_{1} + a_{2} \cdot x_{2} + a_{12} \cdot x_{1} \cdot x_{2} + \delta_{12} \cdot x_{1} \cdot x_{2} \cdot (x_{1}-x_{2})\] dengan parameter yang mengacu pada \(a_{1}=50\), \(a_{2}=100\), \(a_{12}=0\) dan \(\delta_{12}=-200\)
Gambar 5.2: Plot respon model kubik Scheffe untuk campuran biner x1, x2 dengan efek sinergis dan antagonis dalam ruang campuran.
Model kubik penuh, dalam persamaan (5.3), sangat fleksibel karena banyaknya parameter regresi, sehingga cenderung mudah menyesuaikan data dengan mudah. Seringkali model kubik tereduksi (atau khusus), dalam persamaan (5.4) dibawah ini, cukup kompleks untuk menggambarkan situasi eksperimental secara memadai. \[\begin{equation} y = \sum_{i=1}^{K}a_{i}\cdot x_{i} + \sum_{j>i}^{K} a_{ij} \cdot x_{i} \cdot x_{j} + \sum_{k>j>i}^{K} a_{ijk} \cdot x_{i} \cdot x_{j} \cdot x_{k} + \epsilon \tag{5.4} \end{equation}\]
Model campuran yang dijelaskan di atas - model kubik linier, kuadrat, penuh dan tereduksi (khusus) - cukup untuk memodelkan sebagian besar sistem campuran dengan tepat dan dapat digunakan sebagai blok bangunan untuk memodelkan desain proses campuran, di mana komponen campuran K \(x_{1}, x_{2} ,… x_{K}\) digabungkan dengan faktor proses P \(z_{1}, z_{2},… z_{P}\).
Model proses campuran pertama dihasilkan dari perpotongan penuh model Scheffe kuadrat dengan model proses faktorial (bilinear), yaitu; \[y =\Bigg ( \sum_{i=1}^{K}a_{i}\cdot x_{i} + \sum_{j>i}^{K} a_{ij} \cdot x_{i} \cdot x_{j} \Bigg ) \cdot \Bigg (\alpha_{0} + \sum_{l=1}^P \alpha_{l} \cdot z_{l} + \sum_{m>l}^R \alpha_{lm} \cdot z_{l}\cdot z_{m}\Bigg) + \epsilon\] dan dari ekspansi produk mengikuti langsung model parametrik, dalam persamaan (5.5) \[\begin{align} y = & \sum_{i=1}^{K}b_{i}\cdot x_{i} + \sum_{j>i}^{K} b_{ij} \cdot x_{i} \cdot x_{j} + \sum_{i=1}^K \sum_{l=1}^P b_{il}^{'} \cdot x_{i} \cdot z_{l} + \sum_{j>i}^K \sum_{l=1}^P b_{ijl}^{'} \cdot x_{i} \cdot x_{j} \cdot z_{l} \ + \notag \\ & \sum_{i}^K \sum_{m>l}^P b_{ilm}^{'} \cdot x_{i} \cdot z_{l} \cdot z_{m} + \sum_{j>i}^K \sum_{m>l}^P b_{ijlm}^{'} \cdot x_{i} \cdot x_{j} \cdot z_{l} \cdot z_{m} + \epsilon \tag{5.5} \end{align}\] Jumlah koefisien regresi dalam rumus (5.5) sebagai fungsi K dan P tumbuh besar dengan cepat membuat desain silang menjadi penggunaan terbatas untuk sebagian besar aplikasi. Model yang lebih pelit diperoleh dengan secara khusus menghubungkan model kuadrat Scheffe dengan model RSM yang mengarah ke persamaan (5.6) \[\begin{equation} y = \sum_{i=1}^{K}b_{i}\cdot x_{i} + \sum_{j>i}^{K} b_{ij} \cdot x_{i} \cdot x_{j} + \sum_{i=1}^K \sum_{l=1}^P b_{il}^{'} \cdot x_{i} \cdot z_{l} + \sum_{m>l}^P b_{lm}^{'} \cdot z_{l} \cdot z_{m}+ \sum_{l=1}^P b_{ll} \cdot z_{l}^2 + \epsilon \tag{5.6} \end{equation}\]
Empat model campuran, dengan persamaan: # 1: (5.1), # 2: (5.2), # 3: (5.3), # 4: (5.4) dan dua desain proses campuran, desain proses campuran bersilangan penuh dan sebagian, persamaan; # 5: (5.5) dan # 6: (5.6), dapat dengan mudah direferensikan oleh argumen model dalam R-function mixexp :: MixModel (…, model = #) dengan # yang menunjukkan nomor model dalam daftar di atas.
Saat menganalisis eksperimen campuran, semua prinsip validasi model OLS berlaku seperti pengujian statistik dan analisis residual seperti yang akan ditunjukkan dengan beberapa contoh data yang disediakan dengan packages mixexp di R dari (R.H. Myers, D.C. Montgomery 2002). Tujuan dari percobaan ini adalah untuk menemukan perpaduan yang optimal dari komponen x1, x2, x3 yang memaksimalkan laju etch wafer silikon (respon erate).
Rancangan percobaan ini adalah Rancangan Simplex-Centroid yang direplikasi sebagian (akan dijelaskan nanti) ditambah dengan beberapa titik interior. Desainnya secara maksimal mendukung model kubik yang dikurangi dan akan berfungsi di sini sebagai contoh untuk mendemonstrasikan beberapa fungsi dari mixexp. Berbeda dari model ortogonal, dalam model campuran biasanya tidak mungkin untuk menghilangkan satu per satu istilah model yang tidak signifikan karena hal ini akan melanggar batasan campuran dan membuat model menjadi tunggal. Biasanya cukup untuk memilih di antara model #1 sampai model #4 untuk mendeskripsikan data dengan tepat dalam domain model.
Data contoh dibawah ini (ditunjukkan pada gambar 5.3) sebagai plot terner dan mendukung model #1, #2 atau #4. Replikasi ditandai dengan lingkaran dan diseimbangkan dengan baik di atas ruang desain agar tidak mengganggu ortogonalitas desain.
Gambar 5.3: Plot terner dari data etch dengan ulangan yang ditandai dengan lingkaran.
Mengikuti prinsip kesederhanaan, data pertama-tama dianalisis dengan model kuadrat Scheffe dan selanjutnya dengan model kubik tereduksi, dan residu dibandingkan pada gambar 5.4 (perhatikan bahwa \(R^2\) dalam keluaran lm dari model tanpa intersep bias dan harus diabaikan. Nilai \(R^2\) yang dilaporkan oleh MixModel sudah benar dan sebagai gantinya harus digunakan).
library(mixexp)
data("etch")
summary(lm(erate ~ -1 + x1 + x2 + x3 +
x1:x2 + x1:x3 + x2:x3,
data = etch)) ##
## Call:
## lm(formula = erate ~ -1 + x1 + x2 + x3 + x1:x2 + x1:x3 + x2:x3,
## data = etch)
##
## Residuals:
## Min 1Q Median 3Q Max
## -157.68 -30.13 11.62 38.04 177.90
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## x1 534.64 83.35 6.414 0.000206 ***
## x2 329.16 83.35 3.949 0.004241 **
## x3 252.73 83.35 3.032 0.016254 *
## x1:x2 1343.11 469.58 2.860 0.021145 *
## x1:x3 644.53 469.58 1.373 0.207133
## x2:x3 711.68 469.58 1.516 0.168101
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 120.3 on 8 degrees of freedom
## Multiple R-squared: 0.9721, Adjusted R-squared: 0.9511
## F-statistic: 46.38 on 6 and 8 DF, p-value: 8.741e-06# Perhatikan: R2 salah karena itu tidak ada intersep
# maka SS.model tidak disesuaikan dengan rata-rata
# normalnya, intersep a0 memastikan penyesuaian rata-rata
(res1 <- MixModel(etch, "erate",
mixcomps = c("x1", "x2", "x3"), model = 2))##
## coefficients Std.err t.value Prob
## x1 534.6383 83.34866 6.414481 0.000205933
## x2 329.1622 83.34866 3.949220 0.004240658
## x3 252.7336 83.34866 3.032245 0.016253810
## x2:x1 1343.1061 469.58404 2.860204 0.021144784
## x3:x1 644.5346 469.58404 1.372565 0.207132664
## x2:x3 711.6775 469.58404 1.515549 0.168101200
##
## Residual standard error: 120.261 on 8 degrees of freedom
## Corrected Multiple R-squared: 0.7583617##
## Call:
## lm(formula = mixmodnI, data = frame)
##
## Coefficients:
## x1 x2 x3 x1:x2 x1:x3 x2:x3
## 534.6 329.2 252.7 1343.1 644.5 711.7# kuadrat Scheffe R2 benar
(res2 <- MixModel(etch, "erate",
mixcomps = c("x1", "x2", "x3"), model = 4)) ##
## coefficients Std.err t.value Prob
## x1 550.199515 23.22446 23.69051468 6.067419e-08
## x2 344.723325 23.22446 14.84311192 1.509476e-06
## x3 268.294753 23.22446 11.55224716 8.203511e-06
## x2:x1 689.537037 146.51489 4.70625916 2.192427e-03
## x3:x1 -9.034392 146.51489 -0.06166193 9.525557e-01
## x2:x3 58.108466 146.51489 0.39660451 7.034720e-01
## x2:x3:x1 9243.333333 940.85336 9.82441444 2.404146e-05
##
## Residual standard error: 33.43177 on 7 degrees of freedom
## Corrected Multiple R-squared: 0.9836603##
## Call:
## lm(formula = mixmodnI, data = frame)
##
## Coefficients:
## x1 x2 x3 x1:x2 x1:x3 x2:x3 x1:x2:x3
## 550.200 344.723 268.295 689.537 -9.034 58.108 9243.333## [1] 9.353397## [1] 57.53584par(mfrow=c(1,2))
colcol <- rep("black")
colcol[duplicated(etch[,1:3])]<- "red" # warna label direplikasi
par(mfrow=c(1,2))
plot(predict(res1), rstudent(res1),ylim=c(-4,3),
xlab=expression( paste("predicted response ", hat(y)) ),
ylab="studentized residuals",type="n",
main=expression(paste("(",frac(paste("y - ",
hat(y)),sigma), ") ~ ", hat(y) )) )
text(predict(res1), rstudent(res1),label=1:nrow(etch), col=colcol)
abline(h=c(-2,0,2),lty=c(2,1,2))
grid()
# periksa sisa model kubik tereduksi
plot(predict(res2), rstudent(res2),ylim=c(-4,3),
xlab=expression( paste("predicted response ", hat(y)) ),
ylab="studentized residuals",type="n",
main=expression(paste("(",frac(paste("y - ",
hat(y)),sigma), ") ~ ", hat(y) )) )
text(predict(res2), rstudent(res2),label=1:nrow(etch),col=colcol)
abline(h=c(-2,0,2),lty=c(2,1,2))
grid()
Gambar 5.4: Perbandingan Residual mahasiswa untuk prediksi model dari kuadratik (panel kiri) dengan model kubik Scheffe tereduksi dari data etch. Pasangan ulangan (3,10), (2,9), (8,1) dan (7,11) lebih baik dijelaskan dalam kesalahan relikasi oleh model kubik (panel kanan) dibandingkan dengan model kuadrat (panel kiri).
Nilai \(R^2\) yang dikoreksi dari model kuadrat Scheffe, \(R^2\) = 0,76 lebih kecil dibandingkan dengan model kubik, \(R^2\) = 0,98, menunjukkan bahwa istilah x1: x2: x3 berkontribusi secara substansial terhadap varians dari respon. Kesimpulan serupa dapat ditarik dari gambar 5.4 yang mengungkapkan model kuadrat Scheffe menjadi bias, misalnya, tidak tepat menggambarkan pasangan replikasi #7, #11. Dalam model kubik Scheffe suku terner x1: x2: x3 ditemukan signifikan dan residual dari model kubik sekarang lebih sesuai dengan asumsi normalitas dari kesalahan \(\epsilon\).
Model kubik dapat divisualisasikan dengan plot kontur terner yang ditunjukkan pada gambar 5.5. Maksimum terletak dekat dengan sentroid \((x_{1} = 0.33, x_{2} = 0.33, x_{3} = 0.33)\).
library(mixexp)
data("etch")
invisible(capture.output(res <- MixModel(etch, "erate",
mixcomps = c("x1", "x2", "x3"), model = 4)) )
# plot model kubik campuran sebagai plot kontur terner
ModelPlot(model = res,
dimensions = list(x1 = "x1", x2 = "x2", x3 = "x3"),
contour = TRUE, cuts = 6, fill = TRUE)
Gambar 5.5: Plot kontur dari model kubik tereduksi dari data etch.
Plot kontur juga dapat digunakan untuk dimensi yang lebih tinggi dari K = 3 dengan mengkondisikan (mengiris) plot pada variabel yang ditetapkan konstan. Namun, ini dapat menjadi membingungkan ketika ada banyak irisan untuk dibandingkan, oleh karena itu cara lain yang populer dan tidak membingungkan untuk memvisualisasikan model campuran berdimensi tinggi adalah plot efek Cox. Pertama, arah Cox didefinisikan sebagai garis yang menghubungkan titik centroid dengan simpul-simpul (lihat gambar 5.6) dari ruang campuran. Saat bergerak ke arah Cox, katakanlah \(x_{1}\) rasio komponen campuran yang tersisa tetap konstan, dengan \(\frac {x_{2}}{x_{3}}=1\), dan ini setara dengan efek “independen” \(x_{1}\) dalam ruang campuran. Efek Cox adalah jejak yang diperoleh dari mengiris permukaan respons di sepanjang arah Cox. Gambar 5.6 menunjukkan arah Cox dalam ruang campuran tiga dimensi bersama dengan plot efek Cox dari model kubik tereduksi dari data etch.
rm(list=ls())
library(mixexp)
library(Ternary)
data(etch)
invisible(capture.output(res <- MixModel(etch, "erate",
mixcomps = c("x1", "x2", "x3"), model = 4)) )
# plot cubic mixture model as ternary contour plot
par(mfrow=c(1,2))
TernaryPlot(atip="x1", btip="x2", ctip="x3", axis.cex=1.2,lab.cex = 1.3,
axis.labels=seq(0,1,0.1))
TernaryArrows(c(0,0.5,0.5),c(1,0,0),col="red",length=0.15, lwd=2)
TernaryArrows(c(0.5,0,0.5),c(0,1,0),col="red",length=0.15, lwd=2)
TernaryArrows(c(0.5,0.5,0),c(0,0,1),col="red",length=0.15, lwd=2)
AddToTernary(points, list(x1=rep(0.33,3),x2=rep(0.33,3),
x3=rep(0.33,3)),pch=16,
col="black",cex=1.5)
title(main="Cox direction",cex.main=1.5)
ModelEff(nfac = 3, mod = 4, dir = 2, ufunc = res,
dimensions = list("x1", "x2", "x3"))
Gambar 5.6: Arah Cox dalam ruang campuran 3D (panel kiri) dan plot efek sepanjang arah Cox dari model etch kubik tereduksi (panel kanan).
1.2 Mixture Design (Desain Campuran)
Bagian sebelumnya memperkenalkan enam model parametrik yang berbeda untuk pemodelan campuran dan sistem proses campuran. Terkait dengan masing-masing model adalah desain proses campuran dan desain campuran yang secara optimal mendukung model yang dipilih. Model dan desain adalah dua sisi dari coin yang sama: Model parametrik menentukan desain dan desain menentukan dan membatasi model yang dapat diperkirakan dari data. Desain campuran dapat secara luas dibagi menjadi desain biasa di mana desain tersebut mencakup seluruh ruang campuran \((0 \leq x_{i} \leq 1)\) dan desain tidak beraturan (atau dibatasi) di mana hanya subruang, \(LB_{i}(>0) \leq x_{i} \leq UB_{i}(<1)\) dari seluruh ruang campuran terisi. Gambar 5.7 menunjukkan contoh desain campuran beraturan dan tidak beraturan dalam tiga dimensi.
Gambar 5.7: Contoh desain biasa (panel kiri) dan tidak beraturan (dibatasi) (panel kanan)
Model Scheffe linier dan kuadrat (model #1, #2) dapat diestimasi dengan (regular) Simplex Lattice Designs (SLD) yang ditunjukkan untuk K = 3 pada gambar 5.8 dan dibuat dengan mixexp R-function SLD (fac, lev). Fac argumen menunjukkan jumlah faktor dan lev menentukan jumlah level dengan \(x_{i} > 0\).
Gambar 5.8: Desain Simplex Lattice untuk memperkirakan model linear (SLD (3,1), panel kiri) dan kuadratik Scheffe (SLD (3,2), panel kanan).
Komponen murni (titik puncak pada gambar 5.8) cukup untuk memperkirakan model linier Scheffe. Model kuadrat Scheffe didukung dengan menambah desain sebelumnya dengan semua campuran biner di tengah tiga sisi. Kedua desain jenuh dan oleh karena itu perlu proses tambahan untuk memperkirakan varian kesalahan dan untuk memeriksa adanya efek urutan yang lebih tinggi.
Model kubik tereduksi (model # 4) dapat diperkirakan dengan Simplex Centroid Designs (SCD). Sekali lagi, SCD sudah jenuh dan perlu proses tambahan untuk memperkirakan varian error dan Lack of Fit.
Akhirnya, Desain Simplex Lattice tiga tingkat dapat digunakan untuk memperkirakan model kubik penuh (model #3). Semua desain dapat dibuat dengan mudah dengan fungsi mixexp :: SLD () dan mixexp :: SCD () dan contoh untuk K = 3 dari yang pertama dan yang terakhir digambarkan pada gambar 5.8 dan 5.9.
Gambar 5.9: Simplex Centroid Design atau SCD (3) (panel kiri) untuk model kubik tereduksi dan Simplex Lattice Design atau SLD tiga tingkat (3,3) (panel kanan) untuk model kubik penuh Scheffe.
Model dan desain campuran klasik kaku dalam hal struktur model dan ukuran desain mirip dengan apa yang telah dikatakan tentang desain ortogonal klasik. Simplex-Lattice Design (SLD) orde kedua misalnya akan membuat semua interaksi dua arah dapat diperkirakan, terlepas dari apakah mereka menarik atau mungkin secara fisik.
Struktur model dan ukuran desain dapat ditentukan secara lebih fleksibel dengan desain optimal seperti yang ditunjukkan dengan R-code dibawah ini. Di sini tujuannya adalah untuk secara pelit mengidentifikasi model Scheffe non-standar, dalam notasi-R \(~ -1 + x1 + x2 + x3 + x1: x2 + x1: x2: x3\), dengan enam putaran dari satu set kandidat dari 60 kandidat campuran reguler.
Gambar 5.10 menggambarkan titik desain optimal D dalam ruang campuran yang terdiri dari komponen murni (mengacu pada efek utama x1, x2, x3), campuran biner x1, x2 (untuk istilah interaksi x1: x2) dan dua titik yang dekat dengan pusat massa untuk mengidentifikasi interaksi tiga arah x1: x2: x3.
rm(list=ls())
library(AlgDesign)
library(Ternary)
set.seed(1234)
xx <- expand.grid(x1=seq(0,1,0.1),
x2=seq(0,1,0.1),
x3=seq(0,1,0.1) ) # grid design
x <- xx[apply(xx,1,sum)==1,] # sum_x.i=1 only
mix <- optFederov(~ -1 + x1 + x2 + x3 +
x1:x2 + x1:x2:x3, data=x,nTrials=6,
nRepeats=1000)$design
TernaryPlot(atip="x1", btip="x2", ctip="x3",
axis.labels=seq(0,1,0.1), lab.cex=1)
AddToTernary(points, mix,pch=16,col="red",cex=1)
Gambar 5.10: Desain campuran D-optimal model ~ -1 + x1 + x2 + x3 + x1: x2 + x1: x2: x3 ditulis dalam notasi rumus-R
Desain untuk model proses campuran bersilangan penuh (#5), dalam persamaan (5.5), hasil dari campuran persilangan dan desain proses. Misalkan fac menjadi desain faktorial penuh atau pecahan dengan baris \(N_{1}\) dan kolom P dan campur desain campuran dimensi (\(N_{2}\)) x (K), hasil kali Kartesius dari dimensi (\(N_{1}\) x \(N_{2}\)) x (\(P + K\)) mendukung persilangan penuh desain proses campuran seperti yang ditunjukkan R-code berikut untuk K = 3 dan P = 3, masing-masing.
rm(list=ls())
set.seed(1234)
fac <- FrF2::FrF2(8,3,randomize=F) # 2^3 penuh dengan desain faktorial; N=8
mix <- mixexp::SLD(3,3) # SLD desain N=10
# 3vertex+3x2edges+centroid
fac <- data.frame(sapply(fac,function(x)
as.numeric(as.character(x))))
# convert faktor ke numerikal
colnames(fac) <- paste("z",1:3,sep="")
x <- merge(mix,fac) # Cartesian produk
dim(x)## [1] 80 6x$y <- rnorm(nrow(x))
(res <- mixexp::MixModel(x, "y",
mixcomps = c("x1", "x2", "x3"),
procvars = c("z1", "z2", "z3"),
model = 5)) # model proses campuran (mixture) yang bersilangan sepenuhnya #5##
## coefficients Std.err t.value Prob
## x1 0.05990870 0.3082853 0.19432872 0.84695368
## x2 -0.57818736 0.3082853 -1.87549413 0.06842376
## x3 -0.20611270 0.3082853 -0.66857769 0.50780621
## x2:x1 0.05555502 1.3647251 0.04070785 0.96774194
## x3:x1 -1.00098144 1.3647251 -0.73346746 0.46777458
## x2:x3 0.85963536 1.3647251 0.62989636 0.53253190
## x1:z1 -0.17240678 0.3082853 -0.55924415 0.57927649
## x1:z2 0.48903461 0.3082853 1.58630507 0.12095767
## x1:z3 0.12188272 0.3082853 0.39535684 0.69479036
## x2:z1 0.22176294 0.3082853 0.71934310 0.47633042
## x2:z2 0.26971744 0.3082853 0.87489542 0.38712724
## x2:z3 -0.03665631 0.3082853 -0.11890384 0.90597790
## x3:z1 0.18778266 0.3082853 0.60911964 0.54606984
## x3:z2 0.04506181 0.3082853 0.14616916 0.88456052
## x3:z3 -0.12555989 0.3082853 -0.40728467 0.68608424
## x2:x1:z1 -1.82763803 1.3647251 -1.33919867 0.18846181
## x2:x1:z2 -0.75822632 1.3647251 -0.55558905 0.58174831
## x2:x1:z3 0.06996290 1.3647251 0.05126520 0.95938274
## x3:x1:z1 0.13423785 1.3647251 0.09836256 0.92216144
## x3:x1:z2 -0.90892644 1.3647251 -0.66601431 0.50942517
## x3:x1:z3 1.30367180 1.3647251 0.95526330 0.34548243
## x2:x3:z1 -0.72218119 1.3647251 -0.52917705 0.59976052
## x2:x3:z2 -1.33291497 1.3647251 -0.97669119 0.33489917
## x2:x3:z3 1.88583297 1.3647251 1.38184092 0.17509079
## x1:z1:z2 0.26702234 0.3082853 0.86615321 0.39184221
## x1:z1:z3 -0.58110294 0.3082853 -1.88495153 0.06710090
## x1:z2:z3 -0.14089678 0.3082853 -0.45703367 0.65024895
## x2:z1:z2 -0.13425104 0.3082853 -0.43547654 0.66568022
## x2:z1:z3 0.01472591 0.3082853 0.04776715 0.96215199
## x2:z2:z3 -0.17024018 0.3082853 -0.55221625 0.58403377
## x3:z1:z2 -0.30992865 0.3082853 -1.00533047 0.32109681
## x3:z1:z3 -0.69175568 0.3082853 -2.24388117 0.03074588
## x3:z2:z3 0.45376014 0.3082853 1.47188357 0.14928650
## x2:x1:z1:z2 -0.45968552 1.3647251 -0.33683378 0.73809548
## x2:x1:z1:z3 0.93480706 1.3647251 0.68497829 0.49751482
## x2:x1:z2:z3 2.83909884 1.3647251 2.08034487 0.04428803
## x3:x1:z1:z2 -0.64270562 1.3647251 -0.47094146 0.64037490
## x3:x1:z1:z3 3.11414546 1.3647251 2.28188481 0.02818482
## x3:x1:z2:z3 2.41847375 1.3647251 1.77213255 0.08439084
## x2:x3:z1:z2 -1.50904885 1.3647251 -1.10575299 0.27578515
## x2:x3:z1:z3 3.01676321 1.3647251 2.21052813 0.03316260
## x2:x3:z2:z3 -1.35284331 1.3647251 -0.99129365 0.32781272
##
## Residual standard error: 0.926512 on 38 degrees of freedom
## Corrected Multiple R-squared: 0.5948409##
## Call:
## lm(formula = mixmod, data = frame)
##
## Coefficients:
## x1 x2 x3 x2:x1 x3:x1 x2:x3
## 0.05991 -0.57819 -0.20611 0.05556 -1.00098 0.85964
## x1:z1 x1:z2 x1:z3 x2:z1 x2:z2 x2:z3
## -0.17241 0.48903 0.12188 0.22176 0.26972 -0.03666
## x3:z1 x3:z2 x3:z3 x2:x1:z1 x2:x1:z2 x2:x1:z3
## 0.18778 0.04506 -0.12556 -1.82764 -0.75823 0.06996
## x3:x1:z1 x3:x1:z2 x3:x1:z3 x2:x3:z1 x2:x3:z2 x2:x3:z3
## 0.13424 -0.90893 1.30367 -0.72218 -1.33291 1.88583
## x1:z1:z2 x1:z1:z3 x1:z2:z3 x2:z1:z2 x2:z1:z3 x2:z2:z3
## 0.26702 -0.58110 -0.14090 -0.13425 0.01473 -0.17024
## x3:z1:z2 x3:z1:z3 x3:z2:z3 x2:x1:z1:z2 x2:x1:z1:z3 x2:x1:z2:z3
## -0.30993 -0.69176 0.45376 -0.45969 0.93481 2.83910
## x3:x1:z1:z2 x3:x1:z1:z3 x3:x1:z2:z3 x2:x3:z1:z2 x2:x3:z1:z3 x2:x3:z2:z3
## -0.64271 3.11415 2.41847 -1.50905 3.01676 -1.35284## [1] 15.35277Desain dan model yang saling bersilangan sangat besar dan hampir tidak dapat dilakukan dalam pengaturan laboratorium biasa, dan desain yang sebagian disilang, model #6 (persamaan (5.6)), mungkin menjadi alternatif yang lebih realistis untuk sebagian besar proyek. Desain proses campuran yang disilangkan sebagian dan pelit, dapat dihasilkan dengan memilih titik secara D-optiomally dari kumpulan kandidat yang sepenuhnya bersilang dengan AlgDesign :: optFederov () seperti yang ditunjukkan oleh kode berikut,
rm(list=ls())
set.seed(1234)
rsm <- rsm::bbd(3, n0=1,randomize=F)[,-(1:2)] # Box-Behnken desain;N=13
sld <- mixexp::SLD(3,3) # kubik Simplex-lattice desain;N=10
colnames(rsm) <- paste("z",1:3,sep="")
candidate <- merge(sld,rsm) # Cartesian produk dim=(130 x 6)
mix <- AlgDesign::optFederov(~ -1 + (x1+x2+x3)^3 +
x1:z1 + x1:z2 + x1:z3 +
x2:z1 + x2:z2 + x2:z3 +
x3:z1 + x3:z2 + x3:z3 +
z1:z2 + z1:z3 + z2:z3 +
I(z1^2) + I(z2^2) + I(z3^2) ,
data=candidate,nTrials=30,
criterion="D", nRepeats=1000)$design
mix$y <- rnorm(nrow(mix))
(res <- mixexp::MixModel(mix, "y",
mixcomps = c("x1", "x2", "x3"),
procvars = c("z1", "z2", "z3"),
model = 6))## [1] "Warning, when using Model 6 the design in the process variables allow fitting the full quadratic model"
##
## coefficients Std.err t.value Prob
## x1 1.1568064 1.1661825 0.9919600 0.3471356
## x2 1.4896200 1.0769369 1.3832007 0.1999491
## x3 0.5302152 1.2014390 0.4413168 0.6693956
## I(z1^2) -0.4855809 0.5771805 -0.8412983 0.4219706
## I(z2^2) -0.2406845 0.5947270 -0.4046974 0.6951494
## I(z3^2) -0.2628178 0.6226630 -0.4220868 0.6828656
## x1:x2 -3.0211103 3.3953217 -0.8897862 0.3967412
## x1:x3 0.3875330 3.5732353 0.1084544 0.9160147
## x2:x3 -2.6004735 3.4484189 -0.7541060 0.4700549
## x1:z1 -0.5310433 0.5583275 -0.9511321 0.3663794
## x1:z2 -0.4529325 0.5666571 -0.7993062 0.4446950
## x1:z3 -0.2559378 0.4999943 -0.5118815 0.6210569
## x2:z1 -0.3321763 0.5764581 -0.5762367 0.5785786
## x2:z2 -0.3683906 0.5352710 -0.6882320 0.5086556
## x2:z3 -0.2032059 0.5336694 -0.3807712 0.7122050
## x3:z1 0.4478906 0.6118973 0.7319702 0.4828108
## x3:z2 0.4189375 0.6242892 0.6710632 0.5190298
## x3:z3 -0.8799695 0.5389427 -1.6327700 0.1369515
## z1:z2 0.1701386 0.4439679 0.3832228 0.7104494
## z1:z3 0.4408563 0.4035429 1.0924644 0.3030079
## z2:z3 0.1727396 0.4084885 0.4228750 0.6823111
##
## Residual standard error: 1.188226 on 9 degrees of freedom
## Corrected Multiple R-squared: 0.5440434##
## Call:
## lm(formula = mixmod, data = frame)
##
## Coefficients:
## x1 x2 x3 I(z1^2) I(z2^2) I(z3^2) x1:x2 x1:x3
## 1.1568 1.4896 0.5302 -0.4856 -0.2407 -0.2628 -3.0211 0.3875
## x2:x3 x1:z1 x1:z2 x1:z3 x2:z1 x2:z2 x2:z3 x3:z1
## -2.6005 -0.5310 -0.4529 -0.2559 -0.3322 -0.3684 -0.2032 0.4479
## x3:z2 x3:z3 z1:z2 z1:z3 z2:z3
## 0.4189 -0.8800 0.1701 0.4409 0.1727## [1] 24.00369Desain terbatas atau tidak beraturan dapat dibangun dengan algoritma simpul ekstrim yang diimplementasikan dalam fungsi Xvert () dari packages mixexp dengan menentukan dimensi dan wilayah yang dibatasi \((LB_{i} \ge 0) \leq x_{i} \leq UB_{i} \le 1\) dengan \(LB_{i}\) dan \(UB_{i}\) masing-masing menunjukkan batas bawah dan atas \(x_{i}\). Konsistensi membutuhkan jumlah dari batas atas melebihi 1, yaitu \(\sum_{i} UB_{i}> 1\).
Desain biasa dan tidak beraturan seringkali jarang dalam mendukung istilah tatanan yang lebih tinggi. Titik interior tambahan dapat dibuat dengan fungsi Fillv () dengan membuat rata-rata dan dengan demikian membuat titik tengah di antara semua kemungkinan pasangan.
Sebaliknya, jika desain yang dibatasi memiliki terlalu banyak proses yang diberikan tujuan proyek, hal itu dapat dikurangi dengan memilih subset titik desain secara optimal dengan R-function AlgDesign :: optFederov ().
Augmentasi desain terbatas dengan titik interior ditunjukkan dengan kode-R berikut, dan gambar 5.11 menunjukkan desain sebelum dan sesudah augmentasi.
rm(list=ls())
x1 <- mixexp::Xvert(3,uc=c(.3,.3,1),lc=c(0,0,0),ndm=1,plot=F)[,-4]
x2 <- mixexp::Fillv(3,x1)
par(mfrow=c(1,2))
TernaryPlot(atip="x1", btip="x2", ctip="x3",axis.cex=1.2,lab.cex = 1.5,
axis.labels=seq(0,1,0.1))
AddToTernary(points, x1,pch=16,col="red",cex=1.5)
TernaryPlot(atip="x1", btip="x2", ctip="x3",axis.cex=1.2,lab.cex = 1.5,
axis.labels=seq(0,1,0.1))
AddToTernary(points, x2,pch=16,col="red",cex=1.5)
Gambar 5.11: Batas desain \(x_{1} < 0.3; x_{2} < 0.3; \sum_{i} x_{i} = 1\) (panel kiri) dan desain yang sesuai ditambah dengan titik interior (panel kanan)
Tabel 5.1 dibawah ini memberikan gambaran umum tentang model dan desain campuran yang telah dibahas sejauh ini yang dikombinasikan dengan petunjuk singkat tentang cara membuat desain ini di R.
Tabel 5.1: Model dan desain campuran beserta petunjuk tentang cara membuat desain di R. Angka (1-6) dalam daftar secara langsung merujuk ke argumen model di R-function mixexp :: MixModel (…, model = ( 1-6))
Terakhir dan serupa dengan desain ortogonal, desain campuran dapat diblokir jika terdapat parameter gangguan seperti batch bahan baku, teknisi, atau peralatan laboratorium yang berbeda. Meskipun ada beberapa desain campuran yang diblokir secara optimal yang berasal dari prinsip pertama (untuk detailnya lihat (John A. Cornell 1990)), pemblokiran algoritmik dengan R-function AlgDesign :: optBlock () sudah cukup untuk sebagian besar tujuan praktis. Kode contoh berikut akan memblokir desain Simplex-lattic orde tiga dengan sepuluh jalan menjadi dua blok dan menggambarkan kedua blok pada gambar 5.12. Selain itu, kode contoh menunjukkan beberapa solusi untuk menghindari beberapa masalah yang tidak tercakup oleh fungsionalitas standar.
Pertama, tentukan variabel \(x_{ij} = x_{i}-x_{j}\) untuk menentukan suku kubik \(x_{i}: x_{j}: (x_{i}-x_{j})\) dalam ekspresi rumus. Kedua, tentukan variabel baru sesuai dengan block2 = as.integer (block == "B2") sebagai indikator blok unik. Hal ini diperlukan karena lm () secara default menjatuhkan level pertama dari suatu faktor (di sini B1 dalam blok = {B1, B2}) ketika model memiliki intersep. Namun, ketika tidak ada intersep, baik blockB1 dan blockB2 disimpan dalam matriks model sehingga menjadikan model singular karena blockB1 + blockB2 = 1 dan x1 + x2 + x3 = 1. Masalah ini dielakkan dengan variabel blok baru block2 = as.integer (blok == "B2").
Kemudian parameter lokasi model lm kubik penuh sesuai dengan model yang sebenarnya meskipun tidak ditemukan istilah yang signifikan. Namun, ini tidak mengejutkan karena hanya ada 1 derajat kebebasan yang membuat statistik pengujian sangat konservatif.
rm(list=ls())
set.seed(12345)
x <- mixexp::SLD(3,3)
x$x12 <- x$x1-x$x2 # solusi pertama
x$x13 <- x$x1-x$x3
x$x23 <- x$x2-x$x3
res <- AlgDesign::optBlock(~ -1 +(x1+x2+x3)^3 +
x1:x2:x12 + x1:x3:x13 + x2:x3:x23 , x,
c(6,6),criterion="Dpc", nRepeats=1000)
x <- res$design
x$block <- factor(c(rep("B1",6),rep("B2",6) ))
x$y <- 20*(x$block=="B2") + 3*x$x1 + 2*x$x2 +
5*x$x3 + 10*x$x1*x$x2*x$x12 + rnorm(nrow(x),0,0.1)
x$block2 <- as.integer(x$block=="B2") # second workaround
# itu tidak bekerja sebagai blockB1+blockB2=1
# => dibuat bingung karena x1+x2+x3=1
summary(res <- (lm(y ~ -1 + block +(x1+x2+x3)^3 +
x1:x2:x12 + x1:x3:x13 + x2:x3:x23 ,data=x)))##
## Call:
## lm(formula = y ~ -1 + block + (x1 + x2 + x3)^3 + x1:x2:x12 +
## x1:x3:x13 + x2:x3:x23, data = x)
##
## Residuals:
## 1 2 4 5 6 7 3
## 1.251e-17 -1.431e-17 2.339e-17 1.079e-01 -1.079e-01 -1.002e-17 -7.683e-18
## 51 61 8 9 10
## -1.079e-01 1.079e-01 -2.503e-17 6.674e-18 -1.067e-17
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## blockB1 5.0703 0.3052 16.613 0.0383 *
## blockB2 24.9924 0.2158 115.811 0.0055 **
## x1 -2.0946 0.3738 -5.604 0.1124
## x2 -2.9911 0.3738 -8.002 0.0791 .
## x3 NA NA NA NA
## x1:x2 -0.3166 1.0857 -0.292 0.8194
## x1:x3 0.3885 0.9403 0.413 0.7505
## x2:x3 -0.5591 0.9711 -0.576 0.6674
## x1:x2:x3 1.5672 5.7401 0.273 0.8303
## x1:x2:x12 10.8875 2.6148 4.164 0.1501
## x1:x3:x13 0.3372 1.9270 0.175 0.8897
## x2:x3:x23 -1.3537 2.3787 -0.569 0.6706
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2158 on 1 degrees of freedom
## Multiple R-squared: 1, Adjusted R-squared: 0.9998
## F-statistic: 6666 on 11 and 1 DF, p-value: 0.009553## [1] 1.90344e+17## trik ini akan membantu
x$block2 <- as.integer(x$block=="B2")
summary(res <- (lm(y ~ -1 + block2 +(x1+x2+x3)^3 +
x1:x2:x12 + x1:x3:x13 + x2:x3:x23 ,data=x)))##
## Call:
## lm(formula = y ~ -1 + block2 + (x1 + x2 + x3)^3 + x1:x2:x12 +
## x1:x3:x13 + x2:x3:x23, data = x)
##
## Residuals:
## 1 2 4 5 6 7 3
## -1.739e-17 2.315e-17 -1.655e-17 1.079e-01 -1.079e-01 9.307e-18 -7.520e-18
## 51 61 8 9 10
## -1.079e-01 1.079e-01 -1.215e-18 6.300e-18 -3.478e-17
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## block2 19.9221 0.2158 92.316 0.0069 **
## x1 2.9757 0.2158 13.789 0.0461 *
## x2 2.0792 0.2158 9.635 0.0658 .
## x3 5.0703 0.3052 16.613 0.0383 *
## x1:x2 -0.3166 1.0857 -0.292 0.8194
## x1:x3 0.3885 0.9403 0.413 0.7505
## x2:x3 -0.5591 0.9711 -0.576 0.6674
## x1:x2:x3 1.5672 5.7401 0.273 0.8303
## x1:x2:x12 10.8875 2.6148 4.164 0.1501
## x1:x3:x13 0.3372 1.9270 0.175 0.8897
## x2:x3:x23 -1.3537 2.3787 -0.569 0.6706
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2158 on 1 degrees of freedom
## Multiple R-squared: 1, Adjusted R-squared: 0.9998
## F-statistic: 6666 on 11 and 1 DF, p-value: 0.009553## [1] 47.5957
Gambar 5.12: Desain campuran kubik penuh dalam dua blok B1 (panel kiri) dan B2 (panel kanan)
Reference
John A. Cornell. 1990. Experiments with Mixtures. 2nd ed. Wiley, New York.
John Lawson, Cameron Willden. 2016. “Mixture Experiments in R Using Mixexp.” Journal of Statistical Software 72. https://www.jstatsoft.org/article/view/v072c02.
R.H. Myers, D.C. Montgomery. 2002. Response Surface Methodology: Process and Product Optimization Using Designed Experiments. 2nd ed. John Wiley & Sons.