Optimization-Lab13-Translation

1 Mixture Experiments

Semua desain yang sudah didiskusikan di bagian 4 sebelumnya, berhubungan dengan kelas dari desain ortogonal dimana semua faktornya bervariasi independen. Namun, dalam banyak aplikasi dan terutama di pengetahuan bahan kimia, terdapat kasus dimana faktor-faktornya tidak bisa dan harus bukan bervariasi independen. Jika dikatakan ada tiga komponen yaitu \(x1, x2, x3\), semua didalam jarak 0-1, dengan 1 merupakan komponen yang murni. Desain dari menggabungkan kompinen ini harus selalu mematuhi aturan campuran \(∑_i x_i = 1\).

Sesungguhnya, semua masalah dalam bahan 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 dengan desain yang 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 ada didalam buku (John A. Cornell 1990). (John Lawson, Cameron Willden 2016) memberikan gambaran umum tentang mixexp paket, kumpulan fungsi untuk membuat dan menganalisis desain campuran di R.

1.1 Mixture Models

Mencoba untuk mengidentifikasi parameter dari model OLS \(\hat y = a_0 + a_1 ⋅ x_1 + a_2 ⋅x_2 + a_3 ⋅ x_3\) diberikan dari kendala campuran \(∑_{i=1}^3 x_i =1\) akan gagal karena kendala campuran tersebut membuat satu variabel menjadi berlebihan. Satu jalan yang berurusan dengan singularitas merupakan dengan menjatuhkan salah satu variabel, katakanlah \(x_3\), dan dilanjutkan dengan pengurangan model \(\hat y =a_0^{'} + a_1 \cdot x_1^{'} + a_2 \cdot x_2^{'}\) dan dengan kondisi tambahan bahwa \(x_3 = 1 - x_1 - x_2\). Ini dipanggil dengan pendekatan variabel slack satu sama lain, disini \(x_3\), biasanya dipilih menjadi komponen yang paling lembam.

Ketahuilah bahwa variabel lain dapat dikeluarkan tanpa mempengaruhi statistik model, karena semua model variabel kendur secara statistik merosot dan akan memberikan prediksi yang sama. Cara yang lebih alami untuk menangani kendala linear adalah dengan mengintegrasikan pembatas campuran \(∑_i x_i =1\) ke dalam persamaan regresi, yang mengarah ke apa yang disebut model Scheffe, yang dinamai setelah ilmuwan 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 campuran kontrasnya \[x_1 + x_2 +x_3 =1\] diikuti dengan subtitusi \[\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 Schefe \[\hat y = a_1^{'} \cdot x_1 + a_2^{'} \cdot x_2 + a_3^{'} \cdot x_3\] dimana meluas untuk komponen campuran K ke model efek utama Scheffe, persamaannya \[y=∑_{i=1}^K a_i^{'} \cdot x_i + ϵ\]

Dari persamaan diatas dan batasan campuran mengikuti identitas dari \(\hat y_{i=1} = a_i^{'}\), maka koefisien regresi adalah nilai \(\hat y_i\) yang diharapkan dari komponen murni \(x_i\) dengan campuran antara 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 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 menintegrasi kendala campuran model Scheffe kuadrat, dapat diperoleh \[\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 Scheffe kuadrat yang umum untuk komponen K yang diberikan oleh persamaan \[y= ∑_{i=1}^K a_i^{'} \cdot x_i + ∑_{j>i}^K a_{ij}^{'} \cdot x_i \cdot x_j + \epsilon\]

Untuk setiap campuran bineri \(x_i = 0,5\), \(x_j = 0,5\) mengikuti persamaan diatas, dan kesetaraannya \[\hat y_{x_{i}=0.5;x_j=0.5} = \frac {a_i + a_j} {2} + \frac {a_{ij}} {4}\] yang digambarkan pada rumus efek utama schefe secara grafis besar untuk dua campuran bineri dengan aij> 0 (sinergisme) dan aij <0 (antagonisme). Garis lurus menunjukkan model Scheffe linier dengan aij = 0. Plot respon model Scheffe kuadrat untuk campuran bineri \(x_1\), \(x_2\) dengan efek sinergis dan antagonis berlabel merah dan biru.

Penghitungan untuk orde efek tinggi seperti interaksi terner dapat dicapai dengan menambah model Scheffe kuadrat dengan istilah interaksi orde tinggi seperti pada model Scheffe kubik penuh,

\[y = ∑_{i=1}^{K}a_{i}\cdot x_{i} + ∑_{j>i}^{K} a_{ij} \cdot x_{i} \cdot x_{j} + ∑_{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\]

Model campuran yang dijelaskan di atas sampai model kubik liniar, kuadrat, penuh dan tereduksi (khusus) sampai 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 (∑_{i=1}^Ka_i \cdot x_i + ∑_{j>i}^K a_{ij} \cdot x_i \cdot x_j \Bigg ) \cdot \Bigg (\alpha_0 + ∑_{l=1}^P \alpha_l \cdot z_l + ∑_{m>l}^R \alpha_{lm} \cdot z_l \cdot z_m \Bigg) + \epsilon\] dan dari ekspansi produk mengikuti langsung model parametrik, persamaannya adalah \[y = ∑_{i=1}^Kb_i \cdot x_i + ∑_{j>i}^K b_{ij} \cdot x_i \cdot x_j + ∑_{i=1}^K ∑_{l=1}^P b_{il}^{'} \cdot x_i \cdot z_l + ∑_{j>i}^K ∑_{l=1}^P b_{ijl}^{'} \cdot x_i \cdot x_j \cdot z_l + ∑_{i}^K ∑_{m>l}^P b_{ilm}^{'} \cdot x_{i} \cdot z_l \cdot z_m + ∑_{j>i}^K ∑_{m>l}^P b_{ijlm}^{'} \cdot x_i \cdot x_j \cdot z_l \cdot z_m + \epsilon\] Jumlah koefisien regresi dalam rumus diatas sebagai fungsi K dan P tumbuh secara besar dengan cepat, yang membuat desain silang 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 ini \[y = ∑_{i=1}^Kb_i \cdot x_i + ∑_{j>i}^K b_{ij} \cdot x_i \cdot x_j + ∑_{i=1}^K ∑_{l=1}^P b_{il}^{'} \cdot x_i \cdot z_l + ∑_{m>l}^P b_{lm}^{'} \cdot z_l \cdot z_m + ∑_{l=1}^P b_{ll} \cdot z_l^2 + \epsilon\] 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 mixexp paket 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 Desain Simplex-Centroid yang direplikasi sebagian (dijelaskan kemudian) ditambah dengan beberapa titik interior.

Desainnya secara maksimal mendukung model kubik yang dikurangi dan di sini akan berfungsi sebagai contoh untuk mendemonstrasikan beberapa fungsi dari mixe xp. Berbeda dari model ortogonal, dalam model campuran biasanya tidak mungkin untuk menghapus satu per satu istilah model yang tidak signifikan karena ini akan melanggar batasan campuran dan membuat model menjadi tunggal. Biasanya cukup untuk memilih di antara model di atas untuk mendeskripsikan data dengan tepat dalam domain model.

Data contoh ditunjukkan pada gambar sebagai plot terner dan mendukung model mana pun. Replikasi ditandai dengan lingkaran dan diimbangi dengan baik di atas ruang desain agar tidak mengganggu ortogonalitas desain. Plot terner dari data 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 dalam gambar (perhatikan bahwa \(R^2\) dalam lm-output dari model tanpa-intersep bias dan harus diabaikan. Nilai \(R^2\) yang dilaporkan oleh model campuran sudah 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.design

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

# Catat: R2 salah karena saat tidak ada intersep
# maka SS.model tidak disesuaikan dengan rata-rata
# Biasanya, 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

# pengurangan kubik
kappa(res1)

## [1] 9.353397

kappa(res2)

## [1] 57.53584

par(mfrow=c(1,2))
colcol                    <- rep("black") 
colcol[duplicated(etch[,1:3])]<- "red" #warna replikasi label
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()

# mengecek nlai residual dari model kubik yang dikurangi
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()

Residual pelajar versus prediksi model dari kuadrat (panel kiri) dan 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 Scheffe kuadrat, \(R^2 = 0,76\) lebih kecil dibandingkan dengan model kubik, $ R^2 = 0,98$, menunjukkan bahwa istilah \(x_1:x_2:x_3\) berkontribusi secara substansial dengan varian respons. Kesimpulan serupa dapat ditarik dari gambr yang mengungkapkan model kuadrat Scheffe menjadi bias, misalnya, tidak secara tepat menggambarkan pasangan replikasi # 7, # 11. Dalam model kubik Scheffe istilah terner \(x_1:x_2:x_3\) ditemukan signifikan dan residual dari model kubik sekarang lebih sesuai dengan asumsi normalitas dari istilah kesalahan ϵ. Model kubik dapat divisualisasikan dengan plot kontur terner yang ditunjukkan pada gambar. 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 cubic mixture model as ternary contour plot
ModelPlot(model = res,
          dimensions = list(x1 = "x1", x2 = "x2", x3 = "x3"),
          contour = TRUE, cuts = 6, fill = TRUE)

Plot kontur 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, 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 ruang campuran. Saat bergerak ke arah Cox, katakanlah, \(x_1\) rasio komponen campuran yang tersisa tetap konstan, di sini \(\frac {x_2} {x_3} = 1\), dan ini setara dengan efek “independen” sebesar \(x_1\) di ruang campuran. Efek Cox adalah jejak yang diperoleh dari sclicing permukaan respons sepanjang arah Cox. Gambar 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"))

Arah Cox pada ruang campuran 3D (panel kiri) dan plot efek sepanjang arah Cox model etch kubik tereduksi (panel kanan)

1.2 Mixture Design

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 dapat secara luas dibagi menjadi desain biasa di mana desain tersebut mencakup seluruh ruang campuran \((0≤xi≤1)\) dan desain tidak beraturan (atau dibatasi) di mana hanya subruang, \(LBi (> 0) ≤xi≤UBi (<1)\), dari seluruh ruang campuran terisi. Gambar menunjukkan contoh desain campuran beraturan dan tidak beraturan dalam tiga dimensi.

Contoh desain biasa (panel kiri) dan tidak teratur (dibatasi) (panel kanan).

Komponen murni (titik puncak pada gambar) cukup untuk memperkirakan model Scheffe linier. Model Scheffe kuadrat didukung dengan menambah desain sebelumnya dengan semua campuran bineri 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 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. Semua desain dapat dibuat dengan mudah dengan fungsi mixexp :: SLD () dan mixexp :: SCD () dan contoh untuk $ K = 3 $ yang pertama dan yang terakhir digambarkan dalam gambar.

Simplex Centroid Design SCD (3) (panel kiri) untuk model kubik yang dikurangi dan Simplex Lattice Design SLD tiga tingkat (3,3) (panel kanan) untuk model Scheffe kubik penuh.

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 + x_1 + x_2 + x_3 + x_1: x_2 + x_1: x_2: x_3 $, dengan enam angka dari set kandidat 60 reguler kandidat campuran.

Gambar menggambarkan titik desain optimal-D dalam ruang campuran yang terdiri dari komponen murni (mengacu pada efek utama \(x_1, x_2, x_3\)), campuran biner \(x_1,x_2\) (untuk istilah interaksi $ x_1: x_2 $) dan dua titik dekat dengan sentroid untuk mengidentifikasi interaksi tiga arah \(x_1:x_2:x_3\).

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) ) # desain grid
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)

Desain campuran D-optimal model \(~ -1 + x_1 + x_2 + x_3 + x_1: x_2 + x_1: x_2: x_3\) ditulis dalam notasi rumus R.

Desain untuk model proses campuran bersilangan penuh, persamaan, hasil dari campuran persilangan dan desain proses. Misalkan fac menjadi desain faktorial penuh atau pecahan dengan \(N_1\) baris dan \(P\) kolom dan menggabungkan desain campuran dimensi \((N_2) x (K)\), hasil kali Kartesius dari dimensi \((N_1 * N_2) x (P + K)\) mendukung desain proses campuran yang saling silang, seperti yang ditunjukkan koding R berikut untuk masing-masing\(K=3\) dan \(P=3\).

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)))) 
# memasukan faktor ke numerik
colnames(fac) <- paste("z",1:3,sep="")
x <- merge(mix,fac) # Produksi Cartesian
dim(x)

## [1] 80  6

x$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.35284

kappa(res)

## [1] 15.35277

Desain dan model yang saling bersilangan sangat besar dan hampir tidak dapat dilakukan dalam pengaturan laboratorium biasa, dan desain model yang sebagian bersilangan, mungkin menjadi alternatif yang lebih realistis untuk sebagian besar proyek. Desain proses campuran yang disilangkan sebagian dan dapat dihasilkan dengan memilih titik D secara optimal dari kumpulan kandidat yang sepenuhnya disilangkan 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)] # desain Box-Behnken;N=13
sld  <- mixexp::SLD(3,3) # desain kubik Simplex-lattice ;N=10
colnames(rsm) <- paste("z",1:3,sep="")
candidate <- merge(sld,rsm) # produksi Cartesian 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

kappa(res)

## [1] 24.00369

Desain terbatas atau tidak beraturan dapat dibangun dengan algoritma simpul ekstrim yang diimplementasikan dalam fungsi Xvert () dari paket mixexp dengan menentukan dimensi dan wilayah terbatas \((LBi≥0) ≤xi≤ (UBi≤1)\) dengan LBi dan UBi masing-masing menunjukkan batas atas dan bawah \(x_i\). Konsistensi membutuhkan jumlah dari batas atas melebihi 1, yaitu \(∑iUBi> 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 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)

Desain terbatas \(x_1 <0,3; x_2 <0,3; ∑ixi = 1\) (panel kiri) dan desain yang sesuai ditambah dengan titik-titik interior (panel kanan);

Tabel 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: Model dan desain campuran beserta petunjuk tentang cara membuat desain di R. Angka (1-6) didalam daftar secara langsung mengarah kepada model argumen di fungsi R mixexp :: MixModel (…, model = (1-6))

Akhirnya 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 fungsi-R AlgDesign :: optBlock () akan 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. 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 istilah 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. 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 membuat model menjadi tunggal karena blockB1 + blockB2 = 1 dan \(x_1 + x_2 + x_3 = 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       
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 

# tidak akan bekerja menjadi blockB1+blockB2=1 
# => dikompon dengan 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

kappa(res) # 1.90344e+17 desain singular

## [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

kappa(res) # 47.5957

## [1] 47.5957

Desain campuran kubik penuh didalam dua blok B1 (panel kiri) dan B2 (panel kanan)