statcastデータにおけるlaunch_speed, launch_angleの欠測値について（中編）

↓の記事の続きです。

tsuyupon.hatenablog.com

後編にしようと思って書いたんですが、思いのほか長くなってしまったので、中編として欠測値補完の手法についてまとめました。

前編の振り返り
欠測値補完の方法
まとめ
参考文献

前編の振り返り

前編では、「statcastデータにおける打球速度（launch_speed）、打球角度（launch_angle）の欠測の扱い」をテーマにデータの現状を確認していきました。

結果として、わかったことは以下の通りです。

欠測には打球内容（詳細不明）ごとの平均値が代入されているらしい。
ただどこが欠測なのかは不明。
集計してみると異常に頻度の高い打球速度、打球角度の組み合わせがあったので、今回はそれを欠測として扱うこととした（全体の15％程度）。
欠測はバッテリー間（バント）、一塁側・三塁側ファールゾーン（ファールフライ）、内野後方（ポップフライ）への打球がほとんど。
さらに共変量を探索していくと以下の変数が欠測確率に寄与していそうだとわかりました。
- hc_x, hc_y（打球位置）
- bb_type（打球種類）
- events（打球結果）
- park（球場）
- zone（投球ゾーン）

「球場ごとの機器の違い」「機器が測定しづらい打球」と複数要因が絡んでいそうです。

では、この欠測を補完する方法はどのようなものがあるのでしょうか？

欠測値補完の方法

欠測メカニズムがMCARでない場合、リストワイズ除去したデータを使って分析してしまうと、偏った推定結果が出てきてしまいます。

前編の結果からMCARではないことがわかったので、なんらかの方法でこの欠測を補完しなければなりません。

欠測補完の方法は大きく分けて「単一代入法」と「多重代入法」があります。

単一代入法

１. 平均値代入法

欠測していないデータから算出した平均値を代入する方法。
代入後に平均値を計算しなおしても平均値は変化しないが、分散は過剰に小さくなってしまう。今回のデータのようにグループごとの平均値を代入した方が良い。

２. 確定的回帰代入法

回帰モデルから算出した予測値を代入する方法。
変数 $Y_i$ が欠測のあるデータ、変数 $X_i$ が欠測のないデータとすると、欠測していない組み合わせのデータ（完全データ）のみを用いて下式のような回帰モデルが構築できる。そしてモデル構築後に、そのモデルから代入値 $\hat{Y_i}$ を算出する。

$\hat{Y _ i} = \hat{\beta_0} + \hat{\beta_1}X_i$

誤差は考慮していないため、完全データから構築した回帰モデルの回帰直線上の点が代入値になる。
回帰代入法によるパラメータ推定量は特にガウス・マルコフの仮定とMARの仮定の両方が満たされる場合、最良線形不偏推定量（BLUE）になる。

３. 比率代入法

切片0、分散不均一な単回帰モデルから算出した予測値を代入値とする方法。
従業員数と給与総額など従業員数が0人のとき給与総額も0円となるような関係で、かつ分散が不均一な場合に用いる。

代入値を計算するのは比較的容易で

$\hat{Y_i} = \hat{\beta_1}X_i$ $\hat{\beta_1} = \frac{\sum Y_i / n}{\sum X_i / n} = \frac{\bar{Y}}{\bar{X}}$

このように単純な1次関数で表現できる。こちらも誤差は考慮していないため、回帰直線上の点が代入値になる。

４. ホットデック法

欠測がある行と完全データの行を比較し、補助変数の値が最も似ている完全データの値をそのまま代入値として使用する方法（似ているデータが複数あった場合はランダムに1つのデータを選択）。
何をもって似ているとするかは通常、距離関数を定義して最近隣法で決めることが多い。量的変数・質的変数問わず簡単に使える。

５. 確率的回帰代入法

ここまでは誤差の考慮しない確定的単一代入法であったが、この手法は先ほどの確定的回帰代入法において誤差を考慮した方法。

【手順】

確定的回帰代入法と同様の方法で $\hat{Y_i}$ を算出
残差 $\hat{u_i} = Y_i - \hat{Y_i}$ を算出
残差の分散を利用して、 $N(0, \sigma _ {\hat{u_i}}^2)$ に従う誤差項を発生
先ほどの $\hat{Y_i}$ に誤差項を加えた値を代入値とする

誤差を考慮しても、必ずしも真値と一致するわけではないが、これによって分布の復元を図っている。

ここまでは単一代入ということで欠測に対し、代入値をひとつだけ考えてきました。
最終的に求めたいものが算術平均や分散であれば確率的回帰代入でも大きな問題はないですが、代入した値を使って回帰分析していく場合は代入値のバラつきを考慮できるよう次の多重代入法を使うべきです。

多重代入法

そもそも代入法の目的は個別の値の完全復元ではなく、あくまで母集団パラメータの推定であるため、欠測値代入をしたことによって生じる誤差を適切に評価しなければなりません。その問題を解決する方法としてRubinによって提案された方法が多重代入法です。

多重代入法とは、欠測のあるデータの分布から独立かつ無作為に抽出されたM個（M>1）のシミュレーション値によって欠測値を置き換えるものです。

シミュレーションをM回することで欠測データの不確実性を標準誤差として反映し、その後の統計的推測が出来るようになります。

しかし、実際は欠測データは測定できていないので、分布も正確にわかりません。そこで、代わりにMAR（or MCAR）を仮定し、完全データから欠測データの事後予測分布を構築し、そこから独立かつ無作為な抽出をすることによってシミュレーションを行っていきます。

多重代入法は、先ほどの確率的単一代入法を単純に複数回繰り返したものとは全く異なります。確率的単一代入法は誤差を考慮したひとつのモデルで、回帰係数の安定性を評価できるものではありません。

欠測データの事後予測分布から無作為抽出して都度モデリングするため、シミュレーションごとに異なる回帰係数を持った回帰モデルになります。そのため、コチラの図2.1にあるように、データの補完から目的の分析までシミュレーションごとに実施し、最後の最後に結果を統合します。

多重代入法のアルゴリズム

では、どのようにシミュレーションを実施するのか。
多重代入法のアルゴリズムは大きく分けて３つあります。

DA（data augmentation：データ拡大）アルゴリズム

DAアルゴリズムはマルコフ連鎖モンテカルロ法（MCMC）に基づくアルゴリズム。
代入ステップ（式1）で観測値 $Y_{obs}$ と繰り返し時点 $t$ におけるパラメータ値 $\theta^{(t)}$ を条件とした欠測値の予測分布から代入値 $Y^{(t+1)}$ を生成することで疑似的にデータを拡大し、事後ステップ（式2）で観測値 $Y$ と繰り返し時点 $t+1$ における代入値 $Y^{(t+1)}$ を条件として、事後分布からパラメータ値 $\theta^{(t+1)}$ を生成する。これらのステップを収束するまで繰り返す。
※misは欠測、obsは観測データ

式1： $Y _ {mis}^{(t+1)}\sim Pr(Y _ {mis}\mid Y_{obs},\theta^{(t)})$

式2： $\theta^{(t+1)} \sim Pr(\theta \mid Y _ {obs}, Y _ {mis}^{(t+1)})$

この最終的に収束した値がそのシミュレーションにおける代入値になり、このプロセスを通常100回程度繰り返し、100通りの代入値を生成する。
※なお、初期値 $\theta^{(0)}$ はEMアルゴリズムによって得られる値を推奨している。

このアルゴリズムでは多変量分布の条件付分布から代入値を生成しているが、もし真の同時分布が多変量正規分布で近似できるならば、統計分析も妥当なものになるとのこと(Drechsler, 2009)。

Rではnorm2パッケージを使って実行可能だが、CRANから削除されているので、アーカイブからインストールする必要がある。

FCS（fully conditional specification：完全条件付き指定）アルゴリズム

DAアルゴリズムは多変量分布を仮定する手法だったが、FCSアルゴリズムでは多変量欠測データの補完を変数ごとに行う。つまり、欠測のある変数ごとにモデルを構築し、代入値を作成できる。そのため、適切な多変量分布が存在していなくても補完可能となる。

例えば、 $Y$ を欠測のないデータ、 $X_1,X_2$ を欠測のあるデータとすると、
$X_1$ の欠測を補完するモデルは式3、 $X_2$ の欠測を補完するモデルは式4のように表せる。

式3： $X _ 1^{(t+1)} \sim f(X _ {1,mis} \mid X_{1,obs},X_2^{(t)},Y)$

式4： $X _ 2^{(t+1)} \sim f(X _ {2,mis} \mid X_{1}^{(t+1)},X _ {2,obs}^{(t)},Y)$

そして、この例の場合、欠測補完の手順は式3、式4を用いて以下のようになる。

観測データ $X _ {1,obs}$ から無作為に代入値の初期値 $\tilde{X} _ 1^{(0)}$ を抽出する。
観測データ $X _ {2,obs}$ から無作為に代入値の初期値 $\tilde{X} _ 2^{(0)}$ を抽出することで疑似完全データを作る。
式3をもとに疑似完全データのうちもともと $X _ 1$ で欠測していない行からモデル（例えば重回帰モデル）を構築し $X _ {1,mis}$ の代入値を更新。
式4をもとに疑似完全データのうちもともと $X _ 2$ で欠測していない行からモデル（例えばロジスティック回帰モデル）を構築し $X _ {2,mis}$ の代入値を更新。
3., 4.のプロセスが収束するまで繰り返す。

そして、この手順全体をM回（M＞1）繰り返し、M通りの完全データを生成する。

Rではmiceパッケージを使って実行可能。

今回の欠測は打球速度、打球角度が同時に欠測しているため、このアルゴリズムは適さなそう。

EMB（Expectation-Maximization with Bootstrapping：ノンパラメトリック・ブートストラップを用いた期待値最大化）アルゴリズム

まず、ノンパラメトリック・ブートストラップ法では、観測された標本データを擬似的に母集団として扱い、標本サイズ $n$ と同じサイズの副標本 $S _ {boot}$ を無作為復元抽出（欠測が含まれていても良い）し、これをM回（M＞1）行う。

このM個の副標本のうち欠測の含まれているデータに対して以下のEMアルゴリズムを利用して、欠測を補完する。

EMアルゴリズムでは、まず分布を仮定しパラメータの初期値を仮に設定する。この初期値を使用して、モデル尤度の期待値を計算し、尤度を最大化するモデルパラメータを更新する。そして、この期待値ステップ(Expectation Step: E-Step)と最大化ステップ(Maximization Step: M-Step)を繰り返し値が収束するまで実施する。

E-step： $Q(\theta \mid \theta^{(t)}) = \int l(\theta \mid Y) Pr(Y _ {mis} \mid Y _ {obs},\theta^{(t)})dY_{mis}$ 　※ $l(\theta \mid Y)$ は対数尤度