Optimasi Week13
Theodora Putrina Gea
11/30/2020
5.1 Percobaan Campuran
Semua desain yang dibahas dalam babsebelumnya 4 termasuk dalam kelas desain ortogonal di mana semua faktor bervariasi secara independen. Namun, dalam banyak aplikasi dan terutama dalam ilmu kimia ada kasus dimana faktor tidak dapat dan tidak boleh divariasikan secara independen. Jika, katakanlah, ada tiga senyawa x1, x2, x3, semuanya dalam kisaran 0-1 dengan 1 yang menunjukkan komponen murni, desain dari pencampuran komponen ini harus mematuhi aturan campuran \(∑_ix_i=1\). Sebenarnya, semua masalah dalam kimia adalah masalah campuran, tetapi dalam banyak kasus, pelarut diperlakukan sebagai “pengisi” inert yang tidak berkontribusi pada batasan campuran dan dengan demikian memungkinkan faktor konsentrasi yang tersisa untuk bervariasi 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 tentangpaket mixexp, kumpulan fungsi untuk membuat dan menganalisis desain campuran di R.
5.1. Model Campuran.
Mencoba mengidentifikasi model OLS parametrik \(\hat y=a_0+a_1.x_1+a_2.x_2+a_3.x_3\) diberikan batasan campuran \(\sum_{i=1}^3x_{i}=1\) akan gagal karena batasan campuran membuat satu variabel menjadi mubazir. 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}^{'}\) dan kondisi tambahan \(x_3 = 1 - x_1 - x_2\). Ini disebut pendekatan variabel slack dengan variabel slack, di sini \(x_3\), biasanya dipilih sebagai komponen yang paling inert. Perhatikan bahwa variabel lain dapat dikecualikan tanpa memengaruhi statistik model, karena semua model variabel kendur 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 mengarah ke apa yang disebut model Scheffe, dinamai menurut ilmuwan yang pertama kali menurunkan 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 camupuran
\[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)
\[y=\sum_{i=1}^Ka'_1.x_1+ϵ\]
Dari persamaan (5.1) dan batasan campuran mengikuti identitas \(\hat {y}_{x_{i}=1} = a_{i}^{'}\), jadi koefisien regresi adalah nilai yang diharapkan \(\hat y_{i}\) dari komponen murni \(x_i\) dengan campuran antara yang menjadi penjumlahan 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 Scheffe kuadrat. 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, diperoleh model Scheffe kuadrat \(\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 umum model Scheffe kuadrat untuk komponen K yang diberikan oleh persamaan (5.2)
\(y=\sum_{i=1}^Ka'_i.x_i+\sum_{j>1}^Ka'_{ij}.x_i.x_j+ϵ\)
Untuk semua campuran biner xi= 0,5, xj= 0,5 mengikuti persamaan (5.2) persamaan
\(\hat y_{x_{i}=0.5;x_{j}=0.5} = \frac {a_{i} + a_{j}} {2} + \frac {a_{ij}} {4}\)
yang digambarkan pada gambar 5.1 secara grafis untuk dua campuran biner dengan \(a_{ij}> 0\) (sinergisme) dan \(a_{ij}<0\) (antagonisme). Garis lurus menunjukkan model Scheffe linier dengan \(a_{ij}= 0\).
Gambar 5.1: Plot respons model Scheffe kuadrat untuk campuran biner \(x_1\), \(x_2\) dengan efek sinergis dan antagonis berlabel merah dan biru, masing-masing. Menghitung efek orde tinggi seperti interaksi terner dapat dicapai dengan menambah model Scheffe kuadrat dengan lebih tinggi urutan interaksi istilah seperti dalam model Scheffe kubik penuh, lihat persamaan (5.3).
\[y=\sum_{i=1}^Ka_i.x_i+\sum_{j>1}^Ka_{ij}.x_i.x_j+\sum_{j>1}^Kδ_{ij}.x_i.x_j.(x_i-x_j)+\sum_{k>j>1}^Ka_{ijk}.x_i.x_j.x_k+ϵ\]
Gambar 5.2 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 berbentuk parametrik
\(\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 disetel ke \(a_1= 50, a_2= 100, a_{12}= 0\) dan \(δ_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, rumus (5.3), sangat fleksibel karena banyaknya parameter regresi, oleh karena itu itu cenderung overfit data dengan mudah. Seringkali model kubik tereduksi (atau khusus), persamaan (5.4), cukup kompleks untuk menggambarkan situasi eksperimental secara memadai.
\[y=\sum_{i=1}^Ka_i.x_i+\sum_{j>1}^Ka_{ij}.x_i.x_j+\sum_{k>j>1}^Ka_{ijk}.x_i.x_j.x_k+ϵ\]
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 x1, x2,… xK digabungkan dengan faktor proses P z1,z2,… zP.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, persamaan (5.5)
\[\begin{align}y = &\sum_{i=1}^{K}b_{i}x_{i}+\sum_{j>i}^{K}b_{ij}.x_{i}.x_{j}+\sum_{i=1}^K\sum_{l=1}^Pb_{il}^{'}x_{i}.z_{l}+\sum_{j>i}^K\sum_{l=1}^Pb_{ijl}^{'}.x_{i}.x_{j}.z_{l}\+\notag\\&\sum_{i}^K\sum_{m>l}^Pb_{ilm}^{'}tx_{i} z_{l}z_{m}+\sum_{j>i}^K\sum_{m>l}^Pb_{ijlm}^{'}x_{i} x_{j}z_{l} z_{m}+\epsilon\tag{5.5}\end{align}\]
Jumlah koefisien regresi dalam rumus (5.5) sebagai a fungsi K dan P tumbuh besar dengan cepat membuat desain bersilangan menjadi penggunaan terbatas untuk sebagian besar aplikasi. Model yang lebih pelit diperoleh dengan secara khusus menghubungkan model Scheffe kuadrat 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, 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 denganpaket R mixexp dari (RH Myers, DC Montgomery 2002). Tujuan dari percobaan ini adalah untuk menemukan perpaduan yang optimal dari komponen x1, x2, x3 memaksimalkan laju etsa wafer silikon (respon erate). Rancangan percobaan ini adalah Desain Simplex-Centroid yang direplikasi sebagian (dijelaskan kemudian) 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 menghilangkan istilah model non-signifikan secara individual karena hal ini akan melanggar batasan campuran dan membuat model menjadi tunggal. Biasanya cukup untuk memilih di antara model # 1- # 4 untuk mendeskripsikan data dengan tepat dalam domain model41.
Data contoh 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 etsa dengan ulangan yang ditandai dengan lingkaran Mengikuti prinsip kesederhanaan, data pertama-tama dianalisis dengan model Scheffe kuadrat dan selanjutnya dengan model kubik tereduksi, dan residu dibandingkan pada gambar 5.4 (perhatikan bahwa R2 dalam lm-output dari model tanpa intersep bias dan harus diabaikan. Nilai R2 yang dilaporkan oleh MixModel benar dan sebagai gantinya harus digunakan).
library(mixexp)## Loading required package: gdata## gdata: Unable to locate valid perl interpreter
## gdata:
## gdata: read.xls() will be unable to read Excel XLS and XLSX files
## gdata: unless the 'perl=' argument is used to specify the location of a
## gdata: valid perl intrpreter.
## gdata:
## gdata: (To avoid display of this message in the future, please ensure
## gdata: perl is installed and available on the executable search path.)## gdata: Unable to load perl libaries needed by read.xls()
## gdata: to support 'XLX' (Excel 97-2004) files.## ## gdata: Unable to load perl libaries needed by read.xls()
## gdata: to support 'XLSX' (Excel 2007+) files.## ## gdata: Run the function 'installXLSXsupport()'
## gdata: to automatically download and install the perl
## gdata: libaries needed to support Excel XLS and XLSX formats.##
## Attaching package: 'gdata'## The following object is masked from 'package:stats':
##
## nobs## The following object is masked from 'package:utils':
##
## object.size## The following object is masked from 'package:base':
##
## startsWith## Loading required package: lattice## Loading required package: grid## Loading required package: daewr## Registered S3 method overwritten by 'DoE.base':
## method from
## factorize.factor conf.designdata("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# note: R2 is wrong because when there is no intercept
# then SS.model is not adjusted by the mean
# Normally, the intercept a0 ensures mean adjustment
(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# quadratic Scheffe R2 correct
(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# reduced cubic
kappa(res1)## [1] 9.353397kappa(res2)## [1] 57.53584par(mfrow=c(1,2))
colcol <- rep("black")
colcol[duplicated(etch[,1:3])]<- "red" #color label replicates
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()
# check residuals of reduced cubic model
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: Residual pelajar versus prediksi model dari kuadratik (panel kiri) dan model kubik Scheffe tereduksi dari data etsa. 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).
\(R^2\) yang dikoreksi Nilai dari model Scheffe kuadrat, R2= 0,76 adalah kecil dibandingkan dengan model kubik, R2= 0,98, menunjukkan bahwa x1 istilah: x2: x3 memberikan kontribusi substansial untuk varians dari respon. Kesimpulan serupa dapat ditarik dari gambar 5.4 yang mengungkapkan model kuadrat Scheffe menjadi bias dengan, 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 suku kesalahan
Model kubik dapat divisualisasikan dengan plot kontur terner yang ditunjukkan pada gambar 5.5. Maksimum terletak dekat dengan sentroid (x1= 0,33, x2= 0,33, x3= 0,33)42.
library(mixexp)
data("etch")
invisible(capture.output(res <- MixModel(etch, "erate",
mixcomps = c("x1", "x2", "x3"), model = 4)) )
# plot cubic mixture model as ternary contour plot
ModelPlot(model = res,
dimensions = list(x1 = "x1", x2 = "x2", x3 = "x3"),
contour = TRUE, cuts = 6, fill = TRUE)
Gambar 5.5: Plot kontur model kubik tereduksi dari data etsa 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, dan oleh karena itu, cara lain yang populer dan kurang 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) ruang campuran. Ketika bergerak ke arah Cox, katakanlah, x1 rasio komponen campuran yang tersisa tetap konstan, di sini
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 etsa kubik tereduksi (panel kanan).
5.2 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 campuran yang secara optimal mendukung model yang dipilih. Model dan desain adalah dua sisi dari mata uang yang sama: Model parametrik menentukan desain dan desain menentukan dan membatasi model yang dapat diperkirakan dari data.
Desain campuran secara luas dapat dibagi menjadi desain biasa di mana desain tersebut mencakup seluruh ruang campuran \((0≤xi≤1)\) dan dan desain tidak beraturan (atau dibatasi) di mana hanya ada sub ruang, \(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 diperkirakan dengan (regular) Simplex Lattice Designs (SLD) yang ditunjukkan untuk K = 3 pada gambar 5.8 dan dibuat dengan mixexp R-function SLD (fac, lev). Agument fac menunjukkan jumlah faktor dan lev menentukan jumlah level dengan xi> 0.
Gambar 5.8: Desain Kisi Simplex untuk memperkirakan model linear (SLD (3,1), panel kiri) dan kuadrat Scheffe (SLD (3,2), panel kanan)
Komponen murni (titik puncak pada gambar 5.8) cukup untuk memperkirakan model Scheffe linier. Model Scheffe kuadrat didukung dengan menambah desain sebelumnya dengan semua campuran biner di tengah tiga sisi. Kedua desain jenuh43 dan oleh karena itu perlu proses tambahan untuk memperkirakan varian kesalahan dan untuk memeriksa keberadaan 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 Kisi Simplex 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 SCD (3) (panel kiri) untuk model kubik tereduksi dan Simplex Lattice Design 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. Desain Simplex-Lattice 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 berikut. 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), persamaan (5.5), hasil dari campuran persilangan dan desain proses. Mari fac menjadi desain faktorial penuh atau pecahan dengan N1 baris dan kolom P dan campuran desain campuran dimensi (N2)x (K), produk Cartesian dimensi (N1 * N2)x (P + K) mendukung desain proses campuran yang sepenuhnya bersilangan seperti yang ditunjukkan kode-R berikut untuk K = 3 dan P = 3, masing-masing
rm(list=ls())
set.seed(1234)
fac <- FrF2::FrF2(8,3,randomize=F) # 2^3 full factorial design; N=8## creating full factorial with 8 runs ...mix <- mixexp::SLD(3,3) # SLD design N=10
# 3vertex+3x2edges+centroid
fac <- data.frame(sapply(fac,function(x)
as.numeric(as.character(x))))
# convert factor to numerical
colnames(fac) <- paste("z",1:3,sep="")
x <- merge(mix,fac) # Cartesian product
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)) # Fully crossed mixture-process model #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.35284kappa(res)## [1] 15.35277Desain dan model yang sepenuhnya bersilangan sangat besar dan hampir tidak layak dalam pengaturan laboratorium biasa, dan desain yang sebagian bersilangan, model # 6 (persamaan (5.6))), mungkin merupakan alternatif yang lebih realistis untuk sebagian besar proyek. Desain proses campuran yang bersilangan sebagian dan pelit dapat dihasilkan dengan memilih titik secara optimal D dari kumpulan kandidat yang sepenuhnya bersilangan dengan AlgDesign :: optFederov () seperti yang ditunjukkan oleh kode berikutterbatas
rm(list=ls())
set.seed(1234)
rsm <- rsm::bbd(3, n0=1,randomize=F)[,-(1:2)] # Box-Behnken design;N=13
sld <- mixexp::SLD(3,3) # cubic Simplex-lattice design;N=10
colnames(rsm) <- paste("z",1:3,sep="")
candidate <- merge(sld,rsm) # Cartesian product 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.1727kappa(res)## [1] 24.00369Desain atau tidak beraturan dapat dibangun dengan algoritma simpul ekstrim yang diimplementasikan di fungsi Xvert () daripaket mixexp dengan menentukan dimensi dan wilayah terbatas \((LB_{i} \ge 0) \leq x_{i} \leq UB_{i} \le 1\) dengan \(LB_i\) dan \(UB_i\) menunjukkan lebih rendah dan batas atas \(x_i\). Konsistensi membutuhkan jumlah dari batas atas untuk 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: Desain terbatas \(x_{1} < 0.3; x_{2} < 0.3; \sum_{i} x_{i} = 1\) (panel kiri) dan desain yang sesuai ditambah dengan titik-titik interior (panel kanan);
Tabel 5.1 memberikan gambaran umum tentang model dan desain campuran yang telah dibahas sejauh ini dikombinasikan dengan petunjuk singkat tentang cara membuat desain ini di R.
Table 5.1:
Akhirnya dan serupa dengan desain ortogonal, desain campuran dapat diblokir jika ada parameter gangguan seperti batch bahan baku yang berbeda. , teknisi atau peralatan laboratorium. 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 ke 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 xij= xi-xj untuk menentukan suku kubik xi: xj: (xi-xj) dalam ekspresi rumus. Kedua, tentukan variabel baru sesuai dengan block2 = as.integer (block == “B2”) sebagai indikator blok unik. Ini diperlukan karena lm () secara default menjatuhkan level pertama dari sebuah 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 # first workaround
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
# that won't work as blockB1+blockB2=1
# => confounded with 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.009553kappa(res) # 1.90344e+17 singular design## [1] 1.90344e+17## but this trick will help
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.009553kappa(res) # 47.5957## [1] 47.5957
Gambar 5.12: Desain campuran kubik penuh dalam dua blok B1 (panel kiri) dan B2 (panel kanan).