統計学の基礎をすぐに思い出せるような、まとめチートシートを作成しました。
「確率分布編」「推定・検定編」「線形回帰編」「その他の話題（理工学分野）」の４部構成です。本ページは「線形回帰編」です。
注意
あくまで個人用のメモであり、誤植や不正確な情報が含まれている可能性があります。適宜更新していきます。
9 線形回帰モデル
線形回帰モデルは多くの教科書では単回帰モデルから説明されることが多いが、線形代数に慣れている者であれば、むしろ重回帰モデルの一般論を説明し、単回帰モデルはその特別な場合とみなすほうが見通しがよいであろう。ここでは、重回帰モデルの一般論から説明する。
重回帰モデルの設定
データ数を $N$、説明変数の個数を $k$ とし、各データ $i=1,\ldots,N$ について次の線形モデルを考える。
\begin{align} y_i &\overset{?}{=} \beta_0+\beta_1x_{i1}+\cdots+\beta_kx_{ik} \end{align}ここで $y_i$ は被説明変数、$x_{ij}$ は第 $j$ 番目の説明変数である。 $i$ は観測の添字で、 $j$ は説明変数の添字であることに注意しよう。以下では我々は $\beta_j$ （ $j=0,1,\cdots,k$ ）の本当の値を知らず、これをこれから推定するものとする。
実際の問題はこれだけでは説明できない揺らぎが生じる。そこで誤差項 $u_i$ を入れて、以下によりデータが作成されると考える。
ここで誤差項 $u_i$ は確率変数である。行列で書くと、以下のようになる。
ただし、各記号は次で与えられる。
\begin{align} \boldsymbol{y} &= \begin{pmatrix} y_1\\ \vdots\\ y_N \end{pmatrix}, \quad X = \begin{pmatrix} 1 & x_{11} & \cdots & x_{1k}\\ \vdots & \vdots & & \vdots\\ 1 & x_{N1} & \cdots & x_{Nk} \end{pmatrix}, \quad \boldsymbol{\beta} = \begin{pmatrix} \beta_0\\ \beta_1\\ \vdots\\ \beta_k \end{pmatrix}. \end{align}以下では、$\mathrm{rank}(X)=k+1$、すなわち $X^\top X$ が正則であると仮定する（ $X^\top X$ が正則でない場合の意味は、後に説明する）。
典型的な仮定は、誤差項について次が成り立つというものである。
これは「観測間の誤差が無相関であり、線形回帰で説明できない誤差は、方向性のある偏りを残さない」という仮定である。正規性まで仮定する場合には、さらに $\boldsymbol{u}\sim \mathcal{N}(\boldsymbol{0},\sigma^2I_N)$ と置く。
最小二乗推定量
まずはデータ $X, \boldsymbol{y}$ は既知として、パラメータ $\boldsymbol{\beta}$ を推定する方法を考えよう。パラメータの推定値を $\hat{\boldsymbol{\beta}}$ と書く。もちろんこれは、データ $X, \boldsymbol{y}$ の関数であることが期待される。
最小二乗法では、残差平方和を最小にするように $\boldsymbol{\beta}$ を選ぶ。つまり、次の目的関数を最小化する。
$\boldsymbol{\beta}$ で微分すると、正規方程式が得られる。
\begin{align} X^\top(\boldsymbol{y}-X\hat{\boldsymbol{\beta}})=\boldsymbol{0} \label{eq:normal-equation} \end{align}この式は、残差 $\boldsymbol{y}-X\hat{\boldsymbol{\beta}}$ が、説明変数の列空間と直交することを表している。したがって、最小二乗法は「縦方向の誤差を小さくする計算」であると同時に、「$\boldsymbol{y}$ を説明変数の空間に直交射影する計算」でもある。$X^\top X$ が正則なら、最小二乗推定量は次の形になる。
この式の意味は、「$\boldsymbol{y}$ を $X$ の列ベクトルで最もよく表す係数を選んでいる」ということである。
ちなみに、最小二乗推定量は、$\boldsymbol u$ が正規分布に従うときの最尤推定量と一致する。つまり、正規性を仮定しても、最小二乗推定量は変わらない。このことは、正規分布を仮定したときの対数尤度に $(\ref{eq:objective-function})$ 式と同じ形が出てくることを考えれば明らかであろう。
ハット行列と残差
推定値を $\hat{\boldsymbol y}=X\hat{\boldsymbol{\beta}}$ とおく。上の最小二乗推定量を代入すると、
と書ける。ここで次の行列をハット行列という。
$\boldsymbol{y}$ にハットをつけて推定値 $\hat{\boldsymbol y}$ にするので「ハット行列」である。
ハット行列の性質はこの先さまざまな計算で多用するので、ここでまとめておこう。まず $H$ は回帰分析における射影行列なので、次の性質を持つ。
また、以下の式も大変重要である。
最初の式は、残差が説明変数の列空間と直交していることを表している。$2$ 行目は $1$ 行目の転置を取ったものであり、$(\ref{eq:normal-equation})$ 式の亜種である。最後の式は $2$ 行目の特別な場合で、$X^\top$ の $1$ 列目が $\boldsymbol{1}$ ベクトルになっているものを抜き出したものである。
射影行列（補足）
寄り道になるが、ハット行列 $\hat H$ は射影行列の性質を満たしていることを確認しよう。ここで射影行列とは、次の条件を満たす正方行列のことである。
以外に思うかもしれないが、射影行列の定義はこれで十分である。というのも、$\boldsymbol{v}$ をベクトルとすると、$1=P+(1-P)$ より、$\boldsymbol{v} = P\boldsymbol{v} + (1-P)\boldsymbol{v}$ と分解でき、かつ上の定義より
\begin{align} (1-P)^T P = (1-P) P = 0 \end{align}より、分割された $2$ つのベクトルが互いに直交していることがわかるからである。
また、$P^2=P$ であることから、射影行列の固有値は $0$ または $1$ である。さらに、$P$ の $1$ 固有値の固有ベクトルは、$P$ が射影する空間 $M_P$ に属するベクトルである。特に、$M_P$ の正規直交基底を $\left< \boldsymbol{e}_1, \dots, \boldsymbol{e}_m \right>$ とすると、
\begin{align} P = \sum_{i=1}^{m} \boldsymbol{e}_i \boldsymbol{e}_i^\top \end{align}となる。また、射影行列のトレースは、その射影する空間の次元に等しい。
\begin{align} \mathrm{Tr}\, P = \dim M_P \end{align}ハット行列の場合は、射影する空間の次元は $k+1$ であるから、トレースは $k+1$ になる。
ここで $(\ref{eq:normal-equation})$ 式の意味について改めて考えてみるのも面白い。というのも、任意の $\widetilde{\boldsymbol{\beta}}$ に対して、$(\ref{eq:normal-equation})$ 式より
\begin{align} (X \widetilde{\boldsymbol{\beta}})^\top(\boldsymbol{y}-X\hat{\boldsymbol{\beta}})=\boldsymbol{0} \end{align}が成り立つ。これは、残差ベクトル $\boldsymbol{y}-X\hat{\boldsymbol{\beta}}$ が、 $X \widetilde{\boldsymbol{\beta}}$ の形で説明できる任意のベクトル（部分空間）に直交していることを意味している。要するに、最小二乗法では、残差は「$X \widetilde{\boldsymbol{\beta}}$ の形で説明できる空間」に直交するように決まるのである。
分散推定
次に、確率変数 $\boldsymbol{u}$ の分散 $\sigma^2$ を推定しよう。残差平方和を $\mathrm{RSS}=\|\boldsymbol{y}-X\hat{\boldsymbol{\beta}}\|^2$ とおく。ハット行列の性質を用いると
\begin{align} E[\mathrm{RSS}]&=E\big[\|\boldsymbol{y}-X\hat{\boldsymbol{\beta}}\|^2\big]=E\big[\boldsymbol{y}^\top(1-H)^2\boldsymbol{y}\big]\\ &=E\big[\boldsymbol{u}^\top(1-H)\boldsymbol{u}\big]=(N-k-1)\sigma^2 \end{align}となる。したがって、誤差分散 $\sigma^2$ の不偏推定量は次で与えられる。
分母が $N$ ではなく $N-k-1$ になるのは、切片を含む $k+1$ 個の回帰係数を推定した分だけ自由度を消費しているためである。
推定量の平均・分散
$\boldsymbol{\hat \beta}$ を確率変数 $\boldsymbol{u}$ の関数と見て、その平均と分散を計算しよう。これは、パラメータ推定の誤差評価に対応する。
まず、$(\ref{eq:ols-estimator})$ に $\boldsymbol{y}=X\boldsymbol{\beta}+\boldsymbol{u}$ を代入して、次が得られる。
\begin{align} \hat{\boldsymbol{\beta} }-\boldsymbol{\beta} &= (X^\top X)^{-1}X^\top\boldsymbol{u} \label{eq:beta-error} \end{align}これより、$E[\boldsymbol{u}]=\boldsymbol{0}$ であれば、
である。言い換えれば、最小二乗推定量 $\hat{\boldsymbol{\beta} }$ は不偏量である。
また、$(\ref{eq:beta-error})$ より、$\mathrm{Cov}(\boldsymbol{u})=\sigma^2I_N$ であれば
である。
最後に、以上の結果を用いると、標準誤差 $\widehat{\mathrm{SE}}(\beta_j)$ を求めることができる。標準誤差とは標準偏差の推定値だったので、
\begin{align} \widehat{\mathrm{SE}}(\beta_j) &= \left.\sqrt{\mathrm{Cov}(\hat{\boldsymbol{\beta}})_{jj}}\right|_{\sigma^2 = \hat{\sigma}^2} \end{align}とかける。これに \eqref{eq:sigma-hat} および \eqref{eq:cov-beta} を用いると、標準誤差の推定式が得られる。
Gauss-Markov の定理
最小二乗推定量 $ \hat{\boldsymbol{\beta}}$ は、誤差項 $\boldsymbol{u}$ の正規性を仮定しなくても、（一定の条件のもとで）最も優れた推定量である。これを示すのが Gauss-Markov の定理である。
主張に入る前に、まずは用語を整理しよう。まずは正方行列 $A,B$ の大小関係（これは全順序ではない）を、以下のように定義する。
このとき例えば、任意の分散共分散行列 $\Sigma = E\big[ (\boldsymbol{X} - \overline{\boldsymbol{X} }) (\boldsymbol{X} - \overline{\boldsymbol{X} })^\top \big]$ が $\Sigma \geq 0$ であることが、以下の計算により分かる。
\begin{align} \boldsymbol{c}^\top\Sigma\boldsymbol{c} &= E\big[ \boldsymbol{c}^\top (\boldsymbol{X} - \overline{\boldsymbol{X} }) (\boldsymbol{X} - \overline{\boldsymbol{X} })^\top \boldsymbol{c} \big]\\ &= E\big[ (\boldsymbol{c}^\top (\boldsymbol{X} - \overline{\boldsymbol{X} }))^2 \big]\geq 0 \end{align}次に、線形不偏推定量について説明しよう。これは文字通り「線形」かつ「不偏」である推定量を指す。今回で言えば、推定量が定数行列 $A$ を用いて $\tilde{\boldsymbol{\beta} }=A\boldsymbol{y}$ の形で書け、かつ $E[\tilde{\boldsymbol{\beta} }]=\boldsymbol{\beta}$ が成り立つことである。
最後に、「最良」の意味を説明しよう。これは平たく言えば「分散共分散行列を最小にする」ということであり、より厳密に言えば、任意の推定量 $\tilde{\boldsymbol{\theta} }$ に対して、推定量 $ \hat{\boldsymbol{\theta}}$ が
\begin{align} \mathrm{Cov}(\tilde{\boldsymbol{\theta}}) - \mathrm{Cov}(\hat{\boldsymbol{\theta}}) \geq 0 \end{align}を満たすとき、$\hat{\boldsymbol{\theta}}$ は最良という。
以上の用語の準備のもとで、Gauss-Markov の定理の主張はこうである。$E[\boldsymbol{u}]=\boldsymbol{0}$、$\mathrm{Cov}(\boldsymbol{u})=\sigma^2I_N$、$\mathrm{rank}(X)=k+1$ を仮定する。このとき、最小二乗推定量 $\hat{\boldsymbol{\beta}}$ は、$\boldsymbol{\beta}$ の最良線形不偏推定量（BLUE）である。いいかれば、任意の線形不偏推定量 $\tilde{\boldsymbol{\beta} }$ に対して、最小二乗推定量 $ \hat{\boldsymbol{\beta}}$ は
ということである。つまり、どの方向 $\boldsymbol{c}$ で係数の線形結合 $\boldsymbol{c}^\top\boldsymbol{\beta}$ を見ても、最小二乗推定量より分散の小さい線形不偏推定量は存在しない。
正規誤差のもとでの分布
さらに正規性を仮定して、$\boldsymbol{u}\sim\mathcal{N}(\boldsymbol{0},\sigma^2I_N)$ とする。このとき $\hat{\boldsymbol{\beta}}$ は正規ベクトルの線形変換なので、
である。また、残差平方和については、
が成り立つ。さらに、正規線形モデルでは $\hat{\boldsymbol{\beta}}$ と $\mathrm{RSS}$ は独立である。この独立性によって、回帰係数の $t$ 検定や、複数制約に対する $F$ 検定がきれいな形で導かれる。雑に言えば、係数推定の揺らぎと、誤差分散を測るための残差平方和の揺らぎを、別々の部品として扱えるようになる。
回帰係数の検定
第 $j$ 番目の回帰係数について、$H_0:\beta_j=b$ と $H_1:\beta_j\neq b$ を検定したいとする。$(X^\top X)^{-1}$ の $(j,j)$ 成分を $a_{jj}$ と書くと、正規誤差のもとで
である。しかし通常は $\sigma^2$ が未知なので、$\hat{\sigma}^2=\mathrm{RSS}/(N-k-1)$ で置き換える。このとき帰無仮説のもとで、
が成り立つ。したがって、$|T_j|$ が自由度 $N-k-1$ の $t$ 分布の上側臨界値より大きければ、$H_0$ を棄却する。単に係数が $0$ かどうかを見るだけでなく、任意の基準値 $b$ との比較として書いておくと、検定の意味がはっきりする。
複数の線形制約を同時に検定したい場合には、一般に $F$ 検定を用いる。たとえば $q$ 個の制約をまとめて $R\boldsymbol{\beta}=\boldsymbol{r}$ と書けるなら、制約なしモデルと制約付きモデルの $\mathrm{RSS}$ の差を用いて、制約を入れることで当てはまりがどれほど悪化するかを見る。ここで見ているのは、「制約を入れても残差平方和がほとんど増えないなら、その制約はデータと大きく矛盾していない」ということである。
単回帰モデル
統計検定 $1$ 級などの試験では、単回帰モデルについて詳しく聞かれることが多い。そのため、単回帰モデルを「重回帰の特別な場合」などと馬鹿にせず、慣れ親しんでおくことが重要である。
説明変数が $1$ 個だけの場合、モデルは $y_i=\alpha+\beta x_i+u_i$ である。これは
\begin{gather} X = \begin{pmatrix} 1 & x_1\\ 1 & x_2\\ \vdots & \vdots\\ 1 & x_N \end{pmatrix} = \begin{pmatrix} \boldsymbol{1} & \boldsymbol{x} \end{pmatrix},\quad \boldsymbol{\beta} = \begin{pmatrix} \alpha\\ \beta \end{pmatrix} \end{gather}のときであり、このように見れば重回帰モデルの特殊ケースである。
実際に $\alpha, \beta$ を直接求めてみよう。そのための準備として、$\overline{x} = \frac{1}{N}\sum_{i=1}^{N}x_i, \ \overline{y} = \frac{1}{N}\sum_{i=1}^{N}y_i$ とおき、
\begin{gather} Q_{xx} = \sum_{i=1}^{N}(x_i-\overline{x})^2, \quad Q_{yy} = \sum_{i=1}^{N}(y_i-\overline{y})^2,\\ Q_{xy} = \sum_{i=1}^{N}(x_i-\overline{x})(y_i-\overline{y}) \end{gather}と書くことにする。このときまず
である。これを用いて、 $\boldsymbol{\beta} = (X^\top X)^{-1}X^\top y$ を計算・整理すると
である。傾き $\hat{\beta}$ は、$x$ と $y$ の共分散を $x$ の分散で割ったものである。したがって、$x$ が大きい観測ほど $y$ も大きくなる傾向があれば正になり、逆なら負になる。切片 $\hat{\alpha}$ は、回帰直線が点 $(\overline{x},\overline{y})$ を通るように決まる。
ちなみに、ここまで $\boldsymbol{\beta} = (X^\top X)^{-1}X^\top y$ を計算してきたが、実際に上式を計算すると分かる通り、この計算はやや面倒である。単回帰の場合は、最初から最小二乗法で計算したほうが早い。もちろん、計算しなくても上の結果は暗記しておくべきである。
重回帰のときに示したように、 $\hat{\alpha}, \hat{\beta}$ は不偏推定量である。また、分散共分散行列についても
とすればよい。
決定係数
回帰モデルの当てはまりを表す基本的な指標が決定係数 $R^2$ である。これは、推定値 $\hat{\boldsymbol{y}}$ と観測データ $\boldsymbol{y}$ の「データに関する相関関数の $2$ 乗」として定義される。
\begin{align} R^2 = \frac{ \{ (\hat{\boldsymbol{y}}-\overline{\boldsymbol{y}} )^\top (\boldsymbol{y}-\overline{\boldsymbol{y}}) \}^2 }{ |\hat{\boldsymbol{y}}-\overline{\boldsymbol{y}}|^2 |\boldsymbol{y} -\overline{\boldsymbol{y}}|^2 } \end{align}相関関数の $2$ 乗なので $R^2$ と書かれる。ここから、よく使われる $1-\mathrm{RSS}/\mathrm{TSS}$ の形を導こう。そのためには
\begin{align} (\hat{\boldsymbol{y}}-\overline{\boldsymbol{y}})^\top (\boldsymbol{y}-\hat{\boldsymbol{y}})=0 \end{align}を活用すればよい。上の式は、$\boldsymbol{y}-\hat{\boldsymbol{y}} = (1 - H) \boldsymbol{y}$ としてハット行列の性質を用いればすぐに分かる。したがって、
\begin{align} (\hat{\boldsymbol{y}}-\overline{\boldsymbol{y}} )^\top (\boldsymbol{y}-\overline{\boldsymbol{y}}) &= (\hat{\boldsymbol{y}}-\overline{\boldsymbol{y}} )^\top \{(\hat{\boldsymbol{y}}-\overline{\boldsymbol{y}}) + (\boldsymbol{y}-\hat{\boldsymbol{y}})\}\\ &= ||\hat{\boldsymbol{y}}-\overline{\boldsymbol{y}}||^2 \end{align}である。これを最初の定義に代入すると、
\begin{align} R^2 = \frac{||\hat{\boldsymbol{y}}-\overline{\boldsymbol{y}}||^2} {||\boldsymbol{y}-\overline{\boldsymbol{y}}||^2} \end{align}が得られる。つまり、相関の $2$ 乗として定義した $R^2$ は、全体のばらつきに対して、推定値のばらつきが占める割合としても読める。さらに同じ直交性を用いると、
\begin{align} ||\boldsymbol{y}-\overline{\boldsymbol{y}}||^2 &= ||\hat{\boldsymbol{y}}-\overline{\boldsymbol{y}}||^2 + ||\boldsymbol{y}-\hat{\boldsymbol{y}}||^2 \end{align}がすぐに分かる。これは大変重要な式である。左辺は 全平方和と呼ばれ、$ \mathrm{TSS} $ と書かれる。右辺は回帰平方和と残差平方和の和であり、$ \mathrm{ESS} + \mathrm{RSS} $ と書かれる。つまり、
である。ここで $\mathrm{TSS}$ は全体のばらつき、$\mathrm{ESS}$ は回帰式で説明できたばらつき、$\mathrm{RSS}$ は説明しきれずに残ったばらつきである。上で見たように $R^2=\mathrm{ESS}/\mathrm{TSS}$ であり、さらに平方和分解を $\mathrm{TSS}$ で割ると $1=\mathrm{ESS}/\mathrm{TSS}+\mathrm{RSS}/\mathrm{TSS}$ なので、次の形も得られる。
つまり $R^2$ は、目的変数 $\boldsymbol{y}$ の全体のばらつきのうち、回帰式によって説明できた割合を表している。通常の最小二乗回帰では、$0\leq R^2\leq 1$ が成り立つ。$\mathrm{RSS}$ が小さいほど $R^2$ は $1$ に近づき、回帰式がデータによく当てはまっていると見なせる。一方、$\mathrm{RSS}$ が $\mathrm{TSS}$ に近ければ $R^2$ は $0$ に近く、平均 $\overline{y}$ だけで予測する場合とあまり変わらない。
単回帰の場合には、切片を含むモデルで $R^2$ は標本相関係数の $2$ 乗に一致する。したがって、単回帰では「$x$ と $y$ の直線的な関係の強さ」を表す指標としても読める。ただし重回帰では、複数の説明変数をまとめてどれだけ $y$ の変動を説明したかを見る指標になる。
説明変数の個数を増やすことの問題点
説明変数を増やすと、モデルはより柔軟になる。最小二乗法では、変数を追加した後のモデルは追加前のモデルを特別な場合として含むので、訓練データ上の残差平方和 $\mathrm{RSS}$ は基本的に増えない。そのため、決定係数 $R^2$ も下がりにくい。
ただし、ここで大事なのは、$R^2$ が上がることと、よいモデルになることは同じではないという点である。ここでは「よい回帰モデル」を、以下の $2$ つの条件を満たすものとして考えよう。
- 係数推定量の分散 $\mathrm{Cov}(\hat{\boldsymbol{\beta}}) = \sigma^2(X^\top X)^{-1}$ が小さい。
- 分散推定 $\hat{\sigma}^2 = \dfrac{\mathrm{RSS}}{N-k-1}$ が不自然に大きくない。
前者は「係数推定がどれくらい安定しているか」を表し、後者は「自由度を考慮したときに、残差がどれくらい残っているか」を表す。変数を増やして $R^2$ だけが少し上がっても、係数が不安定になったり、自由度の割に残差が小さくなっていなかったりするなら、モデルとしては改善したとは言いにくい。
（1）多重共線性によって、係数の分散が大きくなる
まず問題になるのが、多重共線性である。これは、説明変数どうしが強く相関している状態をいう。たとえば、地域データで人口と世帯数を同時に入れる場合や、健康データで体重とBMIを同時に入れる場合が典型例である。どちらも似た情報を持つので、モデル全体としては $y$ を説明できているように見えても、「どちらの変数がどれだけ効いたのか」を切り分けにくい。
数式で見ると、この問題は $X^\top X$ の逆行列の発散に現れる。簡単のため、切片を省き、$2$ つの説明変数だけを考える。このとき計画行列は
\begin{align} X = \begin{pmatrix} x_{11} & x_{12}\\ x_{21} & x_{22}\\ \vdots & \vdots\\ x_{N1} & x_{N2} \end{pmatrix} = \begin{pmatrix} \boldsymbol{x}_1 & \boldsymbol{x}_2 \end{pmatrix} \end{align}である。ここで、説明変数は中心化されているとする。つまり、各列の標本平均について $\overline{x}_1 = \frac{1}{N}\sum_{i=1}^{N}x_{i1} =0, \overline{x}_2 = \frac{1}{N}\sum_{i=1}^{N}x_{i2} =0 $ である。このとき
であり、この行列の行列式は次のようになる。
この行列式は、 $x_1$ と $x_2$ の相関が強い場合 $0$ に近づく。よってこの場合は $(X^\top X)^{-1}$ （の各成分）が大きくなり、係数推定量 $\mathrm{Cov}(\hat{\boldsymbol{\beta}})=\sigma^2(X^\top X)^{-1}$ も大きくなる。つまり、個々の係数の推定値は大きく揺れやすくなる。
このとき起きているのは、似た説明変数の効果をデータから分離しにくいということである。たとえば人口と世帯数を同時に入れると、どちらも地域の規模を表しているため、片方の係数を大きくしてもう片方の係数を小さくしても、予測値があまり変わらない場合がある。その結果、モデル全体の当てはまりは悪くないのに、個々の係数の標準誤差が大きくなり、係数の解釈が不安定になる。
（2） 自由度が減り、見かけの当てはまりだけがよくなる
次に問題になるのが、自由度の減少である。説明変数を $1$ つ増やすと、$\mathrm{RSS}$ は下がりやすいが、分母 $N-k-1$ も小さくなる。
\begin{align} \hat{\sigma}^2 &= \frac{\mathrm{RSS}}{N-k-1} \end{align}したがって、$\mathrm{RSS}$ の減り方が十分でなければ、自由度で割った残差の大きさは、むしろ悪化することがある。簡単に言えば、変数を増やしたことでモデルの当てはまりがよくなっても、「その改善は、自由度を消費するだけの価値があったのか」という別のチェックが必要になる。
これは幾何的には、単回帰と重回帰の違いを考えると変わりやすい。単回帰では、点の散らばりに対して一本の「直線」を当てはめる。一方、説明変数が $2$ つになれば、点の散らばりに対して「平面」を当てはめることになる。平面のほうが直線より手元のデータに合わせやすい。しかし、その自由さが本当に意味のある構造を拾っているのかは別問題である。単に偶然の揺らぎを拾っているだけかもしれない。
この考え方を反映した指標が、自由度調整済み決定係数である。これはいわば、$\mathrm{TSS}, \mathrm{RSS}$ をそれぞれの自由度で割って比較するものである。
説明変数の数が異なるモデルを比較するときには、決定係数よりも、むしろ自由度調整済み決定係数を指標に用いるべきである。
要するに、説明変数を増やすときには、単に $R^2$ が上がったかどうかを見るのではなく、多重共線性によって係数が不安定になっていないか、自由度の減少に見合うだけ残差が減っているかを見る必要がある。説明変数は多ければよいわけではなく、目的に照らして入れる理由があるかを確認する必要がある。
赤池情報量規準（AIC）
「当てはまりの良さ」と「モデルの複雑さ」の両方のバランスを取りつつモデル選択する際の指標として、赤池情報量規準（AIC）がよく用いられる。最後にこれについて述べよう。
まずは、モデル精度について異なる考え方を提示しよう。これまで我々は、データ $(X, \boldsymbol{y})$ に対して残差平方和 $\mathrm{RSS}$ などを用いてモデルを評価してきた。しかし、これは機械学習的な立場で言えば、いわば「学習データ」でモデルの精度を評価しているものであり、適切ではないだろう。これを評価するためには、同じ $X$ に対して $2$ 種類のデータ
を考え、$\boldsymbol{y}$ を学習データ、$\tilde{\boldsymbol{y}}$ をテストデータとして扱うことにする。$\boldsymbol{u}$ と $\tilde{\boldsymbol{u}}$ はそれぞれ独立に、同じ確率分布
に従うものとする。データ $\boldsymbol{y}, \tilde{\boldsymbol{y}}$ は同じ確率分布に従うので、これらの確率密度関数を $f = f(\cdot)$ とする。
一方で、線形回帰による推定モデルは、学習データ $\boldsymbol{y}$ を用いて構築される。これは具体的には、これまでの議論に基づいて
\begin{align} \hat{\boldsymbol{y}} &\sim \mathcal{N}(X\hat{\boldsymbol{\beta}},\hat{\sigma}^2 I_N) \end{align}で想定され、この確率密度関数を $\hat f$ とかく。$\hat{\boldsymbol{\beta}}, \hat{\sigma}^2$ は $\boldsymbol{y}$ （つまり確率変数 $\boldsymbol{u}$ ）の関数であることに注意しよう。
さて、AICでは、$2$ つの分布 $f,\hat f$ のズレを測る指標として、Kullback-Leibler (KL) ダイバージェンスを用いる。これは以下で与えられる。
簡単に言えば、未知のテストデータに対して、推定モデル $\hat f$ がどれくらい高い対数尤度を与えられるかを見たいのである。さらに少し変形すると、
\begin{align} \mathrm{KL}(f, \hat f) &= E_{\tilde{\boldsymbol{y}}}[\log f(\tilde{\boldsymbol{y}})] - E_{\tilde{\boldsymbol{y}}}[\log \hat f(\tilde{\boldsymbol{y}})] \end{align}となる。第 $1$ 項は真の分布 $f$ だけで決まり、モデル比較では変わらない。したがって、候補モデルを比べるときに本質的に小さくしたいのは、第 $2$ 項の $-E_{\tilde{\boldsymbol{y}}}[\log \hat f(\tilde{\boldsymbol{y}})]$ である。対数尤度が大きいほど、この値は小さくなる。これを $2$ 倍して、$\boldsymbol{y}$ についても期待値を取ったものを $\mathrm{AI}$ と書く。
ここで
\begin{align} - 2 \log \hat f(\tilde{\boldsymbol{y}}) &= - 2 \log \left( \frac{1}{(2\pi \hat{\sigma}^2)^{N/2}} \exp\left( -\frac{1}{2\hat{\sigma}^2} ||\tilde{\boldsymbol{y}}-X\hat{\boldsymbol{\beta}}||^2 \right) \right) \\ &= N \log(2\pi \hat{\sigma}^2) + \frac{||\tilde{\boldsymbol{y}}-X\hat{\boldsymbol{\beta}}||^2}{\hat{\sigma}^2} \end{align}で、
\begin{align} \frac{||\tilde{\boldsymbol{y}}-X\hat{\boldsymbol{\beta}}||^2}{\hat{\sigma}^2} &= \frac{||X\boldsymbol{\beta} + \tilde{\boldsymbol{u}} - X\hat{\boldsymbol{\beta}}||^2}{\hat{\sigma}^2}\\ &= \frac{||X(\boldsymbol{\beta} - \hat{\boldsymbol{\beta}}) + \tilde{\boldsymbol{u}}||^2}{\hat{\sigma}^2}\\ &= \frac{||X(\boldsymbol{\beta} - \hat{\boldsymbol{\beta}})||^2}{\hat{\sigma}^2} + \frac{||\tilde{\boldsymbol{u}}||^2}{\hat{\sigma}^2} + \frac{2\tilde{\boldsymbol{u}}^\top X(\boldsymbol{\beta} - \hat{\boldsymbol{\beta}})}{\hat{\sigma}^2} \end{align}なので、代入して
\begin{align} \mathrm{AI} &= E_{\tilde{\boldsymbol{y}},\boldsymbol{y} }\left[ N \log(2\pi \hat{\sigma}^2) + \frac{||X(\boldsymbol{\beta} - \hat{\boldsymbol{\beta}})||^2}{\hat{\sigma}^2} + \frac{||\tilde{\boldsymbol{u}}||^2}{\hat{\sigma}^2} + \frac{2\tilde{\boldsymbol{u}}^\top X(\boldsymbol{\beta} - \hat{\boldsymbol{\beta}})}{\hat{\sigma}^2} \right] \end{align}である。そして、
\begin{gather} E_{\tilde{\boldsymbol{u}} }\left[ ||\tilde{\boldsymbol{u}}||^2 \right] = N\sigma^2, \quad E_{\tilde{\boldsymbol{u}} }\left[ \tilde{\boldsymbol{u}}^\top \right] = 0,\\ E_{\boldsymbol{u} }\left[ || X (\boldsymbol{\beta} - \hat{\boldsymbol{\beta}}) ||^2 \right] = (k+1)\sigma^2 \end{gather}および
\begin{align} E_{\boldsymbol{u} }\left[ \frac{1}{\hat{\sigma}^2} \right] &= \frac{N-k-1}{\sigma^2} E_{\boldsymbol{u} }\left[ \frac{1}{\chi^2_{N-k-1}} \right] \\ &=\frac{N-k-1}{\sigma^2(N-k-3)} \end{align}なので、以上をまとめて
となる。右辺第一項 $E_{ \boldsymbol{u} }\left[ N \log(2\pi \hat{\sigma}^2) \right]$ はこのままでは計算できないが、
\begin{align} E_{ \boldsymbol{u} }\left[ -2 \log \hat f(\boldsymbol{y}) \right] &= E_{ \boldsymbol{u} }\left[ N \log(2\pi \hat{\sigma}^2) + \frac{1}{\hat{\sigma}^2} \mathrm{RSS} \right]\\ &= E_{ \boldsymbol{u} }\left[ N \log(2\pi \hat{\sigma}^2) \right] + N-k-1 \end{align}となるので、これを逆に代入することで
となる。最後に、データ数 $N$ が十分多いときの近似式として
\begin{align} \frac{N-k-1}{N-k-3} &\approx 1 \end{align}を用いると、AI は $E_{ \boldsymbol{u} }\left[ -2 \log \hat f(\boldsymbol{y}) \right] + 2(k+2)$ で近似される。ここで、$E_{ \boldsymbol{u} }\left[ -2 \log \hat f(\boldsymbol{y}) \right]$ を、手元の $1$ 回の観測値 $-2 \log \hat f(\boldsymbol{y})$ で推定する。もちろん $1$ 回の観測値なのでばらつきは残るが、大標本のもとで平均的な予測の悪さを比較する基準として、
を用いる。これはもとを辿れば $f,\hat f$ のズレを測るものだった。つまり、$\mathrm{AI}$ は本来評価したい期待損失であり、$\mathrm{AIC}$ はそれを観測データから計算できる形に置き換えた基準である。よって、AICを最小にするモデルを選ぶことが望ましい。
なお、ここで第 $2$ 項の $(k+2)$ は、直感的にはモデル $\hat f$ のパラメータ数である。今回は $\hat \beta $ のパラメータ $k+1$ 個に加えて、$ \hat \sigma^2$ も推定するので、合計 $k+2$ 個のパラメータがある。
実際のモデル比較では、AIC の絶対値そのものよりも、候補モデル間の AIC の差を見る。同じデータに対して計算した AIC が小さいモデルほど、予測分布として相対的によいと判断する。重回帰分析の場合、確率密度関数は
\begin{align} \hat f(\boldsymbol{y}) &= (2\pi \hat{\sigma}^2)^{-\frac{N}{2}} \exp\left( -\frac{1}{2\hat{\sigma}^2} ||\boldsymbol{y}-X\hat{\boldsymbol{\beta}}||^2 \right) \end{align}であるから、AICは $(\ref{eq:aic})$ に代入して
\begin{align} \mathrm{AIC} &= N \log(2\pi \hat{\sigma}^2) + \frac{\mathrm{RSS}}{\hat{\sigma}^2} + 2(k+2) \end{align}である。最尤推定量では $\hat{\sigma}^2 = \mathrm{RSS}/N$ であるから、AICは
\begin{align} \mathrm{AIC} &= N \log\left(2\pi \frac{\mathrm{RSS}}{N}\right) + N + 2(k+2) \end{align}となる。あとは残差平方和 $\mathrm{RSS}$ を計算すれば良い。
10 実験計画法
実験では、ただデータをたくさん集めればよいわけではない。どの条件を試すか、どの順番で何回繰り返すかによって、得られる結論の信頼性は大きく変わる。実験計画法は、限られた実験回数の中で、知りたい効果をできるだけはっきり取り出すための統計的方法論の総称である。
実験計画法の「三種の神器」
実験計画法の「三種の神器」はあくまでここでの俗称であり、一般的ではない。ここでは「三種の神器」を、フィッシャーの3原則、分散分析、直交表の $3$ つと定める。実際はこれらの手法に限るものではないが、これら $3$ つは特に、実験の設計と分析においてきわめて重要な役割を果たす。
前提として、実験条件として変化させる項目を因子（または要因）と呼ぶ。たとえば温度、圧力、処理時間、材料の種類などが因子である。また、それぞれの因子が取りうる具体的な値や条件を水準と呼ぶ。たとえば温度という因子に対して、$100^\circ\mathrm{C}$ と $150^\circ\mathrm{C}$ の $2$ つで実験するなら、この因子は $2$ 水準である。
水準を定めて実験を行うと、興味ある特性値が実験結果として得られる。たとえば、製品の強度、収穫量、反応速度などが特性値である。実験計画法では、因子と水準を組み合わせて実験を行い、「どの因子が特性値に影響を与えるか」「その因子をどのような値（水準）に設定すれば、望ましい特性値が得られるか」を分析することが目的である。
フィッシャーの3原則
フィッシャーの3原則とは、反復、無作為化、局所管理の $3$ つである。これは、実験から得られる差が本当に処理の違いによるものなのか、それとも偶然や外部条件の違いによるものなのかを見分けるための基本原則である。
- 反復：同じ条件の実験を複数回行い、偶然誤差の大きさを見積もれるようにする。
- 無作為化：処理の割り当てをランダムにし、実験者の意図や隠れた条件差が特定の処理に偏らないようにする。
- 局所管理：似た条件の実験単位をまとめて比較し、場所・時間・個体差などの影響を小さくする。
反復がないと、誤差のばらつきを推定できない。無作為化がないと、処理効果と他の因子の影響が混ざる危険がある。局所管理がないと、本来は取り除けるはずのばらつきが誤差に入ってしまい、検出力が落ちる。
たとえば、肥料 $A$ と肥料 $B$ の効果を比較したいとする。畑の東側に肥料 $A$、西側に肥料 $B$ をまとめて撒くと、収穫量の差が肥料の違いなのか、日当たりや土壌の違いなのか分からない。そこで区画を細かく分け、各区画に肥料を無作為に割り当て、さらに似た条件の区画どうしをブロックとして扱う。このようにして、処理効果をできるだけ公平に比較する。
一元配置法
まずは因子が $1$ つだけの場合を考えよう。因子 $A$ について、複数の水準（または処理）$A_1, \cdots , A_I$ で実験を行ったとき、それぞれの水準 $A_i$ における特性値の母平均に有意差があるかどうかを検定しよう。各 $A_i$ の実験数を $J_i$ として、 $A_i$ における $j$ 番目の実験データを $x_{ij}$ とする。全データ数は $N= \sum_{i=1}^I J_i$ である。
| 因子 $A$ の水準 | 観測データ | 標本数 | 標本平均 |
|---|---|---|---|
| $A_1$ | $x_{11}\quad x_{12}\quad\cdots\quad x_{1J_1}$ | $J_1$ | $\overline{x}_{1\cdot}$ |
| $A_2$ | $x_{21}\quad x_{22}\quad\cdots\quad x_{2J_2}$ | $J_2$ | $\overline{x}_{2\cdot}$ |
| $\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ |
| $A_I$ | $x_{I1}\quad x_{I2}\quad\cdots\quad x_{IJ_I}$ | $J_I$ | $\overline{x}_{I\cdot}$ |
上の表で見るとイメージがしやすい。ここで、
\begin{align} \overline{x}_{i\cdot} &= \frac{1}{J_i}\sum_{j} x_{ij},\\ \overline{x}_{\cdot\cdot} &= \frac{1}{\sum_{i} J_i}\sum_{i}\sum_{j} x_{ij} \end{align}と書き、
と書く。$\mathrm{TSS}$ は総平方和（Total Sum of Squares）、$\mathrm{WSS}$ は群内平方和（Within-Group Sum of Squares）、$\mathrm{BSS}$ は群間平方和（Between-Group Sum of Squares）と呼ばれる。「群間」という名前がついているが、「水準間の違い」と考えたほうが分かりやすいだろう。
このとき、平方和は次のように分解できる。
実際、
\begin{align} x_{ij}-\overline{x}_{\cdot\cdot} &= (x_{ij}-\overline{x}_{i\cdot}) + (\overline{x}_{i\cdot}-\overline{x}_{\cdot\cdot}) \end{align}と分解して二乗和を取ると、交差項は
\begin{align} 2\sum_i\sum_j (x_{ij}-\overline{x}_{i\cdot}) (\overline{x}_{i\cdot}-\overline{x}_{\cdot\cdot}) &= 2\sum_i (\overline{x}_{i\cdot}-\overline{x}_{\cdot\cdot}) \sum_j(x_{ij}-\overline{x}_{i\cdot})\\ &=0 \end{align}となる。つまり、全体のばらつき $\mathrm{TSS}$ は、群内平方和 $\mathrm{WSS}$ と、群間平方和 $\mathrm{BSS}$ に分解される。
さて、ここからは統計モデルについて考えよう。実験計画法では通常、各水準ごとに母集団をそれぞれ設定する（今回は水準は $I$ 個なので、$I$ 個の母集団があると考える）。ここでは、実験値 $x_{ij}$ がすべて独立に、以下の確率モデルで記述できると仮定しよう。
ここで $\varepsilon_{ij}$ は誤差項である。ただし
\begin{align} \sum_i J_i\alpha_i=0 \end{align}として、 $\mu$ に全体平均、$\alpha_i$ に水準 $A_i$ の効果としての意味を持たせる。
（1） 分散分析
まずは、因子 $A$ が特性値 $x$ に影響するかを分散分析によって調べよう。分散分析と言っても、基本的な考え方は仮説検定と同じである。帰無仮説、対立仮説は以下のように表される。
帰無仮説のもとでは、水準間のばらつきも誤差によって生じたものにすぎない。そのため、
\begin{align} \frac{\mathrm{BSS}}{\sigma^2} \sim \chi^2_{I-1},\quad \frac{\mathrm{WSS}}{\sigma^2} \sim \chi^2_{N-I} \label{eq:F} \end{align}が成り立ち、さらにこれらは独立である（証明は後に行う）。したがって、
となる。水準間の平均差がない場合、$F$ は $1$ 付近である。一方で、水準間の平均差が大きいほど $\mathrm{BSS}$ が大きくなり、したがって $F$ も大きくなる。よって、有意水準 $\alpha$ の検定では、片側検定により
\begin{align} F \geq F_{I-1,N-I,\alpha} \end{align}のときに $H_0$ を棄却する。ただし $F_{I-1,N-I,\alpha}$ は、自由度 $(I-1,N-I)$ の $F$ 分布の上側 $\alpha$ 点である。
実際のデータを分析するときは、下記に示す分散分析表に計算結果を記入して整理することが多い。しかし今回の（一元配置）程度であれば、わざわざ表に整理するまでもないかもしれない。本当に表での整理が効いてくるのは、次節以降で述べる二元配置法からである。
| 変動因子 | 平方和 | 自由度 | 平均平方 | $F$ 値 |
|---|---|---|---|---|
| 水準間 | $\mathrm{BSS}$ | $I-1$ | $\mathrm{BSS}/(I-1)$ | $\dfrac{\mathrm{BSS}/(I-1)}{\mathrm{WSS}/(N-I)}$ |
| 水準内 | $\mathrm{WSS}$ | $N-I$ | $\mathrm{WSS}/(N-I)$ | |
| 全体 | $\mathrm{TSS}$ | $N-1$ |
（2） 母平均の推定
分散分析により影響があることが判明した場合は、母平均の点推定と信頼区間を求める。
まず、点推定値は
である。次に信頼区間を求めよう。まず、統計モデルより
\begin{align} \frac{ \sqrt{J_i}(\overline{x}_{i\cdot}-\mu-\alpha_i)}{\sigma} &\sim \mathcal{N}(0,1) \end{align}である。ここで $\sigma$ は未知なので、
から推定する。これより、$(\ref{eq:F})$ 式、および $\overline{x}_{i\cdot}$ と $\hat{\sigma}_W^2$ が独立であることから、
\begin{align} \frac{ \sqrt{J_i}(\overline{x}_{i\cdot}-\mu-\alpha_i)}{\hat \sigma_W} &\sim t_{N-I} \end{align}となる。したがって、水準 $A_i$ の母平均 $\mu+\alpha_i$ の $1-\alpha$ 信頼区間は
となる。ただし $t_{N-I,\alpha/2}$ は、自由度 $N-I$ の $t$ 分布の上側 $\alpha/2$ 点である。
（3） データの予測区間
最後に、データの予測区間を求める。水準 $A_i$ の新しい観測値 $x$ （表記を簡単にするため、添字は省略する）の予測区間を求めよう。
新しい観測値 $x$ と、すでに得られている標本平均 $\overline{x}_{i\cdot}$ との差を考える。新しい観測値に含まれる誤差と標本平均に含まれる誤差は独立なので、
\begin{align} x-\overline{x}_{i\cdot} &\sim \mathcal{N}\left(0,\sigma^2\left(1+\frac{1}{J_i}\right)\right) \end{align}となる。ここで、母平均の信頼区間を求めたときと同じく、$\sigma^2$ を $\hat\sigma_W^2=\mathrm{WSS}/(N-I)$ で推定する。$x-\overline{x}_{i\cdot}$ と $\hat\sigma_W^2$ は独立なので、
\begin{align} \frac{x-\overline{x}_{i\cdot}} {\hat\sigma_W\sqrt{1+1/J_i}} &\sim t_{N-I} \end{align}が成り立つ。したがって、水準 $A_i$ における新しい観測値 $x$ の $1-\alpha$ 予測区間は
となる。母平均の信頼区間と比べると、新しい観測値そのものの誤差を表す $1$ が平方根の中に加わるため、予測区間のほうが広くなる。
（4）補足：理論的背景について
基本的な作業の手順は以上の通りである。最後に分散分析の理論的背景について補足しておこう。具体的には、帰無仮説のもとで $F = \frac{\mathrm{BSS}/(I-1)}{\mathrm{WSS}/(N-I)}$ が $F$ 分布に従うことを示す。これは、$\mathrm{BSS}$ と $\mathrm{WSS}$ が独立であり、かつ帰無仮説のもとで $\mathrm{BSS}/\sigma^2$ と $\mathrm{WSS}/\sigma^2$ がそれぞれカイ二乗分布に従うことによって示される。順に示していこう。
まず、$\mathrm{BSS}$ と $\mathrm{WSS}$ の独立性を確認する。これは直感的には、「群平均を決める成分」と「群内残差を決める成分」が独立であることに伴う。具体的には、下記の式によって直接示される。
\begin{align} \operatorname{Cov} \left( \overline{x}_{i\cdot},\, x_{ij}-\overline{x}_{i\cdot} \right) &= \operatorname{Cov} \left( \overline{\varepsilon}_{i\cdot},\, \varepsilon_{ij}-\overline{\varepsilon}_{i\cdot} \right)\\ &= \operatorname{Cov} \left( \overline{\varepsilon}_{i\cdot},\varepsilon_{ij} \right) - \operatorname{Var} \left( \overline{\varepsilon}_{i\cdot} \right)\\ &= \frac{\sigma^2}{J_i} - \frac{\sigma^2}{J_i} = 0 \end{align}ただし $\overline{\varepsilon}_{i\cdot}=J_i^{-1}\sum_j \varepsilon_{ij}$ である。さらに、これらはいずれも、互いに独立な正規誤差 $(\varepsilon_{ij})_{i,j}$ の線形変換なので、まとめて多変量正規分布に従う。多変量正規分布では「無相関 $\Longrightarrow$ 独立」であるから、この二つのベクトルは独立である。
補足だが、同様の議論によって $\overline{x}_{i\cdot}$ と $\hat\sigma_W^2=\mathrm{WSS}/(N-I)$ が独立であることも簡単に示せる。
次に、帰無仮説のもとで $\mathrm{BSS}/\sigma^2$ と $\mathrm{WSS}/\sigma^2$ がそれぞれカイ二乗分布に従うことを示そう。ここで、全観測に対する誤差の平均を $\overline{\varepsilon}_{\cdot\cdot}=N^{-1}\sum_i\sum_j\varepsilon_{ij}=N^{-1}\sum_i J_i\overline{\varepsilon}_{i\cdot}$ と定義する。まず、$X_1, \ldots, X_n$ が独立な正規分布 $\mathcal{N}(\mu, \sigma^2)$ に従うとき、標本平均 $\overline{X}$ と不偏分散 $V^2$ が独立であり、かつ $\overline{X} \sim \mathcal{N}\left(\mu, \frac{\sigma^2}{n}\right)$ かつ $(n-1)V^2/\sigma^2 \sim \chi^2_{n-1}$ であることを思い出す。これにより、
\begin{align} \overline{\varepsilon}_{i\cdot} \sim \mathcal{N}\left(0, \frac{\sigma^2}{J_i}\right)&, \quad \frac{1}{\sigma^2}\sum_{i,j} (\varepsilon_{ij}-\overline{\varepsilon}_{i\cdot})^2 \sim \chi^2_{N-I},\\ \overline{\varepsilon}_{\cdot\cdot} \sim \mathcal{N}\left(0, \frac{\sigma^2}{N}\right)&, \quad \frac{1}{\sigma^2}\sum_{i} J_i (\overline{\varepsilon}_{i\cdot}-\overline{\varepsilon}_{\cdot\cdot})^2 \sim \chi^2_{I-1} \label{eq:BSS_chi2} \end{align}である。ただし \eqref{eq:BSS_chi2} の後半の式については、$J=J_1=\cdots=J_I$ のときのみ容易に示すことができる。一方で各 $J_i$ が異なる場合も、適切な線形変換を用いることで同様の結果を得ることができる。
次に、効果の加重平均を $\overline{\alpha}=N^{-1}\sum_i J_i\alpha_i$ と定義する。すでに課した制約 $\sum_i J_i\alpha_i=0$ より、$\overline{\alpha}=0$ である。この記号を用いると、
\begin{align} \mathrm{BSS} &= \sum_i J_i (\overline{x}_{i\cdot}-\overline{x}_{\cdot\cdot})^2\\ &= \sum_i J_i (\alpha_i - \overline{\alpha} + \overline{\varepsilon}_{i\cdot} - \overline{\varepsilon}_{\cdot\cdot})^2\\ &= \sum_i J_i (\alpha_i - \overline{\alpha})^2 + \sum_i J_i (\overline{\varepsilon}_{i\cdot} - \overline{\varepsilon}_{\cdot\cdot})^2 + 2 \sum_i J_i (\alpha_i - \overline{\alpha})(\overline{\varepsilon}_{i\cdot} - \overline{\varepsilon}_{\cdot\cdot}) \label{eq:BSS_decomposition} \end{align}であり、同様に
\begin{align} \mathrm{WSS} &= \sum_{i,j} (x_{ij}-\overline{x}_{i\cdot})^2 = \sum_{i,j} (\varepsilon_{ij} - \overline{\varepsilon}_{i\cdot})^2 \label{eq:WSS_decomposition} \end{align}である。これより、帰無仮説 $H_0: \alpha_1 = \cdots = \alpha_I$ の下で、独立に
であることが示される。上記２つを割ることにより、$F$ 検定統計量が得られる。
これ以降、簡単のため $J_1 = \cdots = J_I = J$ として議論を進める。
一元配置法（変量効果モデル）
前項述べた一元配置の固定効果モデルでは、$\alpha_i$ が母数（固定値）であることを仮定した。一方で、$\alpha_i$ を確率変数とみなすモデルを一元配置の変量効果モデルという。具体的には、
というモデルを考える。ただし今回の場合は $\alpha_i$ は平均 $0$ の正規分布に従うことから、 $\alpha_i$ の総和が $0$ であるという条件付けは不要である。これまで考えてきた帰無仮説 $H_0: \alpha_1 = \cdots = \alpha_I$ は、$\sigma_A^2 = 0$ に対応する。
帰無仮説のもとでは、今回の変量効果モデルは固定効果モデルと全く同じであり、分散分析の手順はまったく同様である。一方で、分散分析の結果 $\sigma_A^2$ が $0$ でないと認められる場合には、各 $\alpha_i$ ではなくその分散 $\sigma_A^2$ を推定することになる。
これは、具体的には $\mathrm{BSS}, \mathrm{WSS}$ によって表現できる。\eqref{eq:BSS_decomposition} および \eqref{eq:WSS_decomposition} 式により、
\begin{align} E[\mathrm{BSS}] &= (I-1)J E\left[\frac{1}{I-1} \sum_i (\alpha_i - \overline{\alpha}_{\cdot})^2 \right] + J \sum_i E\left[ (\overline{\varepsilon}_{i\cdot} - \overline{\varepsilon}_{\cdot\cdot})^2\right] = (I-1)(J \sigma_A^2 + \sigma^2),\\ E[\mathrm{WSS}] &= (N-I) \sigma^2. \end{align}である。これにより、$\sigma_A^2$ の不偏推定量として、
が得られる。以上の結果は、この後述べる乱塊法の理解の助けにもなるだろう。
二元配置法
次に、$2$ つの因子 $A,B$ がある場合を考える。因子 $A$ の水準を $A_1,\cdots,A_I$、因子 $B$ の水準を $B_1,\cdots,B_J$ とし、水準の組 $(A_i,B_j)$ における $k$ 番目の実験で得られた特性値を $x_{ijk}$ と書く。各組み合わせで $K$ 回ずつ観測するとして、全データ数は $N=IJK$ である。
下図のようにデータが揃っていると考えるとイメージしやすい。
| $B_1$ | $B_2$ | $\cdots$ | $B_J$ | |
|---|---|---|---|---|
| $A_1$ | $x_{111} \quad x_{112} \quad \cdots \quad x_{11K}$ | $x_{121} \quad x_{122} \quad \cdots \quad x_{12K}$ | $\cdots$ | $x_{1J1} \quad x_{1J2} \quad \cdots \quad x_{1JK}$ |
| $A_2$ | $x_{211} \quad x_{212} \quad \cdots \quad x_{21K}$ | $x_{221} \quad x_{222} \quad \cdots \quad x_{22K}$ | $\cdots$ | $x_{2J1} \quad x_{2J2} \quad \cdots \quad x_{2JK}$ |
| $\vdots$ | $\vdots$ | $\vdots$ | $\ddots$ | $\vdots$ |
| $A_I$ | $x_{I11} \quad x_{I12} \quad \cdots \quad x_{I1K}$ | $x_{I21} \quad x_{I22} \quad \cdots \quad x_{I2K}$ | $\cdots$ | $x_{IJ1} \quad x_{IJ2} \quad \cdots \quad x_{IJK}$ |
先程と同様に、平均を次のように書く。
\begin{align} \overline{x}_{ij\cdot} = \frac{1}{K}\sum_k x_{ijk},\quad \overline{x}_{i\cdot\cdot} = \frac{1}{JK}\sum_{j,k} x_{ijk},\\ \overline{x}_{\cdot j\cdot} = \frac{1}{IK}\sum_{i,k} x_{ijk}, \quad \overline{x}_{\cdot\cdot\cdot} = \frac{1}{IJK}\sum_{i,j,k} x_{ijk} \end{align}ここで、前節と同じように
\begin{align} \mathrm{TSS} &= \sum_{i,j,k} (x_{ijk}-\overline{x}_{\cdot\cdot\cdot})^2 \end{align}の平方和分解を考えよう。まず、前と同じように
\begin{align} \mathrm{TSS} &= \sum_{i,j,k} (x_{ijk}-\overline{x}_{ij\cdot})^2 + \sum_{i,j,k} (\overline{x}_{ij\cdot} - \overline{x}_{\cdot\cdot\cdot})^2\\ &=\mathrm{WSS} + \mathrm{BSS} \end{align}と分解できる。二元配置の場合、さらに $\mathrm{BSS}$ は、下記を代入することにより $3$ つの項に分解できる。
実際に計算すると、
\begin{align} \mathrm{BSS} &= \sum_{i,j,k} \left\{ (\overline{x}_{i\cdot\cdot}-\overline{x}_{\cdot\cdot\cdot}) + (\overline{x}_{\cdot j\cdot}-\overline{x}_{\cdot\cdot\cdot}) + (x_{ij\cdot}-\overline{x}_{i\cdot\cdot}-\overline{x}_{\cdot j\cdot}+\overline{x}_{\cdot\cdot\cdot}) \right\}^2\\ &=\sum_{i} JK (\overline{x}_{i\cdot\cdot}-\overline{x}_{\cdot\cdot\cdot})^2 + \sum_{j} IK (\overline{x}_{\cdot j\cdot}-\overline{x}_{\cdot\cdot\cdot})^2 + \sum_{i,j} K (\overline{x}_{ij\cdot}-\overline{x}_{i\cdot\cdot}-\overline{x}_{\cdot j\cdot}+\overline{x}_{\cdot\cdot\cdot})^2 \end{align}となる。なお、上の計算において、交差項の和がすべて $0$ になることを用いた。このことは、例えば以下のように簡単に確認できる。
\begin{align} & \sum_{i,j,k} (\overline{x}_{i\cdot\cdot}-\overline{x}_{\cdot\cdot\cdot}) (x_{ij\cdot}-\overline{x}_{i\cdot\cdot}-\overline{x}_{\cdot j\cdot}+\overline{x}_{\cdot\cdot\cdot})\\ =& \sum_{i,k} J (\overline{x}_{i\cdot\cdot}-\overline{x}_{\cdot\cdot\cdot}) (x_{i\cdot\cdot}-\overline{x}_{i\cdot\cdot}-\overline{x}_{\cdot \cdot\cdot}+\overline{x}_{\cdot\cdot\cdot}) =0 \end{align}以上の議論をまとめよう。平方和 $\mathrm{TSS}$ は
\begin{align} \mathrm{SSA} &= JK\sum_i (\overline{x}_{i\cdot\cdot}-\overline{x}_{\cdot\cdot\cdot})^2,\\ \mathrm{SSB} &= IK\sum_j (\overline{x}_{\cdot j\cdot}-\overline{x}_{\cdot\cdot\cdot})^2,\\ \mathrm{SSAB} &= K\sum_i\sum_j (\overline{x}_{ij\cdot}-\overline{x}_{i\cdot\cdot}-\overline{x}_{\cdot j\cdot}+\overline{x}_{\cdot\cdot\cdot})^2,\\ \mathrm{SSE} &=\mathrm{WSS}= \sum_i\sum_j\sum_k (x_{ijk}-\overline{x}_{ij\cdot})^2 \end{align}を用いて
と分解される。ここで、$\mathrm{SSA}$ は因子 $A$ の主効果、$\mathrm{SSB}$ は因子 $B$ の主効果、$\mathrm{SSAB}$ は交互作用、$\mathrm{SSE}$ は誤差に対応する平方和と解釈できる。なお、実際のデータを計算する場合は、$\mathrm{SSAB}$ については一度 $\mathrm{BSS}$ を計算したあと、$\mathrm{SSAB} = \mathrm{BSS} - \mathrm{SSA} - \mathrm{SSB}$ で計算してしまったほうが早い。
さて、確率モデルを考えよう。二元配置の分散分析では、水準の組み合わせ $A_iB_j$ ごとに母集団を想定する（全部で $IJ$ 個の母集団を想定）。そして下記のモデルを考える。
ここで $\alpha_i$ は因子 $A$ の効果、$\beta_j$ は因子 $B$ の効果、$(\alpha\beta)_{ij}$ は交互作用である。交互作用とは、いわば因子 $A$ と因子 $B$ の相乗効果である。制約条件としては
\begin{align} \sum_i J \alpha_i = 0,\quad \sum_j \beta_j = 0,\quad \sum_i (\alpha\beta)_{ij} = 0,\quad \sum_j (\alpha\beta)_{ij} = 0 \end{align}を課す。このようにおけば、例えば $(\ref{eq:two-way-model})$ 式において両辺を $i$ で和を取ると、
\begin{align} \frac{1}{I}\sum_i x_{ijk} &\approx \mu + \beta_j \end{align}となる。このようにして、$\beta_j$ は、因子 $A$ の効果を平均化したときの因子 $B$ の効果として解釈できる。
（1） 相互作用の有無の検定
まず、分散分析では興味ある検定が $3$ つあることに注意しよう。それぞれの検定の帰無仮説を以下に述べる（対立仮説は、母数空間において帰無仮説の主張を反転して作成すればよい）。
ここで重要なのは、相互作用が認められる場合とない場合で、検定・推定の手順が異なるということである。つまり、帰無仮説が $H_0^{AB}$ の検定結果によって、次のステップでの作業が異なるということである。
このように説明すると、「まずは早く、相互作用を検定しよう！」となるかもしれない。しかし、まずは焦らず、得られた実験データに基づいて、下記の分散分析表を作成しよう。この表はこの先のステップでも活用することになる。
| 因子 | 平方和 | 自由度 | 分散 | $F$ 値 |
|---|---|---|---|---|
| $A$ | $\mathrm{SSA}$ | $\nu_A=I-1$ | $V_A^2=\mathrm{SSA}/\nu_A$ | $F_A = V_A^2/V_E^2$ |
| $B$ | $\mathrm{SSB}$ | $\nu_B=J-1$ | $V_B^2=\mathrm{SSB}/\nu_B$ | $F_B = V_B^2/V_E^2$ |
| $AB$ | $\mathrm{SSAB}$ | $\nu_{AB}=(I-1)(J-1)$ | $V_{AB}^2=\mathrm{SSAB}/\nu_{AB}$ | $F_{AB} = V_{AB}^2/V_E^2$ |
| 誤差 | $\mathrm{SSE}$ | $\nu_E=IJ(K-1)$ | $V_E^2=\mathrm{SSE}/\nu_E$ | |
| 全体 | $\mathrm{TSS}$ | $N-1$ |
相互作用の有無の判断方法は、基本的には一元配置のときと同様である。具体的には、上の分散分析表において、$F_{AB} = V_{AB}^2/V_E^2$ を検定統計量として、自由度 $((I-1)(J-1),IJ(K-1))$ の $F$ 分布の上側 $\alpha$ 点と比較すればよい。そこで $F$ 値が上側 $\alpha$ 点より大きい場合は、相互作用があると判断し、このまま次節へ文章を読み進めよう。有意でない場合も相互作用がないとは断定できないが、実験の背景から相互作用を無視できると判断する場合は、次々節のステップへと進む。
（2） 相互作用がある場合
相互作用がある場合の検定・推定は、基本的に一元配置の分散分析と同じである。まず分散分析については、これまでと同様に帰無仮説が $H_0^A, H_0^B$ で考えられるそれぞれの検定に対して、検定統計量として $F_A, F_B$ を用いて $F$ 検定を行えばよい。帰無仮説のもとでは、これらはいずれも自由度 $(\nu_A,\nu_E), (\nu_B,\nu_E)$ の $F$ 分布に従う。
また推定については、水準の組 $A_iB_j$ を、一元配置における一つの水準だと思ってしまえばよい。例えば、母数の推定については
とすればよい。また、母数の信頼区間についても下記のようにおけばよい。
（3） 相互作用がない場合
相互作用を $0$ と仮定する場合、相互作用に対応する平方和を誤差へプールする（つけかえる）ことができる。すなわち、平方和 $\mathrm{SSAB}$ を、誤差の平方和 $\mathrm{SSE}$ に加えて、分散分析を再計算する。この作業はプーリングと呼ばれる。ただし、検定で有意でないことだけを理由にプーリングを行うと、モデル選択の影響により、その後の検定の有意水準や信頼区間の被覆確率が名目通りになるとは限らない。以下の分布と信頼区間は、相互作用を $0$ としたモデルのもとでの結果である。
具体的には、前のステップで作った分散分析表において、「因子 $AB$」の行を、「誤差」の行に足し上げる。
| 因子 | 平方和 | 自由度 | 分散 | $F$ 値 |
|---|---|---|---|---|
| $A$ | $\mathrm{SSA}$ | $\nu_A=I-1$ | $V_A^2=\mathrm{SSA}/\nu_A$ | $F_A = V_A^2/V_E^2$ |
| $B$ | $\mathrm{SSB}$ | $\nu_B=J-1$ | $V_B^2=\mathrm{SSB}/\nu_B$ | $F_B = V_B^2/V_E^2$ |
| $AB$ | $\mathrm{SSAB}$ | $\nu_{AB}=(I-1)(J-1)$ | $V_{AB}^2=\mathrm{SSAB}/\nu_{AB}$ | $F_{AB} = V_{AB}^2/V_E^2$ |
| 誤差 | $\mathrm{SSE}$ | $\nu_E=IJ(K-1)$ | $V_E^2=\mathrm{SSE}/\nu_E$ | |
| 全体 | $\mathrm{TSS}$ | $N-1$ |
上の赤枠内の「平方和」「自由度」列を足し上げて、「誤差」として一行にまとめ上げよう。
| 因子 | 平方和 | 自由度 | 分散 | $F$ 値 |
|---|---|---|---|---|
| $A$ | $\mathrm{SSA}$ | $\nu_A=I-1$ | $V_A^2=\mathrm{SSA}/\nu_A$ | $F_A = V_A^2/V_{E'}^2$ |
| $B$ | $\mathrm{SSB}$ | $\nu_B=J-1$ | $V_B^2=\mathrm{SSB}/\nu_B$ | $F_B = V_B^2/V_{E'}^2$ |
| 誤差 | $\mathrm{SSE}'=\mathrm{SSAB}+\mathrm{SSE}$ | $\nu_{E'}=\nu_{AB}+\nu_{E}$ | $V_{E'}^2=\mathrm{SSE}'/\nu_{E'}$ | |
| 全体 | $\mathrm{TSS}$ | $N-1$ |
誤差分散 $V_{E'}^2$ 、誤差自由度 $\nu_{E'}$ が全面的に書き換わったことに注意しよう。このようにして新しく作成した分散分析表を用いて、$A,B$ の主効果を検定する。より具体的には、新しく書き換えた分散分析表を用いて、検定統計量として $F_A, F_B$ を用いて $F$ 検定を行えばよい。帰無仮説のもとでは、これらはいずれも自由度 $(\nu_A,\nu_{E'}), (\nu_B,\nu_{E'})$ の $F$ 分布に従う。
推定についても、相互作用がある場合と大きく異なる。相互作用がある場合は $\overline{x}_{ij\cdot}$ によって推定を行ったが、今回の場合は
にて推定を行うべきである。イメージとしては、推定 $\widehat{\mu+\alpha_i}$ と $\widehat{\mu+\beta_j}$ をそれぞれ行い、そこから推定 $\widehat{\mu}$ を引くというイメージである。
このようにして一見少々ややこしい推定を行う理由は、こちらのほうがより信頼区間が狭い（ブレが小さい）からである。上のように推定した場合の信頼区間は、
である。ここで $n_e$ は有効反復数であり、下記で与えられる。
\begin{align} n_e = \frac{N}{I+J-1} \end{align}ただし統計検定の試験では、$n_e$ の値まで暗記を要求されることはないであろう。
（4） 繰り返しがない場合の注意
以上が二元配置法についての基本的な説明になるが、ここで繰り返しがない場合について補足しておこう。
ここで「繰り返しがない」というのは、$K=1$ の場合である。ただし注意すべきは、「水準設定を含めた実験自体の繰り返し」のみが、繰り返しとして認められるということである。例えば因子として「炉の温度」を選んだ場合、一回の実験が終わったら炉の温度をベース温度（室温）まで戻して、再び水準まで設定して実験を行うことで、初めて水準設定の繰り返しが行われる。同じ水準設定のまま何回も測定だけ行うのは、水準設定の繰り返しとしては認められない。
さて、$K=1$ の場合、自動的に
\begin{align} \mathrm{SSE} = 0 \end{align}となることに注意しよう。これは病的な状態に思えるが、実際に意味するのは「繰り返しがない場合は、相互作用の有無が検定・推定できない」ということである。というのも、$A,B$ がともに異なる水準間で誤差が生じた場合、それが誤差による影響なのか、あるいは相互作用による影響なのか、全く分からないからである。より数学的に言えば、
\begin{align} \mathrm{SSAB} = \sum_{i,j} K (\overline{x}_{ij\cdot}-\overline{x}_{i\cdot\cdot}-\overline{x}_{\cdot j\cdot}+\overline{x}_{\cdot\cdot\cdot})^2 \end{align}を、相互作用とみなすこともできるし、誤差とみなすこともできる。つまり、相互作用があるのかないのか、全く判断できないのである。よって大前提として、「繰り返しのない二元配置法は、相互作用が存在しないことがはっきりしている場合にのみ用いる」ことが原則であるといえよう。
このような場合は、初めから相互作用がないことが前提とされるため、前節で述べた「相互作用がない場合」の検定のステップに従って検定・推定を進めればよい。
乱塊法
乱塊法は、実験の条件を揃えることが難しい場合に用いられる。例えば、因子 $A$ において毎日 $4$ 水準の実験を行う際に、実験日によって天候や気温、作業担当者が異なることにより実験結果にばらつきが生じる場合がある。このような場合、実験日 $j$ を新しい因子 $R$ として扱い、各水準 $R_j$ のもとで因子 $A$ のすべての水準の実験を（ランダムな順序で）行おうと思うのは自然なことだろう。
この場合、形式的にはこれは、繰り返しのない二元配置法に等しい。ここで「繰り返しのない」と述べたのは、因子 $R$ の水準設定には再現性がない（実験者が水準を意図的に設定できない）ことを意味している。このように、「誤差からその効果を取り除く」ことを目的として導入される因子はブロック因子と呼ばれる。
考えるモデルは下記である。
ここで $r_j$ は変量効果として導入されていることに注意しよう。これは、ブロック因子 $R$ の水準に再現性はなく、ランダムに $r_j$ が選ばれたものとみなすことを意味する。一方で、因子 $A$ は水準設定の再現性がある（実験者が水準を意図的に設定できる）ことから、制御因子として（つまり $\alpha_i$ は固定値として）扱われる。
なお、各ブロック $R_j$ 内では制御因子 $A$ の水準が一通り全部実験でき、かつ各ブロック内は均質な条件とみなせるときのみ、乱塊法が適用できることに注意しよう。また、各ブロック $R_j$ 内ではランダムな順序で実験が行わねばならない。
（1） 分散分析
変量効果モデルのところで述べたように、分散分析については固定効果モデルのときと同じである。よって、前節の二元配置法（相互作用なし）と同じように行えばよい。
（2） ブロック因子が無視できる場合
ブロック因子 $R$ の分散 $\sigma_R^2$ を $0$ と仮定できる場合、下記のモデルに基づいて母数の推定を行う。ただし、ブロック効果の検定で有意でないことだけから $\sigma_R^2=0$ とは断定できない。検定結果だけでブロック効果を除くモデルを選ぶと、その後の信頼区間の被覆確率が名目通りになるとは限らない。
これは一元配置法と同じであり、通常の一元配置法の手順に従って行えばよい。
（3） ブロック因子が無視できない場合
分散分析の結果、ブロック因子 $R$ が無視できない場合、点推定値 $\widehat{\mu+\alpha_i}$ については、これまでと同様に
から推定すればよい。しかし、区間推定がこれまでの方法と異なってくる。信頼区間については、
となる。ここで $s_i^2$ は、水準 $A_i$ の観測値をブロック間で比較したときの標本分散であり、
\begin{align} s_i^2 &= \frac{1}{J-1}\sum_{j=1}^J (x_{ij}-\overline{x}_{i\cdot})^2 \end{align}である。ただし、$J\geq 2$ とし、ブロック効果 $r_j$ と誤差 $\epsilon_{ij}$ は互いに独立であるとする。先程までの誤差分散だけを用いた信頼区間と異なるのは、母平均 $\mu+\alpha_i$ の推定には、ブロック間のばらつきも含まれるからである。
このことを確率モデルから確認しよう。水準 $A_i$ を固定すると、異なるブロックの観測値は独立に
\begin{align} x_{ij} &\sim \mathcal{N}(\mu+\alpha_i,\sigma_R^2+\sigma^2) \end{align}に従う。したがって、正規標本の平均と標本分散の性質より、
\begin{align} \frac{\sqrt{J}(\overline{x}_{i\cdot}-\mu-\alpha_i)} {\sqrt{\sigma_R^2+\sigma^2}} &\sim \mathcal{N}(0,1),\\ \frac{(J-1)s_i^2}{\sigma_R^2+\sigma^2} &\sim \chi^2_{J-1} \end{align}が成り立ち、さらに両者は独立である。よって、これらの比を取ることで
\begin{align} \frac{\overline{x}_{i\cdot}-\mu-\alpha_i} {\sqrt{s_i^2/J}} &\sim t_{J-1} \end{align}となり、冒頭の $1-\alpha$ 信頼区間が得られる。ただし $t_{J-1,\alpha/2}$ は、自由度 $J-1$ の $t$ 分布の上側 $\alpha/2$ 点である。この方法は、水準 $A_i$ のデータだけを用いて、$\sigma_R^2+\sigma^2$ をまとめて推定する方法である。
直交表実験
最初に実験を始めたての段階では、考えられる因子（温度、圧力、材料、など）が数多くあり、まずは主要な（注目すべき）因子を絞り込みたいことが多い。このような場合、まずは考えうる多くの因子を $2$ 水準で実験し、主要な因子を検定で数個に絞り込みたいと考えるのが自然だろう。
たとえば、調べたい因子が $7$ つある場合、全ての組み合わせで $2$ 水準実験を行うと、最低でも $2^7=128$ 回の実験が必要となる。このように、これまで述べてきた要因配置実験では、因子数が増えると実験回数は指数関数的に増加する。そこで、なるべく少ない実験回数で、主要な因子を最大限に特定できるよう実験を工夫したい。
ちなみに、実験回数を $2^7=128$ 回より減らしたいと考えることには、一定の理論的な正当性もある。一般に、$3$ 因子以上の相互作用の影響は実際には無視できることが多い。そのため、$2$ 因子までの相互作用を全て勘定した場合でも、自由度は
\begin{align} 7 + \frac{7(7-1)}{2} = 7 + 21 = 28 \end{align}程度である。これは、$128$ 回よりも少ない実験回数で、主要な因子を特定できる可能性があることを示唆している。
このような場合には、直交表を用いた実験方法・回数の工夫（直交表実験）が有効である。本節では、直交表を使った具体的な実験計画および理論的な背景について説明しよう。
なお、以下では因子間の相互作用は（ $2$ 因子を含め）すべて無視できるものとする。実のところ、今回述べる方法を少し改良すれば、 $2$ 因子間相互作用がある場合でも分析可能である（ただし、この場合の解説は専門的な教科書に説明を委ねる）。
| 実験 | [1] | [2] | [3] | [4] | [5] | [6] | [7] |
|---|---|---|---|---|---|---|---|
| 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 |
| 2 | 1 | 1 | 1 | 2 | 2 | 2 | 2 |
| 3 | 1 | 2 | 2 | 1 | 1 | 2 | 2 |
| 4 | 1 | 2 | 2 | 2 | 2 | 1 | 1 |
| 5 | 2 | 1 | 2 | 1 | 2 | 1 | 2 |
| 6 | 2 | 1 | 2 | 2 | 1 | 2 | 1 |
| 7 | 2 | 2 | 1 | 1 | 2 | 2 | 1 |
| 8 | 2 | 2 | 1 | 2 | 1 | 1 | 2 |
以上は $L_8=L_8(2^7)$ 直交表であり、これを例として説明を続ける。
（0） 導入（直交表の見方）
まず、直交表の各行（ $1 \sim 8$ ）は実験番号であり、各列（ $[1] \sim [7]$ ）は因子を表す。つまり、今回の $L_8$ 直交表を用いると「 $7$ 個未満の因子数を $8$ 回の実験で判定可能」というわけである。これより多くの因子数で実験をしたい場合は、より大きな $L_{16}$ や $L_{32}$ などの直交表を用いればよい。
ここで、この方法で検定できる因子数は、直交表に用意された列数未満でなくてはならない点に注意しよう。というのも、少なくとも一つの因子（例えば因子 $[7]$ ）は、「誤差」の自由度を表すために残す必要があるからである。たとえば今回は列数が $7$ なので、検定できる因子数は $6$ まで（$A,B,C,D,F,G$）である。なお、検定したい因子が $6$ 未満であることは問題ない。その場合は、空けた列（つまり誤差自由度）が $2$ 以上になる。
また、各行の実験 $i$ において、$1$（または $2$ ）の数字が羅列してあるが、これは「因子 $[k]$ については、水準 $1$（または $2$ ）で実験を行いなさい」ということを示している。例えば実験 $3$ では、因子 $[1], [2], [3], [4], [5], [6]$ について、それぞれ水準 $1,2,2,1,1,2$ で実験を行うことになる（ $[7]$ は誤差項として空けておく）。
考える確率モデルは、これまでの要因配置実験と同様である。具体的には、全体平均を $\mu$ 、因子 $[1], \cdots, [6]$ の主効果をそれぞれ $\tau_{[1]}, \cdots, \tau_{[6]}$ で表し、
\begin{align} x_{i} = \mu \pm \tau_{[1]} \pm \tau_{[2]} \pm \tau_{[3]} \pm \tau_{[4]} \pm \tau_{[5]} \pm \tau_{[6]} + \varepsilon_i, \qquad \varepsilon_{i} \overset{\text{i.i.d.}}{\sim}\mathcal{N}(0,\sigma^2) \end{align}とかく。ただし、$\pm$ の符号は直交表において水準 $1$ が $+$、水準 $2$ が $-$ に対応するものとする。なお、以下の説明では因子 $[7]$ を誤差項として空けておくことにする。
たとえば、実験 $3$ における測定値は、$L_8$ 直交表より
\begin{align} x_{3} = \mu + \tau_{[1]} - \tau_{[2]} - \tau_{[3]} + \tau_{[4]} + \tau_{[5]} - \tau_{[6]} + \varepsilon_3, \qquad \varepsilon_{3}\sim\mathcal{N}(0,\sigma^2) \label{eq:experiment3} \end{align}と書けば良い。
（1） 直交表の数学的意味
直交表についての理論的説明は省略されることが多いが、これは直交法の理解に大変重要である。ここでは $L_8$ 直交表を用いて説明するが、その他の $2$ 水準直交表への一般化は難しくないだろう。以下では実験回数を $n=8$ と書く。
まずは、$L_8$ 直交表の $k$ 列（ $[k]$ 因子の列）において、成分の $2$ を $-1$ に置き換えたベクトルをそれぞれ $\ \boldsymbol{\ell}_k$（ $k=1, \cdots, 7$ ）とかく。たとえば、
\begin{align} \boldsymbol{\ell}_1 = \begin{pmatrix} 1 \\ 1 \\ 1 \\ 1 \\ -1 \\ -1 \\ -1 \\ -1 \end{pmatrix}, \quad \boldsymbol{\ell}_3 = \begin{pmatrix} 1 \\ 1 \\ -1 \\ -1 \\ -1 \\ -1 \\ 1 \\ 1 \end{pmatrix} \end{align}である。このとき $\boldsymbol{\ell}_k \in \mathbb{R}^8$ であり、各成分は $\pm 1$ のいずれかである。各列では $1$ と $-1$ が同数ずつ現れ、異なる列では符号の組み合わせが均等に現れるため、
\begin{align} \boldsymbol{1}^\top\boldsymbol{\ell}_k=0,\qquad \boldsymbol{\ell}_{k}^\top \boldsymbol{\ell}_{k'} = \begin{cases} n & (k = k'), \\ 0 & (k \neq k') \end{cases} \end{align}が成り立つ。さらに $\boldsymbol{\ell}_0=\boldsymbol{1}$ とおくと、$\boldsymbol{\ell}_0,\boldsymbol{\ell}_1,\ldots,\boldsymbol{\ell}_7$ は $\mathbb{R}^8$ の互いに直交する基底になる。したがって、これらを長さ $1$ に正規化した $\boldsymbol{\ell}_k/\sqrt{n}$ は正規直交基底であり、次の射影の恒等式が成り立つ。
よって、任意の $\boldsymbol{y}\in\mathbb{R}^8$ は各直交列への射影の和として
\begin{align} \boldsymbol{y} &= \frac{\boldsymbol{1}^\top\boldsymbol{y}}{n}\boldsymbol{1} + \sum_{k=1}^7 \frac{\boldsymbol{\ell}_k^\top\boldsymbol{y}}{n} \boldsymbol{\ell}_k\\ &= \overline{y} \boldsymbol{1} + \sum_{k=1}^7 \frac{\boldsymbol{\ell}_k^\top\boldsymbol{y}}{n} \boldsymbol{\ell}_k\\ \end{align}と一意に分解できる。ただし最後の等号では $\overline{y} = \frac{1}{n}\boldsymbol{1}^\top\boldsymbol{y}$ を用いた。この式は、全体平均に対応する方向と、直交表の各列に対応する方向が、データの変動を重複なく受け持つことを表している。
また、誤差ベクトルが $\boldsymbol{\epsilon}\sim\mathcal{N}(\boldsymbol{0},\sigma^2\boldsymbol{I}_8)$ に従うとき、異なる列への射影は独立である。というのも、$k\neq k'$ ならば
\begin{align} \mathrm{Cov} \left( \boldsymbol{\ell}_k^\top\boldsymbol{\epsilon}, \boldsymbol{\ell}_{k'}^\top\boldsymbol{\epsilon} \right) &= \sigma^2\boldsymbol{\ell}_k^\top\boldsymbol{\ell}_{k'} =0 \end{align}であり、これらの射影は同時正規分布に従うので、「無相関 $\Longrightarrow$ 独立性」が従うためである。この直交分解と独立性が、以下の平方和分解と分散分析の基礎になる。
（2） 点推定
ここで、実験 $i$ （ $i=1, \cdots, 8$ ）に対応する測定値を $x_i$ として、測定ベクトル $\boldsymbol{x}$ を $\boldsymbol{x} = \begin{pmatrix} x_1 & x_2 & \cdots & x_8 \end{pmatrix}^\top \in \mathbb{R}^8$ と定める。 $x_3$ の具体的表式については $(\ref{eq:experiment3})$ 式で説明したが、より一般に $\boldsymbol{x}$ は
とまとめられる。これに例えば、左から $\boldsymbol{\ell}_1^\top$ を掛けると、直交性により
\begin{align} \boldsymbol{\ell}_1^\top \boldsymbol{x} = n \tau_{[1]} + \boldsymbol{\ell}_1^\top \boldsymbol{\epsilon} \end{align}となる。これはそのまま $\tau_{[1]} = \alpha_{[1],1}$ の点推定になっていることに注意しよう。実際、点推定量は次のように書ける。
\begin{align} \hat{\tau}_{[1]} = \frac{1}{n} \boldsymbol{\ell}_1^\top \boldsymbol{x} \end{align}ここで、因子 $[k]$ が水準 $i$ である実験値の合計を $T_{[k],i}$ と書くことにする。たとえば、
\begin{align} T_{[1],1} &= x_1 + x_2 + x_3 + x_4,\\ \qquad T_{[1],2} &= x_5 + x_6 + x_7 + x_8 \end{align}である。$\boldsymbol{\ell}_k$ の成分は、水準 $1$ なら $1$ 、水準 $2$ なら $-1$ なので、$\boldsymbol{\ell}_k^\top \boldsymbol{x}$ は結局、水準 $1$ における合計と水準 $2$ における合計の差を計算している。すなわち、
である。これを代入すると、点推定値は（一般の $k$ について）
である。$L_8$ 直交表では $n=8$ なので、水準ごとの合計差を $8$ で割ることになる。直交表を用いると、他の因子の水準が均等に混ざるため、このような単純な差であっても、因子 $[k]$ の主効果だけをうまく取り出すことができる。
（3） 直交表における平方和の分解
次に、全体のばらつきが各列の平方和へ分解されることを確認しよう。まず、データから全体平均を取り除いた中心化ベクトルを
\begin{align} \boldsymbol{x}_c &= \boldsymbol{x}-\overline{x}\boldsymbol{1}, \qquad \overline{x}=\frac{1}{n}\boldsymbol{1}^\top\boldsymbol{x} \end{align}とおく。$\boldsymbol{1}^\top\boldsymbol{x}_c=0$ であり、$\boldsymbol{\ell}_k^\top\boldsymbol{1}=0$ なので、$\boldsymbol{\ell}_k^\top\boldsymbol{x}_c=\boldsymbol{\ell}_k^\top\boldsymbol{x}$ である。前節の射影の恒等式を $\boldsymbol{x}_c$ に用いると、全体平均の方向 $\boldsymbol{\ell}_0=\boldsymbol{1}$ に対応する項は消えて、
\begin{align} \boldsymbol{x}_c &= \sum_{k=1}^7 \frac{(\boldsymbol{\ell}_k^\top\boldsymbol{x})}{n} \boldsymbol{\ell}_k \end{align}となる。したがって、全平方和は $\mathrm{SST}=\lVert\boldsymbol{x}_c\rVert^2$ と書けること、および直交性 $\boldsymbol{\ell}_k^\top\boldsymbol{\ell}_{k'}=0$（$k\neq k'$）を使うと、
\begin{align} \mathrm{SST} &= \left\lVert \sum_{k=1}^7 \frac{\boldsymbol{\ell}_k^\top\boldsymbol{x}}{n} \boldsymbol{\ell}_k \right\rVert^2 = \sum_{k=1}^7\sum_{k'=1}^7 \frac{(\boldsymbol{\ell}_k^\top\boldsymbol{x}) (\boldsymbol{\ell}_{k'}^\top\boldsymbol{x})}{n^2} \boldsymbol{\ell}_k^\top\boldsymbol{\ell}_{k'}\\ &= \sum_{k=1}^7 \frac{(\boldsymbol{\ell}_k^\top\boldsymbol{x})^2}{n} \end{align}となる。途中の交差項がすべて $0$ になるため、各列に対応する変動を重複なく足し合わせることができる。そこで、列 $[k]$ に対応する平方和を
と定義すれば、全平方和は次のように分解される。
各列は一次元の方向への射影なので、各 $\mathrm{SS}_{[k]}$ の自由度は $1$ である。今回のように列 $[1],\ldots,[6]$ を因子に割り当て、列 $[7]$ を空けておく場合は、残った列が表す変動を誤差として
\begin{align} \mathrm{SSE}=\mathrm{SS}_{[7]}, \qquad \nu_E=1 \end{align}とおく。より一般に、因子を割り当てず誤差に用いる列の集合を $\mathcal{E}$ とすれば、直交する各列の平方和と自由度を足し合わせて
とする。誤差列を多く残すほど誤差自由度が増えるため、後述する誤差分散の推定は安定する。
（4） 分散分析による検定
直交表実験でも、因子の影響があるかどうかは分散分析によって検定する。因子 $[k]$ について調べたい帰無仮説と対立仮説は
である。
ここで、各列に現れる確率変動を調べよう。ここでも出発点となるのは、これまでの議論でも大活躍した下記の式
である。また、$\boldsymbol{\epsilon}\sim\mathcal{N}(\boldsymbol{0},\sigma^2\boldsymbol{I}_8)$ であることと、 $\boldsymbol{\ell}$ の各成分は $\pm 1$ であることより、
なので、これより
となる。帰無仮説 $H_0^{[k]}$ のもとでは $\tau_{[k]}=0$ なので、因子 $[k]$ の平方和について
が成り立つ。因子 $[k]$ の自由度は $1$ なので、その平均平方を
\begin{align} V_{[k]} &= \frac{\mathrm{SS}_{[k]}}{1} = \mathrm{SS}_{[k]} \end{align}と書けば、帰無仮説のもとで $V_{[k]}/\sigma^2\sim\chi_1^2$ である。
同様の議論を誤差自由度についても考えよう。空けておいた列 $[e]\in\mathcal{E}$ には因子を割り当てていないため、仮定したモデルのもとでは $\boldsymbol{\ell}_e^\top\boldsymbol{x}=\boldsymbol{\ell}_e^\top\boldsymbol{\epsilon}$ である。したがって、各誤差列については、（帰無仮説を仮定するまでもなく）
となる。異なる列への誤差射影は互いに独立なので、これらを $\nu_E=|\mathcal{E}|$ 列について足し合わせると、カイ二乗分布の再生性より
を得る。誤差平均平方を
\begin{align} V_E &= \frac{\mathrm{SSE}}{\nu_E} \end{align}とおけば、$\nu_E V_E/\sigma^2\sim\chi_{\nu_E}^2$ である。また、これまで見てきたように因子列 $[k]$ と誤差列は直交しているため、$V_{[k]}$ と $V_E$ は独立である。したがって、帰無仮説のもとでは、
となる。よって、有意水準 $\alpha$ で検定する場合は、上の $F$ 値が $F_{1,\nu_E}(\alpha)$ を超える場合に、因子 $[k]$ は特性値に有意な影響を与えていると判断する。ただし、ここで $F_{1,\nu_E}(\alpha)$ は自由度 $(1,\nu_E)$ の $F$ 分布の上側 $\alpha$ 点を表す。
実際の計算では、次のような分散分析表に整理するとよい。
| 因子 | 平方和 | 自由度 | 平均平方 | $F$ 値 |
|---|---|---|---|---|
| $[1]$ | $\mathrm{SS}_{[1]}$ | $1$ | $V_{[1]}=\mathrm{SS}_{[1]}$ | $F_{[1]}=V_{[1]}/V_E$ |
| $[2]$ | $\mathrm{SS}_{[2]}$ | $1$ | $V_{[2]}=\mathrm{SS}_{[2]}$ | $F_{[2]}=V_{[2]}/V_E$ |
| $\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ |
| 誤差 | $\mathrm{SSE}$ | $\nu_E$ | $V_E=\mathrm{SSE}/\nu_E$ | |
| 合計 | $\mathrm{SST}$ | $n-1$ |
なお、$L_8$ 直交表では実験回数が少ないため、誤差自由度が小さくなりやすい。誤差自由度が小さいと、$V_E$ の推定が不安定になり、検出力も弱くなる。したがって、場合によっては次項で紹介する点推定と信頼区間はあくまで参考程度にとどめ、影響因子の検定までにとどめておく事も考えたほうがよい。
（5） 点推定と信頼区間
最後に、直交表実験における点推定と信頼区間を考えよう。まず、因子 $[k]$ の効果 $\tau_{[k]}$ の点推定量は、すでに述べたように
である。確率モデルを代入すると、
\begin{align} \hat{\tau}_{[k]} &= \tau_{[k]} + \frac{1}{n}\boldsymbol{\ell}_k^\top\boldsymbol{\epsilon} \end{align}となる。ここで $\boldsymbol{\ell}_k^\top\boldsymbol{\epsilon}\sim\mathcal{N}(0,n\sigma^2)$ なので、
\begin{align} \hat{\tau}_{[k]} &\sim \mathcal{N}\left(\tau_{[k]},\frac{\sigma^2}{n}\right) \end{align}つまり
\begin{align} \frac{\sqrt{n}(\hat{\tau}_{[k]}-\tau_{[k]})}{\sigma} &\sim \mathcal{N}(0,1) \end{align}が成り立つ。$\sigma^2$ が未知なので、これまで通り $V_E$ に置き換えて推定する。 これより、$\hat{\tau}_{[k]}$ と $V_E$ は独立であることに注意して、
となる。よって、因子効果 $\tau_{[k]}$ の $1-\alpha$ 信頼区間は
となる。ここで $t_{\nu_E,\alpha/2}$ は、自由度 $\nu_E$ の $t$ 分布の上側 $\alpha/2$ 点である。
次に、因子 $[k]$ の各水準における母平均を推定しよう。直交表では他の因子の水準が均等に現れるため、ここでいう母平均は他の因子の効果を平均化した周辺平均である。全体平均 $\mu$ の点推定量は
\begin{align} \hat{\mu} &= \overline{x} = \frac{1}{n}\boldsymbol{1}^\top\boldsymbol{x} \end{align}であり、確率モデルと $\boldsymbol{1}^\top\boldsymbol{\ell}_k=0$ より
\begin{align} \overline{x} &= \mu+\frac{1}{n}\boldsymbol{1}^\top\boldsymbol{\epsilon} \sim \mathcal{N}\left(\mu,\frac{\sigma^2}{n}\right) \end{align}となる。また、
\begin{align} \mathrm{Cov} \left( \frac{1}{n}\boldsymbol{1}^\top\boldsymbol{\epsilon}, \frac{1}{n}\boldsymbol{\ell}_k^\top\boldsymbol{\epsilon} \right) &= \frac{\sigma^2}{n^2} \boldsymbol{1}^\top\boldsymbol{\ell}_k =0 \end{align}であり、両者は同時正規分布に従うので、$\overline{x}$ と $\hat{\tau}_{[k]}$ は独立である。したがって、水準 $1$ と水準 $2$ の母平均 $\mu+\tau_{[k]}$、$\mu-\tau_{[k]}$ の点推定量はそれぞれ
\begin{align} \widehat{\mu+\tau_{[k]}} &= \overline{x}+\hat{\tau}_{[k]},\\ \widehat{\mu-\tau_{[k]}} &= \overline{x}-\hat{\tau}_{[k]} \end{align}であり、独立な二つの推定誤差の分散を足し合わせることにより、
\begin{align} \overline{x}\pm\hat{\tau}_{[k]} &\sim \mathcal{N}\left( \mu\pm\tau_{[k]}, \frac{2\sigma^2}{n} \right) \end{align}となる。これらの推定量も $V_E$ と独立なので、
\begin{align} T_{[k],\pm} &= \frac{ (\overline{x}\pm\hat{\tau}_{[k]}) - (\mu\pm\tau_{[k]}) }{ \sqrt{2V_E/n} } \sim t_{\nu_E} \end{align}が成り立つ。したがって、水準 $1$ の母平均 $\mu+\tau_{[k]}$ の $1-\alpha$ 信頼区間は
であり、水準 $2$ の母平均 $\mu-\tau_{[k]}$ の $1-\alpha$ 信頼区間は
となる。平方根の中に $2$ が現れるのは、全体平均方向と因子方向が直交しており、それぞれの推定誤差の分散 $\sigma^2/n$ が加算されるからである。
以上より、直交表実験の分析は次のようにまとめられる。まず、直交表の各列に対応する $\boldsymbol{\ell}_k^\top\boldsymbol{x}$ を計算して、各因子の効果を点推定する。次に、その値から平方和を作り、誤差列またはプーリングによって誤差分散 $V_E$ を推定する。そして最後に、$F$ 検定によって有意な因子を判断し、必要に応じて効果や水準平均の信頼区間を求める。つまり直交表とは、少ない実験回数で多くの因子を調べるための単なる表ではなく、平方和の分解・独立性・分散分析を同時に成立させるための設計なのである。