統計学の基礎をすぐに思い出せるような、まとめチートシートを作成しました。
「確率分布編」「推定・検定編」「線形回帰編」「その他の話題（理工学分野）」の４部構成です。本ページは「線形回帰編」です。
注意
あくまで個人用のメモであり、誤植や不正確な情報が含まれている可能性があります。適宜更新していきます。
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}$ ベクトルになっているものを抜き出したものである。
推定量の平均・分散
$\boldsymbol{\hat \beta}$ を確率変数 $\boldsymbol{u}$ の関数と見て、その平均と分散を計算しよう。これは、パラメータ推定の誤差評価に対応する。
まず、$(\ref{eq:ols-estimator})$ に $\boldsymbol{y}=X\boldsymbol{\beta}+\boldsymbol{u}$ を代入して、次が得られる。
これより、$E[\boldsymbol{u}]=\boldsymbol{0}$ であれば、
である。言い換えれば、最小二乗推定量 $\hat{\boldsymbol{\beta} }$ は不偏量である。
また、$(\ref{eq:beta-error})$ より、$\mathrm{Cov}(\boldsymbol{u})=\sigma^2I_N$ であれば
である。
分散推定
次に、確率変数 $\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$ 個の回帰係数を推定した分だけ自由度を消費しているためである。
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原則とは、反復、無作為化、局所管理の $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},x_{12},\ldots,x_{1J_1}$ | $J_1$ | $\overline{x}_{1\cdot}$ |
| $A_2$ | $x_{21},x_{22},\ldots,x_{2J_2}$ | $J_2$ | $\overline{x}_{2\cdot}$ |
| $\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ |
| $A_I$ | $x_{I1},x_{I2},\ldots,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}$ に分解される。
ここで、実験値が以下の確率モデルで記述できると仮定しよう。
ここで $\mu$ は全体平均、$\alpha_i$ は水準 $A_i$ の効果、$\varepsilon_{ij}$ は誤差項である。ただし
\begin{align} \sum_i J_i\alpha_i=0 \end{align}である。このとき、水準による差がないという帰無仮説は
で表される。対立仮説は、少なくとも $1$ つの水準効果 $\alpha_i$ が $0$ ではない、というものである。
帰無仮説のもとでは、水準間のばらつきも誤差によって生じたものにすぎない。そのため、
が成り立ち、さらにこれらは独立である。したがって、
となる。$\mathrm{BSS}/(I-1)$ は群間平均平方、$\mathrm{WSS}/(N-I)$ は群内平均平方である。水準間の平均差が大きいほど $\mathrm{BSS}$ が大きくなり、したがって $F$ も大きくなる。よって、有意水準 $\alpha$ の検定では、
のときに $H_0$ を棄却する。ただし $F_{I-1,N-I}(1-\alpha)$ は、自由度 $(I-1,N-I)$ の $F$ 分布の $1-\alpha$ 分位点である。
| 変動要因 | 平方和 | 自由度 | 平均平方 | $F$ 値 |
|---|---|---|---|---|
| 水準間 | $\mathrm{RSS}$ | $I-1$ | $\mathrm{RSS}/(I-1)$ | $\dfrac{\mathrm{RSS}/(I-1)}{\mathrm{WSS}/(N-I)}$ |
| 水準内 | $\mathrm{WSS}$ | $N-I$ | $\mathrm{WSS}/(N-I)$ | |
| 全体 | $\mathrm{TSS}$ | $N-1$ |
この表を分散分析表という。一元配置分散分析では、群間平均平方が群内平均平方に比べて十分大きいかどうかを見ることで、水準間に有意な差があるかを判定する。
乱塊法
一元配置分散分析では、処理の違い以外のばらつきはすべて誤差に入れてしまう。しかし実験では、「処理」とは別に、結果に影響しそうな条件があらかじめ分かっていることがある。たとえば、同じ処理を試すとしても、実験日、材料のロット、測定者、圃場の場所、被験者の個人差などによって結果が変わることがある。
このようなばらつきを放っておくと、処理の効果が見えにくくなる。そこで、似た条件のものをあらかじめひとまとまりにして、そのまとまりの中で各処理を比較する。これが乱塊法である。ここで「塊」はブロックのことであり、乱塊法はフィッシャーの3原則のうち、特に局所管理を具体化した方法と見ればよい。
処理を $A_1,\cdots,A_I$、ブロックを $B_1,\cdots,B_J$ とし、各ブロックの中で各処理を $1$ 回ずつ行うとする。このとき、各ブロック内で処理の割り付け順や割り付け位置は無作為に決める。処理 $A_i$ をブロック $B_j$ で行ったときの観測値を $x_{ij}$ と書くと、基本的なモデルは
である。ここで $\alpha_i$ は処理 $A_i$ の効果、$\beta_j$ はブロック $B_j$ の効果である。通常は $\sum_i\alpha_i=0,\ \sum_j\beta_j=0$ と制約を置いて、効果の基準をそろえる。誤差項は互いに独立に平均 $0$、分散 $\sigma^2$ の正規分布に従うと仮定する。
このモデルでは、全体のばらつきを、処理によるばらつき、ブロックによるばらつき、残りの誤差に分解する。平均を
\begin{align} \overline{x}_{i\cdot} &= \frac{1}{J}\sum_j x_{ij}, & \overline{x}_{\cdot j} &= \frac{1}{I}\sum_i x_{ij}, & \overline{x}_{\cdot\cdot} &= \frac{1}{IJ}\sum_i\sum_j x_{ij} \end{align}と書くと、平方和は
\begin{align} \mathrm{TSS} &= \sum_i\sum_j (x_{ij}-\overline{x}_{\cdot\cdot})^2,\\ \mathrm{SSA} &= J\sum_i (\overline{x}_{i\cdot}-\overline{x}_{\cdot\cdot})^2,\\ \mathrm{SSB} &= I\sum_j (\overline{x}_{\cdot j}-\overline{x}_{\cdot\cdot})^2,\\ \mathrm{SSE} &= \sum_i\sum_j (x_{ij}-\overline{x}_{i\cdot}-\overline{x}_{\cdot j}+\overline{x}_{\cdot\cdot})^2 \end{align}であり、
と分解できる。ブロックを入れることで、ブロック間の差 $\mathrm{SSB}$ を誤差から取り除けるので、処理効果 $\alpha_i$ の検出力が上がることがある。乱塊法の狙いは、ブロック効果そのものを主役として調べることではなく、処理効果を見やすくするために、邪魔なばらつきを先に分けておくことにある。
| 変動要因 | 平方和 | 自由度 | 平均平方 | $F$ 値 |
|---|---|---|---|---|
| 処理 | $\mathrm{SSA}$ | $I-1$ | $\mathrm{MSA}=\mathrm{SSA}/(I-1)$ | $F_A=\mathrm{MSA}/\mathrm{MSE}$ |
| ブロック | $\mathrm{SSB}$ | $J-1$ | $\mathrm{MSB}=\mathrm{SSB}/(J-1)$ | $F_B=\mathrm{MSB}/\mathrm{MSE}$ |
| 誤差 | $\mathrm{SSE}$ | $(I-1)(J-1)$ | $\mathrm{MSE}=\mathrm{SSE}/((I-1)(J-1))$ | |
| 全体 | $\mathrm{TSS}$ | $IJ-1$ |
処理効果がないという帰無仮説 $H_0^A:\alpha_1=\cdots=\alpha_I=0$ のもとでは、
となる。したがって、乱塊法では、ブロックによるばらつきを取り除いた後の誤差平均平方と処理の平均平方を比べることで、処理間に有意な差があるかを判定する。ただし、各ブロック内で各処理を $1$ 回ずつしか行わない基本的な乱塊法では、処理とブロックの交互作用を誤差と分けて推定できない。そのため、処理効果がブロックによって大きく変わらない、という加法的な見方が前提になる。
二元配置分散分析
次に、$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_{11k}$ | $x_{12k}$ | $\cdots$ | $x_{1Jk}$ |
| $A_2$ | $x_{21k}$ | $x_{22k}$ | $\cdots$ | $x_{2Jk}$ |
| $\vdots$ | $\vdots$ | $\vdots$ | $\ddots$ | $\vdots$ |
| $A_I$ | $x_{I1k}$ | $x_{I2k}$ | $\cdots$ | $x_{IJk}$ |
ここで、平均を次のように書く。
\begin{align} \overline{x}_{ij\cdot} &= \frac{1}{K}\sum_k x_{ijk},\\ \overline{x}_{i\cdot\cdot} &= \frac{1}{JK}\sum_j\sum_k x_{ijk},\\ \overline{x}_{\cdot j\cdot} &= \frac{1}{IK}\sum_i\sum_k x_{ijk},\\ \overline{x}_{\cdot\cdot\cdot} &= \frac{1}{IJK}\sum_i\sum_j\sum_k x_{ijk} \end{align}二元配置分散分析では、全体のばらつきを、要因 $A$ によるばらつき、要因 $B$ によるばらつき、交互作用によるばらつき、誤差によるばらつきに分解する。
このとき、平方和は
と分解される。$\mathrm{SSA}$ は要因 $A$ の主効果、$\mathrm{SSB}$ は要因 $B$ の主効果、$\mathrm{SSAB}$ は交互作用、$\mathrm{SSE}$ は誤差に対応する平方和である。
さて、確率モデルとして
を考えよう。ここで $\alpha_i$ は要因 $A$ の効果、$\beta_j$ は要因 $B$ の効果、$(\alpha\beta)_{ij}$ は交互作用である。交互作用とは、いわば要因 $A$ と要因 $B$ の相乗効果である。
| 変動要因 | 平方和 | 自由度 | 平均平方 | $F$ 値 |
|---|---|---|---|---|
| 要因 $A$ | $\mathrm{SSA}$ | $I-1$ | $\mathrm{SSA}/(I-1)$ | $F_A$ |
| 要因 $B$ | $\mathrm{SSB}$ | $J-1$ | $\mathrm{SSB}/(J-1)$ | $F_B$ |
| 交互作用 | $\mathrm{SSAB}$ | $(I-1)(J-1)$ | $\mathrm{SSAB}/((I-1)(J-1))$ | $F_{AB}$ |
| 誤差 | $\mathrm{SSE}$ | $IJ(K-1)$ | $\mathrm{SSE}/(IJ(K-1))$ | |
| 全体 | $\mathrm{TSS}$ | $IJK-1$ |
検定では、次の $3$ つの帰無仮説を考える。ただし、実際にはまず交互作用の有無を調べるのが自然である。
- 交互作用なし：$H_0^{AB}:(\alpha\beta)_{ij}=0\ \style{font-family:inherit;}{\text{ for all }}i,j$
- 要因 $A$ の主効果なし：$H_0^A:\alpha_1=\cdots=\alpha_I=0$
- 要因 $B$ の主効果なし：$H_0^B:\beta_1=\cdots=\beta_J=0$
順番としては、最初に交互作用 $H_0^{AB}$ について検定するべきである。交互作用が有意であれば、要因 $A$ の効果は要因 $B$ の水準によって変わるので、$A$ や $B$ の主効果だけを単独で解釈するのは危険である。一方、交互作用が有意でないと判断できるなら、その上で $H_0^A$ と $H_0^B$ のもとで、要因 $A,B$ の主効果を検定する。
基本的な方針としては、各平方和を自由度で割った平均平方を作り、誤差平均平方 $\mathrm{MSE}=\mathrm{SSE}/(IJ(K-1))$ と比べる。なお、例えば $H_0^A$ を検定するときに、要因 $B$ の効果 $\beta_j$ が $0$ である必要はないという点に注意しよう。二元配置分散分析では、要因 $A$ による変動 $\mathrm{SSA}$ と要因 $B$ による変動 $\mathrm{SSB}$ を分けているので、$B$ に主効果があっても、交互作用がなければ $F_A$ によって $A$ の主効果を検定できる。同様に、$A$ に主効果があっても $F_B$ によって $B$ の主効果を検定できる。
順番に説明しよう。
（1）相互作用の検定
交互作用なしの帰無仮説 $H_0^{AB}$ のもとでは、主効果 $\alpha_i,\beta_j$ が $0$ である必要はない。バランスした二元配置では、主効果の平方和と交互作用の平方和が直交して分解されるので、$\alpha_i,\beta_j$ が値を持っていても、交互作用の平方和 $\mathrm{SSAB}$ は誤差分散と比較できる。よって帰無仮説 $H_0^{AB}$ のもとで、
この検定で $H_0^{AB}$ を棄却するなら、交互作用があると判断する。この場合、要因 $A$ の効果は要因 $B$ の水準によって変わるので、主効果だけを単独で読むのは慎重にする。
（2）各要因の検定
交互作用が有意でないと判断できる場合には、次に要因 $A$ と要因 $B$ の主効果を調べる。要因 $A$ の主効果なしの帰無仮説 $H_0^A$ のもとでは、要因 $B$ の主効果 $\beta_j$ は $0$ でなくてもよい。よって、$H_0^A$ のもとで以下が成り立つ。
$F_A$ が大きいほど、要因 $A$ の主効果が有意である可能性が高い。
$H_0^A$ を検定するときに、要因 $B$ の効果 $\beta_j$ が $0$ である必要はないという点に注意しよう。二元配置分散分析では、要因 $A$ による変動 $\mathrm{SSA}$ と要因 $B$ による変動 $\mathrm{SSB}$ を分けているので、$B$ に主効果があっても、交互作用がなければ $F_A$ によって $A$ の主効果を検定できる。同様に、$A$ に主効果があっても $F_B$ によって $B$ の主効果を検定できる。
同様に、要因 $B$ の主効果なしの帰無仮説 $H_0^B$ のもとでは、要因 $A$ の主効果 $\alpha_i$ は $0$ でなくてもよい。よって、$H_0^B$ のもとで
である。このように、二元配置分散分析では、まず $F_{AB}$ で交互作用を検定し、その後に必要に応じて $F_A,F_B$ で主効果を検定する、という順番で読むとよい。
直交表
実際の実験では、要因はたくさん存在することが普通である。例えば毎日の料理を実験と見たとき、温度、圧力、材料、その人の腕前や愛情などの要因が影響を与えるだろう。これらの要因をそれぞれ $A,B,C,D$ としよう。各要因について $2$ つの水準を用意すると、すべての組み合わせを試すには $2^4=16$ 通りの実験が必要になる。
これは要因や水準が増えると爆発的に増えていく。全ての組み合わせを試すことができない場合、どの組み合わせを選ぶべきだろうか？ 直交表は、この問題に対する解決策の一つである。
ここでは、それぞれの要因に相互作用がないと仮定し、以下の
\begin{align} x_{ijk\ell} &= \mu+\alpha_i+\beta_j+\gamma_k+\delta_\ell+\varepsilon_{ijk\ell} \end{align}というモデルを考える。ここで $\alpha_i,\beta_j,\gamma_k,\delta_\ell$ は、それぞれ要因 $A,B,C,D$ の主効果である。
このように主効果だけを見たいなら、すべての組み合わせを試さなくても、各要因の水準がバランスよく現れるように実験を選べばよい。これを実現するための表が直交表である。
具体的には、次の $L_8$ 直交表を使うことで、$16$ 通りの候補から $8$ 通りだけを選んで実験できる。$L_8$ 直交表は、$2$ 水準の列を $7$ 本もつ直交表である。今回はその中からランダムに $4$ 列を選び、列 $1,3,6,7$ に要因 $A,B,C,D$ を割り付けることにする。残りの列 $2,4,5$ は今回は使わない。
| 実験番号 | 要因 $A$ | 未使用 | 要因 $B$ | 未使用 | 未使用 | 要因 $C$ | 要因 $D$ |
|---|---|---|---|---|---|---|---|
| 1 | $1$ | $1$ | $1$ | $1$ | $1$ | $1$ | $1$ |
| 2 | $1$ | $1$ | $2$ | $1$ | $2$ | $2$ | $2$ |
| 3 | $1$ | $2$ | $1$ | $2$ | $1$ | $2$ | $2$ |
| 4 | $1$ | $2$ | $2$ | $2$ | $2$ | $1$ | $1$ |
| 5 | $2$ | $1$ | $1$ | $2$ | $2$ | $1$ | $2$ |
| 6 | $2$ | $1$ | $2$ | $2$ | $1$ | $2$ | $1$ |
| 7 | $2$ | $2$ | $1$ | $1$ | $2$ | $2$ | $1$ |
| 8 | $2$ | $2$ | $2$ | $1$ | $1$ | $1$ | $2$ |
実際に、上の $L_8$ 直交表で主効果がどのように取り出せるかを見てみよう。上の $L_8$ 直交表に基づいて $8$ 回実験を行い、それぞれの実験番号 $r$ の観測値を $Y_r$ と書く。要因 $A$ が水準「$1$」のときの実験（上の表だと実験番号 $1,2,3,4$）で得られた観測値の平均を $\overline{Y}_{A1}$ と書き、水準「$2$」のときの実験（実験番号 $5,6,7,8$）で得られた観測値の平均を $\overline{Y}_{A2}$ とする。つまり、
\begin{align} \overline{Y}_{A1} &= \frac{1}{4}(Y_1+Y_2+Y_3+Y_4),\\ \overline{Y}_{A2} &= \frac{1}{4}(Y_5+Y_6+Y_7+Y_8) \end{align}とする。相互作用がないモデル
\begin{align} Y &= \mu+\alpha_i+\beta_j+\gamma_k+\delta_\ell+\varepsilon \end{align}を考え、識別のために
\begin{align} \alpha_1+\alpha_2 = \beta_1+\beta_2 = \gamma_1+\gamma_2 = \delta_1+\delta_2 =0 \end{align}と置く。このとき、$L_8$ 直交表では、要因 $A$ の水準 $1$ の行だけを見ても、他の各要因の水準 $1,2$ は同じ回数だけ現れる。したがって、
\begin{align} E[\overline{Y}_{A1}] &= \mu+\alpha_1 + \frac{\beta_1+\beta_2}{2} + \frac{\gamma_1+\gamma_2}{2} + \frac{\delta_1+\delta_2}{2}\\ &= \mu+\alpha_1,\\ E[\overline{Y}_{A2}] &= \mu+\alpha_2 \end{align}となる。同様に、
\begin{align} E[\overline{Y}_{B1}]&=\mu+\beta_1,& E[\overline{Y}_{B2}]&=\mu+\beta_2,\\ E[\overline{Y}_{C1}]&=\mu+\gamma_1,& E[\overline{Y}_{C2}]&=\mu+\gamma_2,\\ E[\overline{Y}_{D1}]&=\mu+\delta_1,& E[\overline{Y}_{D2}]&=\mu+\delta_2 \end{align}が成り立つ。よって、たとえば要因 $A$ の主効果は
\begin{align} E[\overline{Y}_{A1}-\overline{Y}_{A2}] &= \alpha_1-\alpha_2 = 2\alpha_1 \end{align}で推定できる。$\beta,\gamma,\delta$ についても同様である。つまり直交表を使うと、他の要因の影響が平均の中で打ち消され、各要因の主効果をきれいに取り出せる。
最後に、直交表の数学的な意味をもう少し明確にしておこう。$L_8$ 直交表における $2$ を $-1$ と書き換え、各列を $8$ 次元の列ベクトル $\boldsymbol a,\cdots ,\boldsymbol g$ で書く。すると、これらの列ベクトルは内積 $0$ で互いに直交している。
この直交性が「効果が混ざらない」ことに対応している。符号化したモデルを
\begin{align} \boldsymbol y &= \mu\boldsymbol 1 + \tau_A\boldsymbol a + \tau_B\boldsymbol c + \tau_C\boldsymbol f + \tau_D\boldsymbol g + \boldsymbol\varepsilon \end{align}と書く。ここで $\tau_A,\tau_B,\tau_C,\tau_D$ は、$+1$ 水準と $-1$ 水準の差を表す係数である。要因 $A$ の効果を取り出すには、両辺に $\boldsymbol a^\top$ を掛ければよい。
\begin{align} \boldsymbol a^\top\boldsymbol y &= \mu\boldsymbol a^\top\boldsymbol 1 + \tau_A\boldsymbol a^\top\boldsymbol a + \tau_B\boldsymbol a^\top\boldsymbol c + \tau_C\boldsymbol a^\top\boldsymbol f + \tau_D\boldsymbol a^\top\boldsymbol g + \boldsymbol a^\top\boldsymbol\varepsilon \end{align}直交表では、$\boldsymbol a^\top\boldsymbol 1=0$ であり、さらに $\boldsymbol a^\top\boldsymbol c=\boldsymbol a^\top\boldsymbol f=\boldsymbol a^\top\boldsymbol g=0$ である。したがって、上式は
\begin{align} \boldsymbol a^\top\boldsymbol y &= \tau_A\boldsymbol a^\top\boldsymbol a + \boldsymbol a^\top\boldsymbol\varepsilon \end{align}となる。つまり、$\boldsymbol a^\top\boldsymbol y$ には要因 $B,C,D$ の主効果が入ってこない。これが、列ベクトルが直交していることと「他の要因の効果が混ざらない」ことの対応である。未使用列はモデルに入れていないので、効果推定には直接使わない。ただし、どの列を要因に割り付けても列同士は直交しているため、必要に応じて別の要因を割り付ける余地として残せる。
直交表を使うと、実験回数を減らしながら、各要因の主効果を公平に比較できる。ただし、実験回数を減らしている以上、すべての情報が得られるわけではない。特に、交互作用が強い場合には、主効果と交互作用が混ざって解釈しにくくなることがある。したがって、直交表は「どの交互作用を無視してよいか」という実験上の判断とセットで使う必要がある。