実験計画法 (Design of Experiment) とは何か？ (7)

RSM, Response Surface Methodology（応答曲面法）を利用した実験計画と解析例を、まず、今までの紹介した実験例の延長で扱ってみたいと思います。というのは RSM における実験計画には、回帰分析をすることを前提とした計画があるためですが、それについては次回以降に説明をする予定です。

3 因子 3 水準の実験（繰り返しなし）

今回取り上げる例は、「高」「低」や「大」「小」といった離散的な水準あるいは「A 社品」「B 社品」といった質的な量ではなく、連続量を設定できる因子の水準を扱います。

条件	温度	流量	回転数	収率
1	80	30	250	81.64
2	80	30	300	80.14
3	80	30	350	81.08
4	80	40	250	80.57
5	80	40	300	78.77
6	80	40	350	78.95
7	80	50	250	80.53
8	80	50	300	77.88
9	80	50	350	76.18
10	100	30	250	81.68
11	100	30	300	81.14
12	100	30	350	82.61
13	100	40	250	82.60
14	100	40	300	80.19
15	100	40	350	80.04
16	100	50	250	81.85
17	100	50	300	79.65
18	100	50	350	77.89
19	120	30	250	83.53
20	120	30	300	83.53
21	120	30	350	81.95
22	120	40	250	83.04
23	120	40	300	82.36
24	120	40	350	81.92
25	120	50	250	83.60
26	120	50	300	81.79
27	120	50	350	79.80

例）

ある製品の製造工程において、「温度」、「流量」、「回転数」の 3 つの因子が収量に影響を与えるかどうかを調べるため、3 水準の実験を計画し、各条件の収量を測定した結果を右に示した。この 3 因子が収量に及ぼす影響を分析せよ。ただし、過去の知見から「流量」と「回転数」には交互作用のあることが判っているものとする。

右のテーブル全体を選択して（クリップボードへ）コピーしてから、R 上で次のコマンドを実行してください。クリップボードの値が R のデータフレーム（テーブル）へ読み込まれます。

> tbl <- read.table("clipboard", header = T)
>

> tbl
   条件 温度 流量 回転数  収率
1     1   80   30    250 81.64
2     2   80   30    300 80.14
3     3   80   30    350 81.08
4     4   80   40    250 80.57
5     5   80   40    300 78.77
6     6   80   40    350 78.95
7     7   80   50    250 80.53
8     8   80   50    300 77.88
9     9   80   50    350 76.18
10   10  100   30    250 81.68
11   11  100   30    300 81.14
12   12  100   30    350 82.61
13   13  100   40    250 82.60
14   14  100   40    300 80.19
15   15  100   40    350 80.04
16   16  100   50    250 81.85
17   17  100   50    300 79.65
18   18  100   50    350 77.89
19   19  120   30    250 83.53
20   20  120   30    300 83.53
21   21  120   30    350 81.95
22   22  120   40    250 83.04
23   23  120   40    300 82.36
24   24  120   40    350 81.92
25   25  120   50    250 83.60
26   26  120   50    300 81.79
27   27  120   50    350 79.80
> 

パッケージ rsm を読み込んで、回帰分析をするために、tbl の因子水準を -1, 0, 1 にコード化します。

> library("rsm")
> (tbl.coded <- coded.data(tbl, x1 ~ (温度 - 100)/20, x2 ~ (流量 - 40) / 10, x3 ~ (回転数 - 300) / 50))
   温度 流量 回転数  収率
1    80   30    250 81.64
2    80   30    300 80.14
3    80   30    350 81.08
4    80   40    250 80.57
5    80   40    300 78.77
6    80   40    350 78.95
7    80   50    250 80.53
8    80   50    300 77.88
9    80   50    350 76.18
10  100   30    250 81.68
11  100   30    300 81.14
12  100   30    350 82.61
13  100   40    250 82.60
14  100   40    300 80.19
15  100   40    350 80.04
16  100   50    250 81.85
17  100   50    300 79.65
18  100   50    350 77.89
19  120   30    250 83.53
20  120   30    300 83.53
21  120   30    350 81.95
22  120   40    250 83.04
23  120   40    300 82.36
24  120   40    350 81.92
25  120   50    250 83.60
26  120   50    300 81.79
27  120   50    350 79.80

Data are stored in coded form using these coding formulas ...
x1 ~ (温度 - 100)/20
x2 ~ (流量 - 40)/10
x3 ~ (回転数 - 300)/50
> print(tbl.coded, decode = FALSE)
   条件 x1 x2 x3  収率
1     1 -1 -1 -1 81.64
2     2 -1 -1  0 80.14
3     3 -1 -1  1 81.08
4     4 -1  0 -1 80.57
5     5 -1  0  0 78.77
6     6 -1  0  1 78.95
7     7 -1  1 -1 80.53
8     8 -1  1  0 77.88
9     9 -1  1  1 76.18
10   10  0 -1 -1 81.68
11   11  0 -1  0 81.14
12   12  0 -1  1 82.61
13   13  0  0 -1 82.60
14   14  0  0  0 80.19
15   15  0  0  1 80.04
16   16  0  1 -1 81.85
17   17  0  1  0 79.65
18   18  0  1  1 77.89
19   19  1 -1 -1 83.53
20   20  1 -1  0 83.53
21   21  1 -1  1 81.95
22   22  1  0 -1 83.04
23   23  1  0  0 82.36
24   24  1  0  1 81.92
25   25  1  1 -1 83.60
26   26  1  1  0 81.79
27   27  1  1  1 79.80

Variable codings ...
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  
            80.54815               1.43222              -1.00722  
    FO(x1, x2, x3)x3  TWI(x1, x2, x3)x1:x2  TWI(x1, x2, x3)x1:x3  
            -1.03444               0.37083               0.00250  
TWI(x1, x2, x3)x2:x3    PQ(x1, x2, x3)x1^2    PQ(x1, x2, x3)x2^2  
            -0.90833               0.10889              -0.02278  
  PQ(x1, x2, x3)x3^2  
             0.47556  

>

重回帰分析の統計値の算出は、通常の場合と同じように summary を使います。

> summary(tbl.rsm)

Call:
rsm(formula = 収率 ~ SO(x1, x2, x3), data = tbl.coded)

             Estimate Std. Error  t value  Pr(>|t|)    
(Intercept) 80.548148   0.274184 293.7743 < 2.2e-16 ***
x1           1.432222   0.126922  11.2842 2.563e-09 ***
x2          -1.007222   0.126922  -7.9357 4.073e-07 ***
x3          -1.034444   0.126922  -8.1502 2.831e-07 ***
x1:x2        0.370833   0.155448   2.3856   0.02896 *  
x1:x3        0.002500   0.155448   0.0161   0.98736    
x2:x3       -0.908333   0.155448  -5.8433 1.956e-05 ***
x1^2         0.108889   0.219836   0.4953   0.62672    
x2^2        -0.022778   0.219836  -0.1036   0.91869    
x3^2         0.475556   0.219836   2.1632   0.04506 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Multiple R-squared:  0.9466,    Adjusted R-squared:  0.9184 
F-statistic:  33.5 on 9 and 17 DF,  p-value: 4.426e-09

Analysis of Variance Table

Response: 収率
                Df Sum Sq Mean Sq F value    Pr(>F)
FO(x1, x2, x3)   3 74.445 24.8150 85.5786 1.839e-10
TWI(x1, x2, x3)  3 11.551  3.8504 13.2786 0.0001029
PQ(x1, x2, x3)   3  1.431  0.4771  1.6452 0.2163396
Residuals       17  4.929  0.2900                  
Lack of fit     17  4.929  0.2900                  
Pure error       0  0.000                          

Stationary point of response surface:
       x1        x2        x3 
-1.672520 -2.868849 -1.647805 

Stationary point in original units:
     温度      流量    回転数 
 66.54960  11.31151 217.60974 

Eigenanalysis:
$values
[1]  0.7583657  0.1520116 -0.3487107

$vectors
         [,1]      [,2]       [,3]
x1 -0.1477192 0.9303944  0.3354777
x2 -0.5230889 0.2143790 -0.8248755
x3  0.8393789 0.2973346 -0.4550111

上の赤字の部分は Windows 版 (R 3.1.0) の場合に表示されるのですが、Linux 版 (Fedora 20, R 3.0.2) では表示されません。表示されない理由は不明ですが、重相関係数 (R²) などは算出されています。

> names(summary(tbl.rsm))
 [1] "call"          "terms"         "residuals"     "coefficients"  "aliased"       "sigma"         "df"            "r.squared"    
 [9] "adj.r.squared" "fstatistic"    "cov.unscaled"  "canonical"     "lof"           "coding"       
> summary(tbl.rsm)$r.squared
[1] 0.946626

通常の解析では、ここで変数選択（ステップワイズ）をしてフィッティングを良くするのですが、その場合、他のパッケージを利用する必要がありますので、今回はこのまま続けます。なお、回帰式の変数 x1, x2, x3 は、それぞれ因子 温度, 流量, 回転数の水準が -1, 0, 1 にコード化されていますので、それぞれの係数の推定値の大きさが、応答に及ぼす大きさであると評価することができます。主効果であれば、x1, x2, x3、交互効果であれば x1:x2, x1:x3, x2:x3 の係数の推定値の大きさになります。

残渣などからを回帰分析の善し悪しを診断します。なお、この例は意図的に作ってあるため、極端な外れ値は存在しません。

> par (mfrow = c (2,2))
> plot(tbl.rsm)
>

応答の等高線図 (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)) 
>

応答のターゲットあるいは極小値/極大値を見つけ出すコマンドは、rsm パッケージに用意されていませんが、応答曲面の傾斜に沿って、最も降下あるいは上昇するラインの応答の予測値 (yhat) を算出できますので、求める応答の値を実現する因子の組合せを探索する出発点として利用できます。

> steepest(tbl.rsm)
Path of steepest ascent from ridge analysis:
   dist    x1     x2     x3 |   温度  流量 回転数 |    yhat
1   0.0 0.000  0.000  0.000 | 100.00 40.00 300.00 |  80.548
2   0.5 0.370 -0.150 -0.301 | 107.40 38.50 284.95 |  81.536
3   1.0 0.677 -0.038 -0.735 | 113.54 39.62 263.25 |  82.587
4   1.5 0.864  0.205 -1.209 | 117.28 42.05 239.55 |  83.893
5   2.0 0.997  0.469 -1.669 | 119.94 44.69 216.55 |  85.538
6   2.5 1.106  0.737 -2.118 | 122.12 47.37 194.10 |  87.549
7   3.0 1.204  1.004 -2.558 | 124.08 50.04 172.10 |  89.927
8   3.5 1.296  1.270 -2.993 | 125.92 52.70 150.35 |  92.681
9   4.0 1.383  1.535 -3.426 | 127.66 55.35 128.70 |  95.815
10  4.5 1.467  1.800 -3.855 | 129.34 58.00 107.25 |  99.320
11  5.0 1.549  2.064 -4.282 | 130.98 60.64  85.90 | 103.198
> 

RSM の目的

RSM の目的は、真の理論式を探し出すことではありません。実験範囲において多項式などの回帰モデルを用いて、できるだけ精度よく現象を近似すること、さらには欲しい最適な水準の組合せ（最適条件）を探し出すことです。