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

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)) 
>

参考サイト

1. Central composite design - Wikipedia, the free encyclopedia
2. 5.3.3.6.1. Central Composite Designs (CCD)
3. 応答曲面法のための実験計画

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

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

RSM, Response Surface Methodology（応答曲面法）とは回帰分析の一種で、目的とする変数に対する複数の因子の影響について、有限のデータから連続的な表面として近似させたものです。 複数変数の最適組み合わせを導くのに使われる実験計画法 (DOE) です。今までの分散分析中心の実験計画法から、ひとまず RSM へ話題を移します。

パッケージ rsm

CRAN からパッケージ rsm をインストールして利用できるようにします。

• CRAN - Package rsm

R を起動して、以下のようにして rsm パケージをインストールします。

> install.packages("rsm")
Installing package into ‘/home/bitwalk/R/x86_64-redhat-linux-gnu-library/3.0’
(as ‘lib’ is unspecified)
 --- このセッションで使うために、CRAN のミラーサイトを選んでください --- 

右のように CRAN のミラーサイト一覧が GUI で表示されるので、近いと思われる適当なサイトを選択して、ダブルクリックあるいは下の「OK」ボタンをクリックします。。

 URL 'http://cran.ism.ac.jp/src/contrib/rsm_2.04.tar.gz' を試しています 
Content type 'application/x-gzip' length 1013015 bytes (989 Kb)
 開かれた URL 
==================================================
downloaded 989 Kb

* installing *source* package ‘rsm’ ...
**  パッケージ ‘rsm’ の解凍および MD5 サムの検証に成功しました 
** R
** inst
** preparing package for lazy loading
** help
*** installing help indices
  converting help for package ‘rsm’
    finding HTML links ...  完了 
    ChemReact                               html  
    FO                                      html  
    bbd                                     html  
    ccd                                     html  
    ccd.pick                                html  
    codata                                  html  
    coded.data                              html  
    contour.lm                              html  
    djoin                                   html  
    heli                                    html  
    model.data                              html  
    rsm-package                             html  
    rsm                                     html  
    steepest                                html  
    varfcn                                  html  
** building package indices
** installing vignettes
** testing if installed package can be loaded
* DONE (rsm)

 ダウンロードされたパッケージは、以下にあります 
  ‘/tmp/RtmpSQBlvV/downloaded_packages’ 

以上で無事、パッケージ rsm がインストールされましたので、早速、rsm のマニュアルで紹介されているサンプルから、今後使うことになるプロットを紹介しましょう。

> library('rsm')
> heli.rsm <- rsm (ave ~ block + SO(x1, x2, x3, x4), data = heli)
> heli.rsm

Call:
rsm(formula = ave ~ block + SO(x1, x2, x3, x4), data = heli)

Coefficients:
             (Intercept)                    block2      FO(x1, x2, x3, x4)x1  
               372.80000                  -2.95000                  -0.08333  
    FO(x1, x2, x3, x4)x2      FO(x1, x2, x3, x4)x3      FO(x1, x2, x3, x4)x4  
                 5.08333                   0.25000                  -6.08333  
TWI(x1, x2, x3, x4)x1:x2  TWI(x1, x2, x3, x4)x1:x3  TWI(x1, x2, x3, x4)x1:x4  
                -2.87500                  -3.75000                   4.37500  
TWI(x1, x2, x3, x4)x2:x3  TWI(x1, x2, x3, x4)x2:x4  TWI(x1, x2, x3, x4)x3:x4  
                 4.62500                  -1.50000                  -2.12500  
  PQ(x1, x2, x3, x4)x1^2    PQ(x1, x2, x3, x4)x2^2    PQ(x1, x2, x3, x4)x3^2  
                -2.03750                  -1.66250                  -2.53750  
  PQ(x1, x2, x3, x4)x4^2  
                -0.16250  

> par (mfrow = c (2,3))
> contour (heli.rsm, ~x1+x2+x3+x4, at = xs(heli.rsm))

> contour (heli.rsm, ~x1+x2+x3+x4, at = list(block="2"),
+ atpos = 0, image = TRUE)

> persp (heli.rsm, ~x1+x2+x3+x4, at = xs(heli.rsm))

> persp (heli.rsm, ~x1+x2+x3+x4, at = xs(heli.rsm),
+ contours = "col", col = rainbow(40), zlab = "Flight time",
+ xlabs = c("Wing area", "Wing length", "Body width", "Body length"))
> 

詳細な使い方については、今後、サンプルを使って説明していく予定です。

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

今回は 3 因子の実験です。因子数と水準数が増えると実験数は飛躍的に大きくなりますので、まずは繰り返し実験がなく、水準数を 2 にしたシンプルな例を取り上げます。

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

条件	温度	流量	回転数	収率
1	1低	1小	1低	81.08
2	1低	1小	2高	80.66
3	1低	2大	1低	80.73
4	1低	2大	2高	76.14
5	2高	1小	1低	83.39
6	2高	1小	2高	81.98
7	2高	2大	1低	83.83
8	2高	2大	2高	78.81

例）

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

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

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

因子「流量」と「回転数」の間には交互作用が存在するとして、aov() 関数で分散分析をしてみましょう。

> result.aov <- aov(収率 ~ 温度 + 流量 * 回転数, data = tbl)
> summary(result.aov)
            Df Sum Sq Mean Sq F value  Pr(>F)   
温度         1 11.045  11.045   38.36 0.00848 **
流量         1  7.220   7.220   25.08 0.01533 * 
回転数       1 16.359  16.359   56.82 0.00484 **
流量:回転数  1  7.566   7.566   26.28 0.01437 * 
Residuals    3  0.864   0.288                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>

「収率」に対する各主因子と交互作用因子による影響の差は、有意水準 5% で有意であるという結果になりました。

各水準の平均

model.tables() で、各水準の "平均値" を算出します。

> mm <- model.tables(result.aov, type = "mean")
> mm
Tables of means
Grand mean
        
80.8275 

 温度 
温度
  1低   2高 
79.65 82.00 

 流量 
流量
  1小   2大 
81.78 79.88 

 回転数 
回転数
  1低   2高 
82.26 79.40 

 流量:回転数 
     回転数
流量 1低   2高  
  1小 82.23 81.32
  2大 82.28 77.47
>

上記の数字をグラフにしてみます。

> ylist <- c(mm$tables$温度, mm$tables$流量, mm$tables$回転数, mm$tables$"流量: 回転数")
> plot(0, 0, xlim = c(1, 8), ylim = c(min(ylist), max(ylist)), xaxt = "n", xlab = "", ylab="収 率")
> points(1:6, c(mm$tables$温度, mm$tables$流量, mm$tables$回転数), cex = 2, pch = 16)
> lines(1:2, mm$tables$温度)
> axis(1, at = 1:2, c(names(mm$tables$温度)))
> axis(1, at = 1.5, "温度", tick = F, padj = 1.5) 
> lines(3:4, mm$tables$流量)
> axis(1, at = 3:4, c(names(mm$tables$流量)))
> axis(1, at = 3.5, "流量", tick = F, padj = 1.5)
> lines(5:6, mm$tables$回転数)
> axis(1, at = 5:6, c(names(mm$tables$回転数)))
> axis(1, at = 5.5, "回転数", tick = F, padj = 1.5)
> points(7:8, c(mm$tables$"流量:回転数"[1, ]), cex = 2, pch = 0)
> lines(7:8, mm$tables$"流量:回転数"[1, ], lty = 2)
> points(7:8, c(mm$tables$"流量:回転数"[2, ]), cex = 2, pch = 2)
> lines(7:8, mm$tables$"流量:回転数"[2, ], lty = 2)
> axis(1, at = 7:8, c(names(mm$tables$"流量:回転数"[1, ])))
> axis(1, at = 7.5, "（回転数）", tick = F, padj = 1.5)
> axis(1, at = 7.5, "流量×回転数", tick = F, padj = 3)
> text(7.5, mm$tables$"流量:回転数"[1, 2] + 0.5, labels = "（流量）1小")
> text(7.5, mm$tables$"流量:回転数"[2, 2] + 0.5, labels = "（流量）1大")
>

ここで各水準の "平均値" について、データフレーム tbl を使って、いくつか検算をしておきましょう。

まず Grand mean ですが、これは単純に、全収率の平均値（全平均）になります。

> gm <-mean(tbl$収率)
[1] 80.8275

次に、3 因子の各水準の平均ですが、例えば、因子「温度」の場合、次のように算出します。他の因子も同じように算出します。

> mean(tbl[tbl$温度 == "1低", "収率"])
[1] 79.6525
> mean(tbl[tbl$温度 == "2高", "収率"])
[1] 82.0025

因子「流量」と「回転数」の交互作用については、次のような計算になります。

> mean(tbl[(tbl$流量 == "1小") & (tbl$回転数 == "1低"), "収率"])
[1] 82.235
> mean(tbl[(tbl$流量 == "1小") & (tbl$回転数 == "2高"), "収率"])
[1] 81.32
> mean(tbl[(tbl$流量 == "2大") & (tbl$回転数 == "1低"), "収率"])
[1] 82.28
> mean(tbl[(tbl$流量 == "2大") & (tbl$回転数 == "2高"), "収率"])
[1] 77.475
>

効果の算出

model.tables() で、各水準の "効果" を算出します。

> ee <- model.tables(result.aov, type = "effects")
> ee
Tables of effects

 温度 
温度
   1低    2高 
-1.175  1.175 

 流量 
流量
  1小   2大 
 0.95 -0.95 

 回転数 
回転数
  1低   2高 
 1.43 -1.43 

 流量:回転数 
     回転数
流量 1低     2高    
  1小 -0.9725  0.9725
  2大  0.9725 -0.9725
>

効果の方も、収率の平均と同じようにグラフにしてみます。

> ylist <- c(ee$tables$温度, ee$tables$流量, ee$tables$回転数, ee$tables$"流量:回転数")
> plot(0, 0, xlim = c(1, 8), ylim = c(min(ylist), max(ylist)), xaxt = "n", xlab = "", ylab="要 因 効 果")
> points(1:6, c(ee$tables$温度, ee$tables$流量, ee$tables$回転数), cex = 2, pch = 16)
> lines(1:2, ee$tables$温度)
> axis(1, at = 1:2, c(names(ee$tables$温度)))
> axis(1, at = 1.5, "温度", tick = F, padj = 1.5)
> lines(3:4, ee$tables$流量)
> axis(1, at = 3:4, c(names(ee$tables$流量)))
> axis(1, at = 3.5, "流量", tick = F, padj = 1.5)
> lines(5:6, ee$tables$回転数)
> axis(1, at = 5:6, c(names(ee$tables$回転数)))
> axis(1, at = 5.5, "回転数", tick = F, padj = 1.5)
> points(7:8, c(ee$tables$"流量:回転数"[1, ]), cex = 2, pch = 0)
> lines(7:8, ee$tables$"流量:回転数"[1, ], lty = 2)
> points(7:8, c(ee$tables$"流量:回転数"[2, ]), cex = 2, pch = 2)
> lines(7:8, ee$tables$"流量:回転数"[2, ], lty = 2)
> axis(1, at = 7:8, c(names(ee$tables$"流量:回転数"[1, ])))
> axis(1, at = 7.5, "（回転数）", tick = F, padj = 1.5)
> axis(1, at = 7.5, "流量×回転数", tick = F, padj = 3)
> text(7, ee$tables$"流量:回転数"[1, 1] * 1.1, labels = "（流量）1小")
> text(7, ee$tables$"流量:回転数"[2, 1] * 1.1, labels = "（流量）1大")
> abline(h = 0)
>

ここで各水準の "効果" についても "平均値" のときと同様に、データフレーム tbl を使って、いくつか検算をしておきましょう。

主効果については、各水準の平均から全平均を引いた値を効果としています。以下は、因子「流量」と「回転数」における、各水準の効果の算出例です。

> eff.flow_s <- mean(tbl[tbl$流量 == "1小", "収率"]) - gm
> eff.flow_s
[1] 0.95
> eff.flow_l <- mean(tbl[tbl$流量 == "2大", "収率"]) - gm
> eff.flow_l
[1] -0.95
> eff.rpm_l <- mean(tbl[tbl$回転数 == "1低", "収率"]) - gm
> eff.rpm_l
[1] 1.43
> eff.rpm_h <- mean(tbl[tbl$回転数 == "2高", "収率"]) - gm
> eff.rpm_h
[1] -1.43

交互効果については、各水準の平均から全平均を引いた値から、さらに、各主効果の該当する水準の効果を差し引いた値を、効果としています。因子「流量」と「回転数」の交互効果の算出例です。

> mean(tbl[(tbl$流量 == "1小") & (tbl$回転数 == "1低"), "収率"]) - gm - eff.flow_s - eff.rpm_l
[1] -0.9725
> mean(tbl[(tbl$流量 == "1小") & (tbl$回転数 == "2高"), "収率"]) - gm - eff.flow_s - eff.rpm_h
[1] 0.9725
> mean(tbl[(tbl$流量 == "2大") & (tbl$回転数 == "1低"), "収率"]) - gm - eff.flow_l - eff.rpm_l
[1] 0.9725
> mean(tbl[(tbl$流量 == "2大") & (tbl$回転数 == "2高"), "収率"]) - gm - eff.flow_l - eff.rpm_h
[1] -0.9725
> 

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

今回は二因子の実験です。この辺りから、いろいろと議論したいことが出てきますが、議論を展開するために記事の推敲を重ねていると、いつまで経っても公開できないので、とりあえず、サッと流して説明し、公開後に少しずつ書き足していくことにします。

二元配置分散分析

2 因子の効果を複数の水準で確認する実験には、二元配置分散分析 (two-way ANOVA) を利用することができます。2 因子の実験計画では、因子単独の効果（主効果）に加えて、 2 因子の組み合わせの効果（交互作用）が生じる場合を考慮する必要があります。

条件	触媒	温度	収率
1	A社品	1低	80.8
2	A社品	1低	79.0
3	A社品	1低	81.1
4	A社品	1低	85.0
5	A社品	2中	88.9
6	A社品	2中	81.8
7	A社品	2中	85.0
8	A社品	2中	88.8
9	A社品	3高	85.6
10	A社品	3高	88.5
11	A社品	3高	91.3
12	A社品	3高	81.7
13	B社品	1低	80.5
14	B社品	1低	80.2
15	B社品	1低	80.1
16	B社品	1低	79.2
17	B社品	2中	84.0
18	B社品	2中	81.1
19	B社品	2中	85.7
20	B社品	2中	82.9
21	B社品	3高	81.6
22	B社品	3高	82.0
23	B社品	3高	81.1
24	B社品	3高	85.0

例）

ある化学製品の合成工程において、触媒の種類と温度を 2 つの因子が、その製品の収量に影響を与えるかどうかを調べるため、触媒を 2 種類 (A社品, B社品) と（反応）温度を 3 水準 (低, 中, 高) に設定して実験を計画し、各条件の収量を測定した結果を右に示した。2 因子の収量に及ぼす影響を分析せよ。

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

> tbl <- read.table("clipboard", header = T)
> tbl
   条件  触媒 温度 収率
1     1 A社品  1低 80.8
2     2 A社品  1低 79.0
3     3 A社品  1低 81.1
4     4 A社品  1低 85.0
5     5 A社品  2中 88.9
6     6 A社品  2中 81.8
7     7 A社品  2中 85.0
8     8 A社品  2中 88.8
9     9 A社品  3高 85.6
10   10 A社品  3高 88.5
11   11 A社品  3高 91.3
12   12 A社品  3高 81.7
13   13 B社品  1低 80.5
14   14 B社品  1低 80.2
15   15 B社品  1低 80.1
16   16 B社品  1低 79.2
17   17 B社品  2中 84.0
18   18 B社品  2中 81.1
19   19 B社品  2中 85.7
20   20 B社品  2中 82.9
21   21 B社品  3高 81.6
22   22 B社品  3高 82.0
23   23 B社品  3高 81.1
24   24 B社品  3高 85.0
> 

2 つの因子（触媒と温度）で、因子間の交互作用を考慮しない場合には、因子の単純な線形結合の構造モデル 収率 ~ 触媒 + 温度 に対して分散分析をおこないます。

> result.aov <- aov(収率 ~ 触媒 + 温度, data = tbl)
> summary(result.aov)
            Df Sum Sq Mean Sq F value Pr(>F)   
触媒         1  48.45   48.45   7.217 0.0142 * 
温度         2  83.34   41.67   6.207 0.0080 **
Residuals   20 134.26    6.71                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> 

触媒の水準が有意水準 5% で、温度の水準が有意水準 1% で、収率に差を及ぼすと言えます。

因子間の交互作用を考慮する場合は、以下のように + の代わりに * を使います。

> result.aov2 <- aov(収率 ~ 触媒 * 温度, data = tbl)
> summary(result.aov2)
            Df Sum Sq Mean Sq F value Pr(>F)  
触媒         1  48.45   48.45   6.925 0.0169 *
温度         2  83.34   41.67   5.956 0.0104 *
触媒:温度    2   8.33    4.16   0.595 0.5621  
Residuals   18 125.94    7.00                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> 

触媒と温度との相互作用について有意性が認められませんでしたので、今回は、因子間の交互作用を考慮しない場合の分散分析の結果から解析を進めることにします。

model.tables() は、分散分析につかわれた統計モデルから、各因子、各水準の平均値および効果の大きさを算出します。

> mm <- model.tables(result.aov, type = "mean")
> mm
Tables of means
Grand mean
         
83.37083 

 触媒 
触媒
A社品 B社品 
84.79 81.95 

 温度 
温度
  1低   2中   3高 
80.74 84.78 84.60 
> ee <- model.tables(result.aov, type = "effects")
> ee
Tables of effects

 触媒 
触媒
  A社品   B社品 
 1.4208 -1.4208 

 温度 
温度
    1低     2中     3高 
-2.6333  1.4042  1.2292 

> 

効果 ee の値を用いて、いわゆる要因効果図に相当するプロットを作成してみましょう。

> plot(1:5, c(ee$tables$触媒, ee$tables$温度), xaxt = "n", xlab = "", ylab="Main Effects", cex = 2, pch = 16)
> lines(1:2, ee$tables$触媒)
> lines(3:5, ee$tables$温度)
> abline(h = 0)
> axis(1, at = 1:5, c(names(ee$tables$触媒), names(ee$tables$温度)))
> 

上記のプロットからは、高い収率を得るには、触媒に「A社品」を、温度「中」で、製造すれば良いと結論できそうなのですが、果たしてそれがベストと言えるでしょうか。

念の為、ボックスプロットも作成してみましょう。ggplot2 パッケージを利用すると、2 因子の場合でも簡単に綺麗なボックスプロットを描画できます。

> library(ggplot2)
> bp <- ggplot(tbl, aes(x = 温度, y = 収率, colour = 触媒)) + geom_boxplot()
> plot(bp)
>

上記ボックスプロットによると、確かに A社品の触媒の方が B社品のものより高い収率が得られていますが、見た目だけでも B社品の方が収率のばらつきの少ないことが判ります。

もしこの実験が A社品、B社品どちらかを採用することが目的だとすれば、結論を出すには慎重になります。特性（応答）の平均値の大小だけでなく、ばらつきが少ない条件を選ぶことも重要であるからです。2 因子の実験であれば、この例のように特性（応答）に関するいくつかのプロットを作成して検討することも可能でしょう。しかし、因子数が増えた場合に、合理的にばらつきを評価するにはどのようにすれば良いでしょうか。ここでは、品質工学の SN 比について紹介します。

品質工学の SN 比

品質工学 (Quality Engineering) とは、技術開発・新製品開発を効率的に行う開発技法です。考案者の田口玄一氏 (1924 - 2012) の名を冠してタグチメソッドとも呼ばれます。実験計画法では、実験の効率向上を目指しますが、さらに特性（応答）に SN 比の概念を導入して、特性のばらつきも評価します。これを品質工学では「パラメータ設計」と読んでいます。

ここでは、その中から望大特性に使われる SN 比を、収率（の平均値）に替る特性として評価します。望大特性とは読んで字の如く、（収率のように）値が大きければ大きいほど良い特性です。望大特性に使われる SN 比は以下のように定義されています。単位はデジベルです。

η は SN 比、y は収率、n は各条件の組み合わせの繰り返し回数で、この例の場合 4 になります。

R で SN 比を計算します。まず、tbl に「平方数逆数」という列を追加し、この列で収率の平方数の逆数を計算します。

> tbl$平方数逆数 <- 1 / (tbl$収率 * tbl$収率)
> tbl
   条件  触媒 温度 収率   平方数逆数
1     1 A社品  1低 80.8 0.0001531713
2     2 A社品  1低 79.0 0.0001602307
3     3 A社品  1低 81.1 0.0001520402
...
（以下省略）
...

各因子の水準組み合わせに対して、aggregate コマンドで「平方数逆数」の平均値 (mean) を集計します（データフレーム tbl2）。

> tbl2 <- aggregate(平方数逆数 ~ 触媒 + 温度, data = tbl, mean)
> tbl2
   触媒 温度   平方数逆数
1 A社品  1低 0.0001509626
2 B社品  1低 0.0001562673
3 A社品  2中 0.0001353011
4 B社品  2中 0.0001438573
5 A社品  3高 0.0001334833
6 B社品  3高 0.0001473380
>

tbl2 の「平方数逆数」から、望大特性の定義に従い SN 比を算出します（「SN比」列）。

> tbl2$SN比 <- -10 * log10(tbl2$平方数逆数)
> tbl2
   触媒 温度   平方数逆数     SN比
1 A社品  1低 0.0001509626 38.21131
2 B社品  1低 0.0001562673 38.06132
3 A社品  2中 0.0001353011 38.68699
4 B社品  2中 0.0001438573 38.42068
5 A社品  3高 0.0001334833 38.74573
6 B社品  3高 0.0001473380 38.31685
>

算出した SN 比に対して、分散分析をします。

> result.sn <- aov(SN比 ~ 触媒 + 温度, data = tbl2)
> summary(result.sn)
            Df  Sum Sq Mean Sq F value Pr(>F)  
触媒         1 0.11905 0.11905   12.13 0.0734 .
温度         2 0.22056 0.11028   11.24 0.0817 .
Residuals    2 0.01962 0.00981                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>

SN 比に対しては、どちらの因子も 5% 有意水準で帰無仮説を棄却できません。したがって、ばらつきまで考慮すると、A 社品と B 社品の触媒の優劣は付けられないことになります（まあ、有意水準 10% では帰無仮説を棄却できるのですがから、十中八九のレベルで良ければ、有意差があると言えなくもありません）。

SN 比に対しては有意差がある分析結果を得られませんでしたが、前と同じように特性要因図を作成してみましょう。

> mm2 <- model.tables(result.sn, type = "mean")
> mm2
Tables of means
Grand mean
         
38.40715 

 触媒 
触媒
A社品 B社品 
38.55 38.27 

 温度 
温度
  1低   2中   3高 
38.14 38.55 38.53 
> plot(1:5, c(mm2$tables$触媒, mm2$tables$温度), xaxt = "n", xlab = "", ylab="SN Ratio [dB]", cex = 2, pch = 16)
> lines(1:2, mm2$tables$触媒)
> lines(3:5, mm2$tables$温度)
> axis(1, at = 1:5, c(names(mm2$tables$触媒), names(mm2$tables$温度)))
> 

特性要因図のプロットを比較する限り、SN 比が特性を加工した値だからといって、（ばらつきの大きさに依存しますが）元の特性の平均値と較べて予想外の大きな傾向変化があるわけではありません。

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

今回は一因子の実験で、水準が 3 以上の場合について、分散分析という手法を導入します。分散分析は、実験計画法において切っても切れない関係にあります。ここでは、概念的なことと R による計算例を紹介するに留め、分散分析の原理については詳しく立ち入っていません。

分散分析の必要性

1 因子の 2 水準の平均の違いを調べる方法が t 検定でした。それでは 3 水準（例えば A，B，C）になったらどうしたらよいでしょう。A と B、B と C、そして C と A のペアでそれぞれ t 検定を行ない、いずれかで帰無仮説が棄却されたならば、3 つの水準の平均は等しくない、と結論づけることができます。しかし、このやり方は水準 n の場合、 n * (n - 1) / 2 回の検定をしなければならず、効率的ではありません。

水準が 3 以上の場合には、因子の各水準の母平均に違いがあるかどうかを「分散」の大きさの違いで検定を行ないます。これを「分散分析 (Analysis of Variance, ANOVA)」といいます。複数の水準の平均が互いに等しいという帰無仮設を検定するもので，水準の間で少なくとも 1 つの平均が他と異なっているかどうかをこの分析では調べることになります。

一元配置分散分析

1 因子 1 つの効果を確認しようとする実験において、3 つ以上の水準で見る場合、上述のように分散分析を用いますが、これを一元配置分散分析 (one-way ANOVA) と呼びます。

#	A社品	B社品	C社品
1	75.5	74.3	76.3
2	75.1	74.2	76.4
3	79.5	82.0	75.8
4	84.0	76.2	82.3
5	83.1	74.2	81.1
6	79.0	78.5	81.5
7	82.6	84.2	76.5
8	80.0	80.7	79.2
9	77.9	77.7	75.0
10	81.1	74.3	77.5
11	84.6	77.5	81.7
12	81.2	79.6	75.7
13	78.0	82.4	82.0
14	79.8	73.9	75.4
15	79.0	73.0	82.0
16	84.4	74.5	75.3
17	79.0	79.1	80.9
18	80.3	80.6	78.3
19	79.6	83.5	80.2
20	81.3	78.6	75.8
21	83.4	82.8	79.4
22	77.5	78.5	79.7
23	78.9	79.7	78.8
24	77.3	83.1	77.9
25	83.5	75.6	77.4

例）

ある製品の原料 W について、仕様が同じ 3 社の「A 社品」、「B 社品」、「C 社品」それぞれを使用した場合における収率の違いを評価したい。それぞれの原料を用いた 1 ロット、25 サンプルの収率の結果（右表）を用いて有意水準 5% で検定せよ。

帰無仮説 H₀ は、『原料 W の水準「A 社品」、「B 社品」、「C 社品」では、収率に差は無い』になります。

まずは、右表のデータからデータフレームを作成します。

> a <- c(75.5, 75.1, 79.5, 84.0, 83.1, 79.0, 82.6, 80.0, 77.9, 81.1, 84.6, 81.2, 78.0, 79.8, 79.0, 84.4, 79.0, 80.3, 79.6, 81.3, 83.4, 77.5, 78.9, 77.3, 83.5)
> b <- c(74.3, 74.2, 82.0, 76.2, 74.2, 78.5, 84.2, 80.7, 77.7, 74.3, 77.5, 79.6, 82.4, 73.9, 73.0, 74.5, 79.1, 80.6, 83.5, 78.6, 82.8, 78.5, 79.7, 83.1, 75.6)
> c <- c(76.3, 76.4, 75.8, 82.3, 81.1, 81.5, 76.5, 79.2, 75.0, 77.5, 81.7, 75.7, 82.0, 75.4, 82.0, 75.3, 80.9, 78.3, 80.2, 75.8, 79.4, 79.7, 78.8, 77.9, 77.4)
> y <- c(rep("A", 25), rep("B", 25), rep("C", 25))
> x <- c(rep("A", 25), rep("B", 25), rep("C", 25))
> tbl <- data.frame(原料W = x, 収率 = c(a, b, c))
> tbl
   原料W 収率
1     A  75.5
2     A  75.1
3     A  79.5
...
（省略）
...
73    C  78.8
74    C  77.9
75    C  77.4
> 

R の oneway.test() 関数を利用して一元配置分散分析をしてみましょう。このコマンドも t.test() 同様 、デフォルトではサンプルの等分散を仮定していない Welch の検定をするようになっていますので、var.euqal = T を指定して、サンプルが等分散であることを仮定します。

> oneway.test(収率 ~ 原料W, data = tbl, var.equal = TRUE)

 One-way analysis of means

data:  収率 and 原料W
F = 3.2351, num df = 2, denom df = 72, p-value = 0.04515

有意水準 α = 0.05 の検定において、P-値 0.04515 は 0.05 よりも小さいので有意となり、材料W が収率に与える影響には差があると判断されます。

1 因子 3 水準実験

上の例を『ある製品の原料 W について、仕様が同じ 3 社の「A 社品」、「B 社品」、「C 社品」それぞれを使用した場合における収率の違いを評価するため、原料 W を因子とした 1 因子 3 水準の実験を計画し、右表の結果を得た。』と読みかえれば、（結果は変わらないものの）実験を計画／実施し、実験の解析をする問題に置き換わります。

分散分析モデルのあてはめをする aov() 関数（コマンド）と関連するコマンドを使って統計量を計算すると以下のようになります。

> result.aov <- aov(収率 ~ 原料W, data = tbl)
> coefficients(result.aov)
(Intercept)      原料WB      原料WC 
     80.224      -1.876      -1.740
> summary(result.aov)
            Df Sum Sq Mean Sq F value Pr(>F)  
原料W        2   54.7  27.356   3.235 0.0451 *
Residuals   72  608.8   8.456                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> model.tables(result.aov)
Tables of effects

 原料W 
原料W
      A       B       C 
 1.2053 -0.6707 -0.5347
> residuals(result.aov)
     1      2      3      4      5      6      7      8      9     10     11     12     13     14     15     16     17     18     19     20     21     22     23     24     25     26     27     28     29     30     31     32     33 
-4.724 -5.124 -0.724  3.776  2.876 -1.224  2.376 -0.224 -2.324  0.876  4.376  0.976 -2.224 -0.424 -1.224  4.176 -1.224  0.076 -0.624  1.076  3.176 -2.724 -1.324 -2.924  3.276 -4.048 -4.148  3.652 -2.148 -4.148  0.152  5.852  2.352 
    34     35     36     37     38     39     40     41     42     43     44     45     46     47     48     49     50     51     52     53     54     55     56     57     58     59     60     61     62     63     64     65     66 
-0.648 -4.048 -0.848  1.252  4.052 -4.448 -5.348 -3.848  0.752  2.252  5.152  0.252  4.452  0.152  1.352  4.752 -2.748 -2.184 -2.084 -2.684  3.816  2.616  3.016 -1.984  0.716 -3.484 -0.984  3.216 -2.784  3.516 -3.084  3.516 -3.184 
    67     68     69     70     71     72     73     74     75 
 2.416 -0.184  1.716 -2.684  0.916  1.216  0.316 -0.584 -1.084 
>

分散分析では、水準の間で少なくとも 1 つの平均が他と異なっているかどうかを調べているにも関わらず、各水準の平均値が分析結果に出てこないので、tapply() 関数（コマンド）で、平均、分散、標準偏差を集計してみました。

> tapply(tbl$収率, tbl$原料W, mean)
     A      B      C 
80.224 78.348 78.484
> tapply(tbl$収率, tbl$原料W, var)
      A       B       C 
 7.1144 12.1376  6.1164 
> tapply(tbl$収率, tbl$原料W, sd)
       A        B        C 
2.667283 3.483906 2.473136 
>

この手の分析では、データのばらつきは正規分布に従っていることを前提としていますが、現実のデータでは外れ値 (outlier) が平均値を引っ張っていることがありますので、ボックスプロット（箱ひげ図）で確認します。ボックスプロットは順序統計量を扱っています。

なお、この説明のために準備しているサンプルデータは一様乱数を使って生成していますので、正規分布に従わないばかりか、（十分なサンプル数があれば）外れ値もありません。

> bp <- boxplot(収率 ~ 原料W, data = tbl)
> bp
$stats
     [,1] [,2] [,3]
[1,] 75.1 73.0 75.0
[2,] 78.9 74.5 76.3
[3,] 79.8 78.5 78.3
[4,] 82.6 80.7 80.9
[5,] 84.6 84.2 82.3

$n
[1] 25 25 25

$conf
        [,1]    [,2]    [,3]
[1,] 78.6308 76.5408 76.8464
[2,] 80.9692 80.4592 79.7536

$out
numeric(0)

$group
numeric(0)

$names
[1] "A" "B" "C"

> 

参考文献

1. 岩崎学 統計的データ解析入門 実験計画法 東京図書 (2006)

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

DOE（実験計画法）とは因果関係を調べる手法（学問）である、と言葉で表現するのは簡単ですが、これだけでは、具体的にどのような手法であるのか判然としません。実験計画法のコアな話題を展開する前に、簡単な例を用いて実験計画法を紹介していきます。

実験計画法におけるデータ解析は、統計学の手法がベースとなっています。そこで、出発点として t 検定の例から実験計画法の世界に入っていくことにしましょう。なお、ここでは教科書的な数式による説明を極力省き、R による計算にフォーカスしています。

t 検定

t 検定は、帰無仮説が正しいと仮定した場合に、統計量が t 分布に従うことを利用する統計学的検定法の総称です。二組の標本について平均に有意差があるかどうかの検定などに用いられます。

原料	収率
現行品	87.3
現行品	90.0
現行品	86.7
現行品	90.7
現行品	91.3
現行品	85.8
現行品	84.3
現行品	89.7
現行品	90.4
現行品	86.8
現行品	86.9
現行品	88.3
現行品	88.6
現行品	88.5
現行品	84.2
現行品	84.2
現行品	82.4
現行品	84.6
現行品	88.4
現行品	84.3
現行品	83.5
現行品	84.1
現行品	83.9
廉価品	87.9
廉価品	84.6
廉価品	84.4
廉価品	86.9
廉価品	82.3
廉価品	83.8
廉価品	88.6
廉価品	84.6
廉価品	85.4
廉価品	80.6
廉価品	85.5
廉価品	84.6
廉価品	90.0
廉価品	83.6
廉価品	87.1
廉価品	83.3
廉価品	85.9
廉価品	82.1
廉価品	86.9
廉価品	85.3

例）

ある製品の主原料を、「現行品」から、同じ仕様で仕入れ価格が単位量当たり 5% 安い「廉価品」に切り替えたい。「廉価品」を使用した 20 ロットの試作品の収率（単位 %）と、直近に製造した「現行品」を使用した製品 23 ロットの収率（右表）を、有意水準 5 % で検定せよ。

帰無仮説 H₀ は、『主原料「現行品」と「廉価品」では、収率に差は無い』になります。

まずは、サンプルのデータフレームを作成します。

> x <- c(rep("現行品", 23), rep("廉価品", 20))
> x
 [1] "現行品" "現行品" "現行品" "現行品" "現行品" "現行品" "現行品" "現行品"
 [9] "現行品" "現行品" "現行品" "現行品" "現行品" "現行品" "現行品" "現行品"
[17] "現行品" "現行品" "現行品" "現行品" "現行品" "現行品" "現行品" "廉価品"
[25] "廉価品" "廉価品" "廉価品" "廉価品" "廉価品" "廉価品" "廉価品" "廉価品"
[33] "廉価品" "廉価品" "廉価品" "廉価品" "廉価品" "廉価品" "廉価品" "廉価品"
[41] "廉価品" "廉価品" "廉価品"
> y <- c(87.3, 90.0, 86.7, 90.7, 91.3, 85.8, 84.3, 89.7, 90.4, 86.8, 86.9, 88.3, 88.6, 88.5, 84.2, 84.2, 82.4, 84.6, 88.4, 84.3, 83.5, 84.1, 83.9,87.9, 84.6, 84.4, 86.9, 82.3, 83.8, 88.6, 84.6,85.4, 80.6, 85.5, 84.6, 90.0, 83.6, 87.1, 83.3, 85.9, 82.1, 86.9, 85.3)
> y
 [1] 87.3 90.0 86.7 90.7 91.3 85.8 84.3 89.7 90.4 86.8 86.9 88.3 88.6 88.5 84.2
[16] 84.2 82.4 84.6 88.4 84.3 83.5 84.1 83.9 87.9 84.6 84.4 86.9 82.3 83.8 88.6
[31] 84.6 85.4 80.6 85.5 84.6 90.0 83.6 87.1 83.3 85.9 82.1 86.9 85.3
> tbl <- data.frame(原料 = x, 収率 = y)
> tbl
     原料 収率
1  現行品 87.3
2  現行品 90.0
3  現行品 86.7
...
（省略）
...
41 廉価品 82.1
42 廉価品 86.9
43 廉価品 85.3

R の t.test() 関数（コマンド）を利用して t 検定をしてみましょう。このコマンドは、デフォルトではサンプルの等分散を仮定していない Welch の検定をするようになっていますので、var.euqal = T を指定して、2 つのサンプルが等分散であることを仮定しています。

> t.test(収率 ~ 原料, data = tbl, var.equal = T)

 Two Sample t-test

data:  収率 by 原料
t = 2.0445, df = 41, p-value = 0.04736
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.01913643 3.11042879
sample estimates:
mean in group 現行品 mean in group 廉価品 
            86.73478             85.17000 

>

この例では検定の結果、p-value（p値）が有意水準 5 % (= 0.05) より小さい値になりますので帰無仮説を棄却し、「現行品」と「廉価品」の収率の差は有意であり、「廉価品」を用いた製品の収率は「現行品」より 1.6% 程度低いという結論になります。

1 因子 2 水準実験

上述の例について、主原料「現行品」と「廉価品」の、製品収率に及ぼす影響を調べた「実験」であると捉えることもできます。

実験の主たる目的は、原因と結果の因果関係を明らかにすることです。

因果関係の原因系は、要因あるいは因子 (factor) と呼びます。先の例では、因子は一つで「原料」です。因子の設定条件を水準 (level) と呼びます。因子「原料」の水準は「現行品」と「廉価品」です。因果関係の結果系は、特性値あるいは応答 (response) と呼びます。先の例では、特性値は「収率」になります。先の例を実験として計画 (design) したものであれば、これを「1 因子 2 水準実験」と呼びます。

R で結果だけを手っ取り早く知りたい場合は下記のようにします。aov() 関数は分散分析のために統計モデルのあてはめ（フィッティング）を行うコマンドです。詳細は次の分散分析で説明します。

因子「原料」は 5% の有意水準で有意に応答「収率」に影響を及ぼしており、それぞれの因子の水準のインパクトはボックスプロットの通りであると評価できます。

> summary(aov(収率 ~ 原料, data = tbl))
            Df Sum Sq Mean Sq F value Pr(>F)  
原料         1  26.19  26.194    4.18 0.0474 *
Residuals   41 256.91   6.266                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> bp <- boxplot(収率 ~ 原料, data = tbl)
> bp
$stats
      [,1]  [,2]
[1,] 82.40 80.60
[2,] 84.25 83.70
[3,] 86.80 84.95
[4,] 88.55 86.90
[5,] 91.30 90.00

$n
[1] 23 20

$conf
         [,1]     [,2]
[1,] 85.38335 83.81944
[2,] 88.21665 86.08056

$out
numeric(0)

$group
numeric(0)

$names
[1] "現行品" "廉価品"

>

現実の実験では、5% 安い主原料の価格と、収率が 1.6% 程度低下することを天秤にかけ、結果として「廉価品」を採用する方向で検討するのであれば、次のステップとして製品の品質や信頼性の評価へ進むことになるでしょう。

参考文献

1. 岩崎学 統計的データ解析入門 実験計画法 東京図書 (2006)

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

R の拡張パッケージ AlgDesign を紹介した際に、「最適計画を RSM（応答曲面法）に活用した、エンジニアリングのための実験計画/解析ツールを作りたい、それが自分がコンピュータのプログラミングに取り組む最も基本的なモチベーションになっているからです。」と書きましたが、実験計画法は、誰でも知っているようなエンジニアリング・ツールとして市民権を得ているほどポピュラーな手法だとは考えていません。そこで、自分の知識を整理するためにも、実験計画法とは何かを何回かに分けて紹介していきたいと思います。

はじめに

まずは、実験計画法についてインターネット上でどのように説明がされているか、いくつか例を紹介します。

Wikipedia では、実験計画法について、次のような説明がされています^1.。

実験計画法（じっけんけいかくほう、Experimental design、Design of experiments）は、効率のよい実験方法をデザインし、結果を適切に解析することを目的とする統計学の応用分野である。R.A.フィッシャーが 1920 年代に農学試験から着想して発展させた。特に 1950 年 G.M.コックスと W.G.コクランが標準的教科書を出版し、以後医学、工学、実験心理学や社会調査へ広く応用された。またこれを基にして田口玄一による品質工学という新たな分野も生まれた。

上記は、「実験計画法」を説明した冒頭部分で、実験計画法が発展した歴史的な経緯を端的に説明しています。

Tech On！のものづくり用語では、実験計画法について次のような説明がされています^2.。

寸法や形状，材質，部品の固定方法など，設計パラメータとして検討したい因子が，性能や機能などの評価項目にどのような影響を与えるのかを調べる場合，いい加減な方法で実験しても正しい傾向を知ることはできない。精度良い結果を効率的に（安く，早く）得られるような実験を設計し，その実験で得られた結果を解析して結論を出すことを実験計画法と呼ぶ。実験を設計する部分だけを指す場合は，実験計画となる。

JMP では、実験計画法について次のように説明がされています^3.。

実験計画（DOE）は、複数の因子からなる実験空間を探索するための、実用的かつ普遍的なアプローチです。JMPには、使いやすい、世界標準の実験計画（DOE）の計画と分析の機能が搭載されています。

実験計画法は、効率的かつ効果的な情報収集が必要なさまざまな場面で応用できます。入力（因子）と出力（応答）の関係を明らかにし、モデル化する上で、最も優れたアプローチは、因子を意図的に変化させたときに、応答も変化するかどうかを観察することです。あらかじめ用意した計画に従って、因子を意図的に変化させることが、有意義な新しい発見につながります。

しかし、実際の世界では、複数の因子が存在することがほとんどです。このため、因子を1つずつ変更する方法は、効率的とは言えません。複数の因子が応答に与える影響を正しく把握するためには、実験計画（DOE）が必要になります。

JMP は SAS 社が開発・販売する統計解析ソフトウェアです。

実験計画法は、Wikipedia の説明にあるように Ｒ.Ａ.フィッシャーが農事試験でその基礎を築いたもので、複数の実験因子が結果に及ぼす影響を合理的に分離して評価できるように、言い換えれば、因果関係が判るように、予め、実験因子の組み合わせを選ぶ（デザインする）手法です。農事試験において、例えば収率がテーマであれば、収穫して初めて結果がわかります。このように時間がかかる「実験」において、実験因子の影響を合理的に評価する手法を確立することには、時間的、費用的経済性を考えると、実に切実なるニーズがあったことでしょう。

1. 実験計画法 - Wikipedia
2. 実験計画法とは - 設計・生産とIT - Tech-On！
3. JMPによる実験計画（DOE）

Rで実験計画法 前編 from itoyan110

Rで実験計画法 後編 from itoyan110