実験計画法 多因子-２水準-L16直交表 に関するメモ （Design of experiments & Taguchi method）

まだまだ、実験計画法について復習したいと思います。
ようやく３因子以上の多因子について学びたいと思います。
多因子で水準数が多い場合は、単純には実験数が多くなってしまう問題があり現実的でない場合が多々あります。

実験数を減らしながら、因子の影響をしっかり把握するために、ようやく直交表の出番となります。

それではさっそく例題を使って直交表を使った分散分析と区間推定を実践していきたいと思います。

今回は各因子２水準という制約のある場合について例題を取り上げたいと思います。

■５因子、５因子交互作用の例を実践

以下のような５因子と５因子交互作用を考えたい場合の分散分析を考えたいと思います。

因子	列	成分
A	1	a
B	2	b
C	4	c
D	8	d
F	14	bcd
G	15	abcd
A✕B	3	ab
A✕C	5	ac
A✕D	9	ad
B✕C	6	bc
B✕D	10	bd
C✕D	12	cd
誤差	7	abc
	11	abd
	13	acd

因子には、Eの記号は使っていません（Eは誤差と混乱するため）。

上記を２水準のL16の直交表にあてはめて、出力を得た結果が以下になります。

成分	A	B	A✕B	C	A✕C	B✕C		D	A✕D	B✕D		C✕D		F	G
	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15	出力y
1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	23
2	1	1	1	1	1	1	1	2	2	2	2	2	2	2	2	16
3	1	1	1	2	2	2	2	1	1	1	1	2	2	2	2	26
4	1	1	1	2	2	2	2	2	2	2	2	1	1	1	1	23
5	1	2	2	1	1	2	2	1	1	2	2	1	1	2	2	21
6	1	2	2	1	1	2	2	2	2	1	1	2	2	1	1	18
7	1	2	2	2	2	1	1	1	1	2	2	2	2	1	1	29
8	1	2	2	2	2	1	1	2	2	1	1	1	1	2	2	28
9	2	1	2	1	2	1	2	1	2	1	2	1	2	1	2	28
10	2	1	2	1	2	1	2	2	1	2	1	2	1	2	1	22
11	2	1	2	2	1	2	1	1	2	1	2	2	1	2	1	22
12	2	1	2	2	1	2	1	2	1	2	1	1	2	1	1	23
13	2	2	1	1	2	2	1	1	2	2	1	1	2	2	2	38
14	2	2	1	1	2	2	1	2	1	1	2	2	1	1	2	41
15	2	2	1	2	1	1	2	1	2	2	1	2	1	1	2	38
16	2	2	1	2	1	1	2	2	1	1	2	1	2	2	1	25
成分	a		a		a		a		a		a		a		a
		b	b			b	b			b	b			b	b
				c	c	c	c					c	c	c	c
								d	d	d	d	d	d	d	d

16回の試験を行うことで、5つの因子と、５つの交互作用に関する分散分析表を作成します。

■分散分析表の計算について

分散分析表の計算をまずはしていきましょう。
2水準であることから、各要素の平方和は以下のように計算できます。

■Scilabでの分散分析表の計算

実際にscilabで計算してみました。上記の式からも想像できるように、計算が単純です。

MATLABだfracfactという一部実施要因計画法の関数があるようです。

scilabにはそういった便利な関数はないので、自分で定義しています。

（また、あえて別の統計用のソフトでなくscilabを使っているのは、計算過程を確認しやすく、自習する目的のひとつです。）

//実験データ入力[データ]
Y = [
23
16
26
23
21
18
29
28
28
22
22
23
38
41
38
25
]

N=16//直交表の数

dfF = [
1    1    1    1    1    1    1    1    1    1    1    1    1    1    1
1    1    1    1    1    1    1    -1    -1    -1    -1    -1    -1    -1    -1
1    1    1    -1    -1    -1    -1    1    1    1    1    -1    -1    -1    -1
1    1    1    -1    -1    -1    -1    -1    -1    -1    -1    1    1    1    1
1    -1    -1    1    1    -1    -1    1    1    -1    -1    1    1    -1    -1
1    -1    -1    1    1    -1    -1    -1    -1    1    1    -1    -1    1    1
1    -1    -1    -1    -1    1    1    1    1    -1    -1    -1    -1    1    1
1    -1    -1    -1    -1    1    1    -1    -1    1    1    1    1    -1    -1
-1    1    -1    1    -1    1    -1    1    -1    1    -1    1    -1    1    -1
-1    1    -1    1    -1    1    -1    -1    1    -1    1    -1    1    -1    1
-1    1    -1    -1    1    -1    1    1    -1    1    -1    -1    1    -1    1
-1    1    -1    -1    1    -1    1    -1    1    -1    1    1    -1    1    -1
-1    -1    1    1    -1    -1    1    1    -1    -1    1    1    -1    -1    1
-1    -1    1    1    -1    -1    1    -1    1    1    -1    -1    1    1    -1
-1    -1    1    -1    1    1    -1    1    -1    -1    1    -1    1    1    -1
-1    -1    1    -1    1    1    -1    -1    1    1    -1    1    -1    -1    1
]
//直交表

//平方和の計算
S = dfF'*Y/(N/2)
S = S.*S*(N/4)

上記の平方和から、分散分析表を作成すると以下のようになります。

因子	S	自由度	V	F	p
A	175.5625	1	175.56	11.92	0.041
B	189.0625	1	189.06	12.84	0.037
C	3.0625	1	3.06	0.21	0.679
D	52.5625	1	52.56	3.57	0.155
F	39.0625	1	39.06	2.65	0.202
G	27.5625	1	27.56	1.87	0.265
A✕B	95.0625	1	95.06	6.45	0.085
A✕C	150.0625	1	150.06	10.19	0.050
A✕D	0.0625	1	0.06	0.00	0.952
B✕C	0.5625	1	0.56	0.04	0.858
B✕D	0.0625	1	0.06	0.00	0.952
C✕D	0.5625	1	0.56	0.04	0.858
誤差	22.5625	3	14.73
	7.5625
	14.0625

オレンジ色のセルにあるように、A✕D,B✕C,B✕D,C✕Dに関しては交互作用がまったくないので、誤差の項に入れるのが望ましいと考えるのが一般的です。
いわゆるプーリングの作業です。
Cも省きたいところですが、A✕Cに効果がみられるので、残しておきます。

プーリング作業をした結果が以下になります。

因子	S	自由度	V	F	p
A	175.5625	1	175.56	27.05	0.014
B	189.0625	1	189.06	29.13	0.012
C	3.0625	1	3.06	0.47	0.541
D	52.5625	1	52.56	8.10	0.065
F	39.0625	1	39.06	6.02	0.091
G	27.5625	1	27.56	4.25	0.131
A✕B	95.0625	1	95.06	14.65	0.031
A✕C	150.0625	1	150.06	23.12	0.017
誤差	0.0625	7	6.49
	0.5625
	0.0625
	0.5625
	22.5625
	7.5625
	14.0625

上記の結果からA✕B,A✕Cの2つの交互作用があることが分かりました。

■交互作用が一部でもない場合の区間推定

続いて、各組み合わせにおけるt分布による区間推定を実施しましょう。

今回は、前回の2因子の場合とは違い、交互作用が存在しないものがあるので、t分布の区間推定の方法の算出が少し変わります。

交互作用がある場合は、その組み合わせだけから平均値を求めていましたが、交互作用がないのであれば、それは実験数が少ないため効率的でないといった問題が生じるので、以下のように計算していきます。

まず、各因子の効果をα〜ηで表すとし、２つの交互因子を考えると以下のようにモデルを表現することができます。

各因子の要因効果は、以下のように個別に計算していきます。

　　　　・

　　　　・

　　　　・

　

上記の結果から、μは以下のように表すことができます。

t分布を使った区間推定は以下の式より求めることができます。

ここで、n_eは有効反復数（有効繰り返し数）と呼ばれます。
田口の式よりn_eは以下のように計算できます。

今回は、n_e = 16/9となります。

以上の式を参考に区間推定ができます。
うまくscilabで計算する方法が思いつかなかったので、単純にスプレッドシートで各水準の平均を計算した結果が以下です。

出力y	All	26.3125
	A1B1	22
	A1B2	24
	A2B1	23.75
	A2B2	35.5
	A1C1	19.5
	A1C2	26.5
	A2C1	32.25
	A2C2	27
	A1	23
	A2	29.625
	D1	28.125
	D2	24.5
	F1	27.875
	F2	24.75
	G1	25
	G2	27.625

これを上記の式にそれぞれ64通りあてはめて計算します。
ちなみに、64通りは以下のような感じに並べられます。

No.	A	B	C	D	F	G	AB	AC
1	1	1	1	1	1	1	11	11
2	1	1	1	1	1	2	11	11
3	1	1	1	1	2	1	11	11
4	1	1	1	1	2	2	11	11
5	1	1	1	2	1	1	11	11
6	1	1	1	2	1	2	11	11
7	1	1	1	2	2	1	11	11
8	1	1	1	2	2	2	11	11
9	1	1	2	1	1	1	11	12
10	1	1	2	1	1	2	11	12
11	1	1	2	1	2	1	11	12
12	1	1	2	1	2	2	11	12
・
・
・
62	2	2	2	2	1	2	22	22
63	2	2	2	2	2	1	22	22
64	2	2	2	2	2	2	22	22

上記を計算し、例えば出力yを最大にするものを選択するとすると、A=2,B=2,C=1,D=1,F=1,G=2の時が選択できます。
区間推定も実施すると以下のような結果になります。

平均	42.8125
上側信頼区間	38.2941
下側信頼区間	47.3309

■Rで解く

追加でRで分散分析を実施します。
データを読み込み分散分析を実施するのはすぐできるので、Rを使うのを激しくオススメします。
ただし、そこから区間推定などを各因子で実施するのは、スプレッドシートなどを用いた方がやはり良さそうな印象を持ちます。

#作業ディレクトリの設定
setwd("作業ディレクトリを記入")

#データの読み出し
d<-read.table("data.csv",header = TRUE,sep = ",")
str(d)#読み出しデータ構造の確認

#モデルあてはめ
d.lm <- lm(y~A+B+C+D+F+G+A:B+A:C+A:D+B:C+B:D+C:D,data = d)
summary(d.lm)
anova(d.lm)

#モデルあてはめ(Pooling)
d.pool.lm <- lm(y~A+B+C+D+F+G+A:B+A:C,data = d)
summary(d.pool.lm)
anova(d.pool.lm)

#グラフの描画（推定モデルの確認）
par(mfrow=c(2,2))
plot(d.pool.lm)

本ページは下の書籍を参考に自習用にまとめています。



他の参考書籍です



関連ページ
２因子ー質的因子の場合
１因子ー量的因子ー二次モデルの場合
１因子ー量的因子ー一次モデルの場合
１因子ー質的因子の場合