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

前々回は１因子の量的因子の一次モデルについて
前回は　１因子の量的因子の二次モデルについて　扱いました。

続いては、２因子/質的因子について考えていきたいと思います。

２因子になると交互作用の考慮も必要になるので、それについての議論が重要な要素になります。
しかし、今回は比較しやすいのと、個人的な学習も含めて１因子の分散分析や区間推定について、まずはまとめます。

■一因子の質的因子をまずは確認

最も基本的な一因子の質的因子についてまずは確認しておきたいと思います。
以下のようなデータがあるとします。

入力A	出力y
A	9.7
A	8.7
A	10.2
A	11.3
A	11.2
A	11.7
B	9.8
B	11.8
B	13.1
B	10.9
B	11.3
B	10.3
C	9.2
C	10
C	10.2
C	8.9
C	10.4
C	10.6
D	13.1
D	12.6
D	12.7
D	12.6
D	14.3
D	12.9
E	10.8
E	10.5
E	13
E	11.9
E	13.4
E	10.3

入力Aは質的因子です。

今までと同様に、分散分析表とt分布を使った区間推定を計算したいと思います。

区間推定について。
一次モデルや二次モデルを扱った際は量的因子に関する項が√の中にありましたが、それがなくなり、以下のようにかけます。

95%信頼区間

95%予測区間

■scilabで解く

では、実際にscilabで解いていきたいと思います。

//実験データ入力[水準番号　データ]
X = [
1    9.7
1    8.7
1    10.2
1    11.3
1    11.2
1    11.7
2    9.8
2    11.8
2    13.1
2    10.9
2    11.3
2    10.3
3    9.2
3    10
3    10.2
3    8.9
3    10.4
3    10.6
4    13.1
4    12.6
4    12.7
4    12.6
4    14.3
4    12.9
5    10.8
5    10.5
5    13
5    11.9
5    13.4
5    10.3
]

n=6;//１因子の実験繰り返し数
a=5;//因子の水準数
phiA = a-1;//自由度
phiE = a*(n-1)//自由度
Y = X(:,2);//出力Yだけ分離

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

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

//分散分析表の実際の計算（結果その①）
SA = (L*Y - J*Y)'*(L*Y - J*Y)//各因子平均-全平均
SE = (Y-L*Y)'*(Y-L*Y)//各データー各因子平均
VA = SA/phiA
VE = SE/phiE
f  = VA/VE
p = 1 - cdff("PQ",f,phiA,phiE)

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

//グラフ描画
x = [1:a]';
plot(X(:,1),Y,'*r');
plot(X(:,1),yAbar(:,1),'+b');
plot(x(:,1),t(:,1),'+g');
plot(x(:,1),t(:,2),'+g');
plot(x(:,1),t(:,3),'+k');
plot(x(:,1),t(:,4),'+k');
legend(["data","平均","上側信頼区間","下側信頼区間","上側予測区間","下側予測区間"],2);

上記を解くと以下のような結果が得られます。

分散分析表は以下のようにかけます

要因	S	φ	V	F	p
因子	34.944667	4	8.7361667	8.2014332	0.0002279
誤差	26.63	25	1.0652
	ST	φT

因子に対して優位な差異があることが分かります。

■個人的メモ

複数水準から２水準を選び検定を行うことを繰り返すと、それぞれの検定での有意水準は保証できても全体としてみると有意水準が保証できなくなることになる。

水準間の差を精確に比較するには、多重比較（multiple comparison）が必要になり、そのためにボンフェロニ法やテューキー法などがある

多重比較で検索すると参考になるページが多く出てきます。

■追記 Rを使って解いてみました

今後の自分のことを考えて、Rを使って解くことにしました。

初学者としては、以下が混乱しました。

・パッケージの理解

・行列やベクトル演算がMATLABと違う点

・=と<- or -<の使い方のち外

・データフレームの使い方

ネットで強くおすすめされた情報を信じて、RStudioを開発環境に選んでいます。

分散分析を調べると、Rコマンダーを使う場合の事例が多かったのですが、RコマンダーのパッケージやGUIを極力使わないでトライしてみました。実際のところ自分にはどちらが良いのか判断に迷っています。

データはcsvファイルから読み込む形になっています。

因子の違いを数字にすると、量的因子とみなされてしまったようなので、A,B,C,D,Eという形にしました。（この辺りもうまくできそうですが。）

記念すべき、初めてのRのプログラムなので、おかしな点あれば是非教えて欲しいです。

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

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

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

#区間推定の実行
VE <- aovT$`Sum Sq`[2]/aovT$Df[2]
M <- tapply(y,A,mean)
ue <- M + qt(0.975,aovT$Df[2]) *sqrt(VE/6)
shita <- M - qt(0.975,aovT$Df[2]) *sqrt(VE/6)

#結果の描画
xaxis <- 1:length(table(A))
plot(0, 0, type = "n", xlim = range(xaxis), ylim = range(y),
     xlab = "Input_A", ylab = "Output_y")
points(A, y,pch="+", col = "red")
points(xaxis, ue,pch="*", col = "green")
points(xaxis, shita,pch="*", col = "green")

legend("topleft", legend = c("data","interval estimation"),pch=c("+","*"), col = c("red","green"))

上記を実行すると以下のような出力が得られます。

Analysis of Variance Table

Response: y
          Df Sum Sq Mean Sq F value    Pr(>F)    
A          4 34.945  8.7362  8.2014 0.0002279 ***
Residuals 25 26.630  1.0652                      
---

Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

今後もRを使った計算結果を増やしていければと思っています。

■追記その２ Rのプログラムを変えてみました。

上記のプログラムは区間推定の算出をベクトル計算にて行っていました。

 今回は豊富なRの関数を使って算出したいと思います。

実際には、t.testを用いて信頼区間を求めました。 

上記のプログラムとどちらが良いのか？と問われるとどちらでも良い気もしますが、Rっぽいプログラムは今回のほうな気がしています。 

（予測区間も何か適してた関数がないかと探していますが、よく分かっていません。）

※上記のように考えてプログラムを構築したのですが、よくよく考えてみると下記のプログラムは最初のプログラムと同じになりません。

　何故かというと、各因子毎に誤差を計算するのと、全ての試験のトータルで誤差を計算することに違いがあるからです。

　私のようなミスをしないように気をつけてください。

#作業ディレクトリの設定
setwd("作業ディレクトリを記入してください")

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

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

#区間推定の実行
ue <- c(NULL)
shita <-c(NULL)
for(i in levels(d$Input_A)){
  tmp <- subset(d, A==i)
  results <- t.test(tmp[,2])
  ue <- c(ue,results$conf.int[2])
  shita <- c(shita,results$conf.int[1])
}

#結果の描画
xaxis <- 1:length(table(A))
plot(0, 0, type = "n", xlim = range(xaxis), ylim = range(y),
     xlab = "Input_A", ylab = "Output_y")
points(A, y,pch="+", col = "red")
points(xaxis, ue,pch="*", col = "green")
points(xaxis, shita,pch="*", col = "green")
legend("topleft", legend = c("data","interval estimation"),pch=c("+","*"), col = c("red","green"))

このページは以下の書籍を参考に自習用にまとめました

他関連書籍は以下です。