実験計画法 ２因子-質的因子 に関するメモ （Design of experiments & Taguchi method）

今回はやっと交互作用が議論される２因子の質的因子について、考えたいと思います。

今までは、１因子の質的因子、１因子の量的因子の１次モデル、１因子の量的因子の２次モデルについて考えてきました。

さて、ここからは２因子の質的因子についてまとめます。

■２因子の質的因子について

今までとは異なり、２つ目の因子（因子Bとします）が登場します。
そのため、平方和や分散もBに関する計算が出てきます。
それだけでなく、交互作用に関する平方和と分散の計算も登場してくるので注意が必要です。

因子Aを第i水準（トータル数a）、因子Bを第j水準（トータル数b）、繰り返しをk（トータル数n）とおくと、各平方和は以下のように書くことができます。

交互作用をしめすS_(AxB)は以下のように求めることができます

■２因子の質的因子について

それでは、実際にscilabで解いてみたいと思います

//実験データ入力[水準番号(A,B)　データ]
X = [
1    1    7.56
1    1    7.99
1    2    8.33
1    2    8.57
1    3    9.09
1    3    8.06
2    1    8.49
2    1    9.91
2    2    9.74
2    2    9.92
2    3    9.15
2    3    8.56
3    1    9.26
3    1    10.69
3    2    9.88
3    2    10.64
3    3    7.53
3    3    7.79
4    1    8.29
4    1    9.71
4    2    9.33
4    2    9.7
4    3    6.59
4    3    6.33
]

n=2;//１因子の実験繰り返し数
a=4;//A因子の水準数
b=3;//B因子の水準数
phiA = a-1;//自由度
phiB = b-1;//自由度
phiAxB = (a-1)*(b-1)//自由度
phiE = a*b*(n-1)//自由度
phiT = a*b*n-1//自由度
Y = X(:,3);//出力Yだけ分離

//分散分析表の計算
//下準備①
J = ones(a*b*n,a*b*n)*(1/(a*b*n));//全体平均を出すために利用(Jn)

//下準備②
A = eye(a,a);
B = (1/(b*n))*ones(b*n,b*n);
LA = kron(A,B);//クロネッカー テンソル積　各A水準での平均算出に利用
YAbar=LA*Y;

A = eye(b,b);
B = (1/(a*n))*ones(n,n);
LB = kron(A,B);//クロネッカー テンソル積　各B水準での平均算出に利用
LB = repmat(LB,a);
YBbar=LB*Y;

A = eye(a*b,a*b);
B = (1/(n))*ones(n,n);
LAB = kron(A,B);//クロネッカー テンソル積　各AxB水準での平均算出に利用
YABbar=LAB*Y;

////分散分析表の実際の計算（結果その①）
Ybar = J*Y;
SA = (YAbar - Ybar)'*(YAbar - Ybar)//各因子平均-全平均
SB = (YBbar - Ybar)'*(YBbar - Ybar)//各因子平均-全平均
SAB = (YABbar - Ybar)'*(YABbar - Ybar);//各AB因子平均-全平均
SAxB = SAB - SA - SB//SAｘBの計算
SE = (Y-YABbar)'*(Y-YABbar)//各データー各因子平均

VA = SA/phiA
VB = SB/phiB
VAxB = SAxB/phiAxB
VE = SE/phiE

fA = VA/VE
fB = VB/VE
fAxB = VAxB/VE

pA = 1 - cdff("PQ",fA,phiA,phiE)
pB = 1 - cdff("PQ",fB,phiB,phiE)
pAxB = 1 - cdff("PQ",fAxB,phiAxB,phiE)

//区間推定（結果その②）
for i=1:a*b
    t(i,1)=YABbar(i*n,1)+cdft("T",phiE,0.975,0.025)*sqrt(VE/n);
    t(i,2)=YABbar(i*n,1)-cdft("T",phiE,0.975,0.025)*sqrt(VE/n);
    t(i,3)=YABbar(i*n,1)+cdft("T",phiE,0.975,0.025)*sqrt(VE*(1+1/n));
    t(i,4)=YABbar(i*n,1)-cdft("T",phiE,0.975,0.025)*sqrt(VE*(1+1/n));
end

t//結果表示

上記を計算すると以下のような結果を得ることができます。

分散分析表

要因	S	φ	V	F	p
因子A	6.0202458	3	2.0067486	5.5928799	0.0123468
因子B	11.017675	2	5.5088375	15.353326	0.0004922
交互作用A✕B	9.5343917	6	1.5890653	4.4287816	0.0136891
誤差S_E	4.30565	12	0.3588042
	30.8779625	23

区間推定結果

A	B	上側信頼区間	下側信頼区間	上側予測区間	下側予測区間
1	1	8.70	6.85	9.37	6.18
1	2	9.37	7.53	10.05	6.85
1	3	9.50	7.65	10.17	6.98
2	1	10.12	8.28	10.80	7.60
2	2	10.75	8.91	11.43	8.23
2	3	9.78	7.93	10.45	7.26
3	1	10.90	9.05	11.57	8.38
3	2	11.18	9.34	11.86	8.66
3	3	8.58	6.74	9.26	6.06
4	1	9.92	8.08	10.60	7.40
4	2	10.44	8.59	11.11	7.92
4	3	7.38	5.54	8.06	4.86

一般的にF値2以上を因子の効果がある閾値として使われます。分散分析表を確認すると、因子A、因子B、交互作用A✕Bにおいて2以上のF値を有しています。そのため、それぞれが因子として、出力Yに対して効果があることが分かります。

もし、交互作用の要因がないと判断される場合は、S_（A✕B）を、S_Eに加えてモデルを考えるのが妥当です。その際は自由度も足し合わされます。
上記のように、変動が認められない要因を誤差にあわせることを、プーリング(pooling)と呼ぶそうです。

また、信頼区間の推定も実施しているため、所望の結果が欲しいものを水準に選択することができます。

さて、次回は、もっと因子が増える場合について考えたいと思います。直交表の出番ですね。

■追記　Rで解く

上記の例題をRで解いてみたいと思います。
データはcsvファイルから読み込む形で、因子は質的因子が分かるように、A,Bなどの文字で違いを表現するようにします。

分散分析を以下の１行で実現できるのは感動的です。
aov(formula = y~A*B)
ちなみに、AxBの作用が無視でき、プーリングしたい場合は、以下で可能です。
aov(formula = y~A+B)

各因子の平均を求めるスマートな方法がわからず、そこが煩雑になってしまった印象を抱きます。

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

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

#分散分析の実行
s <- aov(formula = y~A*B)
aovT <- anova(s)
aovT#分散分析表の表示
interaction.plot(B,A,y)
interaction.plot(A,B,y)

#区間推定の実行
VE <- aovT$`Sum Sq`[4]/aovT$Df[4]

A_1 = subset(d,Input_A=='A')
A_2 = subset(d,Input_A=='B')
A_3 = subset(d,Input_A=='C')
A_4 = subset(d,Input_A=='D')

M_a <- tapply(A_1$Output_Y,A_1$Input_B,mean)
M_b <- tapply(A_2$Output_Y,A_2$Input_B,mean)
M_c <- tapply(A_3$Output_Y,A_3$Input_B,mean)
M_d <- tapply(A_4$Output_Y,A_4$Input_B,mean)

M =c(M_a,M_b,M_c,M_d)#各因子の平均

ueA <- M + qt(0.975,aovT$Df[4]) *sqrt(VE/2)
shitaA <- M - qt(0.975,aovT$Df[4]) *sqrt(VE/2)


#結果のまとめ
Est = data.frame(
  InputA = c('A','A','A','B','B','B','C','C','C','D','D','D'),
  InputB = c('a','b','c','a','b','c','a','b','c','a','b','c'),
  Uegawa = ueA,
  Shitagawa = shitaA
)

上記の結果は以下のような形です。

Analysis of Variance Table

Response: y
          Df  Sum Sq Mean Sq F value    Pr(>F)  
A          3  6.0202  2.0067  5.5929 0.0123468 *
B          2 11.0177  5.5088 15.3533 0.0004922 ***
A:B        6  9.5344  1.5891  4.4288 0.0136891 *
Residuals 12  4.3057  0.3588  

グラフとデータフレームも以下に掲載します。

■追記2　Rのプログラムを少し改良

各A✕B因子の平均を求める部分を少し改良しました。
Rに慣れてないとかえってわかりづらいかもしれません。

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

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

#分散分析の実行
s <- aov(formula = y~A*B)
aovT <- anova(s)
aovT#分散分析表の表示

#区間推定の実行
VE <- aovT$`Sum Sq`[4]/aovT$Df[4]

M <-c(NULL)
for(i in levels(d$Input_A)){
  tmp_A <- subset(d,d$Input_A == i)
  tmp_mean <- tapply(tmp_A$Output_Y,tmp_A$Input_B,mean)
  M <- c(M,tmp_mean)
}

ue <- M + qt(0.975,aovT$Df[4]) *sqrt(VE/2)
shita <- M - qt(0.975,aovT$Df[4]) *sqrt(VE/2)

#結果のまとめ
Est = data.frame(
  InputA = c('A','A','A','B','B','B','C','C','C','D','D','D'),
  InputB = c('a','b','c','a','b','c','a','b','c','a','b','c'),
  Uegawa = ue,
  Shitagawa = shita
)

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



他参考書籍です