一元配置分散分析(one-way ANOVA)をわかりやすく解説

この記事では、一元配置分散分析(one-way ANOVA)をわかりやすく解説します。
適当にネットで調べてもあまりしっくりくる説明がなかったので、自分用の備忘録も兼ねています。
前提知識としては、重回帰分析の線形代数による記述スタイルを掴んでいたらよいと思います。

問題意識(簡単のために3群で)

経済学専門の予備校で100点満点の経済学の試験をしました。この予備校にはA・B・Cの3クラスがあります。
それぞれのクラスから2人ずつ試験結果をみてみると、次の表のような結果となっていました。
（話を簡単にするために便宜上2人ずつとしましたが、同じ人数なら何人でもよいです）。

クラス	試験の点数
A	13, 21
B	14, 37
C	21, 51

ここで、次のような仮定を置くことにします。
Aクラスの学生の試験結果は\(N(\mu_A, \sigma^2)\)に従って発生
Bクラスの学生の試験結果は\(N(\mu_B, \sigma^2)\)に従って発生
Cクラスの学生の試験結果は\(N(\mu_C, \sigma^2)\)に従って発生
つまり、適当な平均で、分散が等しい正規分布から発生していると思うことにします。

つまりこれは、
\begin{align} \begin{pmatrix} y_1 \\ y_2 \\ y_3 \\ y_4 \\ y_5 \\ y_6 \end{pmatrix} =\begin{pmatrix} 1 & 0 & 0 \\ 1 & 0 & 0 \\0 & 1 & 0 \\ 0 & 1 & 0 \\0 & 0 & 1 \\ 0 & 0 & 1\end{pmatrix} \begin{pmatrix} \beta_1 \\ \beta_2 \\ \beta_3 \end{pmatrix} + \begin{pmatrix} \varepsilon_1 \\ \varepsilon_2 \\ \varepsilon_3 \\ \varepsilon_4 \\ \varepsilon_5 \\ \varepsilon_6 \end{pmatrix}, \quad \varepsilon_i \sim N(0, \sigma^2), \text{i.i.d.} \end{align}
という線形回帰モデル(正規誤差項つき)を立てていることと同じだと思えます。
（ただし、\(y_1, y_2\)がAクラスの学生の試験結果で、\(y_3, y_4\)がBクラスの学生の試験結果で、\(y_5, y_6\)がCクラスの学生の試験結果ということにします）。
つまり、ここではA・B・Cクラスの平均は
\begin{align*} \beta_1, \beta_2, \beta_3 \end{align*}
であるというモデルになっています。ここで、次の分析をかけることにします。

「A・B・Cクラスの平均は同じであるか？」。つまり、
\begin{align*} \beta_1 = \beta_2 = \beta_3 \end{align*}
かをチェックしてみたいとします。

そこで、帰無仮説を
\begin{align*} h_0 : \beta_1 = \beta_2 = \beta_3 = \beta \end{align*}
として、何かうまくいきそうな検定量を探すことにします。
\begin{align*} X_f = \begin{pmatrix} 1 & 0 & 0 \\ 1 & 0 & 0 \\0 & 1 & 0 \\ 0 & 1 & 0 \\0 & 0 & 1 \\ 0 & 0 & 1\end{pmatrix} \end{align*}
として、
\begin{align*} \varepsilon = \begin{pmatrix} \varepsilon_1 \\ \varepsilon_2 \\ \varepsilon_3 \\ \varepsilon_4 \\ \varepsilon_5 \\ \varepsilon_6 \end{pmatrix} \end{align*}
と表記することにすると、
モデルは、
\begin{align*} y = X_f \begin{pmatrix} \beta \\ \beta \\ \beta \end{pmatrix} + \varepsilon \end{align*}
と思っているということになります。
結局はランダムな要素というのは誤差項に集約されているので、うまく誤差だけになるような形を作りにいくことにします。
として、ハット行列を
\begin{align*} H_f = X_f \left(X_f ^t X_f \right)^{-1} X_f^t \end{align*}
と定めると、OLS推定量による予測値は
\begin{align*} \hat y = H_f y\end{align*}
とかけたことを思い出しておきます。また、
\begin{align*} X_r = \begin{pmatrix} 1 \\ 1 \\ 1 \\ 1 \\ 1 \\1 \end{pmatrix} \end{align*}
と定めると、モデルは
\begin{align*} y = X_r \beta + \varepsilon \end{align*}
とも書くことができます。この場合のハット行列は、
\begin{align*} H_r = X_r \left(X_r ^t X_r \right)^{-1} X_r^t\end{align*}
と書くことにします。
（オリジナルのモデルのOLS推定量を求めた後に\(\beta_1 = \beta_2 = \beta_3 = \beta\)とするのと、退化したモデルでのOLS推定量は同じではない(はず)であることを補足しておきます）。
天下り的にですが（そもそも本当はそうする必要もないと思うのですが）、それぞれの表記での残差平方和をとります。行列の表記を用いて、
\begin{align*} RSS_f = y^t(I- H_f ) y, \quad RSS_r = y^t (I – H_r) y \end{align*}
を考えることにします。
\begin{align*} \beta_r = \begin{pmatrix}\beta \\ \beta \\ \beta \end{pmatrix} \end{align*}
と表記することにすると、
\begin{align*} y^t H_f y &= \beta_f X_f ^t H_f X_f \beta _f + \varepsilon^t H_f X_f \beta_f + \beta_f^t X_f^t H_f \varepsilon + \varepsilon^t H_f \varepsilon
\\&= \beta_f X_f ^t X_f \beta _f + \varepsilon^t X_f \beta_f + \beta_f^t X_f^t \varepsilon + \varepsilon^t H_f \varepsilon
\\&= \beta_r X_r ^t X_r \beta _r + \varepsilon^t X_r \beta_r + \beta_r^t X_r^t \varepsilon + \varepsilon^t H_f \varepsilon \end{align*}
です。ただし\(H_f X_f = X_f\)や\(X_f \beta_f = X_r \beta_r\)をところどころ用いています。
また同様に、
\begin{align*} y^t H_r y &= \beta_r X_r ^t H_r X_r \beta _r + \varepsilon^t H_r X_r \beta_r + \beta_r^t X_r^t H_r \varepsilon + \varepsilon^t H_r \varepsilon
\\&= \beta_r X_r ^t H_r X_r \beta _r + \varepsilon^t H_r X_r \beta_r + \beta_r^t X_r^t H_r \varepsilon + \varepsilon^t H_r \varepsilon
\\&= \beta_r X_r ^t X_r \beta _r + \varepsilon^t X_r \beta_r + \beta_r^t X_r^t \varepsilon + \varepsilon^t H_r \varepsilon \end{align*}
となります。というわけで、
\begin{align*}y^t \left(H_f – H_r \right) y
&= \varepsilon^t \left(H_f – H_r \right) \varepsilon \end{align*}
となり、誤差項の部分を切り出すことができました（ハット行列と共にですが）。
\begin{align} \frac{1}{\sigma} \varepsilon \end{align}
を考えることにします。つまり、
\begin{align} \frac{1}{\sigma} \varepsilon \sim N(0,I)\end{align}
ですので、多変量標準正規分布です。
\(\left(H_f – H_r \right)\)は階数3 – 1 = 2の実対称冪等行列であることがわかります。
適当に直交行列\(S \in O_n \)で対角化することで、
\begin{align} \left(H_f – H_r \right) = S^t \Lambda_{2}S\end{align}
とすることができます。ただし、
\begin{align} \Lambda_2 = \begin{bmatrix}I_{2} & 0 \\0 & 0 \end{bmatrix} \end{align}
という表記を用いています。
\begin{align} \frac{1}{\sigma} \varepsilon^t \left(H_f – H_r \right) \frac{1}{\sigma} \varepsilon &= \frac{1}{\sigma} \varepsilon^t S^t \Lambda_{2} S \frac{1}{\sigma} \varepsilon \end{align}
とできます。
\begin{align} S \frac{1}{\sigma} \varepsilon \end{align}
は多変量標準正規分布の直交行列による変換なので多変量標準正規分布に従います。
というわけで、
\begin{align}\frac{1}{\sigma} \varepsilon^t S^t \Lambda_{2} S \frac{1}{\sigma} \varepsilon \end{align}
は標準正規分布の2乗の\(2\)個の和なので自由度\(2\)のカイ二乗分布に従います。
というわけで、
\begin{align*} \frac{1}{\sigma^2}RSS_f – \frac{1}{\sigma^2}RSS_r = \frac{1}{\sigma^2}y^t (H_f – H_r) y \sim \chi_{2}^2\end{align*}
がわかります。
また、同様の議論から、
\begin{align*} \frac{1}{\sigma^2}RSS_f \sim \chi_{6-3}^2 \end{align*}
であることもわかります。

続いて、
\(H_f\)の具体的な形を見ることにします。
\begin{align*} \begin{pmatrix} 1 & 0 & 0 \\ 1 & 0 & 0 \\0 & 1 & 0 \\ 0 & 1 & 0 \\0 & 0 & 1 \\ 0 & 0 & 1\end{pmatrix} ^t = \begin{pmatrix}
1 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 1
\end{pmatrix}\end{align*}
なので、
\begin{align*} X_f ^t X_ f &= \begin{pmatrix}
1 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 1
\end{pmatrix} \begin{pmatrix} 1 & 0 & 0 \\ 1 & 0 & 0 \\0 & 1 & 0 \\ 0 & 1 & 0 \\0 & 0 & 1 \\ 0 & 0 & 1\end{pmatrix} \\&= \begin{pmatrix} 2& 0 & 0 \\ 0 & 2 & 0 \\ 0 & 0 & 2 \end{pmatrix}\end{align*}
で、
\begin{align*} X_f X_f^t &= \begin{pmatrix} 1 & 0 & 0 \\ 1 & 0 & 0 \\0 & 1 & 0 \\ 0 & 1 & 0 \\0 & 0 & 1 \\ 0 & 0 & 1\end{pmatrix} \begin{pmatrix}
1 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 1
\end{pmatrix} \\&= \begin{pmatrix} 1 & 1 & 0 & 0& 0& 0 \\ 1 & 1 & 0 & 0& 0& 0 \\ 0 & 0 & 1 & 1 & 0& 0 \\ 0 & 0 & 1 & 1 & 0& 0 \\ 0 & 0 & 0 & 0 & 1& 1 \\ 0 & 0 & 0 & 0 & 1& 1 \end{pmatrix} \end{align*}
なので、
\begin{align*} H_f = X_f \left(X_f^t X_f \right)^{-1} X_f^t = \frac{1}{2} \begin{pmatrix} 1 & 1 & 0 & 0& 0& 0 \\ 1 & 1 & 0 & 0& 0& 0 \\ 0 & 0 & 1 & 1 & 0& 0 \\ 0 & 0 & 1 & 1 & 0& 0 \\ 0 & 0 & 0 & 0 & 1& 1 \\ 0 & 0 & 0 & 0 & 1& 1 \end{pmatrix} \end{align*}
となります。また、
\begin{align*} H_r = X_r \left(X_r^t X_r \right)^{-1} X_r^t = \frac{1}{6} \begin{pmatrix} 1& 1 & 1 & 1 & 1 \\ 1& 1 & 1 & 1 & 1 \\ 1& 1 & 1 & 1 & 1 \\ 1& 1 & 1 & 1 & 1 \\ 1& 1 & 1 & 1 & 1 \end{pmatrix} \end{align*}
となります。

ですので、それぞれのクラスの試験結果の標本平均を
\begin{align*} \bar y_A = \frac{1}{2} (y_1 + y_2), \quad \bar y_B = \frac{1}{2} (y_3 + y_4), \quad \bar y_C = \frac{1}{2} (y_5 + y_6) \end{align*}
と表記することにすると、
\begin{align*} H_f y = \begin{pmatrix} \bar y_A \\ \bar y_A \\ \bar y_B \\ \bar y_B \\ \bar y_C \\ \bar y_C \end{pmatrix} \end{align*}
ですし、
\begin{align*} \bar y = \frac{1}{6}(y_1 + y_2 + y _3 + y_4 + y _5 + y_6)\end{align*}
とすると、
\begin{align*} H_r y = \begin{pmatrix} \bar y \\ \bar y \\ \bar y \\ \bar y \\ \bar y \\ \bar y \end{pmatrix} \end{align*}
であることがわかります。
従って、
\begin{align*} RSS_f &= (y_1 – \bar y)^2 + (y_2 – \bar y)^2 + (y_3 – \bar y)^2 + (y_4 – \bar y)^2 + (y_5 – \bar y)^2 + (y_6 – \bar y)^2 \\ RSS_r &= (y_1 – \bar y_A)^2 + (y_2 – \bar y_A)^2 + (y_3 – \bar y_B)^2 + (y_4 – \bar y_B)^2 + (y_5 – \bar y_C)^2 + (y_6 – \bar y_C)^2 \end{align*}
であることがわかります。

未知かもしれない\(\sigma^2\)を適当に消すことにします。
\begin{align*} RSS_f = y^t (I – H_f) y, \quad RSS_f – RSS_r = y^t (H_f – H_r) y \end{align*}
は違いに独立であることに注意すると、
\begin{align*} \frac{1}{\sigma^2} (RSS_f – RSS_r) \sim \chi_{3-1}^2 , \frac{1}{\sigma^2}RSS_f \sim \chi_{6-3}^2 \end{align*}
なので、F分布\(F_{6-3}^{3-1}\)に従う量
\begin{align*} \frac{RSS_f – RSS_r}{RSS_f} \sim F_{6-3}^{3-1}\end{align*}
が得られます。

つまり結論としては、