機械の開発で分散分析した話【要因実験】

#分散分析を使ったハードウエアの検証

機械のハードの開発をした時、どんな因子が影響を与えているのか、どの程度の影響を与えるのかなどのモデルを分散分析で作成した。
交互作用も見られるし使い勝手が良かった。

要因実験といえば実験法などは本やサイトで乗っているので割愛。

ノイズはN1とN2の二種類。
水準はデフォルトを含めて3つの水準を振った。
因子はA,B,Cとして、デフォルトの状態で測定したデータをA0と名前を付けている。
仮にこんな感じでテストデータを作った。

平均50という出力をする機械に対して、ノイズと、いろんな水準を振った時どんな影響が観察できるのか？

```R

N1A0 <- rnorm(100, 50, 0.1)
N2A0 <- rnorm(100, 50.5, 0.1)
N1A1 <- rnorm(100, 50, 0.1)
N1A2 <- rnorm(100, 50, 0.1)
N1B1 <- rnorm(100, 50, 0.1)
N1B2 <- rnorm(100, 50, 0.1)
N1C1 <- rnorm(100, 50, 0.1)
N1C2 <- rnorm(100, 50, 0.1)
N2A1 <- rnorm(100, 50.5, 0.1)
N2A2 <- rnorm(100, 50.5, 0.1)
N2B1 <- rnorm(100, 50.7, 0.1)
N2B2 <- rnorm(100, 50.9, 0.1)
N2C1 <- rnorm(100, 52, 0.1)
N2C2 <- rnorm(100, 54, 0.1)

value <- c(N1A0, N2A0, N1A1, N2A1, N1A2, N2A2, N1A0, N2A0, N1B1, N2B1, N1B2, N2B2, N1A0, N2A0, N1C1, N2C1, N1C2, N2C2)
noise <- factor(rep(rep(c("N1","N2"),each=100),9))
exp_no <- factor(
c(rep("A0",200),
rep("A1",200),
rep("A2",200),
rep("A0",200),
rep("B1",200),
rep("B2",200),
rep("A0",200),
rep("C1",200),
rep("C2",200))
)
df <- data.frame(value = value, noise=noise, exp_no=exp_no)
```
まずデータをとってきた後はボックスプロット。

```R

par(las=2)
boxplot(value ~  nois+exp_no, data=df, col=c(3,4), 
  ylim=c(49, 54))
abline(v=seq(2.5,12.5,2))
legend(-0.5, 54, "default",box.lty = 0)
legend(3, 54, "   A      ",box.lty = 0, bg="red")
legend(7, 54, "   B       ",box.lty = 0, bg="red")
legend(11, 54, "   C       ",box.lty = 0, bg="red")
```
![box.png](https://qiita-image-store.s3.amazonaws.com/0/275572/29fde97b-5579-5427-152b-6d54a6ca2e92.png)


視覚的にこれで分散分析になっている。

最初のA0を見ると、N1とN2の間で0.5ほど差がある。
N2というノイズは平均値を0.5上げる影響を持っていることが分かる。
A1,A2を見てみるとノイズ以外の影響は受けていない(A0のN1,N2間と同じ関係性)からAという水準は変化させても影響がないと分かった。

Bを見てみるとN2のグループはノイズの影響以上に影響を与えているように見える。
しかしN1では変化は見られないので、N2とBの水準の交互作用による影響が出ているとわかる。

Cでは、N1でも変化がみられる。
N1とN2の関係性自体は0.5の差の関係を維持している。
こうして観察するとC1とC2自体に平均値を上げる影響があることが分かる。

では、どの程度の影響なのかモデル化することを考える。
AOV関数を使って、その回帰係数を見て影響度を調べる。

```R
fit <- lm(aov(value~noise+exp_no, data=df))
result <- summary(fit)
res_coe <- result$coefficient#回帰係数を抽出
res_coe
```
![coe.png](https://qiita-image-store.s3.amazonaws.com/0/275572/ed8a45c5-4eb0-4c83-6f0e-efeaa363f798.png)

こんな結果になりました。

Interceptがデフォルトの状態の機械の出力値で、50に近いので問題ない。
noiseのN2による影響は平均値を1.118上げるということが分かる。(←実際は0.5なのに)
exp_noのA1,2を見てみるとほとんど0であり、影響は与えていない水準ということが分かる。t値も小さい。
水準BではB2が平均値を0.19上げている。（←実際はN2の時もっと影響がある）
水準CではC1の時0.7、 C2の時1.7ほど平均値を上げている。

ほとんど影響のない因子は機械の出力値のモデルに対して説明を持たないので省きます。
いかに大変な苦労をして取ったデータでも無理に複雑なモデルにするのは解釈を複雑にするだけです。

ということで、水準Aというデータをデフォルトの中に食わせます。行を消すと係数が変化してしまうの、A0に含ませます。

その前にr二乗値は今回の回帰式がどの程度データを説明できているかと解釈できるので、このr二乗値の変化が、すなわち影響を与える程度といえるでしょう。
r二乗値を確認。

```R
res_adj <- result$adj.r.squared
res_adj
> res_adj
[1] 0.6565302
```
残差も確認しておきましょう

```R
resi <- aov(value~noise+exp_no, data=df)
resi_hist <- residuals(resi)
hist(resi_hist)
```
![resi.png](https://qiita-image-store.s3.amazonaws.com/0/275572/1b850f74-c673-d455-6fae-0540e0c7224e.png)

残差はほとんど0の周辺にあるようですが、ばらついてますね。

では水準AをA0に食わせます。
そのためにも少しファクター列を変形して回帰の計算をわかりやすいものにします

```R

exp_A <- factor(
c(rep("A0",200),
rep("A1",200),
rep("A2",200),
rep("A0",1200)))
exp_B <- factor(
c(rep("A0",600),
rep("A0",200),
rep("B1",200),
rep("B2",200),
rep("A0",600)))
exp_C <- factor(
c(rep("A0",1200),
rep("A0",200),
rep("C1",200),
rep("C2",200)))
df2 <- data.frame(value = value, noise=noise, 
  exp_A = exp_A, exp_B = exp_B, exp_C=exp_C)
fit <- lm(aov(value~noise+exp_A+exp_B+exp_C, data=df2))
result <- summary(fit)
res_coe <- result$coefficient
res_coe
```
![df2.png](https://qiita-image-store.s3.amazonaws.com/0/275572/5cac7347-0a62-fdf3-c9e2-255a872688f3.png)

係数は変わってないことが確認できたと思います。
水準Aを消しましょう。

![delA.png](https://qiita-image-store.s3.amazonaws.com/0/275572/dc1a7f57-5786-091c-68dc-40157c942a83.png)

```R
res_adj <- result$adj.r.squared
res_adj
> res_adj
[1] 0.6569123
```

r二乗値も大きく変化していないので"影響のない水準"と言えそうです。

さらにBも消してみると、

![bb.png](https://qiita-image-store.s3.amazonaws.com/0/275572/ef9b691c-89cd-d5bc-d8c1-8bb1cfc9c7ac.png)

```R
> res_adj
[1] 0.6530927
```

これでも大きく変化はなさそうです。

こうして影響のない因子を確かめていくと、

水準Cを消したlm(aov(value~noise+exp_A+exp_B, data=df2))のとき

```R
> res_adj
[1] 0.3801139
```
ノイズを消したlm(aov(value~exp_A+exp_B+exp_C, data=df2))のとき

```R
> res_adj
[1] 0.3237663
```
の二つがr二乗値を引き下げるものと分かりました。

よって、この機会はノイズと水準Cで説明されるモデルで説明できます。

機械の出力値 = 水準Cの状態 + ノイズ

水準Bはr二乗値としては変化しませんでしたが、微小の変化でも許されない出力であるのならば、より細かい水準の振り方や、相互作用について検証する必要があると思います。そこは実用の人が考えることということで。

こんなやり方もあるという紹介でした。