はじめに

今回は、若干マニアックな内容ですが、分散分析における推定可能関数(Estimable Function)について取り上げます。この推定可能関数は、分散分析の教科書で出てくる内容ですが、正直言って、分散分析をするうえでは知らなくても特に問題はありません。なぜなら分散分析では、因子の水準間に違いがあるかどうかにのみ興味があり、各水準の具体的な値(パラメータの値)自体にあまり興味がないからです。問題になるのは、混合モデルのように、分散分析と同じようなモデルを使い、パラメータの値を推定するような場面です(単回帰なども似たモデルを使いますが、興味深いことに単回帰では推定可能関数は問題にはなりません)。しかし、混合モデルは、推定可能関数が問題にならない単回帰の延長として習うことが多いため、重要な割に知名度が低いというのが現状でした。この記事を通して推定可能関数の重要さが伝われば幸いです。

分散分析におけるパラメータ推定

具体的に話を進めていくために、以下では次のような分散分析のモデルを考えます。
$$
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$に対して線形であり、単回帰モデルと同じような見た目になります。これが分散分析は本質的に単回帰と同じだと言われる所以です。一方、単回帰モデルと決定的に異なるのが行列$X$の構造です。

例えば、単回帰モデル$y_{i}=\mu+\alpha x_{i}+e_{i}$の場合、先ほどの行列$X$に対応する行列は次のような形を取ります。

X=\begin{pmatrix}1 & x_1 \\ \vdots & \vdots \\ 1 &x_n \end{pmatrix}

ここで$x_1,...x_n$には実数が入りますので、基本的に行列$X$のランクは列数に一致します。一方、分散分析の場合、データがどこかの水準に属する関係で、ランクはほぼ確実に列数より小さくなってしまいます。これが単回帰の行列$X$と分散分析の行列$X$の決定的な違いです。

この結果、分散分析ではパラメータが推定できなくなるという問題が生じます。なぜなら、パラメータ推定の過程で出現する$(X^TX)^{-1}$は、ランクが$X$の列数に一致しない場合には、求めることができないためです。つまり何も考えない場合、パラメータは"推定不可能"であるわけです。

パラメータの線形結合から成る推定可能な値

このように、分散分析ではパラメータ推定ができないわけですが、興味深いことに$\beta$の要素同士を線形結合した結果生じる値の中には推定可能なものが含まれます。例えば次のように定義される値は推定可能です。

\mu+\alpha_1 =\begin{pmatrix}1 & 1 & \cdots &0\end{pmatrix} \begin{pmatrix} \mu \\ \alpha_1 \\ \vdots \\ \alpha_a \end{pmatrix}

なぜ推定可能なのか、についての証明は省きますが、このようにパラメータ$\beta$の線形結合から成り、推定可能な値のことを推定可能関数と呼んでいます。繰り返しになりますが、これは行列$X$のランクが基本的に列数と一致する単回帰においては考える必要がない現象です。

推定可能関数の条件

ここで気になるのが、どのような値が推定可能関数となるか、です。
先ほどの$\mu+\alpha_1$は$\beta$に$(1 1 ...0)$を作用させることで得られましたが、この作用させたベクトルが、推定可能関数になるかどうかを決めています。
先に結論だけ言ってしまうと、作用させるベクトルが$X^T$の各列のベクトルが張る空間内にある場合(作用させるベクトルが$X^T$の各列のベクトルの線形結合で表せる場合)は推定可能関数となり、そうでない場合は、推定可能な関数にはなりません。例えば先ほど使った$(1 1 ...0)$というベクトルは次のように、$X^T$の各列のベクトルの線形結合で表すことが可能です。

\begin{pmatrix}1 \\ 1 \\ \vdots \\ 0\end{pmatrix} = X^T \begin{pmatrix}1 \\ 0 \\ \vdots \\ 0\end{pmatrix} = 1 \begin{pmatrix}1 \\ 1 \\ \vdots \\ 0\end{pmatrix} + 0 \begin{pmatrix}1 \\ 1 \\ \vdots \\ 0\end{pmatrix} +  ... + 0 \begin{pmatrix}1 \\ 0 \\ \vdots \\ 1\end{pmatrix}

同様に、$\mu+\alpha_2$や$\mu+\alpha_a$などが推定可能関数であることはすぐ確かめられます。逆に$\mu$や$\alpha_1$単体は、対応する$(1 0 ...0)$や$(0 1 ...0)$は、どう頑張っても$X^T$の各列のベクトルの線形結合で作り出すことができませんので、推定可能関数にはなりません。

今回は推定可能関数の必要性を書くのが目的ですので、なぜ推定可能関数をこのように判断できるかについては述べません。これについては「統計解析スタンダード実験計画法と分散分析」が詳しいですので、興味がある方はぜひこちらの9章あたりをご覧ください。

補足

先ほどの説明において、ランクが行列$X$の列数に一致しない場合、パラメータを推定できない、と書きましたが、パラメータが一意に定まらない、というのが正確な表現です。今回の場合は$\beta$の値が一意に定まらないわけですが、非常に興味深いことに、推定可能関数はどんな$\beta$の値を使っても値が一意に定まることが分かっています。
また、パラメータ推定を可能にする別の方法としてパラメータに適当な制限をかける方法があります。例えば、$\Sigma \alpha_i=0$とすることでパラメータを推定できるようになるのですが、これは推定可能関数の話の別表現でしかないことが分かっています。
これについても、「統計解析スタンダード実験計画法と分散分析」に詳しく書かれていますので、ご興味のある方はぜひ。