今回は二因子の実験です。この辺りから、いろいろと議論したいことが出てきますが、議論を展開するために記事の推敲を重ねていると、いつまで経っても公開できないので、とりあえず、サッと流して説明し、公開後に少しずつ書き足していくことにします。
二元配置分散分析
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 比が特性を加工した値だからといって、（ばらつきの大きさに依存しますが）元の特性の平均値と較べて予想外の大きな傾向変化があるわけではありません。
