分散分析(ANOVA)のF検定の導出について１——制約つき最小二乗法を用いたアプローチ

前置き
準備
ガウス・マルコフモデル
線形制約条件つき最小二乗法からのF検定の導出
分散分析のF検定の導出
参考文献

前置き

分散分析(analysis of variance, ANOVA)は複数のグループの平均に統計学的に優位な差があるかを調べる方法である. この文章では分散分析のF検定の検定統計量が数学的にどのように導出されるかを説明する. ここでは, 検定統計量の数学的な導出にフォーカスし, 分散分析自体についての説明はしない. これは, 分散分析を学び始める人にとってはこの内容は過度に詳細だろうし, 一定の理解がある人にとっては一からの説明は恐らく冗長だからだ.

このような古典的な内容をなぜ改めて書くのかを簡単に説明する. 私は最近, 数理統計学を勉強している. 私の読んでいるテキスト[久保川]では¹, 分散分析の概要と、F検定で等平均性をテストする方法が説明されているが, その統計量が何故F分布に従うかは書かれていなかった. そこで, 他の手持ちのテキストやウェブ上で色々と調べてみたが, 検定統計量の導出を解説した文章はあまりなく, あったとしても長大なテキストの一部だったり, 線形回帰モデル(もしくは, ガウス・マルコフモデル)の応用として軽く触れられているだけで詳細は省略されていたりした. 最終的には, [吉田], [ラオ]の中に厳密な導出を見つけることができたが, 前者は洗練されているがやや直感的に理解しにくく, 後者は一部証明が理解できなかった. ともあれ, これらを通じて導出の概要はわかり, [吉田]の方法をやや素朴にした(つまり, 洗練度を減じた)方法で証明を与えることができたので, その内容をまとめておくことにした.

また, ここで説明するアプローチは少し天下り的だが, 尤度比検定を出発点として自然にこのアプローチに至る道筋を以下の記事で説明したので, 興味があれば, 追って読んでほしい.

kihanto.com

今回の内容は, self-contained に説明するのは骨が折れるので, いくつかの命題は証明せずに用いる. といっても, それらはどれも適切な予備知識があれば, 容易に証明は理解できるものばかりである. 具体的には非心カイ二乗分布と分布の特性関数についての基本的な知識があれば十分である. それらを除けば, 必要な道具は基本的には, 学部一回生レベルの線形代数, 微分積分であり, これらは前提知識として断りなく用いる.

ここからの内容のアウトラインを説明する. まず準備として, 線形代数から単純ながらも鍵となる命題を2つ思い出し, 今回の議論で中心的な役割を果たすフィッシャー・コクランの定理を証明なしで述べる. 次に, ガウス・マルコフモデルを導入する. これは分散分析モデルを含む, より一般的なモデルであり, これに対し, 制約つき最小二乗値の従う分布を導出する. 最後に, ガウス・マルコフモデルに対して得た結果を分散分析モデルに適用することで, よく知られているF検定の検定統計量が導出されることを示す.

準備

まず, 線形代数から2つの命題を準備する. いずれもほぼ常識と言える命題である.

以下では, 線形代数の標準的と思われる表記法は断りなく用いる. ベクトル空間は常に適切な次元のユークリッド空間とし, 内積は常に標準内積を考えるものとする. また, 行列とそれの定める線形変換を同一視する. つまり, 行列 $A$ から定まる線形変換 $x \mapsto Ax$ のことを線形変換 $A$ と書いたり, その変換の像のことを $\mathrm{Im} \ A$ と書いたりする. 行列 $A$ の転置を $A'$ とかく.

1つ目の命題は内積と行列の転置の基本的な関係についてである.

命題1 行列 $A$ に対し, $\mathrm{Ker} \ A' = (\mathrm{Im} \ A)^{\perp}$ .

証明ベクトル $x,y$ に対し, $x \in \mathrm{Ker} \ A'$ なら, $x'Ay = (A'x)'y = 0$ . $y$ は任意に取れるので, $x \in (\mathrm{Im} \ A)^{\perp}$ . 逆に, $x \in (\mathrm{Im} \ A)^{\perp}$ なら, $(A'x)'y = x'Ay = 0$ . $y$ は任意より, $A'x = 0$ . つまり, $x \in \mathrm{Ker} \ A'$ . (証明終わり)

線形変換 $P$ が $P^{2} = P$ , $P' = P$ を満たすとき, 直交射影という. 部分空間 $U$ が $P$ の像であるとき, $P$ を $U$ への射影という. 2つ目の命題は, 直交射影が任意のベクトルを射影の像に最も近い点に移すことを示すものである.

命題2

$P$ を部分空間 $U$ への直交射影とするとき, 任意の $y \in \mathbb{R}^{n}$ に対し, $\min_{z \in U}|y - Pz|^{2} = |y - Py|^{2} = y'(I - P)y.$

証明

任意の $z \in U$ に対し, $|y - Pz|^{2} = |y - Py|^{2} + |P(y-z)|^{2} + 2(y-Py)'P(y-z).$ ここで, $(I-P)P = 0$ より, $(y-Py)'P(y-z) = y'(I-P)P(y-z) = 0$ . よって, $|y - Pz|^{2} \geq |y - Py|^{2}$ . (証明終わり)

次に, フィッシャー・コクランの定理と呼ばれる定理を紹介する. これは分散分析の検定統計量の導出において, 最も基本的な役割を果たす定理で, 単位行列の分解とガウス確率変数の二次形式の関係を示している. 自由度 $n$ , 非心度 $\delta$ の非心カイ二乗分布を $\chi^{2}(n,\delta)$ と書く. 非心カイ二乗分布については例えば, [吉田], Wikipedia 参照. ただし, Wikipedia の日本語版記事はかなり内容が薄いので, 英語版を勧める.

注意: 非心カイ二乗分布はややマイナーであり, また, その確率密度関数には無限和が表れ, 厄介そうに見えるかもしれない. しかし, カイ二乗分布と同様の再生性をもち, 特性関数(もしくはモーメント母関数)もかなり単純で, 意外に扱いやすい対象である.

$Y = (Y_1, \dots, Y_n)'$ を $n$ 次元のガウス確率変数で平均ベクトル $\mu = (\mu_1, \dots,\mu_n)'$ , 分散共分散行列 $\sigma^2 I$ ( $I$ は n次単位行列) をもつものとする.

定理1 (フィッシャー・コクランの定理)

$A_1, \dots, A_k$ をn次対称行列で, $A_1 + \cdots + A_k = I$ を満たすものとする. $Q_i := Y'A_i Y$ , $\mathrm{rank} \ A_i := n_i$ とする. このとき, 以下は同値:

各 $1 \leq i \leq k$ に対し, ある $\delta_i \geq 0$ が存在し, $Q_i \sim \chi^2(n_i,\delta_i)$ かつ $Q_1, \dots, Q_k$ は独立.
$n_1 + \cdots + n_k = n$ .
各 $1 \leq i \leq k$ に対し, $A_i^2 = A_i$ .
各 $1 \leq i,j \leq k$ , $i \neq j$ に対し, $A_i A_j = O$ ( $O$ はゼロ行列).

またこのとき, $\delta_i = \mu' A_i \mu$ , $\mu_1^2 + \cdots + \mu_k^2 = \delta_1 + \cdots + \delta_k$ .

証明は省略する. 例えば, [吉田]参照. 以下で用いるのは, 次の系で与えられる, この定理の特別な場合のみである.

系1

$P_1$ , $P_2$ を $\mathbb{R}^n$ 上の直交射影とする. $Q_i := Y'P_iY \ (i = 1,2)$ , $n_i := \mathrm{rank} \ P_i$ , $\delta_i := \mu'P_i\mu$ とする. このとき, 以下が成立:

$Q_i \sim \chi^2(n_i,\delta_i)$ .
$P_1$ と $P_2$ が互いに直交するとき, つまり, $P_1P_2 = P_2P_1 = O$ のとき, $Q_1$ と $Q_2$ は独立.

証明

(i): 分解 $I = P_i + (I - P_i)$ は明らかに $P_i(I - P_i) = O$ を満たすので, フィッシャー・コクランの定理から $Q_i := Y'P_iY$ がわかる.

(ii): 直交性の仮定から, $(I - P_1 - P_2)^2 = I - P_1 - P_2$ . よって, 分解 $I = P_1 + P_2 + (I-P_1-P_2)$ にフィッシャー・コクランの定理を適用することでわかる. (証明終わり)

この系も含め, 以下では単にフィッシャー・コクランの定理と呼ぶ.

ガウス・マルコフモデル

ガウス・マルコフモデルとは, 観測値がいくつかの説明変数の線型結合と誤差との和によって与えられるとするモデルである. より具体的には, $n$ 個の観測値 $Y_1, \dots, Y_n$ が与えられたとき, それぞれの観測値は $p$ 個の説明変数 $x_{i1},\dots, x_{ip} \ (1 \leq i \leq n)$ の線形関係によって生成されると考え, そして, 観測には誤差が伴うとしたものである. また, 誤差は独立同分布で分散 $\sigma^2$ の正規分布に従うとする. つまり, 以下の観測値ベクトル $Y$ と計画行列 $X$ は, あるパラメータベクトル $\beta$ と誤差ベクトル $\varepsilon$ に対し, 以下のような関係を満たすとする:

$Y = X \beta + \varepsilon. \label{}$

ここで,

で, $(\varepsilon_i)_{1 \leq i \leq n}$ は独立同分布で平均 $0$ , 分散 $\sigma^2$ の正規分布に従うとする. $Y$ と $X$ は既知の値として扱い, パラメータ $\beta$ と分散 $\sigma^2$ は未知の値である. $r := \mathrm{rank} \ X \leq n-1$ を仮定する. これは全く制限的な仮定ではない. またいくつかのテキストでは仮定されているようだが, ここでは $X$ がフルランクであること, つまり, $\mathrm{rank}\ X = p$ は仮定しない. これを仮定すると分散分析モデルに適用できないためである.

ガウス・マルコフモデルは分散分析モデルよりも一般的なモデルである. 例えば, 1つの因子が3つの水準をもつ一元配置モデル $Y_{ij} = \alpha_i + \varepsilon_{ij} \quad (1 \leq i \leq 3, \ 1 \leq j \leq 2)$ とするとき,

とすればよい.

$p \times k$ 行列 $H$ 行列は $\mathrm{rank}\ H = k$ かつ $\mathrm{Im}\ H \subset \mathrm{Im}\ X'$ を満たすとし, また $d \in \mathrm{Im}\ H' \subset \mathbb{R}^k$ とする. このとき, 以下の帰無仮説 $H_0$ , 対立仮説 $H_1$ の検定を考える:

$H_0: \ H'\beta = d, \quad H_1: \ H'\beta \neq d \label{}$

線形制約条件 $H'\beta = d$ はパラメータ間に制約を与える. 例えば上の一元配置の例で,

とすると, $H'\beta = d$ は $\alpha_1 = \alpha_2 = \alpha_3$ と同値であり, 3つのパラメータが等しいという制約条件である. よって, $H$ をこのようにとれば, 帰無仮説 $H_0$ , 対立仮説 $H_1$ による検定は一元配置の等平均性の検定と一致する.

線形制約条件つき最小二乗法からのF検定の導出

$H'\beta = d$ を満たす $\beta$ を $1$ つ固定し, $\gamma$ とする. $\beta$ , $\sigma^2$ の真値を $\beta_0$ , $\sigma_0^2$ とする. 次の定理が最も重要である. 制約条件つきの最小二乗値が非心カイ二乗分布に従うことを示している.

定理2

以下を定める:

$R_1^2 := \min_{\beta:H'\beta = d} |Y-X\beta|^2. \label{}$

このとき,

$R_1^2 = (Y-X\gamma)'(I - P_{X|H})(Y-X\gamma). \label{}$

ここで, $P_{X|H}$ は直交射影で, その像は

$\{ X\beta \mid H'\beta = 0 \} \cdots (*)$

で与えられ, $\mathrm{rank} \ P_{X|H} = r-k$ . よって, $R_1^2 / \sigma_0^2$ は非心度 $(X(\beta_0-\gamma))'(I-P_{X|H})(X(\beta_0-\gamma)) / \sigma_0^2$ , 自由度 $n-(r-k)$ の非心カイ二乗分布に従う.

証明

まず, $\mathrm{rank} \ H' = k$ より, $\mathrm{dim} \ (\mathrm{Ker} \ H') = p-k$ である. $H'\gamma = d$ より, $\mathrm{Ker} \ H'$ の基底を $1$ つとり, $g_1,\dots,g_{p-k}$ とし, この基底を並べて作った $p \times (p-k)$ 行列 $G := (g_1, \dots,g_{p-k})$ とすると, $H'\beta = d$ を満たす $\beta$ の集合は $\{ \gamma + G\delta \mid \delta \in \mathbb{R}^{p-k} \}$ である. よって, $R_1^2 = \min_{\delta \in \mathbb{R}^{p-k}}|(Y - X\gamma) - XG\delta|^2$ . 命題2より, $R_1^2 = (Y-X\gamma)'(I - P)(Y - X\gamma)$ . ここで, $P$ は $\mathrm{Im} \ XG$ への直交射影. よって, $\mathrm{Im} \ XG$ が(*)で与えられ, $\mathrm{rank} \ P = r-k$ であることを示せばよい.

$G$ は $\mathrm{Ker} \ H'$ の基底を並べた行列なので, $\mathrm{Ker} \ H' = \mathrm{Im} \ G$ . よって, $\mathrm{Im} \ XG = \mathrm{Im} \ X|_{\mathrm{Ker} \ H'}. \label{}$ ただし, $X|_{\mathrm{Ker} \ H'}$ は線型写像 $X$ の定義域を部分空間 $\mathrm{Ker} \ H'$ に制限したものである. これより, $\mathrm{Im} \ XG$ が(*)に一致することがわかる. $\mathrm{Im} \ H \subset \mathrm{Im} \ X'$ であることを思い出すと,

$\mathrm{Ker} \ X = (\mathrm{Im} \ X')^{\perp} \subset (\mathrm{Im} \ H)^{\perp} = \mathrm{Ker} \ H' = \mathrm{Im} \ G. \label{}$

よって, $\mathrm{rank} \ XG = \mathrm{rank} \ G - \dim (\mathrm{Ker} \ X) = (p-k) - (p-r) = r - k$ . ゆえに, $\mathrm{rank} \ P = r -k$ . 最後の主張はフィッシャー・コクランの定理の直接的な帰結である. (証明終わり)

次に, $R_0^2 := \min_{\beta} |Y-X\beta|^2 \label{}$ と定めると, 命題2より, $R_0^2 = Y'(I-P_X)Y. \label{}$ ただし, $P_X$ は $\mathrm{Im} \ X$ への直交射影. よって, フィッシャー・コクランの定理より, $R_0^2/ \sigma_0^2$ は自由度 $n-r$ のカイ二乗分布に従う. ここで非心度が $0$ となるのは, $(X\beta_0)'(I - P_X)X\beta_0 = 0$ だからである. また明らかに, $R_0^2 = (Y-X\gamma)'(I-P_X)(Y-X\gamma)$ である. これより,

$R_1^2 - R_0^2 = (Y-X\gamma)'(P_X - P_{X|H})(Y-X\gamma) \label{}$

である. 定義から明らかに, $\mathrm{Im} \ P_{X|H} \subset \mathrm{Im} \ P_X$ であるから, $(R_1^2 - R_0^2) / \sigma_0^2$ は自由度 $r - (r-k) = k$ , 非心度 $(X(\beta_0-\gamma))'(P_X-P_{X|H})(X(\beta_0-\gamma)) / \sigma_0^2$ の非心カイ二乗分布に従う. また, $P_X - P_{X|H}$ と $P_{X|H}$ は互いに直交する直交射影であるから, フィッシャー・コクランの定理より, $R_1^2 - R_0^2$ と $R_0^2$ は独立である. 以上より,

$F := \frac{n-r}k\frac{R_1^2 - R_0^2}{R_0^2} \label{}$

は自由度 $k$ , $n-r$ , 非心度 $(X(\beta_0-\gamma))'(P_X-P_{X|H})(X(\beta_0-\gamma)) / \sigma_0^2$ の非心 $F$ 分布に従うことがわかった. 非心F分布については例えば, 英語版のWikipedia 参照 (日本語版記事は無いようである).

これで帰無仮説 $H_0$ , 対立仮説 $H_1$ の場合の検定統計量 $F$ が構成できたように見えるが, $F$ の非心度には知ることのできない分散の真値 $\sigma_0^2$ が含まれている. しかし, これは帰無仮説の下では $0$ になるので問題はない.

命題3

$\beta$ の真値 $\beta_0$ が帰無仮説 $H_0$ を満たすとする. すなわち, $H'\beta_0 = d$ とする. このとき,

$(X(\beta_0-\gamma))'(P_X-P_{X|H})(X(\beta_0-\gamma)) = 0. \label{}$

証明

$H'\beta_0 = H'\gamma = d$ より, $\beta_0 - \gamma \in \mathrm{Ker} \ H'$ . よって, $X(\beta_0-\gamma) \in \mathrm{Im} \ P_{X|H}$ . (証明終わり)

よって, 帰無仮説の下では $F$ は非心度 $0$ の通常の $F$ 分布に従うことがわかり, $F$ は $X$ と $Y$ のみから計算できるので, これで帰無仮説 $H_0$ , 対立仮説 $H_1$ の検定を行うことができる.

分散分析のF検定の導出

前節までの結果を用いて, 二元配置の分散分析モデルで交互作用のない場合の検定量を導出する. より複雑な分散分析モデルでも, 計算が多少複雑になるだけで同様に導出できる. モデルは以下で与えられる:

$Y_{ij} = \alpha_i + \beta_i + \varepsilon_{ij} \quad (1 \leq i \leq I, 1 \leq j \leq J). \label{}$

これがどのようなガウス・マルコフモデルに対応するかを念のため述べておく. 例えば, $I = 2$ , $J = 3$ のときは,

と定めればよい. このように定めるとき, 掃き出し法により, $\mathrm{rank} \ X = I+J-1$ がわかる(やや計算と考察を要する).²

このとき, 以下の帰無仮説 $H_0$ , 対立仮説 $H_1$ の検定を考える:

$H_0: \ \alpha_1 = \cdots = \alpha_I =:\alpha, \quad H_1: \ \alpha_i \neq \alpha_k \quad (\exists \ i,k). \label{}$

既にみたように帰無仮説 $H_0$ は適当な行列 $H$ を用いた線形制約条件で表現することができる. また, 容易に $\mathrm{rank} \ H = I-1$ がわかる.

前節の結果より, 検定を行うには,

$R_0^2 = \min_{\alpha_1,\dots,\alpha_I,\beta_1,\dots,\beta_J}\sum_{i,j}(Y_{ij} - \alpha_i - \beta_j)^2, \quad R_1^2 = \min_{\alpha,\beta_1,\dots,\beta_J}\sum_{i,j}(Y_{ij} - \alpha - \beta_j)^2 \label{}$

を求め, 自由度 $I-1$ , $IJ - I - J - 1 = (I-1)(J-1)$ の $F$ 分布に従う確率変数

$F = \frac{(I-1)(J-1)}{I-1}\frac{R_0^2 - R_1^2}{R_0^2}$

の値が棄却域に入るか否かで検定を行えばよい. $R_0^2$ , $R_1^2$ の値は, 例えば偏微分するなどの通常の方法で求めることができ,

$R_0^2 = \sum_{i,j}(Y_{ij} - \bar{Y}_{\cdot j} - \bar{Y}_{i \cdot} + \bar{Y}_{\cdot \cdot})^2, \quad R_1^2 = \sum_{i,j}(Y_{ij} - \bar{Y}_{\cdot j})^2 \label{}$

がわかる. ただしここで,

$\bar{Y}_{\cdot j} := \frac{\sum_iY_{ij}}I, \quad \bar{Y}_{i \cdot} := \frac{\sum_jY_{ij}}J, \quad \bar{Y}_{\cdot \cdot} := \frac{\sum_{i,j}Y_{ij}}{IJ}.$

また容易に以下がわかる:

$R_0^2 = \sum_{i,j}( (Y_{ij} - \bar{Y}_{\cdot j}) - (\bar{Y}_{i \cdot} - \bar{Y}_{\cdot \cdot} ) )^2 = R_1^2 - J\sum_i(\bar{Y}_{i \cdot} - \bar{Y}_{\cdot \cdot})^2. \label{}$

よって,

$R_1^2 - R_0^2 = J \sum_i(\bar{Y}_{i \cdot} - \bar{Y}_{\cdot \cdot})^2. \label{}$

以上より,

$F = \frac{(I-1)(J-1)}{I-1}\frac{J\sum_i(\bar{Y}_{i \cdot} - \bar{Y}_{\cdot \cdot})^2}{\sum_{i,j}(Y_{ij} - \bar{Y}_{\cdot j} - \bar{Y}_{i \cdot} + \bar{Y}_{\cdot \cdot})^2}$

これは分散分析において通常用いられている検定量であるから, これで分散分析の検定統計量が $F$ 分布に従うことが証明できた.

最後にややどうでもいいコメントをしておくと, 交互作用のない二元配置ぐらいの比較的単純なモデルならば, $R_0^2$ と $R_1^2$ の式を展開し, $R_0^2 = Y'PY$ , $R_1^2 = Y'QY$ なる直交射影 $P,Q$ を具体的に求め, それらに対してフィッシャー・コクランの定理を適用することで, 確率変数 $F$ が $F$ 分布に従うことを示すことができる. ただし, このような力技はモデルが複雑になるにつれて破綻するように思える.³ 一方で, 今回の方法はモデルが複雑化してもそこまで難しくはならないし, 検定統計量の具体系が予めわかっていなくても, 最小値問題を解くことで機械的に検定量を導出できるという点で利点がある.