今回は 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 >
0 件のコメント:
コメントを投稿