直交表は必要ない！？一部実施要因配置計画の作り方

一部実施要因計画の作り方
ジェネレータを定め、交絡関係を知る
L8直交表を作ってみよう
なぜ"一部実施"でよいのか？
交絡行列（Aliase Matrix）
直交表や線点図は必要ないのか？

一部実施要因計画の作り方

直交表や線点図は必要ありません。本記事では、一部実施要因配置計画の一番簡単な作成方法について説明します。

一部実施要因計画とは、要因配置計画（要因全ての組合わせを実験する計画）の一部だけを実施する計画です。試行回数は2^k-pになります。ここで、kは要因の数、pはk要因の完全実施の要因配置計画の2p分の1であることを表しています。

一部実施要因計画は次の2つのステップで作成できます。

要因がk-p個の場合の要因配置計画を書き出す。これが基本の計画。
他に必要な要因を、基本の計画の列を掛け合わせて作成する。

これだけです。以下にk=4, p=1の場合を例として示します。

まず、要因がk-p(=4-1=3)個の場合の完全実施要因配置を書き出します（表1）。これは単純に3つの要因に対する2水準の全組合せです。この試行回数8回の計画が基本の計画となります。

f:id:Sturgeon:20200510104907p:plain — 表1．3要因の完全実施要因配置計画（2^4-1計画の基本計画）

要因の数は4（k=4）なのでA、B、Cとは別に新たな要因Dを加えます。このとき、

$D＝ABC$

の関係を満たすようにDの列を作ります。すると2^4-1の一部実施要因計画は表2のようになります。

f:id:Sturgeon:20200510111154p:plain — 表2．2^4-1の一部実施要因計画

たしかに、D=ABCの関係になっていますね。これで一部実施要因配置計画の作成は完了です。

ジェネレータを定め、交絡関係を知る

表2ではD=ABCの関係でD列を作成しましたが、このような関係を決定関係（Defining Relation）と呼びます。今回のような2水準実験において、決定関係の式は同じ文字を2度かけると必ず＋になります。どいうことかというと、

$D=ABC\rightarrow D*D=I=ABCD$

のように変形できます（Iはすべての列が＋という意味）。このABCDを計画のジェネレータと呼びます。よって、一部実施要因計画は基本計画とジェネレータを決めることで作成できると言えます。

また、ジェネレータを決めることで、効果同士の交絡関係を知ることが出来ます。たとえば、D=ABCの両辺にCをかけて

$CD=ABC*C=AB$

したがってCDとABは交絡することが分かります。

ところで、ジェネレータI=ABCDは適当に選んだわけではなくて、できるだけ多くの要因を掛け合わせたジェネレータを使っています。こうすることで、低次の効果が交絡しないようにするのです。

仮に、ジェネレータをI=ABDに選んだとすれば、決定関係はD=ABとなりますから1次の項Dと交互作用ABが交絡してしまいます。また両辺にBをかければB=ADなのでBにも交互作用ADとの交絡関係ができてしまいますね。

このように、ジェネレータの選び方で交絡関係が決まるため、注意が必要です。

※RなどのプログラムやMinitabのようなソフトウェアを使えば、低次の項同士の交絡がもっとも少なるような計画（Minimum aberration designという計画）を自動で生成してくれますから、ジェネレータについてはあまり気にしないで大丈夫です。

L8直交表を作ってみよう

L8直交表もジェネレータの考えで作ることができます。

A, B, C,は表2と同じように2×2×2の全パターンを用意しておきます。そこから、D=±AB, E=±BC, F=±CA, G=±ABCで残りの4因子を追加することで（±はプラスとマイナスどちらでもOK）、L8直交表ができます。下図のようにExcelで作れます。

なぜ"一部実施"でよいのか？

ところで、要因配置計画では要因と水準の全組合せを実施しますが、その一部だけを実施する一部実施要因配置計画は、そもそもなぜ一部だけで良いのでしょうか？

この理解のためには、効果を表すための数式モデルについて説明する必要があります。実験で得る値をyとすると、モデルは

$y_i=\mu+a_i+b_i+c_i+d_i+\varepsilon_i$

と表されます。ここで、iは試行回を表しており、今回はi=1,2,...8。εは平均0・分散1の正規分布に従うばらつき。µは全試行の平均を示しています。a, b, c, dはそれぞれ要因A, B, C, Dの効果で、今回は2水準ですから、aは水準がーのときはa_－、＋のときはa₊とします。効果は平均からのズレを表すため

$a_{-}+a_{+}=0$

の関係が成り立ちます（b, c, dについても同様）。

さて、表3に実験計画（表2）とモデル式を一緒にしたものを示します。

f:id:Sturgeon:20200510133724p:plain — 表3. 2^4-1の一部実施要因計画とモデル式

そして、要因Aの効果を求めたいときは、Aの水準が+の試行の合計と、Aの水準がーの試行の合計を計算し、両者の差を取ります。すわなち、

$(y_5+y_6+y_7+y_8)-(y_1+y_2+y_3+y_4)=\\ \{4\mu+4a_{+}+4(b_{-}+b_{+})+4(c_{-}+c_{+})+4(d_{-}+d_{+})\}\\ -\{4\mu+4a_{-}+4(b_{-}+b_{+})+4(c_{-}+c_{+})+4(d_{-}+d_{+})\} \\ =4(a_{+}-a_{-})$

です。このように要因Aの効果だけが残り、その他の効果は相殺されることが分かります（εも相殺されますが、上式では省略しました）。相殺されるのは、すべての要因の水準（＋とー）が同じ数ずつ、計画に現れている（バランスしている）からです。

以上より、要因と水準のすべての組合せを実施しない一部実施要因計画でも、要因効果の推定ができることがわかりました。

交絡行列（Aliase Matrix）

ジェネレータによって交絡関係が決まることを先ほど紹介しましたが、すべての交絡関係を一斉に確認する方法として、交絡行列の計算が挙げられます。今、前述の計画は忘れて、回帰分析によって推定したいモデル式を

$\boldsymbol{y}=\boldsymbol{X_m\beta_m}+\varepsilon$

とするとき、真のデータが

$\boldsymbol{y}=\boldsymbol{X_m\beta_m}+\boldsymbol{X_t\beta_t}+\varepsilon$

で表されるとしましょう。ここで、βは回帰係数、Xはデザインマトリックス（詳しくはこちら）です。

すると、βmの推定値は

$E(\boldsymbol{\beta_m})=\boldsymbol{\beta_m}+\boldsymbol{(X_m'X_m)}^{-1}\boldsymbol{X_mX_t\beta_2}=\boldsymbol{\beta_m}+\boldsymbol{A\beta_t}$

で表されます。仮に1~4の要因についてベクトル表記すると、上式は

$E\begin{bmatrix}\beta_0 \\\beta_1 \\ \beta_2 \\\beta_3 \end{bmatrix}=\begin{bmatrix}\beta_0 \\ \beta_1+\beta_{23} \\ \beta_2+\beta_{13} \\ \beta_3+\beta_{12} \end{bmatrix}$