2014-03-09

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

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


3 因子 2 水準の実験(繰り返しなし)

条件温度流量回転数収率
11低1小1低81.08
21低1小2高80.66
31低2大1低80.73
41低2大2高76.14
52高1小1低83.39
62高1小2高81.98
72高2大1低83.83
82高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
> 

0 件のコメント: