はじめに
分散分析は、ある処理が対象に与えている影響の有無を分散の観点から明らかにする手法で、アカデミアのみならず様々な場所で必要不可欠な手法となっています。この分散分析の基礎を築いたのは生物統計学の父Fisherであり、私の所属する生物測定学研究室(生物統計学を主なツールとして用いる)においても必要不可欠な方法です。
さて、この分散分析ですが、体系的に習ったのは大学に入ってからだったかと思います。習った当初の印象は「よくわからない」でした。分散分析の授業では、一元配置分散分析や、二元配置分散分析、交互作用を含む分散分析、乱解法/完全無作為法、分割区法、ブロックの導入、など特定のモデルごと分散分析のやり方を習います。分散分析のやり方というのは、結局のところ、モデルに含まれる各要素(処理)における自由度と二乗和を計算し、そこから平均平方和を求め、誤差における平均平方和とともに検定をすることにほかならず、それは特に問題なく理解することができました。
私がよくわからなかったのは、自由度と二乗和の計算方法です。自由度に関しては、(n-1)のように1がなくなったり、k(n-1)のように片方は1がなくなるのにもう片方はなくならなかったり、(k-1)(n-1)のように両方なくなったりするわけですが、その規則性がよくわかりませんでした。その規則性を考えようと、自由度が減る理由を調べると、「$\sum_{i=1}^{n}\beta_i=0$という制限が入るがら」などというふわっとした根拠しか出てこず、数学的になぜそのような制限が必要なのかなどには答えてくれませんでした。二乗和の計算についても、同様にその計算方法の規則性がよくわかりませんでした。確かになんとなくそうなりそうな気にはなるのですが、完全に未習のモデルが出てきて、処理ごとの二乗和を予想して書きなさい、と言われたときに解ける自信はなかなかありません。結局のところ、日本における分散分析の授業では、ケーススタディとして個別のモデルの解き方しか教えてくれないため、一般化するのが非常に難しく、私はそのことに非常に苦しめられました。
私の、これらのモヤモヤをすっきりさせてくれたのが、たまたま受け始めたアメリカの教授による分散分析の授業でした。その授業では分散分析の各モデルを紹介しつつも、全モデルを線形モデル($\mathbf{y}=X\mathbf{\beta}+\mathbf{e}$)として統一的に扱い、分散分析の一般的化した解法を教えてくれました。これを受けたことにより、私がそれまで分散分析に対して抱いてきたモヤモヤがすっきり晴れ、かなり分散分析への分解能もあがりましたので、まとめも兼ねて今回ご紹介できればと思います。なお、以下では行列の基本的な知識は前提とします。
目次
※ 以下のサイトは分散分析をちゃんと数学的に説明しようとしているサイトです。非常に良い投稿だと思いましたので紹介させていただきます。
分散分析の数学的理解とpythonでの実装方法を理解していく
残差平方和と水準間平方和の期待値の導出方法(1元配置分散分析)
分散分析のモデルと推定可能関数
はじめに、基礎的ではありますが、どの分散分析のモデルも線形モデル($\mathbf{y}=X\mathbf{\beta}+\mathbf{e}$)として書けることを見ていきます。詳しくは私が以前書いた記事にまとめてありますので参考になさってください( 分散分析における推定可能関数 )。
まず次のような一元配置分散分析のモデルを考えます。
$$
y_{ij}=\mu+\alpha_i+e_{ij}
$$
ここで、$y_{ij}$は$i$番目の水準における$j$番目の計測値、$\mu$は全体平均、$\alpha_i$は$i$番目の水準の効果、$e_{ij}$は誤差を表します。さて、これらは行列を意識して書くと次のように書き替えることが可能です。
\begin{pmatrix} y_{11} \\ \vdots \\ y_{1n_1} \\ \vdots \\ y_{a1} \\ \vdots \\ y_{an_a} \end{pmatrix} = \begin{pmatrix} 1 & 1 & & \\ \vdots & \vdots & O & \\ 1 & 1 & & \\ \vdots & & \ddots & \\ 1 & & & 1\\ \vdots & & O & \vdots \\ 1 & & & 1 \end{pmatrix} \begin{pmatrix} \mu \\ \alpha_1 \\ \vdots \\ \alpha_a \end{pmatrix} + \begin{pmatrix} e_{11} \\ \vdots \\ e_{1n_1} \\ \vdots \\ e_{a1} \\ \vdots \\ e_{an_a} \end{pmatrix}
ここで、ベクトル$\mathbf{y},\beta,\mathbf{e}$、行列$X$を次のように定義します。
\mathbf{y}=\begin{pmatrix} y_{11} \\ \vdots \\ y_{1n_1} \\ \vdots \\ y_{a1} \\ \vdots \\ y_{an_a} \end{pmatrix}, \beta=\begin{pmatrix} \mu \\ \alpha_1 \\ \vdots \\ \alpha_a \end{pmatrix} ,\mathbf{e}=\begin{pmatrix} e_{11} \\ \vdots \\ e_{1n_1} \\ \vdots \\ e_{a1} \\ \vdots \\ e_{an_a} \end{pmatrix} , X=\begin{pmatrix} 1 & 1 & & \\ \vdots & \vdots & O & \\ 1 & 1 & & \\ \vdots & & \ddots & \\ 1 & & & 1\\ \vdots & & O & \vdots \\ 1 & & & 1 \end{pmatrix}
これらを用いると上記のモデルは以下のようにまとめられます。
$$
\mathbf{y}=X\mathbf{\beta}+\mathbf{e}
$$
これは$\beta$に対して線形ですので線形モデルと呼ばれます。
次に、以下のような二元配置分散分析のモデルを考えます。
$$
y_{ijk}=\mu+\alpha_i+\beta_j+\gamma_{ij}+e_{ijk}
$$
\begin{pmatrix} y_{111} \\ \vdots \\ y_{11n} \\ \vdots \\ y_{1m1} \\ \vdots \\ y_{1mn} \\ \vdots \\ y_{lmn} \end{pmatrix} = \begin{pmatrix} 1 & 1 & & \\ \vdots & \vdots & O & \\ 1 & 1 & & \\ \vdots & & \ddots & \\ 1 & & & 1\\ \vdots & & O & \vdots \\ 1 & & & 1 \end{pmatrix} \begin{pmatrix} \mu \\ \alpha_1 \\ \vdots \\ \alpha_l \\ \beta_1 \\ \vdots \\ \beta_m \\ \gamma_{11} \\ \vdots \\ \gamma_{lm} \end{pmatrix} + \begin{pmatrix} e_{111} \\ \vdots \\ e_{11n} \\ \vdots \\ e_{1m1} \\ \vdots \\ e_{1mn} \\ \vdots \\ e_{lmn} \end{pmatrix}
こちらも、$\beta=(\mu, \alpha_1, ...,\alpha_n, \beta_1,...,\beta_m, \gamma_{11},...,\gamma_{1m},\gamma_{21},...,\gamma_{nm})^T$とすることで線形モデル($\mathbf{y}=X\mathbf{\beta}+\mathbf{e}$)として書くことができると分かります。
このように分散分析で扱うモデルは基本的にこのような線形モデルとして書くことができるのです。このような線形モデルは回帰分析の最小二乗法を習ったことがある方ですとなじみが深いのではないでしょうか。一方で、分散分析における線形モデルは、回帰分析とは決定的に違う点があります。それが$X$の構造です。分散分析において$\mathbf{y}$の要素は必ずいずれかの水準に属します。例えば最初に扱った一元配置分散分析$y_{ij}=\mu+\alpha_i+e_{ij}$の場合、$\mathbf{y}$の要素はどれでも$\alpha_1,...\alpha_a$のいづれかに属します。これは$X$の2~(a+1)列目のベクトルの和がXの1列目と一致することに他なりません。回帰分析では$X$の要素には様々な実数が入るためこのようなことは起こりませんから、分散分析特有の現象と言えるでしょう。この現象をもう少し数学的にいうと、$X$の列ベクトルは線形独立になっていないということであり、このことから行列$X$のランクが列数に一致しない(ランク落ち)ことが分かります。
一般に、線形モデル$\mathbf{y}=X\mathbf{\beta}+\mathbf{e}$の解は$(X^TX)\hat{\beta}=X^T\mathbf{y}$を解くことにより求まりますが、$X$がランク落ちする場合、$X^TX$の逆行列が存在せず、したがって解の値が一意に定まりません。つまり分散分析では$\beta$を正しく求めることができないのです。とても興味深いことに、$X\beta$をはじめとする$\beta$の線形結合$\mathbf{k}\beta$で表される一部は、たとえ$\beta$が一意に定まらなくとも、値が一意に定まります。そしてこのような一意に値を定められるものは推定可能関数と呼ばれます。推定可能関数はどんな$\beta$の推定値を取ってきてもよいわけですから、非常に便利なわけです。例えば、正方行列$A$に対し、$AGA=A$を満たす正方行列$G$は一般化逆行列と呼ばれ$G^{-}$で表されますが、$\hat{\beta}=(X^TX)^{-}X^T\mathbf{y}$は$\beta$の解の一つ(正確には解空間の一部)をなすことが知られています。なお、推定可能関数である条件は、$M(X)$を$X$の列ベクトルで張られるベクトル空間としたとき、$k\in M(X)$であることです。詳しくは実験計画法と分散分析(朝倉書店)などをご覧ください。
かなり本題から脱線してしまいましたが、分散分析ではどのモデルも線形モデルとしてかけること、分散分析の線形モデルの特徴として$X$がランク落ちすることなどを覚えていていただければと思います。
分散分析における平方和
次にモデル平方和と残差平方和を導入します。
線形モデル$\mathbf{y}=X\mathbf{\beta}+\mathbf{e}$に対して、モデル平方和$
S(\mathbf{\hat{\beta}})$を次のように定義します。
$$
S(\mathbf{\hat{\beta}})=(X\hat{\beta})^TX\hat{\beta}
$$
$\hat{\beta}$の関数になっているため、先ほどの話からして一意に定まらないんじゃないの?と思う方もいらっしゃると思いますが、推定可能関数$X\hat{\beta}$の関数でもあるため、特に問題にはなりません。そして、残差平方和$R(\hat{\beta})$を次のように定義します。
$$
R(\mathbf{\hat{\beta}})=(y-X\hat{\beta})^T(y-X\hat{\beta})
$$
ここで、先ほど少し出てきた$\hat{\beta}=(X^TX)^{-}X^T\mathbf{y}$を$S(\mathbf{\hat{\beta}}),R(\mathbf{\hat{\beta}})$に代入しまとめると、次が得られます。
\begin{align}
S(\mathbf{\hat{\beta}})&=\mathbf{y}^T(X(X^TX)^{-}X^T)\mathbf{y}\\
R(\mathbf{\hat{\beta}})&=\mathbf{y}^T(I-X(X^TX)^{-}X^T)\mathbf{y}
\end{align}
これはどちらも$\mathbf{y}$の二次形式の形をしており、$\mathbf{y}$が正規分布に従う場合、これはカイ二乗分布に従います。カイ二乗分布は自由度によって特徴づけられますが、その自由度を決定しているのが、二次形式の真ん中の行列のランクです。頭で扱った一元配置分散分析$y_{ij}=\mu+\alpha_i+e_{ij}$の場合、$X$のランクは一つ落ちて$(a+1)-1=a$になりますので、$S(\mathbf{\hat{\beta}})$の自由度は$a$です。(説明すると大変なので少し説明をごまかしましたが、$rank(X(X^TX)^{-}X^T)=rank(X)$であることを使っています)。また、各水準$\alpha_i$の反復数が一律$n$である場合、$\mathbf{y}$は$na$の長さを持ちますので$I$の次元は$na×na$になり、従って$rank(I-X(X^TX)^{-}X^T)=rank(I)-rank(X(X^TX)^{-}X^T)=na-a=a(n-1)$となります。従って$R(\mathbf{\hat{\beta}})$の自由度は$a(n-1)$となります。
鋭い方は$S(\mathbf{\hat{\beta}})$のランクが$a-1$ではなく$a$であったことに違和感を覚えるかもしれません。これについては次の章で説明しています。
ANOVA表との対応
ここからは、先ほど定義したモデル平方和や残差平方和とANOVA表との関連を調べていきます。
一元配置分散分析
まずは一元配置分散分析を考えます。
例えば、一元配置分散$y_{ij}=\mu+\alpha_i+e_{ij}$の場合は以下のようなANOVA表が作られます(本当はこのほかにもF値などがかかれるわけですが今回は簡単のため省略しています)。
| 処理 | 自由度 | 平方和 |
|---|---|---|
| 処理A | a-1 | $\sum \sum (y_{i.}-y_{..})^2$ |
| 誤差 | a(n-1) | $\sum \sum (y_{ij}-y_{i.})^2$ |
| 合計 | an-1 | $\sum \sum (y_{ij}-y_{..})^2$ |
このANOVA表における自由度、および二乗和がなぜこのような形で表されるのか、を上記で定義したモデル平方和と残差平方和を用いながら説明していきます。
誤差の平方和
まず最も簡単なのが、一番下段で定義される誤差に対する二乗和で、これは先ほど定義した$R(\mathbf{\hat{\beta}})$です。$R(\mathbf{\hat{\beta}})$の自由度は先ほど計算した通り$a(n-1)$であり、確かにANOVA表と一致することが分かります。次に$R(\mathbf{\hat{\beta}})=\sum \sum (y_{ij}-y_{i.})^2$となることを示します。まず$R(\mathbf{\hat{\beta}})=y^T(I-X(X^TX)^{-}X^Ty)$ですので、$R(\mathbf{\hat{\beta}})=y^Ty-y^TX(X^TX)^{-}X^Ty$が成り立ちます。ここで右辺第二項に対応する$y^TX(X^TX)^{-}X^Ty$ですが、$X(X^TX)^{-}X^T=H$と置くと$H$は次のようになることがわかります(詳しくはAppendixを参照)。
H_{ij} = \left\{
\begin{array}{ll}
\frac{1}{n} & 観測値iとjが同じ処理 \\
0 & その他
\end{array}
\right.
よって、$y^THy=y^T(y_{1.},y_{1.},..., y_{n.})^T=\sum_i^a ny_{i.}^2$となり、結果的に$R(\mathbf{\hat{\beta}})=\sum_i\sum_j y_{ij}^2-\sum_i^a ny_{i.}^2$が求まります。最後に$\sum \sum (y_{ij}-y_{i.})^2=\sum \sum y_{ij}^2+\sum ny_{i.}^2-2\sum \sum y_{ij}y_{i.}=\sum \sum y_{ij}^2-\sum ny_{i.}^2$が成り立ちますので、$R(\mathbf{\hat{\beta}})=\sum \sum (y_{ij}-y_{i.})^2$であることが確かめられました。
処理の平方和
わかりにくいのが処理$\alpha$に対する二乗和で、これは二つの残差平方和の差分として定義されます。一つ目が先ほど定義した$R(\hat{\beta})$で、もう一つが帰無仮説に基づくモデルの残差平方和です。今回は処理$\alpha$の効果があるかどうかを検定したいため、帰無仮説は$\alpha=0$となります。この時、一元配置分散分析のモデル$y_{ij}=\mu+\alpha_i+e_{ij}$は$y_{ij}=\mu+e_{ij}$に変化します。これに対応する線形モデルを$\mathbf{y}=X'\mathbf{\beta'}+\mathbf{e'}$としたとき、モデル平方和と残差平方和はそれぞれ$S(\hat{\beta'}),R(\hat{\beta'})$と書き表せます。処理$\alpha$に対応する平方和の正体はこれらの残差平方和の差$R(\mathbf{\beta'})-R(\mathbf{\beta})$です。定義より$R(\beta')=y^T(I-X'(X'^TX')^{-}X'^T)y$、$R(\beta)=y^T(I-X(X^TX)^{-}X^T)y$ですので、その差分は$R(\mathbf{\beta'})-R(\mathbf{\beta})=y^T(X(X^TX)^{-}X^T-X'(X'^TX')^{-}X'^T)y=S(\beta)-S(\beta')$となります。
これは$\mathbf{y}$の二次形式になっていますので、カイ二乗分布に従い、その自由度は$\mathbf{y}$で挟まれた行列のランクと一致します。帰無仮説に対応するモデルの計画行列$X'$は$X'=(1,1,...,1)^T$と等しいですからそのモデル平方和$S(\hat{\beta'})$の自由度は$1$、残差平方和$R(\hat{\beta'})$の自由度は$an-1$となります。一方、対立仮説のモデルの$S(\beta)$の自由度は上述したように$an-1$となりますので、処理$\alpha$に対応する平方和の自由度は$(an-1)-a(n-1)=a-1$となり、anova表と一致することが確認できました。
次に平方和$\sum \sum (y_{i.}-y_{..})^2$が$R(\mathbf{\beta'})-R(\mathbf{\beta})$に一致することを確認していきます。定義より$R(\mathbf{\beta'})-R(\mathbf{\beta})$は$S(\beta)-S(\beta')$と一致します。従って、今回示したいことは$S(\beta)-S(\beta')=\sum \sum (y_{i.}-y_{..})^2$と言い換えることができるでしょう。ここでは、$S(\beta)-S(\beta')=\sum\sum y_{i.}^2 - \sum \sum y_{..}^2$となること、$\sum \sum (y_{i.}-y_{..})^2= \sum\sum y_{i.}^2 - \sum \sum y_{..}^2$となることを示して証明を行います。
まず$S(\beta)-S(\beta')$についてです。帰無仮説に対応するモデル平方和$S(\beta')$は$\beta'=\mu$ですので、$S(\beta')=y^T(I-X'(X'^TX')^-X')y$において$X'=\mathbf{1}=(1,1,...,1)^T$となります。よって$S(\beta)-S(\beta')=y^T(X(X^TX)^-X^T-\mathbf{1}(\mathbf{1}^T\mathbf{1})^-\mathbf{1})y$となります。ここで、$y$で挟まれた行列の右辺は$y^T(\mathbf{1}(\mathbf{1}^T\mathbf{1})^-\mathbf{1} )y=\frac{1}{N}y^T\mathbf{1}\mathbf{1}^Ty=\frac{1}{N}(\mathbf{1}^Ty)^T(\mathbf{1}^Ty)=\frac{1}{N}(\sum_i y_i)^2=N\bar{y}^2$となります。次に$y$で挟まれた行列の左辺に対応する$y^TX(X^TX)^{-}X^Ty$ですが、先ほどの議論より$\sum_i^a ny_{i.}^2$となります。以上より$S(\beta)-S(\beta')=\sum_i^a ny_{i.}^2-N\bar{y}^2=\sum_i^a ny_{i.}^2-an\bar{y}^2=\sum \sum y_{i.}^2-\sum\sum y_{..}^2$となることが分かりました。
次に$\sum \sum (y_{i.}-y_{..})^2= \sum\sum y_{i.}^2 - \sum \sum y_{..}^2$となることの証明ですが、これは先ほどと同様に左辺の二乗の項を展開してやればよく、$\sum\sum (y_{i.}-y_{..})^2=\sum\sum y_{i.}^2+\sum\sum y_{..}^2-2\sum\sum y_{i.}y_{..}=\sum\sum y_{i.}^2-\sum\sum y_{..}^2$と求まります。ということで、長くなりましたが以上より$R(\mathbf{\beta'})-R(\mathbf{\beta})=\sum \sum (y_{i.}-y_{..})^2$を確認することができました。
全体の平方和
最後に全体ですが、これは$R(\beta)+R(\beta')-R(\beta)=R(\beta')$となりますが、$R(\beta')=y^Ty-S(\beta')$ですので$\sum_i (y_{ij}-y_{..})^2$が成り立ります(詳しく証明しませんがここまでの計算を振り返ればすぐ証明できるはずです)。
というわけで、anova表に出てきた平方和は、線形モデルにおける平方和、あるはその差であることが分かったかと思います。最後にそれをまとめて表にしたものを添付しておきます。
| 処理 | 自由度 | 平方和 | 線形モデル表現 |
|---|---|---|---|
| 処理A | a-1 | $\sum \sum (y_{i.}-y_{..})^2$ | $R(\beta')-R(\beta)$ |
| 誤差 | a(n-1) | $\sum \sum (y_{ij}-y_{i.})^2$ | $R(\beta)$ |
| 合計 | an-1 | $\sum \sum (y_{ij}-y_{..})^2$ | $R(\beta')$ |
Ｆ検定
検定について最後に少し補足しておきます。まず$R(\beta)-R(\beta')$を考えます。この二次形式はより詳しく書くと$y^T(X(X^TX)^{-}X^T-X'(X'^TX')^{-}X'^T)y$という二次形式から成り、これは自由度が$rank(X(X^TX)^{-}X^T-X'(X'^TX')^{-}X'^T)$のカイ二乗分布に従います。この時、非心度$\lambda$は$E[y]^T(X(X^TX)^{-}X^T-X'(X'^TX')^{-}X'^T)E[y]$になるわけですが、帰無仮説においては$E[y]=X'\beta'$であるため、$\lambda=\beta'^TX'^T(X(X^TX)^{-}X^T-X'(X'^TX')^{-}X'^T)X'\beta'=0$が成り立ちます(これはX'を射影行列X'(X'^TX')^{-}X'^T)X'で射影した結果と射影行列X(X^TX)^{-}X^Tで射影した結果が完全一致することから導けますが、説明すると長いので省きます。説明してほしい！という方はご連絡ください)。$R(\beta)$はもとより非心であるので、帰無仮説のもとで$R(\beta')-R(\beta),R(\beta)$ともに非心のカイ二乗分布に従います。非心のカイ二乗分布の商はF分布に従いますので$\frac{R(\beta')-R(\beta)}{R(\beta)}$は帰無仮説においては自由度$(a-1,a(n-1))$のF分布に従うというわけです。
さて、今までは一因子の単純な分散分析を考えてきましたが、以下ではもう少し複雑なモデルについて考えていきます。それらのモデルも残差平方和の差分で一般的に書き表せることを示すのが目的ですが、残差平方和の差分が一般的な平方和の式と一致することまで確認し始めると大変なことになるので、以下では自由度が一致することのみ示してきます
巣ごもりモデル
まずは果樹などで使われる巣ごもりモデル$y_{ijk}=\mu+a_i+b_{ij}+e_{ijk}$を考えます。これは、処理$b$の値が$a$の水準ごとに異なるモデルであり、$1<=i<=l, 1<=j<=m, 1<=k<n$の場合、ANOVA表は次のように書かれます。
| 処理 | 自由度 | 平方和 |
|---|---|---|
| 処理A | l-1 | $\sum \sum (y_{i.}-y_{..})^2$ |
| 処理B | l(m-1) | $\sum \sum (y_{i.}-y_{..})^2$ |
| 誤差 | lm(n-1) | $\sum \sum (y_{ij}-y_{i.})^2$ |
| 合計 | lmn-1 | $\sum \sum (y_{ij}-y_{..})^2$ |
処理Aの平方和
まず、処理$a$に対する平方和ですが、これはモデル$y_{ijk}=\mu+a_i+e_{ijk}$とモデル$y_{ijk}=\mu+e_{ijk}$の残差平方和の差分です。したがって一元配置分散分析のモデルの場合と同じで、自由度は$l-1$となります。
処理Bの平方和
次に処理$b$に対する平方和ですが、これはモデル$y_{ijk}=\mu+a_i+b_{ij}+e_ijk$と$y_{ijk}=\mu+a_i+e_{ijk}$の残差平方和を比較しています。前者についてはモデル平方和、残差平方和の自由度が分かっていませんので、これを求めていきます。
まずモデル平方和の自由度ですが、先ほどの議論からモデル平方和を考えるには、対応するモデルの計画行列$X$のランクを考えれば良いことが分かるかと思います。$1<=i<=2, 1<=j<=2, 1<=k<=2$の場合、今回のモデルは次のように定式化されます。
\begin{pmatrix} y_{111} \\ y_{112} \\ y_{121} \\ y_{122} \\ y_{211} \\ y_{212} \\ y_{221} \\ y_{222} \end{pmatrix} = \begin{pmatrix} 1 & 1 & 0 & 1 & 0 & 0 & 0 \\ 1 & 1 & 0 & 1 & 0 & 0 & 0 \\ 1 & 1 & 0 & 0 & 1 & 0 & 0 \\ 1 & 1 & 0 & 0 & 1 & 0 & 0 \\ 1 & 0 & 1 & 0 & 0 & 1 & 0 \\ 1 & 0 & 1 & 0 & 0 & 1 & 0 \\ 1 & 0 & 1 & 0 & 0 & 0 & 1 \\ 1 & 0 & 1 & 0 & 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} \mu \\ a_1 \\ a_{2} \\ b_{11} \\ b_{12} \\ b_{21} \\ b_{22} \end{pmatrix} + \begin{pmatrix} e_{111} \\ e_{112} \\ e_{121} \\ e_{122} \\ e_{211} \\ e_{212} \\ e_{221} \\ e_{222} \end{pmatrix}
ここで、$a$に対応する列(2,3列目)は足すと$X$の一列目と一致することから、まずランクが1つ落ちます。次に$b$に関する列ですが、$b_{ij}$の$i$の添え字が一致するものに対応する行を足すと$a_i$の列と一致することがわかります(2列目=4列目＋5列目,3列目=6列目＋7列目)。従って、$a$の水準分ランクが落ちてしまいます。行列$X$の列数は$1+l+lm$ですので、ランクは$1+l+lm-(1+l)=lm$になります。よって、モデル平方和$R(\hat{\beta})$の自由度は$lm$とわかります。次に残差平方和$R(\hat{\beta})$の自由度ですが、$I$に対応するランクが$lmn$ですので、$lmn-lm=lm(n-1)$となります。帰無仮説$y_{ijk}=\mu+a_i+e_{ijk}$の残差平方和の自由度は$lmn-l=l(mn-1)$ですので、差分を取った$lm(n-1)-l(mn-1)=l(m-1)$が処理$b$に対する平方和になります。
誤差の平方和
最後に誤差の平方誤差ですが、これは上を引き継ぎ、$lmn-lm=lm(n-1)$となります。
ということで、以下の表が以上のまとめになります。
| 処理 | 自由度 | 平方和 | 線形モデル表現 |
|---|---|---|---|
| 処理A | l-1 | $\sum \sum (y_{i.}-y_{..})^2$ | $R(\beta_{\mu})-R(\beta_a)$ |
| 処理B | l(m-1) | $\sum \sum (y_{i.}-y_{..})^2$ | $R(\beta_{a})-R(\beta_{ab})$ |
| 誤差 | lm(n-1) | $\sum \sum (y_{ij}-y_{i.})^2$ | $R(\beta_{ab})$ |
| 合計 | lmn-1 | $\sum \sum (y_{ij}-y_{..})^2$ | $R(\beta_{\mu})$ |
交互作用を含む二元配置分散分析
次に交互作用を含む二元配置分散分析のモデル$y_{ijk}=\mu+a_{i}+b_{j}+c_{ij}+e_{ijk}$を考えます。$1<=i<=l, 1<=j<=m, 1<=k<n$の場合、この時のanova表は次のようになることが知られています。
| 処理 | 自由度 | 平方和 |
|---|---|---|
| 処理A | l-1 | $\sum \sum (y_{i..}-y_{...})^2$ |
| 処理B | m-1 | $\sum \sum (y_{.j.}-y_{...})^2$ |
| 交互作用 | (l-1)(m-1) | $\sum \sum (y_{ij.}-y_{i..}-y_{.j.}+y_{...})^2$ |
| 誤差 | lm(n-1) | $\sum \sum (y_{ijk}-y_{ij.})^2$ |
| 合計 | lmn-1 | $\sum \sum (y_{ijk}-y_{...})^2$ |
処理A/Bの平方和
処理$a$の平方和については、$y_{ijk}=\mu+a_i+e_{ijk}$と$y_{ijk}=\mu+e_{ijk}$、処理$b$の平方和については$y_{ijk}=\mu+b_i+e_{ijk}$と$y_{ijk}=\mu+e_{ijk}$の残差平方和の差分となりますので、それぞれの処理における平方和の自由度は$l-1$、$m-1$となることが分かります(一元配置分散分析の方法を参考にしてください)。
交互作用の平方和
次に交互作用の平方和ですが、$y_{ijk}=\mu+a_{i}+b_{j}+c_{ij}+e_{ijk}$の残差平方和と$y_{ijk}=\mu+a_{i}+b_{j}+e_{ijk}$の残差平方和の差分となります。この差分に対応する自由度は巣ごもりモデルと同様、対応する$X$のランクを通して考えます。$1<=i<=2, 1<=j<=2, 1<=k<=2$の場合、今回のモデルは次のように定式化されます。
\begin{pmatrix} y_{111} \\ y_{112} \\ y_{121} \\ y_{122} \\ y_{211} \\ y_{212} \\ y_{221} \\ y_{222} \end{pmatrix} = \begin{pmatrix} 1 & 1 & 0 & 1 & 0 & 1 & 0 & 0 & 0 \\ 1 & 1 & 0 & 1 & 0 & 1 & 0 & 0 & 0 \\ 1 & 1 & 0 & 0 & 1 & 0 & 1 & 0 & 0 \\ 1 & 1 & 0 & 0 & 1 & 0 & 1 & 0 & 0 \\ 1 & 0 & 1 & 1 & 0 & 0 & 0 & 1 & 0 \\ 1 & 0 & 1 & 1 & 0 & 0 & 0 & 1 & 0 \\ 1 & 0 & 1 & 0 & 1 & 0 & 0 & 0 & 1 \\ 1 & 0 & 1 & 0 & 1 & 0 & 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} \mu \\ a_1 \\ a_{2} \\ b_{1} \\ b_{2} \\ c_{11} \\ c_{12} \\ c_{21} \\ c_{22} \end{pmatrix} + \begin{pmatrix} e_{111} \\ e_{112} \\ e_{121} \\ e_{122} \\ e_{211} \\ e_{212} \\ e_{221} \\ e_{222} \end{pmatrix}
ここで、処理$a$に対応する列(2,3列目)の和は$X$の一列目と一致するためランクは一つ下がります。また処理$c$に対応する列(6~9列目)について、$a_iの列=\sum_j c_ijの列,b_jの列=\sum_i c_ijの列$が成り立つため、$l+m$個のランクが落ちます。$X$の列数は$1+l+m+n$なので、ランクは$1+l+m+lm-(1+l+n)=lm$となります。従って$y_{ijk}=\mu+a_{i}+b_{j}+c_{ij}+e_{ijk}$のモデル平方和の自由度は$lm$、残差平方和の自由度は$lmn-lm=lm(n-1)$となります。次にモデル$y_{ijk}=\mu+a_{i}+b_{j}+e_{ijk}$に対応する線形モデルは次のように書き表されます。
\begin{pmatrix} y_{111} \\ y_{112} \\ y_{121} \\ y_{122} \\ y_{211} \\ y_{212} \\ y_{221} \\ y_{222} \end{pmatrix} = \begin{pmatrix} 1 & 1 & 0 & 1 & 0 \\ 1 & 1 & 0 & 1 & 0 \\ 1 & 1 & 0 & 0 & 1 \\ 1 & 1 & 0 & 0 & 1 \\ 1 & 0 & 1 & 1 & 0 \\ 1 & 0 & 1 & 1 & 0 \\ 1 & 0 & 1 & 0 & 1 \\ 1 & 0 & 1 & 0 & 1 \end{pmatrix} \begin{pmatrix} \mu \\ a_1 \\ a_{2} \\ b_{1} \\ b_{2} \end{pmatrix} + \begin{pmatrix} e_{111} \\ e_{112} \\ e_{121} \\ e_{122} \\ e_{211} \\ e_{212} \\ e_{221} \\ e_{222} \end{pmatrix}
ここで、$a$に対応する列(2,3列)の和は1列目に一致しますのでランクは一つ落ち、さらに$b$に対応する列(4,5列)の和は$a$に対応する列の和と一致しますのでランクがさらに一つ落ちます。$X$の列数は$l+m+1$ですから、行列$X$のランクは$1+l+m-(1+1)=l+m-1$です。従って、このモデルのモデル平方和の自由度は$l+m-1$、残差平方和の自由度は$lmn-l-m-1$となります。以上より、交互作用に対する平方和の自由度は$lmn-l-m-1-(lmn-lm)=lm-l-m-1=(l-1)(m-1)$となり、ANOVA表のものと一致します。
ということで以上をまとめたものが次の表になります。
Type III
| 処理 | 自由度 | 平方和 | 線形モデル表現 |
|---|---|---|---|
| 処理A | l-1 | $\sum \sum (y_{i..}-y_{...})^2$ | $R(\beta_{B,A:B})-R(\beta_{\beta_{A,B,A:B}})$ |
| 処理B | m-1 | $\sum \sum (y_{.j.}-y_{...})^2$ | $R(\beta_{A,A:B})-R(\beta_{\beta_{A,B,A:B}})$ |
| 交互作用 | (l-1)(m-1) | $\sum \sum (y_{ij.}-y_{i..}-y_{.j.}+y_{...})^2$ | $R(\beta_{A,B})-R(\beta_{\beta_{A,B,A:B}})$ |
| 誤差 | lm(n-1) | $\sum \sum (y_{ijk}-y_{ij.})^2$ | $R(\beta_{\beta_{A,B,A:B}})$ |
| 合計 | lmn-1 | $\sum \sum (y_{ijk}-y_{...})^2$ | --- |
処理A、Bの線形モデル表現ところをよく見るとさっきと言っていることが違うじゃんとなる方がいらっしゃるかと思いますが、今回は説明の都合上このようにさせてもらいました。また今回は合計に対応する線形モデル表現の欄が"---"となっています。ここは本来であれば$R(\mu)$となることが期待されるのですが、線形モデル表現の欄を足しあげてわかるように、実は交互作用を含むモデルでは合計が$R(\mu)$にはなりません。今回ご説明した方法で平方和を構成していく方法はTypeIIIと呼ばれる方法で、処理A,Bを対等に扱えるというメリットがある一方、平方和の合計が合計に一致しないという問題があります。一方TypeIと呼ばれる方法では平方和の合計こそ合計に一致するものの、処理A
,Bを対等に扱えなくなるという問題を抱えています。この問題は反復数が処理の組み合わせごとに違うアンバランスなデータの場合特に問題になるのですが、ここでは本旨から少しずれてしまうので、TypeIのANOVA表を乗せるだけにしておきます。
Type I
| 処理 | 自由度 | 平方和 | 線形モデル表現 |
|---|---|---|---|
| 処理A | l-1 | $\sum \sum (y_{i..}-y_{...})^2$ | $R(\beta_{\mu})-R(\beta_{A})$ |
| 処理B | m-1 | $\sum \sum (y_{.j.}-y_{...})^2$ | $R(\beta_{A})-R(\beta_{A,B})$ |
| 交互作用 | (l-1)(m-1) | $\sum \sum (y_{ij.}-y_{i..}-y_{.j.}+y_{...})^2$ | $R(\beta_{A,B})-R(\beta_{A,B,A:B})$ |
| 誤差 | lm(n-1) | $\sum \sum (y_{ijk}-y_{ij.})^2$ | $R(\beta_{A,B,A:B})$ |
| 合計 | lmn-1 | $\sum \sum (y_{ijk}-y_{...})^2$ | $R(\beta_{\mu})$ |
まとめ
最後にかけてかなり雑になってしまいましたが、今回は線形モデルという観点から統一的に分散分析を説明してみました。実際は今回紹介したデザインのほかにも分割区法といったより複雑なモデルが存在し、その時に線形モデルではどう平方和がかかれるかも書くべきだったのでしょうが、紙面の都合上(qiitaは文章が長くなると重くなる…)今回は割愛させていただきました。いずれこういったものについてもまとめていければと思います。いずれにせよ、今回の説明が分散分析にモヤモヤを抱いている方の一助になれば幸いです。たくさんな違いがあると思いますので見つけ次第ご連絡いただければ幸いです
Appendix： 射影行列Hについて
一元配置分散分析の線形モデル$\mathbf{y}=X\beta+\mathbf{e}$において、射影行列$H=X(X^TX)^{-}X^T$が次の形になることを軽く示します。
H_{ij} = \left\{
\begin{array}{ll}
\frac{1}{n} & 観測値iとjが同じ処理 \\
0 & その他
\end{array}
\right.
まず計画行列$X$は次のように定義されていますが、これは2行目以降を足すと1行目になることから冗長性をはらんでいることが分かります。
X=\begin{pmatrix} 1 & 1 & & \\ \vdots & \vdots & O & \\ 1 & 1 & & \\ \vdots & & \ddots & \\ 1 & & & 1\\ \vdots & & O & \vdots \\ 1 & & & 1 \end{pmatrix}
射影行列$H=X(X^TX)^{-}X^T$は右からベクトルを作用すると$X$の列ベクトルが張る空間に射影する行列となっていますので、$X$の列ベクトルが張る空間を変えない限りにおいて$X$を自由に変えてよいという性質があります。そこで$X$の代わりとなる$X'$を次のように定義します。なお、射影行列の性質より$H=X(X^TX)^{-}X^T=X'(X'^TX')^{-}X'^T$が成り立ちます。
X'=\begin{pmatrix} 1 & & \\ \vdots & O & \\ 1 & & \\ & \ddots & \\ & & 1\\ & O & \vdots \\ & & 1 \end{pmatrix}
この時、$(X'^TX')^{-}=\frac{1}{n}I$となりますので、$H=\frac{1}{n}X'X^T$となり、以下が成立することがわかります。
H_{ij} = \left\{
\begin{array}{ll}
\frac{1}{n} & 観測値iとjが同じ処理 \\
0 & その他
\end{array}
\right.
追加予定の項目
・平方和の期待値を線形モデルの観点からきれいに導出
・分散分析を用いて分散成分を推定する (Hendersonの方法)
・Type I, II, IIIについて線形モデルの観点から説明
・カイ二乗検定
・冪等性