Percobaan Mixture
Widi yantih
December 07, 2020
Semua desain yang dibahas pada bab 4 sebelumnya termasuk dalam kelas desain orthogonal di mana semua faktor bervariasi secara independen. Namun, dalam banyak aplikasi dan terutama dalam ilmu kimia ada kasus di mana faktor-faktor tidak dapat dan tidak boleh bervariasi secara independen Jika, katakanlah, ada tiga senyawa \(x1,x2,x3\) semua dalam rentang 0-1 dengan 1 menunjukkan komponen murni, desain dari pencampuran komponen ini harus selalu mematuhi aturan campuran \(∑_ix_i=1\) Secara ketat, semua masalah dalam kimia adalah masalah campuran, tetapi dalam banyak kasus pelarut diperlakukan sebagai “pengisi” inert tidak berkontribusi pada kendala campuran dan dengan demikian memungkinkan faktor konsentrasi yang tersisa bervariasi secara mandiri. Karena kendala campuran, masalah campuran membutuhkan model parametrik khusus dan sesuai dengan model desain khusus ini. Model campuran dan desain campuran akan diperkenalkan di bagian berikut dan penggunaannya dalam R akan ditunjukkan dengan beberapa contoh yang dipinjam dari literatur. Ikhtisar komprehensif tentang eksperimen campuran dan sumber yang baik untuk membaca lebih lanjut adalah buku (John A. Cornell 1990). (John Lawson, Cameron Willden 2016) memberikan gambaran tentang paket mixexp, kumpulan fungsi untuk membuat dan menganalisis desain campuran di R.
1 Mixture Models
Mencoba mengidentifikasi model OLS parametrik \(\hat{y}=a_0+a_1.x_1+a_2.x_2+a_3+x_3\) mengingat batasan campuran \(∑^3_{i=1} X_i=1\) akan gagal karena batasan campuran membuat satu variabel berlebihan. Salah satu cara mengatasi singularitas adalah dengan menjatuhkan satu variabel, katakan \(x3\) dan lanjutkan dengan model yang berkurang \(\hat{y}=a_0+a_1.x_1+a_2.x_2\) dan kondisi tambahan \(x_3=1-x_1-x_2\).Ini disebut pendekatan variabel kendur dengan variabel slack, di sini x3, biasanya memilih untuk menjadi komponen yang paling inert. Perhatikan bahwa variabel lain dapat dikecualikan tanpa mempengaruhi statistik model, karena semua model variabel slack secara statistik merosot dan akan memberikan prediksi yang sama.Cara yang lebih alami untuk mengatasi kendala linear adalah dengan mengintegrasikan campuran kontraints \(∑_ix_i=1\) ke dalam persamaan regresi, yang mengarah ke apa yang disebut model Scheffe, dinamai sesuai dengan ilmuwan yang pertama kali mendapatkan persamaan ini. Dari model parametrik \[\hat{y}=a_0+a_1.x_1+a_2.x_2+a_3.x_3\] dan Mixture \[x_1+x_2+x_3=1\] diikuti dengan substitusi
model efek utama Scheffe
yang diperluas untuk komponen campuran K ke model efek utama Scheffe, persamaan
Dari persamaan dan batasan campuran mengikuti identitas \(\hat{y_{x=i}}=1=a′\) sehingga koefisien regresi adalah nilai yang diharapkan \(\hat{y}_i\) komponen murni \(X_i\) dengan campuran menengah hanya menjadi jumlah tertimbang dari komponen murni. Dalam model Scheffe linear, komponen tidak berinteraksi dan ini mungkin merupakan anggapan terlalu ketat untuk banyak aplikasi sehingga memotivasi model Scheffe kuadrat. Dimulai dengan model kuadrat
dengan mengintegrasikan batasan campuran model Scheffe kuadrat diperoleh
dengan model Scheffe kuadrat umum untuk komponen K yang diberikan oleh persamaan
Untuk setiap campuran biner \(x_i=0,5\), \(x_j=0,5\) mengikuti dari persamaan kesetaraan
yang digambarkan pada gambar secara grafis untuk dua campuran biner dengan \(aij>0\) (sinergisme) dan \(aij<0\) (antagonisme). Garis lurus menunjukkan model Scheffe linear dengan \(aij=0\).
Plot respons model Scheffe kuadrat untuk campuran biner x1,x2 dengan efek sinergis dan antagonis berlabel merah dan biru, masing-masing
Akuntansi untuk efek pesanan yang lebih tinggi seperti interaksi ternary dapat dicapai dengan menambah model Scheffe kuadrat dengan istilah interaksi pesanan yang lebih tinggi seperti dalam model Scheffe kubik penuh, lihat persamaan
Mengikuti prinsip parsimony, data pertama kali dianalisis dengan model Scheffe kuadrat dan berikutnya dengan model kubik yang berkurang, dan residu sedang dibandingkan pada gambar 5.4 (perhatikan bahwa R2 dalam output lm dari model tanpa pencegatan bias dan harus diabaikan. Nilai R2 yang dilaporkan oleh MixModel sudah benar dan harus digunakan sebagai gantinya).
## 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##
## 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## [1] 9.353397## [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()
Residual siswa versus prediksi model dari quadratic (panel kiri) dan mengurangi model Scheffe kubik data etsa. Pasangan replikasi (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\) kecil dibandingkan dengan model kubik, \(R^2 = 0,98\), menunjukkan bahwa istilah x1:x2:x3 berkontribusi secara substansial terhadap varian respons. Kesimpulan serupa dapat diambil dari gambar yang mengungkapkan model Scheffe kuadrat menjadi bias oleh, misalnya, tidak tepat menggambarkan pasangan replikasi #7.#11. Dalam model Scheffe kubik istilah ternary \(x1:x2:x3\) ditemukan signifikan dan residu model kubik sekarang lebih baik sesuai dengan asumsi normalitas dari istilah kesalahan \(Ε\) . Model kubik dapat divisualisasikan dengan plot kontur ternary yang ditunjukkan pada gambar. Maksimum terletak dekat dengan sentroid \((x_1=0,33,x_2=0,33,x_3=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)
plot kontur juga dapat digunakan untuk dimensi yang lebih tinggi dari K = 3 dengan mengkondisikan (mengiris) plot pada variabel yang ditetapkan konstanta. Namun, ini bisa menjadi membingungkan ketika ada banyak irisan untuk dibandingkan, dan karena itu cara lain yang populer, dan kurang membingungkan untuk memvisualisasikan model campuran dimensi tinggi adalah plot efek Cox.
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 dalam ruang campuran 3D (panel kiri) dan plot efek di sepanjang arah Cox dari model etsa kubik yang berkurang (panel kanan)
2 Mixture designs
Bagian sebelumnya memperkenalkan enam model parametrik yang berbeda untuk pemodelan campuran dan sistem proses campuran. Terkait dengan setiap model adalah desain proses campuran dan campuran secara optimal mendukung model yang dipilih. Model dan desain adalah dua sisi koin yang sama: Model parametrik menentukan desain dan desain menentukan dan membatasi model yang dapat diperkirakan dari data. Desain campuran dapat dibagi secara luas menjadi desain reguler di mana desain mencakup seluruh ruang campuran \((0≤XⅠ≤1)\) dan desain tidak teratur (atau dibatasi) di mana hanya subspace, \(LBⅠ(>0)≤XⅠ≤UBⅠ(<1\)) , dari seluruh ruang campuran dihuni. Gambar menunjukkan contoh desain campuran biasa dan tidak teratur dalam tiga dimensi.
Model dan desain campuran klasik kaku dalam hal struktur model dan ukuran desain yang mirip dengan apa yang telah dikatakan tentang desain ortogonal klasik. Urutan kedua Simplex-Lattice Design misalnya akan membuat semua interaksi dua arah dapat diestimasi, terlepas dari apakah mereka menarik atau mungkin secara fisik. Struktur model dan ukuran desain dapat ditentukan lebih fleksibel dengan desain optimal seperti yang ditunjukkan dengan R-code berikut. Di sini tujuannya adalah untuk secara parsimoniously mengidentifikasi model Scheffe non-standar, dalam R-notation \(~ -1 + x1 + x2 + x3 + x1:x2 + x1:x2:x3\), dengan enam run dari kumpulan kandidat dari 60 kandidat campuran reguler. Gambar menggambarkan titik desain D-optimal 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 poin yang dekat dengan sentroid 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)
Desain untuk model proses campuran yang sepenuhnya disilangkan (#5), persamaan (5,5), hasil dari campuran persimpangan dan desain proses. Biarkan fac menjadi desain faktorial penuh atau fraksional dengan baris N1 dan kolom P dan mencampur desain campuran dimensi (N2)x(K), produk dimensi Cartesian (N1 * N2)x (P + K) mendukung desain proses campuran yang sepenuhnya disilangkan sebagai R-code berikut menunjukkan untuk K = 3 dan P = 3, masing-masing
## 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.35284## [1] 15.35277Desain dan model yang sepenuhnya disilangkan sangat besar dan hampir tidak layak dalam pengaturan laboratorium biasa, dan desain yang disilangkan sebagian, model #6 (persamaan (5,6))), mungkin menjadi alternatif yang lebih realistis untuk sebagian besar proyek. Desain proses campuran yang disilangkan sebagian dan parsimonious dapat dihasilkan dengan memilih titik-titik yang dipilih secara optimal dari set 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)] # 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.1727## [1] 24.00369rm(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)
Akhirnya dan mirip 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-prinsip pertama (untuk detailnya lihat (John A. Cornell 1990)), pemblokiran algoritmik dengan fungsi R AlgDesign::optBlock() akan cukup untuk tujuan yang paling praktis. Contoh kode berikut akan memblokir desain Simplex-lattic urutan ketiga dengan sepuluh berjalan menjadi dua blok dan menggambarkan kedua blok dalam gambar 5.12. Selain itu, kode contoh menunjukkan beberapa work-around untuk menghindari beberapa masalah yang tidak dicakup oleh fungsionalitas standar. Pertama, tentukan variabel xij=xi-xj untuk menentukan istilah kubik xi:xj:(xi-xj) dalam ekspresi rumus. Kedua, tentukan variabel baru menurut block2=as.integer(block==“B2”) sebagai indikator blok unik. Ini diperlukan karena lm() secara default menjatuhkan tingkat pertama faktor (di sini B1 di block={B1,B2}) ketika model memiliki intersep. Namun, ketika tidak ada intersep baik blockB1 maupun blockB2 disimpan dalam matriks model sehingga membuat model tunggal karena blockB1+blockB2=1 dan x1+x2+x3=1. Masalah ini disunat dengan block variable block2=as.integer baru(block==“B2”). Kemudian parameter lokasi dari model lm kubik penuh setuju dengan baik dengan model yang sebenarnya meskipun fakta bahwa tidak ada istilah yang ditemukan signifikan. Namun, ini tidak surpsising karena hanya ada 1 tingkat kebebasan yang membuat statistik tes 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.009553## [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.009553## [1] 47.5957