分散分析(ANOVA)のF検定の導出について２——尤度比検定との等価性

前置き
準備
ガウス・マルコフモデル
尤度比検定とF検定の等価性

前置き

前回の文章では, 分散分析の等平均性の検定の検定統計量がF分布に従うことを示した.

そこでは検定統計量は最小二乗値を用いて表された. 今回の文章の目的は, 尤度比検定を出発点にして検定を構成しても, 同じ最小二乗値を用いた検定統計量が自然に現れることを示すことである.

なぜこのようなことを考えるのか, 目的を明確に述べるために, まず前回の導出の流れを振り返る:

ガウス確率変数の二次形式はフィッシャー・コクランの定理から, 分布, 独立性が容易にわかる.
ガウス・マルコフモデルにおける, 制約条件つきと条件なしの最小二乗値( $R_{1}^{2}$ と $R_{0}^{2}$ は, ガウス確率変数の二次形式で表せる.
フィッシャー・コクランの定理から, $R_{1}^{2} - R_{0}^{2}$ と $R_{0}^{2}$ は独立な二次形式となり, 前者は非心カイ二乗分布, 後者はカイ二乗分布に従う.
制約条件が実際に成り立つとき, $R_{1}^{2} - R_{0}^{2}$ は非心度0のカイ二乗分布に従うので, 比をとることでF検定が導出できる.
分散分析モデルはガウス・マルコフモデルの特別な場合であるため, 以上の内容が適用できる.
等平均性の仮説は線形制約条件で表すことができ, その仮説に対し, $R_{1}^{2} - R_{0}^{2}$ と $R_{0}^{2}$ を具体的に計算すると, よく知られた平方和の式に行き着く.
よって, 分散分析において用いられている検定量は等平均性の仮定の下ではF分布に従う.

この導出は線形代数やガウス確率変数に関する一般的で比較的単純な結果の応用として示すもので, 論理はクリアだが, 不満点もある. それは, 最小二乗値から構成した検定統計量を考える明確な動機づけに乏しく, やや天下り的なアプローチに見えるという点である.

そこでここでは, 尤度比検定を出発点として検定を構成しても, 前回考えた最小二乗値に関する検定統計量が現れることを示す. これにより, 最小二乗値による検定統計量を考えるに至る自然な道筋が存在することを示すのが, この文章の目的である.

ただし, 統計量が「自然に導かれる」ことは, それが良い統計量だということを必ずしも意味しないことを注意しておく. 統計量の良さは検出力や分散の小ささといった目的に沿った具体的な指標で評価されるべきだからだ.

ここからの内容のアウトラインを説明する. まず, 線形代数から, 射影と一般化逆行列に関する幾つかの命題を準備する. 次に, ガウス・マルコフモデルを説明する. ここは前回の引き写しである. そして, ガウス・マルコフモデルに対し, 線形制約条件を帰無仮説とした場合の尤度比検定を考えることで, 最小二乗値を用いた検定統計量を考えることに自然に導かれることを示す.

準備

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

まず前回も用いた, 線形代数から2つの命題を思い出す.

命題 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 $P$ を部分空間 $U$ への直交射影とするとき, 任意の $y \in \mathbb{R}^n$ に対し,

$\min_{z \in U}|y - Pz|^2 = |y - Py|^2 = y'(I - P)y. \label{}$

証明

任意の $z \in U$ に対し,

$|y - Pz|^2 = |y - Py|^2 + |P(y-z)|^2 + 2(y-Py)'P(y-z). \label{}$

ここで, $(I-P)P = 0$ より, $(y-Py)'P(y-z) = y'(I-P)P(y-z) = 0$ . よって, $|y - Pz|^2 \geq |y - Py|^2$ . (証明終わり)

次に一般化逆行列について基本的な命題をいくつか示す. $m \times n$ 行列 $A$ に対し, $n \times m$ 行列 $B$ が $ABA = A$ を満たすとき, $B$ は $A$ の一般化逆行列という.

注意: Wikipedia に一般化逆行列の記事はあるが, やや制限の強い定義を採用しているようである. 英語版では, ここと同じ定義で扱われている.

一般化行列は常に存在する.

命題 3

任意の $m \times n$ 行列 $A$ に対し, その一般化逆行列は存在する.

証明

$\mathrm{rank} \ A = r$ とすると, 適当な正則行列 $P$ , $Q$ によって,

とできる. ただし, $I_r$ は $r \times r$ の単位行列で, それぞれの $O$ は適当なサイズの全ての成分が $0$ の行列である. $A^- := QR'P$ と定めると,

$AA^-A = (P^{-1}RQ^{-1}) (QR'P) (P^{-1}RQ^{-1}) = P^{-1}RR'RQ^{-1} = P^{-1}RQ^{-1} = A. \label{}$

よって, $A^{-}$ は一般化逆行列である. (証明終わり)

一般化逆行列は一意的とは限らない. 例えば,

とするとき, 任意の実数 $a,b,c$ に対し,

は一般化逆行列である. ただし明らかなように, $A$ が正方行列で正則なときには, 通常の逆行列が唯一の一般化逆行列である.

以下では行列 $A$ に対し, 一般化逆行列をどれでもいいので一つ固定し, $A^{-}$ と書く. 特別な一般化逆行列を取る必要があるときはその都度言及する.

命題 4 行列 $A$ に対し, 以下が成立:

$\mathrm{rank} \ A^-A = \mathrm{rank} \ A$ .
$\mathrm{Im}\ A' = \mathrm{Im} \ A'A$ .

特に, $\mathrm{rank} \ A'A = \mathrm{rank} \ A' = \mathrm{rank} \ A$ .

証明

(1): $\mathrm{rank} \ A = \mathrm{rank} \ AA^-A \leq \mathrm{rank} \ A^-A \leq \mathrm{rank} \ A$ よりわかる.

(2): ベクトル $x$ に対し, $A'Ax = 0$ なら, $0 = x A'Ax = |Ax|^2$ より, $x \in \mathrm{Ker} \ A$ . よって, $\mathrm{Ker} \ A = \mathrm{Ker} \ A'A$ なので, $\mathrm{Im} \ A' = (\mathrm{Ker} \ A)^{\perp} = (\mathrm{Ker} \ A'A)^{\perp} = \mathrm{Im} \ A'A. \label{}$
よって, 証明できた. (証明終わり)

一般化逆行列を使うと, 連立一次方程式の解をシンプルに表せる.

命題 5 $A$ を $m \times n$ 行列, $y \in \mathbb{R}^m$ とし, $y \in \mathrm{Im} \ A$ を仮定する. このとき, 連立一次方程式

$Ax = y \label{}$

の解の集合は

$\{ A^-y + (I-A^-A)u \mid u \in \mathbb{R}^n \} \label{}$

である.

証明

$y \in \mathrm{Im} \ A$ より, ある $z \in \mathbb{R}^n$ が存在し, $y = Az$ . よって, 任意の $u \in \mathbb{R}^n$ に対し,

$A(A^-y + (I-A^-A)u) = AA^-Az + (A - AA^-A)u = Az = y. \label{}$

よって, $A^-y + (I-A^-A)u$ は $Ax = y$ の解である. 特に, $u = 0$ とすれば, $A^-y$ は特殊解であるから, あとは $I-A^-A$ が $\mathrm{Ker} \ A$ への射影であることを示せばよい.

$(I-A^-A)^2 = I-A^-A$ , $\mathrm{Im} \ (I-A^-A) \subset \mathrm{Ker} \ A$ は明らかである. よって, $\mathrm{rank} \ (I-A^-A) = \mathrm{dim} \ (\mathrm{Ker} \ A)$ を示せば十分である. $I-A^-A$ は射影なので, $n = \mathrm{rank} \ A^-A + \mathrm{rank} \ (I - A^-A)$ . よって, 命題4より,

$\mathrm{rank} \ (I-A^-A) = n - \mathrm{rank} \ A^-A = n - \mathrm{rank}A = \mathrm{dim} \ (\mathrm{Ker} \ A). \label{}$

よって, 証明できた. (証明終わり)

一般化逆行列は, 射影に関するシンプルな表示も与える.

命題 6 行列 $A$ に対し, $A(A'A)^-A'$ は $\mathrm{Im} \ A$ への射影. また $(A'A)^-$ が対称なら, $A(A'A)^-A$ は直交射影.

証明

$P := A(A'A)^-A'$ とおく. このとき, $\mathrm{Im} \ A' = \mathrm{Im} \ A'A$ より, 任意の $x$ に対し, ある $z$ が存在し, $A'x = A'Az$ . これより, $P^{2}x = A(A'A)^{-}A'A(A'A)^{-}A'Az = A(A'A)^{-}A'Az = Px$ . よって, $P^2 = P$ . 後半の主張は自明. (証明終わり)

ガウス・マルコフモデル

ガウス・マルコフモデルとは, 観測値がいくつかの説明変数の線型結合と誤差との和によって与えられるとするモデルである. より具体的には, $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$ は仮定しない.

$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{}$

尤度比検定とF検定の等価性

帰無仮説 $H_0$ , 対立仮説 $H_1$ の尤度比検定から検定量を導出することを考える. つまり, 以下の検定統計量を考える:

$\lambda(Y) := \frac{\sup_{\beta: H'\beta = d, \sigma^2} L(\beta,\sigma^2|Y)}{\sup_{\beta, \sigma^2} L(\beta,\sigma^2|Y)}. \label{eq05}$

ここで, $L$ は尤度関数. つまり,

$L(\beta,\sigma^2|Y) := \frac{1}{(2\pi\sigma^2)^{n/2}} \mathrm{e}^{-\frac{|Y-X\beta|^2}{2\sigma^2}}. \label{}$

まず, 分母は制約条件の無い尤度の最大化であるから, 例えば, 偏微分して極値を求めることで容易に以下がわかる:

$\sup_{\beta, \sigma^2} L(\beta,\sigma^2|Y) = \left( \frac{n}{2\pi \mathrm{e}} \right)^{n/2} |Y'(I - P_X)Y|^{-n/2}. \label{}$

ここで, $P_X$ は $X$ の像への直交射影. 後の便宜のため, $Y'(I - P_X)Y = (Y - X\beta)'(I-P_X)(Y-X\beta)$ であることを注意しておく. なぜなら, $(I-P_X)X = 0$ だからである.

よって, あとは分子の値を具体的に求めれば, 尤度比が求まる. これは制約条件付きの最大化問題であるからラグランジュの未定乗数法ですぐにできそうである. しかし, ラグランジュ関数の停留点を求める方程式の解き方が分からなかった. そこで, ここではより素朴に, 制約条件を満たす $\beta$ の範囲を先に求め, その中で(制約のない)通常の最大化問題を考える. $H'\beta = d$ を満たす $\beta$ を $1$ つ固定し, $\gamma$ とする.

定理 1 以下が成り立つ:

$\sup_{\beta: H'\beta = d, \sigma^2} L(\beta,\sigma^2|Y) = \left( \frac{n}{2\pi \mathrm{e}} \right)^{n/2} |(Y-X\gamma)'(I - P_{X|H})(Y-X\gamma)|^{-n/2}$

ただし, $P_{X|H}$ は直交射影で, 像は

$\{X\alpha \mid \alpha \in \mathrm{Ker} \ H'\}. \label{eq07}$

また, $\mathrm{rank} \ P_{X|H} = r - k$ が成り立つ.

証明

$\mathrm{rank}\ H' = k$ より, $\mathrm{Ker} \ H'$ の基底 $g_1,\dots,g_{p-k}$ を $1$ つ固定し, それらを並べた行列を $G := (g_1, \dots, g_{p-k})$ とすると, 連立一次方程式 $H'\beta = d$ の任意の解は $\gamma + G\delta \ (\delta \in \mathbb{R}^{p-k})$ と書ける. よって,

$\sup_{\beta: H'\beta = d} L(\beta,\sigma^2|Y) = \sup_{\delta, \sigma^2} \frac{1}{(2\pi\sigma^2)^{n/2}} \exp \left(-\frac{|Y-X\gamma - XG \delta|^2}{2\sigma^2}\right). \label{}$

右辺は制約条件の無い関数の最大値である. 右辺の $\sup$ の中身の関数の対数を $\ell$ と書き, $\frac{\partial}{\partial \delta}\ell = 0$ とすると,

$G'X'XG \delta = G'X'(Y-X\gamma). \label{}$

命題4より, 右辺は $\mathrm{Im} \ G'X'XG$ に含まれるので, これを満たす $\delta$ は,

$(G'X'XG)^-G'X'(Y-X\gamma) + (G'X'XG)^-(G'X'XG)u \quad (u \in \mathbb{R}^{p-k})$

と書ける. ただし, $(G'X'XG)^-$ は $G'X'XG$ の一般化逆行列の $1$ つで対称なものとする. 対称なものが取れることは例えば, $G'X'XG$ のスペクトル分解を $a_1p_1 + \cdots a_mp_m \ (a_i > 0, \ p_i$ は互いに直交する直交射影) とするとき, $(G'X'XG)^- = a_1^{-1}p_1 + \cdots a_m^{-1}p_m$ ととれることからわかる. $P := XG(G'X'XG)^-G'X'$ と定めると, 命題6より, $P$ は $\mathrm{Im} \ XG$ への直交射影であり, $u \in \mathbb{R}^{p-k}$ に依らず,

$XG\delta = P(Y-X\gamma) \label{}$

である.

次に, $\frac{\partial}{\partial \sigma^2}\ell = 0$ を解くと,

$\sigma^2 = \frac{(Y-X\gamma)'(I-P)(Y-X\gamma)}n. \label{}$

よって,

$\sup_{\beta: H'\beta = d, \sigma^2} L(\beta,\sigma^2|Y) = \left( \frac{n}{2\pi \mathrm{e}} \right)^{n/2} |(Y-X\gamma)'(I - P)(Y-X\gamma)|^{-n/2} \label{}$

がわかった.

あとは, $P = P_{X|H}$ , $\mathrm{rank} \ P = r - k$ を示せばよい. $P$ と $(G'X'XG)^-(G'X'XG)$ はべき等であること, 及び, 命題4から,

$P$ の定義から明らかに $\mathrm{Im} \ P \subset \mathrm{Im} \ XG$ であることから, $\mathrm{Im}\ P = \mathrm{Im} \ XG$ がわかる. $G$ は $\mathrm{Ker} \ H'$ の基底を並べた行列なので, $\mathrm{Ker} \ H' = \mathrm{Im} \ G$ で, $\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{}$

これより, $P = P_{X|H}$ であることが分かる. また, $\mathrm{rank} \ XG = \mathrm{rank} \ G - \dim (\mathrm{Ker} \ X) = (p-k) - (p-r) = r - k$ であるから, $\mathrm{rank} \ P = r -k$ . (証明終わり)