2018年統計検定一級(応用理工)略解追加

#はじめに

- 前回の記事で除外した応用理工問4の調査結果のまとめ。
- 公式の略解 2018年11月25日試験(http://www.toukei-kentei.jp/pt/pt-5956/#seikai)

#問4

ある化学物質の合成において、合成物質の収量を最大化する要因を探索し、それらの最適水準を決定したい。
収量に影響のある要因候補をA,B,C,D,Eの5種類とし、スクリーニング実験を行うこととする。

[1]各要因はそれぞれ2水準(-1,1)ずつとし、8回の実験を行うものとする。
以下に示す4因子(A,B,C,D)の実験計画表(S)と5因子(A,B,C,D,E)の実験計画表(T)の空欄部分の水準を埋めて、
適当と思われる実験計画を示せ。

|実験 | A | B | C | D ||実験 | A | B | C | D | E |
|:--:|:--:|:--:|:--:|:--:|:--:|:--:|:--:|:--:|:--:|:--:|:--:|
|1 | 1 | 1 | 1 | 1 ||1 | 1 | 1 | 1 | 1 | 1 |
|2 | 1 | 1 | -1 | -1 ||2 | 1 | 1 | -1 | -1 | 1|
|3 | | | | ||3 | | | | | |
|4 | | | | ||4 | | | | | |
|5 | | | | ||5 | | | | | |
|6 | | | | ||6 | | | | | |
|7 | | | | ||7 | | | | | |
|8 | | | | ||8 | | | | | |

<details><summary>解答例</summary><div>

実験計画表には行や列の選択に任意性があるため、答えは一意に定まらない。
略解から逆算すると以下のようになる。縦の列は互いに直交している。
4因子の実験計画表に要因Eの列を追加したものを5因子の実験計画表とする。

|実験 | A | B | C | D ||実験 | A | B | C | D | E |
|:--:|:--:|:--:|:--:|:--:|:--:|:--:|:--:|:--:|:--:|:--:|:--:|
|1 |  1 |  1 |  1 |  1 ||1 |  1 |  1 |  1 |  1 |  1 |
|2 |  1 |  1 | -1 | -1 ||2 |  1 |  1 | -1 | -1 |  1|
|3 | -1 | -1 |  1 |  1 ||3 | -1 | -1 |  1 |  1 |  1 |
|4 | -1 | -1 | -1 | -1 ||4 | -1 | -1 | -1 | -1 |  1 |
|5 |  1 | -1 |  1 | -1 ||5 |  1 | -1 |  1 | -1 | -1 |
|6 |  1 | -1 | -1 |  1 ||6 |  1 | -1 | -1 |  1 | -1 |
|7 | -1 |  1 |  1 | -1 ||7 | -1 |  1 |  1 |  -1 | -1 |
|8 | -1 |  1 | -1 |  1 ||8 | -1 |  1 | -1 |  1 | -1 |

</div></details>

[2]4因子実験の実験計画につき、実験結果を$y_1,\dots,y_8$とするとき、
要因Aの主効果の計算式、および要因Aと要因Bの2因子交互作用の計算式を示せ。

<details><summary>解答例</summary><div>

要因Aの主効果の計算式の係数は、列の直交性から、要因Aの列の係数となる。

```math
\begin{align*}
\frac{1}{8}(y_1+y_2-y_3-y_4+y_5+y_6-y_7-y_8)
\end{align*}
```

要因Aと要因Bの2因子交互作用の計算式の係数は、それぞれの要因の列の係数の積となる。

```math
\begin{align*}
\frac{1}{8}(y_1+y_2+y_3+y_4-y_5-y_6-y_7-y_8)
\end{align*}
```

ところで、全ての交互作用を考慮した構造モデルは以下のようになる。
ここで$x_{ik}$は実験$k$における要因$i$の水準である。

```math
\begin{align*}
y_k=\mu + \sum_i \beta_i x_{ik} + \sum_{i\ne j} \beta_{ij} x_{ik}x_{jk}+\epsilon_k
\quad(i,j=A,B,C,D;k=1,\dots,8)
\end{align*}
```

この場合、要因Aと要因Bの交互作用項$\beta_{AB}$と要因Cと要因Dの交互作用項$\beta_{CD}$については、
$x_{Ak}x_{Bk}=x_{Ck}x_{Dk}(k=1,\dots,8)$となるため、同時に推定することはできない。

```math
\begin{align*}
\sum_{k=1}^8 x_{ik}&=0\quad(i=A,B,C,D)\\
\text{主効果同士}\quad\sum_{k=1}^8 x_{ik}x_{jk}&=8\delta_{ij}\quad(i,j=A,B,C,D)\\
\text{主効果と交互作用}\quad\sum_{k=1}^8 x_{ik}x_{jk}x_{lk}&=0\quad(i,j,l=A,B,C,D)\\
\text{交互作用同士}\quad\sum_{k=1}^8 x_{ik}x_{jk}x_{lk}x_{mk}&=\mathrm{0\ or\ 8}\quad(i,j,l,m=A,B,C,D)
\end{align*}
```

色々と試行錯誤してみたが、全ての交互作用項が推定可能となる水準の組合せを見つけることはできなかった。
主効果と交互作用の列は互いに直交するようにできるが、交互作用同士の列には直交しない組合せが生じる。
デザイン行列を$X$とすると$X'X$は常に特異となる。交互作用をどこまで考慮すれば良いかが分からない。

</div></details>

[3]5因子実験の実験計画につき、各要因の主効果と各要因間の2因子交互作用との関係を論ぜよ。
具体的に、要因Aの主効果と交絡する2因子交互作用はどれかを示せ。

<details><summary>解答例</summary><div>

要因Aの列と二つの要因の列の積との内積を計算すると、ゼロでないのは要因Bと要因Eの組合せとなる。
但し、追加した要因Eの列の内容により交絡する要因の組合せは変化する。
要因Bと要因Eだけでなく、要因Cと要因Eや要因Dと要因Eもあり得る。
4因子実験の場合とは異なり、主効果と交互作用の列には互いに直交しない組合せが生じる。
全ての交互作用を考慮したデザイン行列を$X$とすると$X'X$は常に特異となる。
要因間の交互作用をどこまで考慮すれば良いかが分からない。

</div></details>

[4]4因子実験について、各要因の主効果の推定を行ったところ、以下の結果を得た。
各要因の主効果に関する検定を行い、5%有意となる要因をすべて示せ。

|要因 | 係数  | 標準偏差  |
|:--:|:--:|:--:|
|切片 | 9.25 | 0.78 |
|A | 2.75 | 0.78 |
|B | 4.00 | 0.78 |
|C | 1.00 | 0.78 |
|D | -0.75 | 0.78 |

<details><summary>解答例</summary><div>

統計検定量Z=係数/標準偏差を計算して、両側検定の上側2.5%点(1.96)と比較する。

```R
df = data.frame(
  v = c("const","A","B","C","D"),
  coef = c(9.25,2.75,4.00,1.00,-0.75),
  sd = c(0.78,0.78,0.78,0.78,0.78)
)
df$Z = df$coef / df$sd
df$sig = abs(df$Z) > qnorm(0.975)
df

      v  coef   sd          Z   sig
1 const  9.25 0.78 11.8589744  TRUE
2     A  2.75 0.78  3.5256410  TRUE
3     B  4.00 0.78  5.1282051  TRUE
4     C  1.00 0.78  1.2820513 FALSE
5     D -0.75 0.78 -0.9615385 FALSE
```

</div></details>

[5]実験結果から、共に連続的な「反応温度」と「反応時間」の影響が大きいことが分かった。
それらの水準の最適な組合せを探索するため、$x_1$を反応温度、$x_2$を反応時間とした
二次多項式$y=\beta_0+\beta_1 x_1+\beta_2 x_2+\beta_{11}x_1^2+\beta_{22}x_2^2+\beta_{12}x_1x_2$の当てはめを行う。
採用する実験計画は複合中心計画とし、実験点は$(0,0),(\pm 1,\pm 1),(0,\pm\sqrt{2}),(\pm\sqrt{2},0)\ $とする。
中心点で4回実験を行い、その他の点では一回ずつの計12回の実験を行うとしたとき、
パラメータ$\beta=(\beta_0,\beta_1,\beta_2,\beta_{11},\beta_{22},\beta_{12})'\ $に関する$12\times6\ $のデザイン行列$X$に対し、
以下となることが示される。

```math
\begin{align*}
(X'X)^{-1}=\frac{1}{32}
\begin{pmatrix}
8 & 0 & 0 & -4 & -4 & 0 \\
0 & 4 & 0 & 0 & 0 & 0 \\
0 & 0 & 4 & 0 & 0 & 0 \\
-4 & 0 & 0 & 5 & 1 & 0 \\
-4 & 0 & 0 & 1 & 5 & 0 \\
0 & 0 & 0 & 0 & 0 & 8
\end{pmatrix}
\end{align*}
```

この計画は回転可能性(rotatability)を満たすこと、すなわち、
実験領域内のある点$x=(x_1,x_2)'$での予測分散が$x_1^2+x_2^2$の関数となることを示せ。

<details><summary>解答例</summary><div>

```math
\begin{align*}
\sigma^2x'(X'X)^{-1}x=\frac{\sigma^2}{32}(8-4(x_1^2+x_2^2)+5(x_1^2+x_2^2)^2)
\end{align*}
```

```R
xy = rbind(
    c(0,0),
    c(0,0),
    c(0,0),
    c(0,0),
    c(-1,-1),
    c(-1, 1),
    c( 1,-1),
    c( 1, 1),
    c(0, sqrt(2)),
    c(0,-sqrt(2)),
    c( sqrt(2),0),
    c(-sqrt(2),0)
)
Xcol = function(x1, x2){
  cbind(c=1,x1=x1,x2=x2,x11=x1^2,x22=x2^2,x12=x1*x2)
}
X = Xcol(xy[,1],xy[,2])
solve(t(X) %*% X)*32

     c x1 x2 x11 x22 x12
c    8  0  0  -4  -4   0
x1   0  4  0   0   0   0
x2   0  0  4   0   0   0
x11 -4  0  0   5   1   0
x22 -4  0  0   1   5   0
x12  0  0  0   0   0   8
```

</div></details>

# 参考文献

- 実験計画法入門 Part 1(https://www.slideshare.net/hajimemizuyama/par-49440366)
- 実験計画法入門 Part 2(https://www.slideshare.net/hajimemizuyama/part-2-51006255)
- 実験計画法入門 Part 3(https://www.slideshare.net/hajimemizuyama/part-3-51768519)
- 実験計画法入門 Part 4(https://www.slideshare.net/hajimemizuyama/part-4-52807120)
- Rで実験計画法 前編(https://www.slideshare.net/itoyan110/r-9139125)
- Rで実験計画法 後編(https://www.slideshare.net/itoyan110/r-14261638)