RSM, Response Surface Methodology（応答曲面法）を利用した実験計画と解析例について、前回は、今までの流れから、因子の水準の全ての組合せの計画（完全実施要因計画, Full Factorial Design）を例に取り上げて解析してみましたが、今回は RSM でよく使われる「中心複合計画」について説明します。※ 説明が不十分ですが、まずは公開して、徐々に書き足していくことをご了承下さい。
中心複合計画（Box-Wilson 法）
中心複合計画 Central Composite Design は、 二次の多項式モデルを用いて最適な条件を探索するのに適した計画の一つです。
中心複合計画の特徴は以下の通りです。
- 各因子について、2 水準の完全実施計画（あるいは部分実施計画）。一次の項と 2 因子の交互作用効果の推定に寄与。
- 軸上の 2k 個の計画点。ある点は、ある因子について水準の α 倍ないしは −α 倍、他の因子について 0 とした水準。二次の項の推定に寄与。
- n 個の中心条件。誤差の推定と二次項の推定に寄与。
3 因子の計画の場合、左図のような水準の組み合わせになります。
実験計画法 (Design of Experiment) とは何か？ (7) で扱った例は、以下のような因子と水準範囲でした。
- 温度［℃］
- 80 - 120
- 流量［リットル毎分］
- 30 - 50
- 回転数［回転毎分］
- 250 - 350
これに中心複合計画を適用すると以下のようになります。中心条件の繰り返し数は 3 としました。
> library("rsm") > (tbl <- ccd(3, n0 = c(3,0), alpha = "rotatable", randomize = FALSE, inscribed = FALSE, coding = list (x1 ~ (温度 - 100)/20, x2 ~ (流量 - 40)/10, x3 ~ (回転数 - 300)/50), oneblock = TRUE)) run.order std.order 温度 流量 回転数 1 1 1 80.00000 30.00000 250.0000 2 2 2 120.00000 30.00000 250.0000 3 3 3 80.00000 50.00000 250.0000 4 4 4 120.00000 50.00000 250.0000 5 5 5 80.00000 30.00000 350.0000 6 6 6 120.00000 30.00000 350.0000 7 7 7 80.00000 50.00000 350.0000 8 8 8 120.00000 50.00000 350.0000 9 9 9 100.00000 40.00000 300.0000 10 10 10 100.00000 40.00000 300.0000 11 11 11 100.00000 40.00000 300.0000 12 1 1 66.36414 40.00000 300.0000 13 2 2 133.63586 40.00000 300.0000 14 3 3 100.00000 23.18207 300.0000 15 4 4 100.00000 56.81793 300.0000 16 5 5 100.00000 40.00000 215.9104 17 6 6 100.00000 40.00000 384.0896 Data are stored in coded form using these coding formulas ... x1 ~ (温度 - 100)/20 x2 ~ (流量 - 40)/10 x3 ~ (回転数 - 300)/50 >
| 条件 | 温度 | 流量 | 回転数 | 収率 |
|---|---|---|---|---|
| 1 | 80 | 30 | 250 | 80.18 |
| 2 | 120 | 30 | 250 | 82.73 |
| 3 | 80 | 50 | 250 | 79.79 |
| 4 | 120 | 50 | 250 | 82.49 |
| 5 | 80 | 30 | 350 | 80.81 |
| 6 | 120 | 30 | 350 | 82.12 |
| 7 | 80 | 50 | 350 | 76.44 |
| 8 | 120 | 50 | 350 | 79.90 |
| 9 | 100 | 40 | 300 | 81.01 |
| 10 | 100 | 40 | 300 | 80.98 |
| 11 | 100 | 40 | 300 | 81.95 |
| 12 | 66 | 40 | 300 | 79.63 |
| 13 | 134 | 40 | 300 | 82.17 |
| 14 | 100 | 23 | 300 | 81.76 |
| 15 | 100 | 57 | 300 | 80.26 |
| 16 | 100 | 40 | 216 | 82.33 |
| 17 | 100 | 40 | 384 | 79.63 |
右のテーブル全体を選択して（クリップボードへ）コピーしてから、R 上で次のコマンドを実行してください。クリップボードの値が R のデータフレーム（テーブル）へ読み込まれます。
> tbl <- read.table("clipboard", header = T)
>
> tbl
条件 温度 流量 回転数 収率
1 1 80 30 250 80.18
2 2 120 30 250 82.73
3 3 80 50 250 79.79
4 4 120 50 250 82.49
5 5 80 30 350 80.81
6 6 120 30 350 82.12
7 7 80 50 350 76.44
8 8 120 50 350 79.90
9 9 100 40 300 81.01
10 10 100 40 300 80.98
11 11 100 40 300 81.95
12 12 66 40 300 79.63
13 13 134 40 300 82.17
14 14 100 23 300 81.76
15 15 100 57 300 80.26
16 16 100 40 216 82.33
17 17 100 40 384 79.63
>
パッケージ rsm を読み込んで、回帰分析をするために、tbl の因子水準を -1, 0, 1 にコード化します。
> (tbl.coded <- coded.data(tbl, x1 ~ (温度 - 100)/20, x2 ~ (流量 - 40) / 10, x3 ~ (回転数 - 300) / 50))
条件 温度 流量 回転数 収率
1 1 80 30 250 80.18
2 2 120 30 250 82.73
3 3 80 50 250 79.79
4 4 120 50 250 82.49
5 5 80 30 350 80.81
6 6 120 30 350 82.12
7 7 80 50 350 76.44
8 8 120 50 350 79.90
9 9 100 40 300 81.01
10 10 100 40 300 80.98
11 11 100 40 300 81.95
12 12 66 40 300 79.63
13 13 134 40 300 82.17
14 14 100 23 300 81.76
15 15 100 57 300 80.26
16 16 100 40 216 82.33
17 17 100 40 384 79.63
Data are stored in coded form using these coding formulas ...
x1 ~ (温度 - 100)/20
x2 ~ (流量 - 40)/10
x3 ~ (回転数 - 300)/50
>
コード化したデータを使って高々二次の線形モデル (SO) にフィッティング（回帰分析）をします。
> (tbl.rsm <- rsm(収率 ~ SO(x1, x2, x3), data = tbl.coded))
Call:
rsm(formula = 収率 ~ SO(x1, x2, x3), data = tbl.coded)
Coefficients:
(Intercept) FO(x1, x2, x3)x1 FO(x1, x2, x3)x2
81.3329 1.0405 -0.7090
FO(x1, x2, x3)x3 TWI(x1, x2, x3)x1:x2 TWI(x1, x2, x3)x1:x3
-0.7663 0.2875 -0.0600
TWI(x1, x2, x3)x2:x3 PQ(x1, x2, x3)x1^2 PQ(x1, x2, x3)x2^2
-0.7450 -0.2255 -0.1874
PQ(x1, x2, x3)x3^2
-0.2044
>
重回帰分析の統計値の算出は、通常の場合と同じように summary を使います。
> summary(tbl.rsm) Call: rsm(formula = 収率 ~ SO(x1, x2, x3), data = tbl.coded) Estimate Std. Error t value Pr(>|t|) (Intercept) 81.33287 0.38906 209.0510 1.513e-14 *** x1 1.04049 0.18174 5.7251 0.0007165 *** x2 -0.70900 0.18174 -3.9011 0.0058903 ** x3 -0.76630 0.18264 -4.1957 0.0040573 ** x1:x2 0.28750 0.23852 1.2053 0.2672386 x1:x3 -0.06000 0.23852 -0.2515 0.8086174 x2:x3 -0.74500 0.23852 -3.1234 0.0167648 * x1^2 -0.22548 0.19782 -1.1398 0.2918503 x2^2 -0.18742 0.19782 -0.9474 0.3749705 x3^2 -0.20440 0.20184 -1.0127 0.3449259 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Multiple R-squared: 0.9185, Adjusted R-squared: 0.8137 F-statistic: 8.766 on 9 and 7 DF, p-value: 0.004582 Analysis of Variance Table Response: 収率 Df Sum Sq Mean Sq F value Pr(>F) FO(x1, x2, x3) 3 29.8579 9.9526 21.8665 0.0006232 TWI(x1, x2, x3) 3 5.1303 1.7101 3.7571 0.0678184 PQ(x1, x2, x3) 3 0.9213 0.3071 0.6747 0.5944141 Residuals 7 3.1861 0.4552 Lack of fit 5 2.5776 0.5155 1.6945 0.4112877 Pure error 2 0.6085 0.3042 Stationary point of response surface: x1 x2 x3 1.453147 -1.286240 0.256269 Stationary point in original units: 温度 流量 回転数 129.0629 27.1376 312.8134 Eigenanalysis: $values [1] 0.2123273 -0.2428227 -0.5868074 $vectors [,1] [,2] [,3] x1 -0.2766053 0.93416169 0.2254583 x2 -0.7064912 -0.03863937 -0.7066663 x3 0.6514290 0.35475192 -0.6706649 >
応答の等高線図 (Contour Map) をプロットしてみましょう。
> par (mfrow = c (2,3)) > contour(tbl.rsm, ~x1+x2+x3, at = xs(tbl.rsm)) > persp (tbl.rsm, ~x1+x2+x3, at = xs(tbl.rsm)) >
0 件のコメント:
コメントを投稿