線形モデルの推定可能関数に対するBLUEとBLUP

# はじめに
本投稿は、線形混合モデルでよく出てくる、最良線形不偏推定量(BLUE: Best Linear Unbiased Estimation)と最良線形不偏予測量(BLUP: Best Linear Unbiased Prediction)について、本来の定義に立ち戻り、導出を試みたものです。特に、BLUEやBLUPにおける最良とは何か、についての誤解が多いように感じますので、本投稿を通して、BLUEやBLUPの正しい理解が広まればと思っています。なお、参考として、佐々木義之さんの「[変量効果の推定とBLUP法](https://www.kyoto-up.or.jp/books/9784876987023.html)」と[こちらの論文](https://ecommons.cornell.edu/server/api/core/bitstreams/295be5f7-3f0b-4d01-86cd-a351f4e4f679/content)を用いました。詳しくはそちらをご覧ください。

### 目次
- [線形モデルと推定可能関数](#線形モデルと推定可能関数)
- [推定値が最良であるとは](#推定値が最良であるとは)
- [推定可能関数に対するBLUE](#推定可能関数に対するblue)
- [推定可能関数に対するBLUP](#推定可能関数におけるblup)
- [ガウス・マルコフの定理](#ガウスマルコフの定理)
- [まとめ](#まとめ)

# 線形モデルと推定可能関数
さて、本投稿で説明をするBLUEとBLUPは、線形モデルにおける推定可能関数に対して定義される推定値です。この推定可能関数については[こちらの記事](https://qiita.com/KimuraSoh/items/c8a49a1ac4fd2c59bb16)に詳しくまとめてありますが、ここで今一度説明をしておきます。

まず、次のような線形モデルを考えます。
```math
\begin{align}
\mathbf{y}=X\mathbf{b}+\mathbf{e}\quad 
E(\mathbf{y})=X\mathbf{b},& \quad Var(\mathbf{y})=V
\end{align}
```
ここで、$\mathbf{b}$と同じ長さを持つベクトル$\mathbf{k}$と、$\mathbf{b}$との内積を取った値である$\mathbf{k}^T\mathbf{b}$を、$\mathbf{y}$から推定することを考えます。様々な推定値の形が考えられますが最も単純な$\mathbf{a}^T\mathbf{y}$から推定してみます($\mathbf{a}$はベクトル)。

さすがにこれだけの情報量からだと推定ができないので、この$\mathbf{a}^T\mathbf{y}$に対して不偏性という性質を加えてみます。この不偏性とは、推定値に対して定義される単語で、推定値の期待値が真値に一致する性質を指します。この$\mathbf{a}^T\mathbf{y}$に不偏性を課すと次のような式が得られます。
```math
\mathbf{k}^T\mathbf{b}=E(\mathbf{a}^t\mathbf{y})
```
ここで、$\mathbf{y}$に対しては$E(\mathbf{y})=X\mathbf{b}$という条件が成り立ちますので上の式はさらに変形できて、次の通りになります。
```math
\begin{align}
\mathbf{k}^T\mathbf{b}=E(\mathbf{a}^t\mathbf{y})&=\mathbf{a}^TE(\mathbf{y})=\mathbf{a}^TX\mathbf{b}\\
\therefore \mathbf{k}^T&=\mathbf{a}^TX
\end{align}
```
この式に対して、転置を取ると$\mathbf{k}=X^T\mathbf{a}$となりますが、これはまず、$\mathbf{k}$が$X$の列ベクトルの線形結合で表されること、つまり、$X$の行ベクトルが張る空間$\mathscr{M}(X^T)$の中に存在していることを示しています。これは不偏性を要求したことにより生じた制限であり、逆に、$X$の行ベクトルが張る空間内に存在している$\mathbf{k}$を用いた$\mathbf{k}^T\mathbf{b}$には、不偏性を持つ推定値$\mathbf{a}^T\mathbf{y}$が存在することが分かります。このような理由から、$\mathbf{k}\in\mathscr{M}(X^T)$を満たす$\mathbf{k}$と内積を取ることで作られる$\mathbf{k}^T\mathbf{b}$は、**推定可能関数**と呼ばれます。

例えば、$\mathbf{b}$のi番目の要素$b_i$は、$\mathbf{k}^T\mathbf{b}$において、$\mathbf{k}=(0,0,...1,...0)$である場合に相当しますが、多くの場合、そのようなベクトルは$\mathscr{M}(X^T)$に含まれません。したがって、$\mathbf{b}$の各要素は推定可能関数ではなく、したがって不偏推定量が存在しません。逆に、$X\mathbf{b}$の各要素は$X$の行ベクトルとの内積になっていますので、すべて推定可能関数です。分散分析で$\mathbf{b}$よりも$X\mathbf{b}$が多く出てくるのはそのような理由に由来しています。

さて、これまでの議論より、$\mathbf{k}\in\mathscr{M}(X^T)$を満たす$\mathbf{k}$と内積を取ることで作られる$\mathbf{k}^T\mathbf{b}$に対して、$\mathbf{a}^T\mathbf{y}$の形を取る不偏推定量が存在することはわかりましたが、そのような不偏推定量は複数あり、これだけの条件では一意に定まりません。そこで、最良性という新しい性質を$\mathbf{a}^T\mathbf{y}$に新しく課すことにより、解を絞っていきます。実は、この最良性には複数の定義があり、ある定義を課した上で定まる解がBLUE、別の定義を課した上で定まる解がBLUPに対応します。以下ではそれをじっくりと見ていきます。

:::note info
**ちょっとしたメモ**
線形モデルに課される$[\mathbf{y}]=X\mathbf{b}$という条件について、私は今まで$\mathbf{y}$の期待値は$X\mathbf{b}$となるという読み方をしていましたが、この推定可能関数を正しく理解するには、$\mathbf{y}$は$X\mathbf{b}$の不偏推定量の一つである、という読み方をした方が良いように思います。つまり、線形モデルで仮定される条件は以下の二つに書き換えられます。
①$\mathbf{y}$は$X\mathbf{b}$の不偏推定量の一つである
②$\mathbf{y}$の分散は$V$である
この条件のみから、$\mathbf{k}^T\mathbf{b}$について言えることを考えようとしているのが、上記の話となります。
:::

# 推定値が最良であるとは
ここで少し脱線し、推定値の最良性について考えてみます。以下のチャプターで書くことをより理解しやすくするための解説ですので、BLUEとBLUPの違いが早く知りたい方は飛ばしていただいてもかまいません。

平均$\mu$、標準誤差$\sigma$従う確率変数Xが10個あり(それぞれ$X_1,X_2,...X_{10}$とする)、ここから平均値を推定するとします。最も一般的な方法は、$\frac{1}{10}\sum X_i$を計算する方法だと思いますが、その妥当性は次のように期待値が真値に一致することで保証されます。
$$
E\Big(\frac{1}{10}\sum X_i\Big)=\frac{1}{10}E\Big(\sum X_i\Big)=\frac{1}{10}\sum E\Big( X_i\Big)=\mu
$$
つまり、この平均による推定値は不偏推定値なのです。しかし、わざわざサンプル10個すべてを使わなくとも、$X_1$だけでもその期待値は真値に一致します。
```math
E(X_1)=\mu
```
また、$max(X_1,X_2,...X_{10})$や$(X_1-X_2+X_3)$なども、期待値が真値することから不偏推定値であることが確かめられます。このように不偏性を持つ推定値はたくさん存在し、不偏性という観点だけではそれぞれの推定値の善し悪しは判断できません。

そこで考えるのが、分散です。10個の確率変数$X_1,X_2,...X_{10}$は、当然、別日に観測すれば別の値が得られます。したがって、推定値は毎回、変動します。この変動が大きい場合は、日ごとに大きく推定値が変わるわけですから、あまり信用なりません。逆に変動が少ない場合は、たとえ1日分の結果でも信用することができるでしょう。このように推定値の変動(=分散)に着目することで、不偏推定量に対して、信頼性のようなものを表すことができ、不偏推定量の善し悪しを判断することが可能になります。そして、特にそのような分散が最小となる不偏推定値は**最良性を持つ**と呼ばれます。これが、私たちが最良性を考えたい直観的な理由です。
なお、不偏推定値の分散には下限があり、クラメール-ラオの不等式より、その下限がフィッシャー情報量の逆数に一致することが分かっています。(フィッシャー情報量については[こちら](https://note.com/kimukimucitrus/n/nace8ba84d129)にその直観的なイメージをまとめています、ご興味のある方はぜひご覧ください)

# 推定可能関数に対するBLUE
それでは、推定可能関数$\mathbf{k}^T\mathbf{b}$の推定量$\mathbf{a}^T\mathbf{y}$に対して、最良性を課して、解を絞っていきます。

一旦、今まで課してきた制限をまとめると、以下の二つが挙げられます。
①線形性($\mathbf{a}^T\mathbf{y}$として$\mathbf{y}$に対して線形な形でかける)
②不偏性(期待値が一致: $\mathbf{k}^T\mathbf{b}=E(\mathbf{a}^T\mathbf{y})$)

ここに次のような最良性に関する条件を付け加えます
③最良性($Var(\mathbf{a}^T\mathbf{y}))$が最小になる)

$\mathbf{a}$はまだ具体的な値が分かっていませんので、これらの条件のもと、計算をすることで$\mathbf{a}$の具体的に求めていきます。しかし、②の不偏性の条件から$\mathbf{k}=X^T\mathbf{a}$という制約が生じますので、正確にはラグランジュの未定乗数を用いた次の値の最小化問題となります。
```math
Var(\mathbf{a}^T\mathbf{y}))+2\mathbf{t}^T(\mathbf{k}-X^T\mathbf{a})
```
ここで、この式はパラメータ$\mathbf{a},\mathbf{t}$に支配されていますので、$f(\mathbf{a},\mathbf{t})$のように書くことにします。
```math
\begin{align}
f(\mathbf{a},\mathbf{t})&=Var(\mathbf{a}^T\mathbf{y}))+2\mathbf{t}^T(\mathbf{k}-X^T\mathbf{a})\\&=\mathbf{a}^TV\mathbf{a}+2\mathbf{t}^T(\mathbf{k}-X^T\mathbf{a})
\end{align}
```

この$f(\mathbf{a},\mathbf{t})$を最小化する$\mathbf{a},\mathbf{t}$を求めるには、微分をして、$=0$と置けば良いので、まず$f(\mathbf{a},\mathbf{t})$を$\mathbf{a},\mathbf{t}$それぞれで微分してみます。

```math
\begin{align}
\frac{\partial}{\partial\mathbf{a}}f(\mathbf{a},\mathbf{t})
&=\frac{\partial}{\partial\mathbf{a}}\Big(\mathbf{a}^TV\mathbf{a}+2\mathbf{t}^T(\mathbf{k}-X^T\mathbf{a})\Big)\\
&=\frac{\partial}{\partial\mathbf{a}}\Big(\mathbf{a}^TV\mathbf{a}+2\mathbf{k}^T\mathbf{t}-2\mathbf{a}^TX\mathbf{t}\Big)\\
&=(V^T+V)\mathbf{a}-2X\mathbf{t}\\
&=2V\mathbf{a}-2X\mathbf{t}
\end{align}
```
```math
\begin{align}
\frac{\partial}{\partial\mathbf{t}}f(\mathbf{a},\mathbf{t})
&=\frac{\partial}{\partial\mathbf{t}}\Big(\mathbf{a}^TV\mathbf{a}+2\mathbf{t}^T(\mathbf{k}-X^T\mathbf{a})\Big)\\
&=2(\mathbf{k}-X^T\mathbf{a})
\end{align}\\
```
ここで、$\mathbf{a}$に対する微分について$=0$と置き、$\mathbf{a}$について解く、
```math
\mathbf{a}=V^{-1}X\mathbf{t}
```
さらに左から$X^T$を作用させることにより、
```math
X^T\mathbf{a}=X^TV^{-1}X\mathbf{t}
```
次に、$\mathbf{t}$に対する微分について$=0$と置くと次を得ます。
```math
\mathbf{k}=X^T\mathbf{a}
```
これをひとつ上の式に代入すると、
```math
\mathbf{k}=X^TV^{-1}X\mathbf{t}
```
これを$\mathbf{t}$について解いて、
```math
\mathbf{t}=(X^TV^{-1}X)^{-1}\mathbf{k}
```
最後に、$\mathbf{a}=V^{-1}X\mathbf{t}$にこの$\mathbf{t}$を代入すると
```math
\mathbf{a}=V^{-1}X(X^TV^{-1}X)^{-1}\mathbf{k}
```
したがって、分散を最小化する推定量$\mathbf{a}^t\mathbf{y}$を次のような形を取ることが分かります。
```math
\mathbf{a}^t\mathbf{y}=\mathbf{k}^TX(X^TV^{-1}X)^{-1}X^TV^{-1}\mathbf{y}
```
このような推定量は、最良性を持つ、線形で不偏性を備えた推定量ということで、最良線形不偏推定量(BLUE: Best Linear Unbiased Estimation)と呼ばれています。

# 推定可能関数におけるBLUP
さて、今まで見てきたように、$\mathbf{a}^T\mathbf{y}$に対して分散が最小という条件を課すことにより推定可能関数$\mathbf{k}^\mathbf{b}$のBLUEが得られたわけですが、次に、少し別の最良性を課し、違う推定量を求めてみます。課す条件は次の通りです
③最良性(真値と推定量の誤差分散($Var(\mathbf{k}^T\mathbf{b}-\mathbf{a}^T\mathbf{y}))$が最小になる)

先ほどは推定値の分散$Var(\mathbf{a}^T\mathbf{y}))$が最小になることを最良として扱っていたので、先ほどとは少し意味合いが異なります。先に答えを述べておくと、この条件のもと導かれる推定量がBLUPの正体です。なお、$\mathbf{k}^T\mathbf{b}$と$\mathbf{y}$との共分散が無い場合はBLUEと同じ解が導かれます。

さて、ではこの最良性のもと、線形性と不偏性を課しながら、推定可能関数$\mathbf{k}^T\mathbf{b}$の推定値を求めてみます。最小化すべき関数$f(\mathbf{a},\mathbf{t})$は先ほどと同様の議論から次のように導かれます。

```math
\begin{align}
f(\mathbf{a},\mathbf{t})&=Var(\mathbf{k}^T\mathbf{b}-\mathbf{a}^T\mathbf{y}))+2\mathbf{t}^T(\mathbf{k}-X^T\mathbf{a})\\
&=Var(\mathbf{k}^T\mathbf{b})+Var(\mathbf{a}^T\mathbf{y})-2Cov(\mathbf{k}^T\mathbf{b},\mathbf{a}^T\mathbf{y})+2\mathbf{t}^T(\mathbf{k}-X^T\mathbf{a})\\
&=Var(\mathbf{k}^T\mathbf{b})+Var(\mathbf{a}^T\mathbf{y})-2Cov(\mathbf{k}^T\mathbf{b},\mathbf{a}^T\mathbf{y})+2\mathbf{t}^T(\mathbf{k}-X^T\mathbf{a})\\
&=Var(\mathbf{k}^T\mathbf{b})+\mathbf{a}^TV\mathbf{a}-2\mathbf{a}^TCov(\mathbf{k}^T\mathbf{b},\mathbf{y})+2\mathbf{t}^T(\mathbf{k}-X^T\mathbf{a})
\end{align}
```
ここで、単純のため、$Cov(\mathbf{k}^T\mathbf{b},\mathbf{y})=C$と置くことにします。したがって、
```math
f(\mathbf{a},\mathbf{t})=Var(\mathbf{k}^T\mathbf{b})+\mathbf{a}^TV\mathbf{a}-2\mathbf{a}^TC+2\mathbf{t}^T(\mathbf{k}-X^T\mathbf{a})
```
ここで、$\mathbf{a},\mathbf{t}$それぞれについて微分をすると
```math
\begin{align}
\frac{\partial}{\partial\mathbf{a}}f(\mathbf{a},\mathbf{t})
&=2V\mathbf{a}-2C+2X\mathbf{t}
\end{align}
```
```math
\begin{align}
\frac{\partial}{\partial\mathbf{t}}f(\mathbf{a},\mathbf{t})
&=2(\mathbf{k}-X^T\mathbf{a})
\end{align}
```
ここで先ほどと同様、$\mathbf{a}$に対する微分について$=0$と置き、$\mathbf{a}$について解くと、
```math
\mathbf{a}=V^{-1}C-V^{-1}X\mathbf{t}
```
さらに左から$X^T$を作用して
```math
X^T\mathbf{a}=X^TV^{-1}C-X^TV^{-1}X\mathbf{t}
```
ここに、$\mathbf{a}$に対する微分について$=0$と置いた結果得られる、$\mathbf{k}=X^T\mathbf{a}$を代入して
```math
\mathbf{k}=X^TV^{-1}C-X^TV^{-1}X\mathbf{t}
```
これを$\mathbf{t}$について解いて、
```math
\mathbf{t}=(X^TV^{-1}X)^{-1}(X^TV^{-1}C-\mathbf{k})
```
最後に、$\mathbf{a}=V^{-1}C-V^{-1}X\mathbf{t}$に対して、この$\mathbf{t}$を代入して、
```math
\begin{align}
\mathbf{a}=V^{-1}C-V^{-1}X(X^TV^{-1}X)^{-1}(X^TV^{-1}C-\mathbf{k})
\end{align}
```
したがって、$\mathbf{k}^T\mathbf{b}$の推定量$\mathbf{a}^T\mathbf{y}$として次が得られます。
```math
\begin{align}
\mathbf{a}^T\mathbf{y}&=C^TV^{-1}\mathbf{y}-(X^TV^{-1}C-\mathbf{k})^T(X^TV^{-1}X)^{-1}X^TV^{-1}\mathbf{y}\\
&=C^TV^{-1}\mathbf{y}-C^TV^{-1}X(X^TV^{-1}X)^{-1}X^TV^{-1}\mathbf{y}+\mathbf{k}^T(X^TV^{-1}X)^{-1}X^TV^{-1}\mathbf{y}\\
&=\mathbf{k}^T(X^TV^{-1}X)^{-1}X^TV^{-1}\mathbf{y}+C^TV^{-1}(\mathbf{y}-X(X^TV^{-1}X)^{-1}X^TV^{-1}\mathbf{y})
\end{align}\\
```
この推定量は、最良性を持ち、線形で不偏性を備えた推定量ではありますが、BLUEと最良性が違うことを明示するために、最良線形不偏予測量(BLUP: Best Linear Unbiased Prediction)と呼ばれています。なお、よく観察をすると、$\mathbf{a}^T\mathbf{y}$の第1項が$\mathbf{k}^T\mathbf{b}$のBLUEになっていることが分かります。また今回は証明していませんが、第2項の括弧内の第2項$X(X^TV^{-1}X)^{-1}X^TV^{-1}\mathbf{y}$は$X\mathbf{b}$のBLUEに対応しています。第2項の$C$は$\mathbf{k}^T\mathbf{b}$と$\mathbf{y}$との共分散を表していますが、これが無い場合、第2項はすべてなくなり、第1項のみが残ります。このことからも、BLUPを特徴づけているのは、$\mathbf{k}^T\mathbf{b}$と$\mathbf{y}$との共分散ということができるでしょう。

# ガウス・マルコフの定理
少し余談ですが、最後にガウス・マルコフの定理について言及しておきます。

本投稿で紹介した導出方法は、どちらかというとマイナーな手法で、特に、単回帰モデル$\mathbf{y}=X\mathbf{b}+\mathbf{e}$における$\mathbf{b}$の推定に関しては、二乗和誤差である$\mathbf{e}^T\mathbf{e}=(\mathbf{y}-X\mathbf{b})^T(\mathbf{y}-X\mathbf{b})$を最小化するように$\mathbf{b}$を求めるのが一般的な方法だと思います(分散分析では$\mathbf{b}$は推定可能関数ではありませんが、単回帰モデルでは基本的に推定可能です)。

本投稿で紹介した方法では$\mathbf{\beta}$を真値と見なし、$\mathbf{a}$や$\mathbf{t}$に対して微分をすることで、その推定値$\mathbf{a}^T\mathbf{y}$を求めたのに対し、この二乗和誤差を最小化する方法では、$\mathbf{b}$そのものをパラメータと見なし、$\mathbf{e}^T\mathbf{e}$を$\mathbf{b}$で微分することを通して$\mathbf{b}$の推定値を求めます。その意味でこの方法と、今回ご紹介した方法は全くスタンスが違うと言えるでしょう。ところがこのように求まる$\mathbf{b}$の推定値は、$\mathbf{b}$のBLUEであることが証明されており、ガウス・マルコフの定理と呼ばれています。

BLUPが、誤差分散の最小化により求められたことから、二乗和誤差の最小化により求まるのが、BLUPではなくBLUEであることに少しの違和感を感じる方もいらっしゃるかと思います。これに関して、私ははっきりとした答えを持ち合わせているわけではないのですが、$\mathbf{b}$と$\mathbf{y}$との共分散が取り入れられない(し、うまく取り入れるように変形もできない)ためだろうと想像しています。

# まとめ
今回は、推定可能関数に対するBLUEとBLUPを、元もとの定義に立ち返り、導出を試みました。両者ともに最良性を持っているものの、前者が予測値の分散の最小化により求まるのに対し、後者は予測誤差分散の最小化により求まる、という違いがある、ということが皆様に伝われば、幸いです。かなり粗い計算をしていますので、間違いがあるかもしれません。何かございましたらご連絡いただければと思います。