Rによる多重代入法の学習記録 - 霞と側杖を食らう

【学習動機】

Kaggleで欠測値を埋めなきゃということがある。他にも埋めなきゃいけない状況がある。
欠測値の理論は、因果推論の一般的な枠組みとして捉えるのに役立つ。
あと、ベイズと関わってくるので、勉強しておきたいと思った。実際に、Rで実行できるように。

【学習記録】

高橋・渡辺『欠測データ処理 Rによる単一代入法と多重代入法』で学んでいくこととする。

＜総括＞

Rで実際にどうやって欠測値に多重代入法を用いて対処するかが分かったのは、とても良い。
今出ている日本語書籍では、理論を学ぶことはできるが、どれをどうやって使うのかがすぐに分からなくて困っていた。
データセットは書籍HPに用意すべき。データは書面に記載しているってのはモダンじゃない。
norm2とmiceとAmelia。代入・分析・統合。

ｃ1：Rによるデータ解析

Rの使い方と欠測あるときの反応を少しチェックしている。

ｃ2：不完全データの統計解析

多くの統計ソフトウェアでは、リストワイズ除去で、欠測値を含む観測対象の行をすべて削除し、擬似的「完全」な状態にして分析することが多い。
欠測パターンと欠測メカニズムで分類される。
欠測パターンには、単変量欠測パターン、単調欠測パターン、計画的欠測パターン、一般的な欠測パターン。多重代入法では、いずれの欠測パターンにも有効で、欠測パターンはそんなに重視されない。
欠測メカニズムには、完全に無作為な欠測(MCAR：Missing Completely At Random)と条件付きで無作為な欠測(MAR：Missing At Random)と無作為ではない欠測(NMAR：Not Missing At Random)がある。 Dをn×pのデータセット。 $D =$ { $D _ {obs}$ , $D _ {miss}$ }。Kは回答指示行列でn×p行列。 $K _ {i j}=1$ で観測、 $K _ {ij}=0$ で欠測
・MCAR
$Pr(K|D) = Pr(K)$ であり、欠測する確率が、データと独立。
・MAR
$Pr(K|D) = Pr(K|D _ {obs})$ で、観測データを条件とし、欠測確率の分布が非観測データから独立。欠測データがMARならば、欠測を除去する分析は偏っている。この偏りは共変量(covariate)または補助変数(auxiliary variable)による代入で是正可能。・NMAR
$Pr(K|D) \neq Pr(K|D _ {obs})$ で、欠測する確率が欠測した変数の値自体に依存している。代入法で偏りを是正できるとは限らず、選択モデルやパターン混合モデルを用いた分析が必要で、感度分析が有用。

分析の目的が母集団の推定なら、多重代入法を使用すべき。

ｃ3：単一代入法

従来、公的統計では、調査データの集計が主目的で、確定的単一代入法を用いることが通例。方法として、よく用いられるのが、回帰代入法、比率代入法、平均値代入法、ホットデック法。
・確定的単一代入法はガウス・マルコフの仮定あるとよい。
・比率代入法は不均一分散のとき使える。
・平均値代入法は百害あって一利なし。
・ホットデック法(hot deck imputation)は、ノンパラな代入法。ホットデックでは、欠測している人の補助変数の情報を利用して、回答者から似通った値のものを選び出し、その回答者の値を代入値として採用。何を基準にするかは議論あり。最近隣法が多い。
HotDeckImputationパッケージのimpute.NN_HD()で実行できる。
・確率的回帰代入法
各々の代入値に誤差項を加える。誤差項は、平均が０、分散が回帰モデルの残差分散の正規分布から無作為抽出するなどのやり方があり、miceのmice()で引数methをnorm.nobとすれば実行できる。

ｃ4：多重代入法の概要

ガウス・マルコフの仮定とMARの過程がすべて満たされた場合、多重代入法によるパラメータ推定量はBLUEになる。
代入、分析、統合の３ステージ
統合方法。 $\tilde{\theta_m}$ をパラメータ $\theta$ のm番目の代入済みデータセットに基づいた推定値とする。
統合した点推定値は
$\bar{\theta_M} = \frac{1}{M} \Sigma \tilde{\theta_m}$
多重代入法の代入値の分散は代入間分散(between)と代入内分散(within)に分解でき、
代入内分散は、
$\bar {W _ m}$ $= \frac{1}{M}$ $\Sigma var(\tilde{\theta_m})$
代入間分散は、
$\bar{B _ M} = \frac{1}{M-1}$ $\Sigma$ $(\tilde {\theta _ m}$ - $\bar{\theta_M}) ^2$
$\bar{\theta_M}$ の分散 $T_M$ は、これらを合算して
$T_M = \bar{W_m} + (1 + \frac{1}{M} ) \bar{B_M}$

この方法をRubin(1987)のルールという。

代入間の繰り返しをどれだけ行えば収束するかは、様々な要素が絡んでくるが、とりわけ欠測情報の比率が重要。
適切な多重代入法(proper multiple imputation)でないと、推測の結果が妥当でなくなるおそれあり。
代入モデルと分析モデルが同一の変数を持ち、同じ数のパラメータを推定するとき、２つのモデルには適合性(congeniality)があるという。適合していない場合、多重代入法のパラメータ推定量の一致性は保証されない。
代入モデルよりも大きな分析モデルを使うと、リストワイズ除去より悪い結果になりうる。
また、多重代入済みデータ数は多い方が望ましいとされる。

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

3つのアルゴリズム。
・DAアルゴリズム。norm2でできる。しかし、norm2が消えてる。
・FCS アルゴリズム。miceでできる。これが主流か？
・EMBアルゴリズム。Ameliaでできる。 (EMアルゴリズムについてはまた別の機会で詳しく学ぶ)

norm2パッケージがCRANから消えてるので

url <- "https://cran.r-project.org/src/contrib/Archive/norm2/norm2_2.0.1.tar.gz"
pkgFile <- "norm2_2.0.1.tar.gz"
download.file(url = url, destfile = pkgFile)
install.packages(pkgs=pkgFile, type="source", repos=NULL)
unlink(pkgFile)

あと、ホットデックのコードを動かした後、Ameliaパッケージの多重代入するとエラー起きるかもしれないらしい。
Amelia、norm2、miceと見比べて、すべての変数が量的な連続変数なら、どれも変わらず、miceが時間的に効率悪いか。
質的データや時系列データになると、アルゴリズムの使い分けが必要。

ｃ6：多重代入モデルの診断

MARが正しいと仮定したとき、代入の結果がMARの仮定と整合性があるかを診断する。
MCARの場合を除くと、観測データの分布と代入データの分布は一致せず、むしろ、異なっていることが期待される。診断で、分布の違いを発見し、MARの仮定と整合性があるか確認する。分布が極端に違う場合も、どの変数に問題があるか見つけるのに診断ツールが役立つ。 Ameliaとmiceの診断ツールを用いる。

ｃ7,8：量的データの多重代入法

・多重代入法を用いた平均値のt検定 AmeliaとMKmiscを使うか、miceのpool.scalar()の返り値から直接計算。
欠測データによって、自由度の計算が大変。分散の比率λと代入の不確実性の分散の相対的な増加rを計算する。

このへん、本のコードが変。本が書かれた後にアップデートされて関数が変わったのか知らないが。 miceのpool.scalar()は引数にmethodはとらないし、返り値にlambdaはない。
rは返ってくるので、 $\lambda = r/(1+r)$ で自分で計算しないとかと。

・重回帰分析
回帰診断で、ジャック・ベラの正規性検定(JB検定)で誤差項の正規性を検証、ブルーシュ・ペイガン検定(BP検定)で不均一分散の検証、分散拡大要因(VIF:Variance Inflation Factor)で多重共線性の検証、ボンフェロー二外れ値検定で外れ値の検証。

ｃ9,10：質的データの多重代入法

・ダミー変数のある重回帰分析質的データの代入法に関して、決定打はないようだ。ここでは、応用研究で評判が高いとされるmiceによる方法とhot.deckの多重ホットデックによる方法が紹介されている。

・ロジスティック回帰
miceでいつも通り代入して、glm.mids()の引数familiyをbinomialとすればよさげ。

ｃ11：時系列データの多重代入法

一般的に、時系列データにおける欠測値は、補完(interpolation)やカルマンフィルタを用いて処理されることが多い。
一方で、多重代入法を用いて欠測値を処理することで不確実性を反映した分析を可能にする。

ｃ12：パネルデータの多重代入法

Ameliaで実行可能。amelia()の引数csに横断面の単位、引数tsに時系列の単位、引数polytimeに3以内の多項式を指定、引数lagsにラグ付き従属変数を指定。未来の値から過去の値を予測することで代入モデルの精度を向上させたいなら引数leadsを指定してもよい。
Ameliaに入っているafricaのデータを用いて、plmでパネルデータ分析をする。
p151のコードを回してみたが、変量効果モデルの推定がうまくいかない。特異行列が出て、solve()で逆行列を求めるところがうまくいってないみたいだが。。。

ｃ13：感度分析：NMARの統計分析

SensMice(SensMiceDA)パッケージを使えば、感度分析は可能。感度分析の目的は、欠測データに関する設定を変えて、分析結果にどう影響するかを調べること。デルタ修正のパラメータデルタをどう設定するかが重要。デルタの設定方法として、変数の標準偏差を基準とする方法がある。標準偏差の-1倍, -0.5倍, 0倍, 0.5倍, 1倍した5つのデルタ値を用意して、どう影響するかを見て、ドメイン知識と照合して検討。

ｃ14：事前分布の導入

Ameliaならamelia()の引数boundsに(変数の位置、最小値、最大値)を行にした行列を指定することができる。リッジ事前分布を使ってモデルの安定性を実現したいなら、amelia()の引数empriで指定。
norm2ならば、emNorm()のprior="ridge",prior.df=0.5などとする。miceはsqueeze()を用いる。