2014-02-28

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

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


二元配置分散分析

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

条件触媒温度収率
1A社品1低80.8
2A社品1低79.0
3A社品1低81.1
4A社品1低85.0
5A社品2中88.9
6A社品2中81.8
7A社品2中85.0
8A社品2中88.8
9A社品3高85.6
10A社品3高88.5
11A社品3高91.3
12A社品3高81.7
13B社品1低80.5
14B社品1低80.2
15B社品1低80.1
16B社品1低79.2
17B社品2中84.0
18B社品2中81.1
19B社品2中85.7
20B社品2中82.9
21B社品3高81.6
22B社品3高82.0
23B社品3高81.1
24B社品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 比が特性を加工した値だからといって、(ばらつきの大きさに依存しますが)元の特性の平均値と較べて予想外の大きな傾向変化があるわけではありません。

2014-02-15

実験計画法 (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社品
175.574.376.3
275.174.276.4
379.582.075.8
484.076.282.3
583.174.281.1
679.078.581.5
782.684.276.5
880.080.779.2
977.977.775.0
1081.174.377.5
1184.677.581.7
1281.279.675.7
1378.082.482.0
1479.873.975.4
1579.073.082.0
1684.474.575.3
1779.079.180.9
1880.380.678.3
1979.683.580.2
2081.378.675.8
2183.482.879.4
2277.578.579.7
2378.979.778.8
2477.383.177.9
2583.575.677.4

例)

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

帰無仮説 H0 は、『原料 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)

2014-02-11

実験計画法 (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 % で検定せよ。

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

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

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