変量効果の推定とBLUP法 ~最良予測量が条件付平均値であることの証明~

## はじめに
本投稿は、佐々木義之さんの「[変量効果の推定とBLUP法](https://www.kyoto-up.or.jp/books/9784876987023.html)」の第1章(p19-20)で解説されている、最良予測量(BP)が条件付平均値であることの証明([Henderson, 1973](https://academic.oup.com/jas/article-abstract/1973/Symposium/10/4697208?redirectedFrom=fulltext))について、その計算をより詳しく記述したものになります。より数学的には、バイアス・バリアンス分解と呼ばれる操作の過程として導かれるものですので、より詳しく学びたい方は「[パターン認識と機械学習 上](https://www.amazon.co.jp/%E3%83%91%E3%82%BF%E3%83%BC%E3%83%B3%E8%AA%8D%E8%AD%98%E3%81%A8%E6%A9%9F%E6%A2%B0%E5%AD%A6%E7%BF%92-%E4%B8%8A-C-M-%E3%83%93%E3%82%B7%E3%83%A7%E3%83%83%E3%83%97/dp/4621061224)」の第3章などをご覧ください。
なお本投稿は、線形混合効果モデルをすでに学習された方を対象としていますので、特に変量効果などの説明は行いません。ご了承ください。

### 目次
- [背景と証明内容](#背景と証明内容)
- [証明](#証明)
- [バイアス・バリアンス分解](#バイアスバリアンス分解)

## 背景と証明内容
まず、変量効果推定の背景と証明内容の確認します。
家畜や作物の品種改良において、変量効果は遺伝的効果に相当します。同効果は優秀な個体を選抜する指標として極めて重要ですが、直接観察することができないため、観察された値をもとに推定する必要があります。つまり、観察値を$\mathbf{y}$、変量効果を$\mathbf{g}$とすると、$\mathbf{y}$をもとに$\mathbf{g}$を推定する必要があるというわけです。
さて、その推定には様々な方法があり、推定方法ごとに異なるの推定値が得られます。本書における証明は、推定された変量効果$\mathbf{\hat{g}}$と真の変量効果$\mathbf{g}$との平均平方誤差$(\mathbf{g}-\mathbf{\hat{g}})^2$を最小にするような推定において、その推定値$\mathbf{\hat{g}}$が、観察値を与えたときの条件付平均、すなわち$E(\mathbf{g}|\mathbf{y})$になることを主張するものです。
なお、変量効果の推定値は、固定効果の推定値と異なり、予測量と呼ばれます。また、平均平方誤差を最小にする推定値は「最良である」という表現されます。したがって、このように推定された推定値は最良予測量(Best Predictor)と呼ばれます。これらの言葉を使って言い直すと、本書における証明は、**最良予測量が条件付平均値になることの証明**となります。

## 証明
さて、ここからは「最良予測量が条件付平均値になること」を証明していきます。簡単のためここではベクトル$\mathbf{g}$のうち、ある要素$g_i$に着目して推定することとします。
スタートポイントは本書と同じく、平均平方誤差の期待値$E(g_i-\hat{g}_i)^2$です。混乱を防ぐため、ここでは$E[(g_i-\hat{g}_i)^2]$と表記することにします。本書に明記はされていませんが、同期待値は正確には以下のように定義されています。
$$E[(g_i-\hat{g}_i)^2]=\int\int (g_i-\hat{g}_i)^2p(g_i, \mathbf{y})dg_id\mathbf{y}$$
さて、これを変形すると次が得られます。
$$E[(g_i-\hat{g}_i)^2]=E[(g_i-E(g_i|\mathbf{y})+E(g_i|\mathbf{y})-\hat{g}_i)^2]$$
ここで、$E(g_i|\mathbf{y})$は、本書において明記されていませんが以下のように定義される値です。
$$E(g_i|\mathbf{y})=\int g_i p(g_i|\mathbf{y})dg_i$$
さらに変形を重ねると、次の式が得られます。
$$E[(g_i-\hat{g}_i)^2]=E[(g_i-E(g_i|\mathbf{y}))^2]+E[(E(g_i|\mathbf{y})-\hat{g}_i)^2]+2E[(g_i-E(g_i|\mathbf{y}))(E(g_i|\mathbf{y})-\hat{g}_i)]$$
ここまでは本書に書いてある内容とほぼ同じで、多くの方が問題なく理解できているのではないでしょうか。さて、本書には続いて次のような一文が書かれています。

```math
\begin{align}
条件付き期待値の定義から第三項は\\
2E[(g_i-E(g_i|\mathbf{y}))(E(g_i|\mathbf{y})-\hat{g}_i)]=0
\end{align}
```
この変形については、本書にこれ以上の説明がないのですが、一文で済ませられるほど単純なものではありませんので、以下では第3項に注目し、同項が0になることを導出します。

まずこの期待値をしっかりと定義通り書き下すと以下の式が得られます。
$$
2E[(g_i-E(g_i|\mathbf{y}))(E(g_i|\mathbf{y})-\hat{g}_i)]=\int\int 2(g_i-E(g_i|\mathbf{y}))(E(g_i|\mathbf{y})-\hat{g}_i)p(g_i,\mathbf{y})dg_id\mathbf{y}
$$
ここで$g_i$に着目すると、被積分関数の中で$g_i$が出てくるのは$(g_i-E(g_i|\mathbf{y}))$の$g_i$だけですので、次のように変形することができます。
$$
2E[(g_i-E(g_i|\mathbf{y}))(E(g_i|\mathbf{y})-\hat{g}_i)]=2\int(E(g_i|\mathbf{y})-\hat{g}_i)\Biggl[\int (g_i-E(g_i|\mathbf{y}))p(g_i,\mathbf{y})dg_i\Biggr]d\mathbf{y}
$$
ここでカギカッコ$[]$の中身に着目します。条件付き確率の定義より次のような変形が可能です。

```math
\begin{align}
\int (g_i-E(g_i|\mathbf{y}))p(g_i,\mathbf{y})dg_i&=\int (g_i-E(g_i|\mathbf{y}))p(g_i|\mathbf{y})p(\mathbf{y})dg_i\\&=p(\mathbf{y})\int (g_i-E(g_i|\mathbf{y}))p(g_i|\mathbf{y})dg_i\\&=p(\mathbf{y})\int g_ip(g_i|\mathbf{y})dg_i-p(\mathbf{y})\int E(g_i|\mathbf{y}))p(g_i|\mathbf{y})dg_i\\&=p(\mathbf{y})\int g_ip(g_i|\mathbf{y})dg_i-p(\mathbf{y})E(g_i|\mathbf{y}))\int p(g_i|\mathbf{y})dg_i\\&=p(\mathbf{y})\int g_ip(g_i|\mathbf{y})dg_i-p(\mathbf{y})E(g_i|\mathbf{y})
\end{align}
```

ところが最終変形した式の第1項は条件付き期待値$E(g_i|\mathbf{y})$の定義そのものですので、第1項と第2項が打消しあって、カギカッコ$[]$の中身は0になります。被積分関数0の積分は0ですので、その結果、$2E[(g_i-E(g_i|\mathbf{y}))(E(g_i|\mathbf{y})-\hat{g}_i)]=0$が導かれます。

第3項が打ち消された結果、平均平方誤差の期待値は次の2項に分解されることが分かります。実はこの分解こそ、最初に少し伏線を張ったバイアスバリアンス分解に相当します。

$$E[(g_i-\hat{g}_i)^2]=E[(g_i-E(g_i|\mathbf{y}))^2]+E[(E(g_i|\mathbf{y})-\hat{g}_i)^2]$$

最良予測量という文脈では、左辺を最小化する推定値が"最良"であるわけですが、そのためには右辺第2項が$0$でなければなりません。ここで、$g_i$が以下を満たせば、それが成り立つことがわかります。
$$\hat{g}_i=E(g_i|\mathbf{y})$$
この事実をもって、「最良予測量は条件付き期待値である」という主張が認められたことになります。

:::note info
ちょっとしたメモ
$E[(g_i-\hat{g}_i)^2]$についてですが、$\hat{g}_i$は考えうるすべての推定値を表しています。したがって、$f(\hat{g}_i)=E[(g_i-\hat{g}_i)^2]$として捉えることができますし、そう考えた方が、今回の話は分かりやすくなるのと思います。
:::

## バイアス・バリアンス分解
最後に、本投稿の冒頭で少し触れたバイアス・バリアンス分解について説明します。
今までは、観測値$\mathbf{y}$から$g_i$を推定することを考えていましたが、ここからは任意の$\mathbf{y}$から$g_i$を予測するモデル$g_i(\mathbf{y})$を作ることを考えます。観測値の$\mathbf{y}$と任意の$\mathbf{y}$を区別するため、前者を$\mathcal{D}$、後者を$\mathbf{y}$とすることにします。
さて、もちろんこの予測関数に対しても平均平方誤差を定義することができますので、これを$E[(g_i(\mathbf{y})-\hat{g}_i(\mathbf{y}))^2]$とすることにします。計算自体は先ほどと同じなので省略しますが、もちろんこの平均平方誤差についても、先ほどと同様の分解ができます。
```math
E[(g_i(\mathbf{y})-\hat{g}_i(\mathbf{y}))^2]=E[(g_i(\mathbf{y})-E(g_i|\mathbf{y}))^2]+E[(E(g_i|\mathbf{y})-\hat{g}_i(\mathbf{y}))^2]
```
バイアス・バリアンス分解が扱うのは、このような予測関数に関する平均平方誤差を分解した第2項です。そこに出てくる$\hat{g}_i(\mathbf{y})$は推定されたモデルであり、通常はデータ$\mathcal{D}$をもとに推定されます。このようなデータ依存性を明示するため$\hat{g}_i(\mathbf{y}|\mathcal{D})$と書くことにします。
```math
E[(g_i(\mathbf{y})-\hat{g}_i(\mathbf{y}|\mathcal{D}))^2]=E[(g_i(\mathbf{y})-E(g_i|\mathbf{y}))^2]+E[(E(g_i|\mathbf{y})-\hat{g}_i(\mathbf{y}|\mathcal{D}))^2]
```
さきほどは、予測ではなく推定を行っていたため、この第二項を0にすることができたのですが、今回のようにデータにない点での予測を行うための関数を扱う場合、基本的に、第二項は0になりません。そこで、この第二項をさらに分解して、さらに性質を調べてみることにします。
さて、やや唐突ですが、$E_\mathcal{D}[\hat{g}_i(\mathbf{y}|\mathcal{D})]$という関数を定義します。これは次のように定義される量となっており、それぞれのデータから推定された予測関数の平均を表しています。
```math
E_\mathcal{D}[\hat{g}_i(\mathbf{y}|\mathcal{D})]=\int \hat{g}_i(\mathbf{y}|\mathcal{D})p(\mathcal{D})d\mathcal{D}
```
この関数$E_\mathcal{D}[\hat{g}_i(\mathbf{y}|\mathcal{D})]$を平均平方誤差の中に入れ、次のような分解をしてみます。
```math
\begin{align}
E[(E(g_i|\mathbf{y})-\hat{g}_i(\mathbf{y}|\mathcal{D}))^2]
&=E[(E(g_i|\mathbf{y})-E_\mathcal{D}[\hat{g}_i(\mathbf{y}|\mathcal{D})]+E_\mathcal{D}[\hat{g}_i(\mathbf{y}|\mathcal{D})]-\hat{g}_i(\mathbf{y}|\mathcal{D}))^2]\\
&=E[(E(g_i|\mathbf{y})-E_\mathcal{D}[\hat{g}_i(\mathbf{y}|\mathcal{D})])^2]+E[(E_\mathcal{D}[\hat{g}_i(\mathbf{y}|\mathcal{D})]-\hat{g}_i(\mathbf{y}|\mathcal{D}))^2]+
2E[(E(g_i|\mathbf{y})-E_\mathcal{D}[\hat{g}_i(\mathbf{y}|\mathcal{D})])(E_\mathcal{D}[\hat{g}_i(\mathbf{y}|\mathcal{D})]-\hat{g}_i(\mathbf{y}|\mathcal{D}))]
\end{align}
```
実は、この第3項目は先ほどと同様の計算をすることで$0$になることが分かりますので、最終的に以下の分解が得られます。
```math
\begin{align}
E[(E(g_i|\mathbf{y})-\hat{g}_i(\mathbf{y}|\mathcal{D}))^2]
&=E[(E(g_i|\mathbf{y})-E_\mathcal{D}[\hat{g}_i(\mathbf{y}|\mathcal{D})])^2]+E[(E_\mathcal{D}[\hat{g}_i(\mathbf{y}|\mathcal{D})]-\hat{g}_i(\mathbf{y}|\mathcal{D}))^2]
\end{align}
```
実はこの分解こそバイアス・バリアンス分解であり、第1項はバイアス、第2項はバリアンスと呼ばれます。ここで、$E(g_i|\mathbf{y})$は予測関数$\hat{g}_i(\mathbf{y}|\mathcal{D})$が近づける最も良い関数の上限であるため、平均平方誤差の第1項は理想と平均の差、第2項は平均と現実の差を表しているととらえることができます。パラメータ数を増やし予測モデルを複雑にすると、期待値の意味で理想とは合うようになりますが(バイアスの減少)、逆に、データに合わせすぎてしまうため(過学習)、データごとの予測モデルのばらつきが増加し、平均的な挙動と大きくずれるようになってしまいます(バリアンスの増加)。逆もしかりで、予測モデルがシンプルであると理想とは合わなくなりますが(バイアスの増加)、データごとの予測モデルのばらつきは小さくなります(バリアンスの減少)。このように、バイアスとバリアンスはモデルの複雑度と予測精度のトレードオフを反映しており、そのバランスが重要であることをここから学べます。

:::note warn
このバイアス・バリアンス分解の内容については少し理解が怪しいです。特にBLUPにおいて$\mathcal{D}$が果たして本当に$\mathbf{y}$にあたるのか、や、$p(\mathcal{D})$がどんな分布になるのか、に関してははっきりとした答えを持ち合わせていません。ご存じの方がいらっしゃいましたら、ぜひご教授ください
:::