Little, R. (2011) Calibrated Bayes, for Statistics in General, and Missing Data in Particular (with comments and a rejoinder). Statistical Science, 26(2), 162-186.
統計学者Little先生があちこちで提案している Calibrated Bayes アプローチについて調べていて、その一環として読んだ奴。
良く引用されるLittle(2011)は2012年にメモをとりながら読んでいたのだが、私には話が大きすぎ、いまいち雲をつかむような感じでよく分からなかった。この論文は「欠損データ」と問題が狭く指定されているので、もう少しわかりやすいかと思ったのだが…
この論文には2人の論者によるコメント、そして著者の返答がついている。そちらも合わせてメモする。
1. イントロダクション
[本論文の構成。省略]
ベイジアン手法は欠損データの扱いに際して特に有用である。[最尤法と比較した短い説明。省略]
2. なぜベイズか? calibrated bayesの哲学
[頻度主義 vs ベイジアン]
統計学の世界は、いまだに頻度主義とベイジアンに大きく分かれている。
頻度主義は未知パラメータ\(\theta\)の推論を、抽出の繰り返しにおける統計量の分布から導出した仮説検定なり信頼区間なりに基づいて行う。
ベイジアンははデータについてのモデルと未知パラメータについての事前分布を形成し、未知パラメータについての推論を事後分布に基づいて行う。ベイジアンは「主観的」だったり「客観的」だったりする(文脈に応じてどちらも有用である)。漸近的最尤推論は大標本ベイズの一形式とみることができる。つまり\(\theta\)の区間を信頼区間というより事後信用区間として解釈できる。
[ここから頻度主義ディス＆ベイジアン推しが始まるんだけど、項目立てせずにダラダラと続くので困ってしまう。きちんとした議論は Little(2006) をみてくれ、ということであろう(イントロでそういう書き方をしていた)。2006年論文を参考に思い切って要約すると、こういうことだと思う。頻度主義の欠点:
- prescriptiveでない。推論システムそのものというより、推論手続きの諸特性を評価する一連の概念に過ぎない。
- 多くの小標本問題に対する解がない。たとえばBehrens-Fisher問題。
- 尤度原理に反する。
- 最尤法はnuisanceパラメータも最大化してしまい、それらのパラメータの不確実性を考慮にいれない。
ベイジアンの欠点:
- 正しいモデルを完全に指定しないといけない。
- 仮説検証のツールが不十分。
って感じかな…]
要約しよう。ベイジアンは仮定されたモデルの下での推論に優れ、頻度主義はモデル開発と評価に優れている。
[Calibrated Bayes]
双方のいいとこ取りをするのがCalibrated Bayes(CB)アプローチである。
- Box(1980)はこう定式化した。$$ p(Y, \theta) | \mathrm{Model}) = p(Y | \mathrm{Model}) p(\theta | Y, \mathrm{Model}) $$ 右辺の\(p(\theta | Y, \mathrm{Model})\)は、データとモデルの下でのパラメータの事後分布である。これが推論の基盤となる。\(p(Y | \mathrm{Model})\)はモデルの下でのデータの周辺分布で、モデルの妥当性検証に用いられる(頻度主義的)。特に、周辺分布と観察データとの乖離の評価を行う。
後者の”posterior predictive checks”は本当に事前分布の選択に対して敏感なのかどうか疑わしい。事前分布の選択が事後の推論に限られたインパクトしか持たないときには敏感でないかもしれない。これは特に、モデルが無情報事前分布を含むときに問題となる。 - Rubin(1984)は、モデルのチェックは\(p(Y | \mathrm{Model})\)ではなくて、モデルとデータの下での将来のデータ\(Y^*\)とパラメータ\(\theta^*\)の予測分布 \(p(Y^*, \theta^* | Y, \mathrm{Model})\)に基づいて行うべきだと考えた。これがposterior predictive checksである[以下PPCと略記]。PPCは観察データを、モデルの下で生成し得る未観察データの系列のなかに埋め込んだとき、観察データが「筋がとおっている」かどうかに関心を向ける。精神としては頻度主義だといえる。
哲学としての頻度主義がベイズとは完全に異なるとしても、頻度主義はMCMCのようなベイズ計算ツールに適用できるし、結果をMLの近似として解釈することができる。
[雲をつかむような話で困るなあ…]
3. 欠損データの統計的分析：小史
- EMアルゴリズム以前(-60年代)。欠損データは欠損行を削ったり、単純に補完されたりしていた。ANOVAで欠損があるとき、釣り合い計画を維持するために欠損を埋めるというのが広く行われていた(いまでは不釣り合い計画でも計算できるので、この方法は歴史の遺物となった)。最尤法では、たとえば二変量正規データで一方に欠損があるときに尤度を分解して反復なしに解くというような方法が提案された。もっと複雑な問題について反復推定する方法も提案された。
- 最尤法の時代(70年代-80年代中盤)。EMアルゴリズムが登場し、最尤法が一般的になった。欠損メカニズムのモデリングと、欠損メカニズムを無視できる十分条件の研究が発展した。
- ベイズと多重代入の時代(80年代中盤-現在)。欠損データ分析は最尤法からベイジアンに移行した。きっかけは、欠損データによるMVNモデルのパラメータの事後分布をdata augmentationにより生成するという提案の登場である。また、Rubinによる多重代入の提案が広まった。
- 頑健性への関心(90年代-現在)。モデル誤指定に対して頑健な手続きへの関心が欠損データ問題に拡張された(たとえば、逆確率ウェイトつき方程式に基づく二重頑健法)。ベイジアンの観点からは、頑健性は構造的過程が比較的に弱いモデルの構築という形をとる(たとえばPenalized Spline of Propensity Prediction)。また、標準的欠損データモデルのためのモデルチェックにも関心がもたれるようになった。
4. 欠損データの尤度ベース手法
尤度法は欠損データモデルに直接適用できる。
\(Y=(y_{ij})\)を\(n\)行\(K\)列のデータ行列とする。\(M = (m_{ij})\)を同じ大きさの欠損インジケータ行列とする。\(Y\)の観察された部分を\(Y_{obs}\), 欠損部分を\(Y_{mis}\)とする。
[(Y, M)の分解]
\((Y, M)\)の分布を分解する完全にパラメトリックなモデルが二種類ある。
- selection model。\((Y, M)\)の分布を、\(f(Y | \theta)\)と\(f(M | Y, \psi)\)に分解する。
- pattern-mixture model。\((Y, M)\)の分布を、\(M\)の周辺分布と、\(M\)の下での\(Y\)の周辺分布に分解する[\(f(M | \psi)\)と\(f(Y | M, \theta)\)ってことかな]。
以下ではselection modelについて考える。
[欠損データへの2つのアプローチ]
もし\(Y\)が完全に観察されていたら、パラメータの事後分布は $$ p_{complete} (\theta, \psi | Y, M) = \mathrm{const.} \times \pi(\theta, \psi) \times L(\theta, \psi | Y) $$ $$ L(\theta, \psi | Y) = f(Y | \theta) \times f(M|Y, \psi)$$ となる。では、\(Y\)が不完全だったらどうなるか。次の2つのアプローチがある。
- パラメータのフルの事後分布について考える。$$ p_{full} (\theta, \psi | Y_{obs}, M) \propto \pi(\theta, \psi) \times L(\theta, \psi | Y_{obs}, M) $$ $$ f(Y_{obs}, M | \theta, \psi) = \int f(Y_{obs}, Y_{mis} | \theta) f(M | Y_{obs}, Y_{mis}, \psi) d Y_{mis} $$
- 欠損データメカニズムを無視する。$$ p_{ign}(\theta | Y_{obs}) \propto \pi(\theta) \times L(\theta | Y_{obs}) $$ $$ L(\theta | Y_{obs}) = \int f(Y_{obs}, Y_{mis} | \theta) d Y_{mis}$$
アプローチ2のほうが扱いやすい。なぜなら、(a)欠損データメカニズムを指定しなくてよい、(b)欠損データの積分がアプローチ1よりも楽, (c)フルモデルは往々にして識別できない(ベイジアンの場合でいうと事前分布の選択に強く依存する)。
というわけで、アプローチ2を正当化できる十分条件が問題となる。それは以下であることが知られている。
- Missing at Random: すべての\(Y_{mis}\)について \(p(M|Y_{obs}, Y_{mis}, \psi) = p(M | Y_{obs}, \psi)\)
- アプリオリ独立: \( \pi(\theta, \phi) = \pi(\theta) \times \pi(\phi) \)
[EMアルゴリズム]
アプローチ1, 2のどちらの場合でも、もっとも重要なのはモデルの選択と計算上の問題である。最尤法の世界では、EMアルゴリズムが、複雑な観察データ尤度と単純な完全データ尤度の橋渡しをしてくれる。
欠損データメカニズムが無視可能であるとしよう。\(\theta\)の現在の推定値を\(\theta^{(t)}\)とする。Eステップで、もし\(\theta\)が\(\theta^{(t)}\)と等しかったら、完全データの期待値がどうなるかを求める。$$ Q(\theta|Y_{obs}, \theta^{(t)}) = \int \log f(Y_{obs}, Y_{mis}; \theta) \cdot f(Y_{mis} | Y_{obs}, \theta = \theta^{(t)} ) d Y_{mis} $$ Mステップで、完全データ対数尤度の期待値を最大化する\(\theta\)を求めこれを\(\theta^{(t+1)}\)とする。つまり、すべての\(\theta\)について $$ Q(\theta^{(t+1)} | Y_{obs}, \theta^{(t)} ) \geq Q(\theta|Y_{obs}, \theta^{(t)}) $$
[data augumentation]
ベイジアンならEMのかわりにdata augumentationが使える。これは一種のGibbsサンプラーである。
[以下、ドローされた値には上添字に\(d\)がついている]
欠損データメカニズムが無視可能であるとしよう。イテレーション\(t\)におけるドローを\(Y^{(dt)}_{mis}, \theta^{(dt)}) \)とする。次のイテレーションでは、まず$$ Y^{(d,t+1)}_{mis} \sim p(Y_{mis} | Y_{obs}, \theta^{(dt)})$$をドローし(Eステップに相当)、次に $$ \theta^{(d,t+1)} \sim p(\theta | Y_{obs}, Y^{(d, t+1)}_{mis}) $$ をドローする(Mステップに相当)。
5. 多重代入
上で \( Y^{(d)}_{mis} \) をドローしたが、これをつかって多重代入データセットが作れる。
引き続き、欠損データメカニズムが無視可能であるとしよう。$$ p(\theta | Y_{obs}, Y_{mis}) \propto \pi(\theta) \times L(\theta | Y_{obs}, Y_{mis}) $$ より、$$ p_{ign}(\theta | Y_{obs}) = \int p(\theta | Y_{obs}, Y_{mis}) p(Y_{mis}|Y_{obs}) d Y_{mis}$$ つまり、まず\(Y_{mis}\)の事後分布\(p(Y_{mis}|Y_{obs})\)から\(Y^{(d)}_{mis}\)をドローして、次に\( p(\theta | Y_{obs}, Y^{(d)}_{mis}) \) から\(\theta\)をドローするのを繰り返し、$$ p_{ign}(\theta | Y_{obs}) \approx \frac{1}{D} \sum_d^D p(\theta | Y_{obs}, Y^{(d)}_{mis}) $$ を得ればいいわけだ。[…事後分布の期待値と分散をどう求めるか。中略…]
頻度主義の観点からは、パラメトリックモデルのためのベイジアン多重代入は、最尤法と類似した大標本特性を持ち、しかも計算は楽である。さらに、代入モデルには分析モデルに含めない変数を含めることもできる。[…後略…]
6. 欠損データ問題へのベイジアン手法の適用
4-5章では欠損データのベイジアン推論の基礎を紹介した。これらの手法の適用例を示そう。
事例1. 欠損値が単調パターンを持つデータ
二変量正規変数\((Y_1, Y_2)\)のデータ(サイズ\(n\))があって、\(Y_1\)は完全で、\(Y_2\)には欠損があり値は\(r\)個だとしよう。欠損メカニズムは無視できるとする。
同時正規分布を、\(Y_1\)の周辺分布\(N(\mu_1, \sigma_{11}^2)\)と、\(Y_2\)の条件付き分布\(N(\beta_{20,1} + \beta_{21,1} Y_1, \sigma^2_{22,1})\)に分解する。
最尤法の場合、尤度は\(\phi_1 = (\mu_1, \sigma_{11})\) (\(n\)件がベース)と\(\phi_2 = (\beta_{20,1}, \beta_{21, 1}, \sigma_{22,1})\) (\(r\)件がベース) に分解できる。ここから最尤推定し、$$\hat{\mu}_2 = \hat{\beta}_{20,1} + \hat{\beta}_{21,1} \hat{\mu}_1 $$とすればよろしい。
ベイジアンならどうするか。\(\phi_1, \phi_2\)に共約事前分布を与え、データを使って事後分布をつくり、そこからパラメータをドローして上の式にいれて\(\mu_2\)の事後分布をつくればよい。
ベイジアンのほうが以下の点で優れている。(1)パラメータについての事前知識があるとき、それを事前分布に組み込める。(2)小標本の場合、事後分布は不確実性をうまく捉えている。(3)不確実性の推定がすぐに手に入る。
今度は、変数が\(Y_1,\ldots,Y_K\)で、\(Y_{k+1}\)が観察されていたら\(Y_k\)はかならず観察されており(欠損パターンは単調)、\((Y_k | Y_1, \ldots, Y_k-1)\)の条件付き分布は未知パラメータ\(\phi_k\)を持っていて、\((\phi_1, \ldots, \phi_K)\)は独立な事前分布を持つとしよう。
ベイジアンなら、観察されたデータに基づいて\(\phi_k\)の事後分布をつくることができる。なお、パラメータが事前に独立でない場合、ないし欠損パターンが単調でない場合は、MCMCのような方法が必要になる。
事例2. 一般的なパターンの欠損を持つMVNモデル
$$ y_i = (y_{i1}, \ldots, y_{ip}) \sim_{ind} N_p (\mu, \Sigma) $$ なのだが\(y_i\)に欠損があるとしよう。非欠損部を\(y_{obs,i}\), 欠損部を\(y_{mis,i}\)とする。
いまパラメータのドロー$$ \theta^{(dt)} = (\mu^{(dt)}, \Sigma^{(dt)}) $$があれば、欠損値のドロー$$ y^{(d, t+1)}_{mis, i} \sim p(y_{mis, i} | y_{obs, i}, \theta^{(dt)}) $$ が可能である。sweep演算子を使えばよい。[勉強不足でなにいってんのかわからん。ドローされた\(\mu, \Sigma\)を、\(i\)さんの観察値ベクトルでベイズ更新して、それをパラメータにしたMVNからドローすればいいじゃんという話かなあ? Little & Rubin(2002, 書籍)をみよとのこと]
で、次は埋めたデータから事後分布(\(\Sigma\)は逆ウィシャート分布, \(\mu\)はMVN)を得てパラメータをドローする。$$ \Sigma^{(d, t+1)} \sim p(\Sigma | Y_{obs}, y^{(d, t+1)_mis}) $$ $$ (\mu^{(d, t+1)} | \Sigma^{(d, t+1)}) \sim p(\mu | \Sigma^{(d, t+1)}, Y_{obs}, Y^{(d, t+1)}_{mis}) $$ 詳細はSchaefer(1997)をみてほしい。なおこれがSASのPROC MIのアルゴリズムの基礎になっている。これはGibbsアルゴリズムなんだけど、EMと似ている。
頻度主義者ならまず最尤推定値を得て次に情報行列からSEを求めるだろう。上記のGibbsアルゴルズムのほうがシンプルである。頻度主義者のみなさんも、Gibbsアルゴリズムのドローを使ってパラメータの検定をしたり信頼区間を求めたりできる。
このやり方の限界として、欠損を持つ変数群がMVNに従うという強い仮定を置いている点がある。つまり欠損値を非欠損値に線形回帰するモデルを考えているわけで、交互作用があるとか二値データとかだと困る。その解決策のひとつとして Sequential regression multivariate imputation がある。事例3をご覧下さい。
事例3. Sequential regression multivariate imputation (SRMI)
\(Y\) (\(K\)列)は欠損あり、\(X\)はなしとする。\(Y_j\)以外の全変数のもとでの\(Y_j\)の条件付き分布のモデル\( g_j(Y_j|X, Y_1, \ldots, Y_{j-1}, Y_{j+1}, \ldots, Y_K, \phi_j) \)を考え、パラメータ\(\phi_j\)に無情報事前分布を与える。
イテレーション\(t\)における\(i\)さんの\(Y_{j}\)の値(観察値または補完値)を\(y^{(t)}_{ji}\)とする。イテレーション\(t\)で\(Y_j, \ldots, Y_K\)を補完し、イテレーション\(t+1\)で\(Y_1,\ldots,Y_{j-1}\)を補完したデータ行列を\(Y^{(jt)}\)とする。[あるイテレーションで列1から列Kまでの欠損を順に埋めていき(どうやって埋めるかはこれから説明する)、イテレーション\(t+1\)で\(j\)まで進んだらその瞬間のスナップショットを\(Y^{(jt)}\)とする、ってことだろうか]
\(j = 1, \ldots, K\)について、\(Y_j\)の欠損を次のように埋める。まず\((X, Y^{(jt)})\)のもとでの\(\psi\)の事後分布からドローする。$$ \psi^{(t+1)}_j \sim p(\psi | X, Y^{(jt)}) $$ つぎに欠損を埋める。\(i=1,\ldots,n\)について、もし\(y_{ji}\)が欠損だったら$$ y^{(t+1)}_{ji} \sim g_j (Y_j | X, y^{(t+s1)}, \ldots, y^{(t+1)}_{j-1,i}, y^{(t)}_{j+1,i}, \ldots, y^{(t)}_{Ki}, \psi^{(t+1)}_j ) $$ [ \(y^{(t+s1)}\)というのが全くわからない。\(y^{(t+1)}\)の誤植ではなかろうか。それにしても、ややこしくてわからん。擬似コードを見せてほしいよ]
補完が安定するまでイテレーションを繰り返す。チェーンを複数走らせて収束を評価する。これを\(D\)回繰り返し\(D\)個の補完データをつくる。
SRMIは、多変量欠損データという問題を、他の全ての変数を所与とした単変量欠損データという問題へと縮小することで、それぞれの変数のモデルに柔軟性を与える。早い話、非線形項や交互作用項を入れたり、たとえば二値変数ならロジスティック回帰にしたりできるわけです。それぞれのモデルは、\(X\)のもとでの\((Y_1, \ldots, Y_K)\)の同時分布とは整合しないが、実務的にはたいした問題でない。ソフト実装としてはMICEがある。
事例4. Penalized spline of propensity prediction (PSPP) によるロバスト・モデリング
単一の変数\(Y\)に欠損があるとしよう。\(i\)さんの観察ベクトルは\((x_{i1}, \ldots, x_{ip}, y_i)\)なんだけど、あいにく\(y_i\)に欠損があるとする(欠損の人は\(i=r+1, \ldots, n\)だとする)。欠損確率は\(x_{i1}, \ldots, x_{ip}\)に依存し\(y_i\)に依存しないとする(つまりMAR)。欠損インジケータを\(m_i\)とする。\(\mu_y = E(Y)\)を推定したい。
観察プロペンシティ\(M\)を、\(X_1, \ldots, X_p\)への回帰で推定したとする(ロジスティック回帰でもプロビット回帰でもなんでもいい)。その逆数をウェイトにして完全ケースのみで平均を求めることができるが、分散を小さくするために、残差を重みづけたカリブレーション項を入れた方が良いことが知られている。
この文脈では、以下のどちらかが成り立っているだけで一致推定量となる推定量のことを二重頑健(DR)推定量という。(a)\(Y\)を\(X_1, \ldots, X_p\)で予測するモデルが正しく指定されている。(b)回答プロペンシティのモデルが正しく指定されている。
DR特性を持つベイジアン欠損データ法を紹介しよう。プロペンシティ予測の罰則付きスプライン(PSPP)である。
観察プロペンシティのロジットを $$ p^*_i(\psi) = \mathrm{logit}(Pr(m_i = 0 | x_{i1}, \ldots, x_{ip}; \psi) $$ と定義する。
MARの仮定の下で、欠損が\(x_{i1}, \ldots, x_{ip}\)だけに、プロペンシティ・スコアだけを通じて依存することをバランシング特性という。この特性のもとで、$$ \mu_y = E[(1-m_i)y_i] + E[m_i \times E(y_i | p^*_i(\psi))] $$ と書ける。
そこで以下のようにモデル化する。$$ (y_i | @^*_i(\psi), x_{i1}, \ldots, x_{ip}; \psi, \beta, \phi, \sigma^2) \sim N(spl(p^*_i(\psi), \beta) + g(p^*_i, x_{i2}, \ldots, x_{ip}; \phi), \sigma^2) $$ 正規分布の平均の第1項 \(spl(p^*_i(\psi), \beta)\)はプロペンシティ・スコアのスプライン関数である。どんな関数でもいいけど、我々は罰則つきスプラインを選ぶ[式省略]。第二項\(g(p^*_i, x_{i2}, \ldots, x_{ip}; \phi)\)はパラメトリック関数である(多重共線性を避けるため\(x_{i1}\)だけ抜いてある)。
推定手順はこうである。まずプロペンシティ・スコアをどうにかして推定する。次に上のモデルで\(Y\)を推定する。[軽くいうねえ… 一気に推定できるんだろうか？]
このモデルで欠損を埋めた\(Y\)の平均はDR特性を持つ。
このモデルでは\(y_i\)に多重代入できるが、その際のパラメータはたくさんある。スプラインモデルのパラメータ、g関数の回帰係数と残差分散、そしてプロペンシティ関数の剰余パラメータ\(\psi\)である。これらのパラメータの不確実性を表現するには、Gibbsサンプラーでパラメータを事後分布からドローするのがよい。
7. 結論
[ううう。この論文、話が逸れまくっていないか？ どういう落ちがついているのか？ 不安に駆られて、ここからほぼ逐語訳]
本論文では、まず統計的推論におけるCBアプローチを支持する議論を要約した。私の思うところ、CBアプローチはベイズ統計学と頻度主義統計学のいいとこ取りである。簡単にいうと、推論はベイジアンでなければならないが、モデル構築とチェックのためには、ベイジアン推論の結果が持つ頻度主義的な諸特性に注意する必要がある。
CBのこの「ロードマップ」は完全な解とはいえない。なぜなら、モデル構築とモデル推論の間には相互作用があるからだ。そこにはたとえば、モデル不確実性のうちどこまでをフォーマルな統計的推論のなかに含めなければならないのかといった、雲をつかむような問いが含まれている。とはいえ、CBアプローチは実務的な統計的推論に接近していくための非常に良い基盤となると思われる。
本論文の残りの部分で、CBアプローチが欠損データ問題を扱う際にとりわけ魅力的だと論じたい。ある意味で、すべての推測統計学は欠損データ問題であり(未知のなにかについての推論なのだから)。その意味では、これから書くことは前のパラグラフの繰り返しである。以下では、欠損データをもっと狭く、データ行列が不完全である状況(ないし一部の変数について部分的情報しか使えない状況)を指すものとしよう。
ベイジアン・パラダイムは概念的にストレートである。なぜなら尤度は完全な矩形データセットを必要としないからである。簡単にいえば、ベイズ統計学では未知の量についてデータのもとでの予測分布を生成する。欠損データにあてはめれば、観察データのもとでの欠損データの予測分布を求める。多重代入は単にこの予測分布からのドローでアリ、その予測分布のモデルが良いモデルなら、他の分析に使うことができる。
4-5章では、欠損データのベイズ推論と多重代入についての基礎理論を概観した(特にMAR仮定の役割を強調した)。この理論の拡張のうち重要なもののひとつとして、coarseされたデータ(一部の変数がroundされたりグループ化されたり打ち切られてりしているデータ)という問題がある。この理論は無情報打ち切りという概念と関係がある(この概念は、打ち切りデータの生存分析の多くの手法の基盤となっている)。多重代入はこうした場面にも適用できるし、主たる分析に含まれていない補足変数で条件付けた補完を行う際に特に有用である。
6章では、欠損データへのベイジアンのアプローチについて概観した。主に正規モデルについて概観したが、非正規な欠損データ問題に対してもベイジアン手法はやはり有用である。SRMIは正規モデルに限らずに用いることができるし、ベイジアンアプローチや多重代入はカテゴリカルデータモデルやカテゴリカル変数と連続変数が混合したデータに対しても適用できるし、共変量に欠損がある一般化線型モデルにも適用できる。階層ベイズモデルは縦断データや小地域推定に対しても自然である。
本論文の事例ではMARモデルに焦点を当てたが、NMARのモデルに対してもベイジアン・アプローチは魅力的である。NMARデータの重要な問題は、モデルの識別性がないということである。ベイジアン手法は、識別されていないパラメータについて、その不確実性を反映した事前分布を形成することで、この問題にフォーマルな解を与えてくれる。結果として得られる推論は、感度分析からみて、頻度主義的手法よりもあきらかに優れている。頻度主義的手法では識別されていないパラメータはおおきくばらつく。
本論文では、私はCBのうちベイズの部分に焦点を当ててきた。これはこの論文のもとになったワークショップがベイズ手法に重点を置いていたからである。[“Bayesian Methods that Frequentists Should Know”というワークショップでのトークが基になっているらしい]
カリブレートのパートについていえば、よい頻度主義的特性を得るためには欠損値の予測分布についての現実的なモデルが必要となる。そのため、モデルのデータへのあてはまりのチェック、そしてMARからの逸脱のインパクトを評価するための感度分析が必要になる。多重代入の際のモデルチェック手法についてはGelman et al.(2005 Biometrics), Abayomi, Gelman, Levy(2008 Appl.Statist.)をみよ。この分野についてのさらなる研究が求められる。
コメント: Larsen
CBアプローチというアイデアを詳述して下さったのはありがたいんですが、ベイジアンの文脈でカリブレーションを実装する方法がどう拡張されたのかを知りたいです。カリブレーションというのは、変数を選択・変換してモデルのデータへのあてはまりを良くすることなんですか？ 欠損データとモデル/変数についての仮定の感度をチェックするためにもっと分析しようということなんですか？ 階層モデルについては？
計算パワーとアルゴリズム・ソフトウェアの進歩のおかげで、ベイジアンアプローチはより多くの人に使えるようになりました。今度は「どのようにカリブレートするか」のガイダンスが必要です。
私が思いつく例を挙げさせて下さい。頻度主義的な分析を報告する手はずなのだが欠損がある、というとき、結果についての信頼性を確かめるために多重代入が有用です。なぜなら、多重代入は欠損についてのさまざまな仮定のもとで行うことができ、その結果を比較できるからです。[事例。中略]
ベイジアンの分析は詳細には報告しないとしても、多重代入手続きは頻度主義的な手続きによる結果に保証を与えてくれます。こういうチェック方法が標準的になれば、統計的実践はよりCBへと近づくでしょう。
上の例で、もし多重代入による分析結果が完全ケースの分析結果と異なっていたら、その理由を知ることが必要になります。そのプロセスでなにか重要なことがみつかるかもしれません。
Littleの返答: 実践においてベイジアン・カリブレーションに到達する方法について、もっと多くの例を示せ、ということですね。仰るとおりです。私の例は上っ面を撫でた程度に過ぎません。
階層ベイズモデルはとても有用ですね。モデルチェックも重要ですが、残念ながら良い予測を保証してくれる訳ではありません。CBというパースペクティブによって、より多くの研究が生まれ、良いモデルの構築方法についてのトレーニングが進むことを願っています。
以下、3点コメントします。
1点目。ベイジアン分析と多重代入の実装は簡単になりましたが、やはり計算の詳細を知っている応用統計家は必要です。たとえば、欠損データからの最尤推定値のSEを得る方法は？ 現在のツールキットで簡単に手に入るでしょうか。さらに、計算時間は大きな壁となるので、効率的な計算・アルゴリズムは依然として重要です。
2点目。Little先生はSRMI(MICE)とPSPPについて紹介し、後者は二重頑健特性を持つと述べています。サーベイ抽出では「ホットデック法」と呼ばれるドナーベースの手法が多く使われています。良いホットデック手法では、類似したドナーをピックアップするために、多変量マッチングやプロペンシティ・マッチングと似たやり方でマッチングを行います。ドナーは現実の、データセットにおける真の連関パターンと一致する値を持っていますが、変数間には依存性があり、そのモデル化が課題となります[ここ、訳に自信がない…]。
よくデザインされた多重代入ホットデックアプローチと、ホットデック法とモデリングの併用によって、多重代入に変わるもう一つのアプローチが得られるでしょう。それはベイジアンにも頻度主義者にも受け入れられるものとなるはずです。
3点目。これは長いです。
Little先生は2章で、頻度主義統計学とベイジアン統計学を分けておられますね。でも、もっと広い見方もできます。自分のことをサーベイ・サンプラーだと思う統計家がいるのだ、という見方です(先生も私もこの世界と関わりを持っています)。サーベイ・サンプラーは、有限母集団の値について推論するために、Cochranの教科書に書いてあるような手続きを踏みます[基本的に確率抽出を考える、ということだろうか]。この手続きにおいて、ランダム性は有限母集団からのランダム選択から生じます。これは頻度主義的推論と関連していますが、「パラメータ」は異なります。たとえば、パラメータは\(\mu\)なんかではなくて\(\bar{y} = \sum_i^N y_i / N\)です。推論の目標は本質的に頻度主義です。「現在の抽出枠から得た確率標本に基づく95%信頼区間は、目標母集団の量を少なくとも95%含んでいなければならない」と考えるわけです[95% confidence interals … shold contain their target population quantity at least 95% of the time…]。この目標が含意しているのは、概念的に類似した抽出枠からの標本であれば、やはり95%のカバレッジとなるべきだ、ということです。
大規模サーベイでは回答のコーディング・エディティング、事後層別、ウェイト調整といったステップがとられることが多いです。これらはモデルパラメータ\(\theta\)を推定するという頻度主義の原理から明確に動機づけられたものではありません[←???]。たとえば事後層別は、分散減少と無回答バイアスの減少を目指しています。ウェイト調整は既知の母集団合計などに合わせるために行われることもあります。
ランダム化による推論の良い点はモデルがいらないという点です。もちろん、すべてのサーベイ推論手続きの良し悪しが文脈から独立に決まるわけではありません。たとえば、サーベイ抽出では比推定量が標準的選択ですが、それがうまくいくのはアウトカムと補足変数の間に強い正の相関があるときだけであり、このアプローチは切片ゼロの線型モデルという制約的なモデルと整合しています。良い推定量を選ぶためにはモデルの限界を知る必要があります。
Little先生のいう意味で「カリブレートされている」ためには、現在のモデルとデータとの適合のチェック、そしてモデリングのより広い選択肢の検討が必要になります。たとえば、Sarndal et al.(1992, 書籍)が述べているカリブレーション推定のようなアイデアです。[…中略…]
サーベイ・リサーチャーがベイジアン手法を受け入れるためには、まずベイジアン手法一般と有限母集団サーベイ抽出との間の関係についてのさらなる精緻化と発展が必要でしょう。ただし、欠損データの多重代入はそれほど高いハードルがありません。代入されたデータに対していつものウェイトとサーベイ推定量を用いればいいわけです。鍵となる問題は、代入が適切かどうかをどうやって判断するか、です。
本当の分断は、データをモデル化する人とモデル化しない人の間の分断です。人はこう問うことができます。そのパラメトリック・モデルの選択は正しいか？ その尤度関数は？ その事前分布は？ ここで、データとの整合性に関心を持つか、モデルを通じたデータからの情報抽出という目標に関心を持つか、そこに真の分断があるのです。多重代入のための柔軟なモデリングは、整合性のほうに関心を持っています。[…後略…]
[だらだらメモしてきたけど、話の要点がつかめない。要するに、Little先生や私は、頻度主義者でもベイジアンでもない、サーベイ・サンプラーでしょう？ 推論の対象は有限母集団特性であり、ランダム性が抽出から生じていると考えている点で本質的に頻度主義的でしょう？ というような主旨かと思ったのだが、最後に出てくる「データへのフィット vs 情報抽出」という話とのつながりがみえない。ううう… きっと読解力が足りないんだろうな…]
Littleの返答 「デザイン・ベース」のサーベイ推論の目標は「本質的に頻度主義的」ですけれど、私にとっては、ベイジアン推論は、モデルパラメータについての推論と同様に有限母集団の量についての推論においても有用であり適切です。
「デザインベース推論にはモデルがない」というのは言い過ぎです。Larsenさんが指摘されているように、デザインベース推論でも、推定量の選択の基板にある暗黙的なモデルについて考慮しないと、Basuの有名な象の例のようになってしまいます。
「カリブレーション」という言葉について、私の使い方と他の人の使い方を区別する必要があると気づきました。まず、推定値を累積統計量に「カリブレート」するのは、ベイジアンモデルではその情報を組み込んだ予測によって(原理的には)自動的に達成されます(ぴったり合わせるのは難しいかも知れませんが)。次に、Larsenさんは「モデリングのより広い選択肢」としてSärndal, et al.(1992)の手法を挙げていますが、これはモデルの残差に基づきモデルの予測をデザイン・ウェイト調整に「カリブレート」する手法で、モデル予測と直接的なデザイン・ベース推定値の折衷案です。こういう手法は突き詰めていえばデザイン・ベースであり、CBアプローチよりも劣っています。PSPPの文脈でそのことを示した研究として、Kang & Schafer(2007 Statist.Sci.), Zhang & Little(2011 J.Statist.Comput.Simul.)を参照して下さい。[どちらも二重頑健法の研究]
コメント: Schenker
1. プラグマティストのことを忘れるな
わたくしは自分が「プラグマティスト」であると捉えたい。異なる哲学の間に橋を架けることなく、手元の問題に役立つ部分をつまみ食いする人である。さらにいえば、私が一緒に仕事している実質的問題の専門家たちは、点推定値とかSEとか信頼区間とかになじんでいるが、事後平均や信用区間を使うことになんら抵抗がない。
かつてLittle先生は、プロとしての信頼を得るためにはプラグマティスト的な「二重人格」を持つべきでないと仰った。先生のCBのご提案は統一的人格の好例であろう。しかし、頻度主義の方法はすでに普及しており、時乳人格から抜け出すのは容易でない。
Littleの返答: 告白しますが、応用研究では私もプラグマティックです。信頼区間もP値も出しちゃいます。でも、そういう日々の現実の話とは別に、我々の方法論を導く原理について考えることも大事だと思います。たとえばサーベイ抽出の文脈で、ある種の問題についてはデザイン・ベース推論は不適切であり、モデル・ベース推論のほうが適切だ、と私は信じています。
2. 頻度主義 vs. ベイジアンの分裂は、おそらくサーベイ抽出において拡大している
私はサーベイ抽出の世界で働いており、モデルパラメータについてというより有限母集団の量について推論したいほうなので、ベイジアンとはなにか、頻度主義者とはなにかがいまいちよくわからない。
有限母集団についての推論ではデザイン・ベースのパラダイム(所与のデザインのもとでの有限母集団からの反復抽出における推定量の分布に基づくパラダイム)が用いられることがいい。頻度主義的推論を定義するとしたら、有限簿手段の値\(Y\)を固定パラメータと捉え、それらのパラメータのある関数\(Q(Y)\)について、抽出された値の関数と反復抽出におけるその分布に基づいて推論すること、とすればいいのかもしれない。またベイジアン推論を定義するとしたら、\(Y\)に対して事前分布\(p(Y|\theta)\)を与え(\(\theta\)はハイパーパラメータで事前分布\(p(\theta)\)を持つ)、抽出された値のもとでの\(Q(Y)\)の事後予測分布に基づいて推論すること、とすればいいのかもしれない。
サーベイ抽出の世界での頻度主義とベイジアンの違いはなんだろうか。量を確率変数とみなすかという点が違うし、事前分布を指定するかという点も違う。でも、最も大事な違いは、頻度主義ないしデザイン・ベースのアプローチでは、推論のための参照分布が有限母集団の値\(Y\)についてのモデルから得られるわけではないのに対し、ベイジアンの事後予測分布はそういうモデルを含んでいるという点である。[??? どういうこと？ よくわからない。デザイン・ベース推論でウェイトは\(Y\)なしで決まるってことをいっているのだろうか]
デザインベース推論とモデルベース推論の対立については、Hansen et al.(1983), Little(2004)を参照してほしい。[…中略…]
Little(2004)はサーベイ抽出の文脈でもベイジアン・パラダイムが有用だと結論しているが、その際、ベイジアンの推論に用いるモデルも標本デザインの特性(ウェイティングとか)を適切に反映していなければならないと指摘している。しかし後述するように、モデル・ベース推論が有用であるような問題の中には、標本デザインの特徴を反映させると複雑になってしまう問題もある。デザイン特性をどう反映させるかが今後の重要な研究領域である。
3. プラグマティストはなぜベイジアンの手法を好むか
ベイジアンの手法が非常に有用であったいくつかのプロジェクトを紹介しよう。[すべてスキップ。見出しのみメモする]
- ときどき観察される共変量をもつ生存分析
- 多重代入を通じて生存分析に補足変数を統合する
- サーベイにおける欠損データの多重代入
- ふたつのサーヴェイから得た情報を結合して小地域推定を拡張する
4. 今後の研究課題
- 柔軟やモデルと手法。SRMIとかPSPPみたいな手法がもっとたくさんほしい。
- モデルの診断。特に欠損データの予測モデルのチェック方法。
- 複雑な標本デザイン特性をモデルに組み込む方法。たとえば、多重代入の文脈では、PSUと層についてのサーベイ・ウェイトとインジケータ変数の導入が推奨されている。しかし倹約のため、PSU選択についてのすべてのインジケータ変数ではなく、その一部だけを含めるということがある。倹約性を高める方法の研究が必要だ。ランダム効果を使うとか。
また、複数のサーヴェイの情報を結合する際には、複雑な標本デザイン特性を組み込むのがさらに難しくなる(調査によって標本デザイン特性がちがうから)。[よくわからんが、複数のサーベイの個票を縦に積むとき、標本デザインが違ってたらどうするんだ、という話だろうか] - 補完モデルに含まれていない変数を使った二次分析のインパクト。欠損値が補完されたデータを二次分析する際、補完モデルに含まれていなかった変数を一緒に使ったらなにが起きるのだろうか。そのとき生じるバイアスは、補完モデルに含まれていた変数が、二次分析で研究されている関係をどの程度よく説明しているかによって決まるだろう。
- CBアプローチの有用性を示す現実場面での事例をもっと下さい。
————————-
いや、勉強になりましたよ？ 欠損データ分析って不勉強でしたし。とても勉強になりましたですが、しかし、読み終えて非常にもやもやしている。
冒頭でCalibrated Bayesアプローチという風呂敷を拡げたきり、そのあとでずーっと欠損データ分析へのベイジアン・アプローチの紹介になってしまい、この話のどこがCalibrated Bayesなのか、よくわからないまま終わってしまった。落ちはどこだ、どこにあったのだ。それとも私の力不足で読み取れていないのか。不安に駆られて最終章をほぼ逐語訳してしまったが、やっぱりわからんかった。
Larsenさんがコメントしているように、Littleさんのいうところの「カリブレート」の例をもっと挙げてほしい… 実例がないとわかんないよ… (泣く)