Little, R.J.A., West, B.T., Boonstra, P.S., & Hu, J. (2020) Measures of the Degree of Departure from Ignorable Sample Selection. Journal of Survey Statistics and Methodology, 8(5), 932–964.

　非確率標本からの母集団特性推定についてあれこれ調べていて、たまたま見つけた論文。たいていの手法は、母集団の個体が標本に包含されないということを欠損と捉えたときにその欠損が調査変数にとってMARであることを想定してバイアスを取り除こうとするのだが(傾向スコアであろうがMRPであろうがそうです)、この論文は珍しく、MNARであることを前提にして、バイアスをどう取り除いたらいいのかわからんけれどどのくらいの大きさのバイアスがありそうかを推定します、という話である。
　なんだかややこしそうな話ではあるが、第一著者はLittle先生。話は難しいけれど書き方はわかりやすい先生だ、というぼんやりした信頼感がある。

1. イントロダクション
　元の標本の一部しか反応しておらず、反応メカニズムが無視可能でないかもしれないとき、その標本をどの程度まで確率標本として扱ってよいのだろうか？ […]
　Rubinいうところの無視不能性とは、観察データで条件づけても欠損確率が欠損データに依存しているということであった。この定義は標本選択にも適応できる。[…]
　この論文では、無視不能な選択バイアスの大きさを示す指標を提案する。こはAndridge & Little(2009, 2011)、West & Little(2013), West, Wagner, Gu, Hubbard(2015)に基づいていて…[簡単な説明がなされているが、ここだけ読んでもわかりにくいのでメモ省略]。

　先行する提案としてR指標がある。調査への回答確率の変動性を、標本全体における補助共変量の関数として測る。値が小さいことは最終的な回答者集合においてバランスがとれていることを示す。Sarndalたちはこの別バージョンを提案していて[…]。いずれにせよ、選択のモデルが必要になるし、選択バイアスというのは調査変数と選択との関係の強さに依存するという事実を反映できていない。また、利用可能な補助変数の値を通じた反応に依存しており、無視不能選択を反映していない。SarndalらのH_1指標は[…これもまた問題がある]。というわけで、本研究ではこれらの指標については検討していない。
　Nishimura, Wagner & Elliott(2016)のシミュレーション実験によれば、R指標は欠損メカニズムが無視不能な時の無回答バイアスの効果的な指標ではない。モデルベース多重代入の副作用としてのアウトカム特定的な指標である「欠損情報の推定割合」(FMI)は、所与の測定値と結びついた無回答率よりも大きい。このことは潜在的な無回答バイアスを示している。これらの結果は、FMIは付加的な検討事項としては有用かもしれないが、潜在的選択バイアスの指標が依然として必要とされていることを示している。[…背景知識がないのでなにがなんだかよくわかんない…]

2. Rubinの欠損データフレームワークの標本選択への適用
[以下、勝手に小見出しをつける]
　RubinのフレームワークはRubin(1976)で、その標本選択への適用はRubin(1978), Rubin(1987), Little(2003)で示された。反応インジケータを標本選択インジケータに差し替えたわけである。

[セッティング]
　母集団の個体\(i (=1,\ldots,N)\)の調査変数の値を\(y_i\)とし\(\mathbf{Y} = (y_1, \ldots, y_N)\)とする。観察されている補助変数ないしデザイン変数を\(\mathbf{Z} = (Z_1, \ldots, Z_n)\)とする。有限母集団特性を\(Q = Q(Y, Z)\)とする[えっ、\(Z\)の関数なの？ どういう状況を考えているんだろうか]。標本包含インジケータのベクトルを\(\mathbf{S} = (S_1, \ldots, S_N)\)とする。\(\mathbf{Y}\)を標本包含されている部分\(\mathbf{Y}_{inc}\), されてない部分\(\mathbf{Y}_{exc}\)にわける。

[調査変数と標本包含の同時モデル]
　\(Y, S\)の同時モデルを考える。$$ f_{Y,S}(\mathbf{Y}, \mathbf{S} | \mathbf{Z}, \theta, \phi) = f_Y(\mathbf{Y} | \mathbf{Z}, \theta) f_{S|Y} (\mathbf{S} | \mathbf{Y}, \mathbf{Z}, \phi) $$ 全尤度は こう書ける。$$ L(\theta, \phi | \mathbf{Z}, \mathbf{Y}_{inc}, \mathbf{S}) \propto f_{Y,S}(\mathbf{Y}_{inc}, \mathbf{S} | \mathbf{Z}, \theta, \phi) = \int f_Y(\mathbf{Y} | \mathbf{Z}, \theta) f_{S|Y} (\mathbf{S} | \mathbf{Y}, \mathbf{Z}, \phi) d \mathbf{Y}_{exc} $$ [落ち着け落ち着け。最初の比例関係は尤度の定義ですね。で、$$ f_{Y,S}(\mathbf{Y}_{inc}, \mathbf{S} | \mathbf{Z}, \theta, \phi) = \int f_{Y,S}(\mathbf{Y}, \mathbf{S} | \mathbf{Z}, \theta, \phi) d \mathbf{Y}_{exc} $$ だよね。で、$$ f_{Y,S}(\mathbf{Y}, \mathbf{S} | \mathbf{Z}, \theta, \phi) = f_{S|Y} (\mathbf{S} | \mathbf{Y}, \mathbf{Z}, \phi) f_Y(\mathbf{Y} | \mathbf{Z}, \theta) $$ ですよね。第1項は選択モデル、第2項はアウトカムモデルですね]

[パラメータの事後分布]
　\(\theta, \phi\)の事前分布と事後分布は$$ p(\theta, \phi | \mathbf{Z}, \mathbf{S}, \mathbf{Y}_{inc}) \propto p(\theta, \phi | \mathbf{Z}) L(\theta | \mathbf{Z}, \mathbf{S}, \mathbf{Y}_{inc}) $$ ここで\(p(\theta, \phi|\mathbf{Z})\)はパラメータの事前分布である。[んんん？ 右辺は\(p(\theta, \phi) L(\theta, \phi | \mathbf{Z}, \mathbf{S}, \mathbf{Y}_{inc})\)じゃないの？ 仮に\(p(\theta, \phi|\mathbf{Z})\)を事前分布として定義できるとしても(どういう状況だろう？)、その場合は\(p(\theta, \phi | \mathbf{Z}) L(\theta, \phi | \mathbf{Z}, \mathbf{S}, \mathbf{Y}_{inc})\)じゃないだろうか。よくわからない…]

[非標本の目的変数の事後分布]
　\(\mathbf{Y}_{exc}\)の事後分布は $$ p(\mathbf{Y}_{exc} | \mathbf{Z}, \mathbf{S}, \mathbf{Y}_{inc}) \propto \int p(\mathbf{Y}_{exc} | \mathbf{Z}, \mathbf{S}, \mathbf{Y}_{inc}, \theta, \phi) p(\theta, \phi | \mathbf{Z}, \mathbf{S}, \mathbf{Y}_{inc}) d\theta d\phi $$ [これはたいしたこと云ってない。パラメータを消しているだけだ] 多くのモデルにおいては\(p(\mathbf{Y}_{exc} | \mathbf{Z}, \mathbf{S}, \mathbf{Y}_{inc}, \theta, \phi) = p(\mathbf{Y}_{exc} | \mathbf{Z}, \mathbf{S}, \theta, \phi)\)だから、\(\mathbf{Y}_{exc}\)の事後分布はパラメータを通じてしか\(\mathbf{S}, \mathbf{Y}_{inc}\)に依存しない。

[もし標本選択が無視可能だったら]
　\(S\)のモデルの指定は難しい。選択メカニズムを無視したパラメータの尤度について考えてみよう。$$ L_{ign}(\theta|\mathbf{Y}_{inc}, \mathbf{Z}) \propto p_Y(\mathbf{Y}_{inc}|\mathbf{Z}, \theta) = \int p_Y(\mathbf{Y}|\mathbf{Z}, \theta) d \mathbf{Y}_{exc}$$ 無視可能だから、\(S\)のモデルは出てこない。パラメータと\(\mathbf{Y}_{exc}\)の事後分布は $$ p(\theta | \mathbf{Y}_{inc}, \mathbf{Z}) \propto p(\theta | \mathbf{Z}) L_{ign}(\theta | \mathbf{Y}_{inc}, \mathbf{Z}) $$ $$ p(\mathbf{Y}_{exc} | \mathbf{Y}_{inc}, \mathbf{Z}) \propto \int p(\mathbf{Y}_{exc} | \mathbf{Y}_{inc}, \mathbf{Z}, \theta) p(\theta | \mathbf{Y}_{inc}, \mathbf{Z}) d\theta $$ このように事後分布が単純化される場合、選択メカニズムは無視可能であるという。

[無視可能性の十分条件]
　無視可能性の一般的で単純な十分条件として以下がある。

• ランダム選択(SAR): すべての\(\mathbf{Y}_{exc}\)について\( f_{S|Y}(\mathbf{S}|\mathbf{Y}, \mathbf{Z}, \phi) = f_{S|Y}(\mathbf{S}|\mathbf{Y}_{inc}, \mathbf{Z}, \phi) \)
• ベイジアンdistinctness: \(p(\theta, \phi | \mathbf{Z}) = p(\theta|\mathbf{Z}) p(\phi|\mathbf{Z}) \)

ふたつめはふつう合理的ですね、標本選択を制御するパラメータと\(Y\)のモデルのパラメータはふつう無関係と想定されるから。
　これらの仮定の下で $$ p(\theta, \mathbf{Y}_{exc} | \mathbf{Y}_{inc}, \mathbf{Z}) = p(\theta, \mathbf{Y}_{exc} | \mathbf{Y}_{inc}, \mathbf{Z}, \mathbf{S}) $$ つまり、有限母集団特性についての推論にデータ収集メカニズムがが影響しない。

　なお、確率標本の場合、すべての\(\mathbf{Y}_{exc}\)について \(f_{S|Y}(\mathbf{S}|\mathbf{Y}, \mathbf{Z}, \phi) = f_{S|Y}(\mathbf{S} | \mathbf{Z}) \)である。SARよりはるかに強い。

3. 連続変数の平均に対する選択バイアスの指標
　\(Y\)が連続変数だとして、その平均における選択バイアスを表す指標を考えよう。

[最良予測子]
　\(Y\)を\(Z\)に回帰したときの最良予測子を\(\mathbf{X}\)としよう[\(\mathbf{Z}\beta\)のことであろう]。で、母集団レベルでの\(X\)の漸近不偏な要約指標が手に入っているとしよ。非確率標本における\(Y\)と\(X\)の分散をそれぞれ\(\sigma^{(1)}_{YY}, \sigma^{(1)}_{XX}\)とする[上添字は\(S=1\)であることを表す]。$$ X^* = X \sqrt{\sigma^{(1)}_{YY}/\sigma^{(1)}_{XX}} $$ とリスケールする (分散を\(Y\)の分散に一致させている)。これを\(Y\)の補助プロキシと呼ぶ。

[プロキシパターン混合モデル]
　これから提案する指標は、\(Y\)を\(X\)と関連付ける正規プロキシパターン混合モデル(PMM)の最尤推定値に基づく。
　次のように仮定する。\(j\)は0か1ね。\(N_2()\)は二変量正規分布。$$ (X, Y | S = j) = N_2 \left( (\mu^{(j)}_X, \mu^{(j)}_Y), \left( \begin{array}{cc} \sigma^{(j)}_{XX} & \sigma^{(j)}_{XY} \\ \sigma^{(j)}_{XY} & \sigma^{(j)}_{YY} \end{array} \right) \right) $$ $$ Pr(S=1 | X, Y) = g(V), \ V=(1-\phi) X^* + \phi Y $$ [しばらく式を眺めていてようやくわかってきた。目的変数\(Y\)と、その最良予測子\(X\)の同時分布に二変量正規性を仮定する(これ、結構強気な仮定ですね)。そのうえで、標本包含確率を、\(Y\)の関数、ないし\(X\)の関数、ないしその中間とする。前者は最悪のケース、後者は(\(X\)で条件づければ標本選択メカニズムを無視できるわけだから)最良のケースである。わざわざ\(X\)をスケーリングしたのは、前者と後者のバランスを表すパラメータ\(\phi\)を0以上1以下にしたいからだ]

　もっと一般的な形式で書いてもよい。利用可能な補足変数\(Z\)を最良予測子\(X\)と(非確率標本において)それと直交な変数\(U\)に分ける。\(Y\)の平均は変換後の\(U\)には依存しない。証明はAppendixに回すが、非標本部分についても\(X\)が\(Y\)の最良予測子であり\(U\)が\(X\)と直交していれば、上のPMMの最尤推定量なりベイズ推定量なりは、$$ Pr(S=1 | X, Y, U) = g(U, V), \ V=(1-\phi) X^* + \phi Y $$ というより一般的なメカニズムにおいても妥当である。
　[このくだり、理解できない。標本における\(Y\)の\(Z\)への回帰モデルが非標本においても成り立っているならば、ということですよね。それって、\(Z\)のなかに選択メカニズムを条件づけている共変量が入っているから成り立つわけですよね。それらの共変量は\(X\)に含まれますよね。だったら、標本包含は\(X\)の下で\(Y\)と独立でしょう？ そのときは\(Pr(S=1 | X, Y, U]) = g(X)\)じゃありませんか？]

[母平均の最尤推定量]
さて、\(\phi\)のもとでの母平均の最尤推定量は、\(g(\cdot)\)がどうであれこうなる。詳しくはAndridge & Little(2011) をみよ。$$ \hat{\mu}_Y(\phi) = \bar{y}^{(1)} + \frac{\phi + (1-\phi) r_{XY}^{(1)}}{\phi r_{XY}^{(1)} + (1-\phi)} \sqrt{ \frac{ s^{(1)}_{YY} }{ s^{(1)}_{XX} }} (\bar{X}-\bar{x}^{(1)}) $$

[提案指標SMUB]
というわけで、ご提案します！ 標本平均の「未調整バイアス指標」です！ $$ MUB(\phi) = \bar{y}^{(1)} – \hat{\mu}_Y(\phi) = \frac{\phi + (1-\phi) r_{XY}^{(1)}}{\phi r_{XY}^{(1)} + (1-\phi)} \sqrt{ \frac{ s^{(1)}_{YY} }{ s^{(1)}_{XX} }} (\bar{X}-\bar{x}^{(1)}) $$ いやしかし、これだと\(Y\)のスケールに依存するので使いにくい。よし、\(Y\)の標本SDで割ろう。
　ご提案します！ 標本平均の「未調整バイアスの標準化指標」です！ $$ SMUB(\phi) = \frac{\phi + (1-\phi) r_{XY}^{(1)}}{\phi r_{XY}^{(1)} + (1-\phi)} \frac{\bar{x}^{(1)} – \bar{X}}{\sqrt{ s^{(1)}_{XX}}} $$ \(\phi = 0\)ならSARである。\(\phi = 1\)ならば、標本選択の\(X, Y\)への依存性は\(Y\)のみを通じてなされる。わからない場合は\(\phi = 0.5\)とすることをお勧めします。結局こうなる。$$ SMUB(0) = r^{(1)}_{XY} \frac{(\bar{x}^{(1)} – \bar{X})}{ \sqrt{s^{(1)}_{XX} }} $$ $$ SMUB(0.5) = \frac{(\bar{x}^{(1)} – \bar{X})}{ \sqrt{s^{(1)}_{XX} }} $$ $$ SMUB(1) = \frac{1}{r^{(1)}_{XY}} \frac{(\bar{x}^{(1)} – \bar{X})}{ \sqrt{s^{(1)}_{XX} }} $$ [うわー。めっさシンプルやんか] なお、Rの関数を配っとるから使ってね。

　提案指標についていくつか注釈したい。

• \(X\)の母平均がわかっていればよいという点に注意。ミクロデータはいらない。
• SMUBはAndridge & Little(2011)で提案した感度分析に対応している。
• 平均の調整済みバイアスの標準化指標 \(SMAB(\phi) = SMUB(\phi) – SMUB(0) \)を使うこともできる。ただし、これはPMMモデルに強く依存している点に注意。
• プロキシ変数\(X\)が\(Y\)の良い予測子でないとき\(r^{(1)}_{XY}\)はゼロに近づきSMUB(1)は不安定になる。
• 標本割合が大きいと指標は小さくなる。
• SMUB(0.5)はBiemer & Peytchev(2011)のバイアス効果サイズと密接に関係している[…]
• \(SMUB(\phi)\)はSARを仮定していない。そのかわりにPMMに依存している。非正規変数に拡張できるけどややこしくなるので、非正規変数であってもこの指標を使うことをお勧めする。たとえ選択バイアスが大きくても\(X\)の標本平均と母平均が近かったらゼロに近くなる。\(X\)には層別抽出の層の変数は入れられない[比例層別抽出の場合はってことだよね？]
• もし\(S=1\)の標本が、母集団の確率標本のなかで反応した個体であり、台帳誤差と無回答のせいでそうなっているのならば、モデルにおける\(S=0\)の個体はより現実的に、無回答者という下位母集団、ないし台帳外の人々に限定するべきである。しかし、このモデルから得られるバイアスの最尤推定値は、少なくとも抽出デザインが等確率な場合は式(7)のモデルから得られる推定値に近いことを示すことができる。[なにいってんだかさっぱりわかんない。式(7)ってたぶん誤植だ]
• この指標に、抽出の不確実性を統合することもできる。[SMUB(0)とSMUB(1)の区間をベイズ推定するという話のようだ]

4. 適用例
[実調査データを母集団にして、スマホユーザに絞った時の真のバイアスとSMUBを比べている。パス]

5. 要約と今後の課題
　SMUBは無視不能な選択バイアスの大きさを変数別に表す指標である。補助変数の母集団での分布と標本での分布に注目している。プロキシ\(X\)と\(Y\)の相関が0.4以上くらいあればよい指標になる。[…]
　今後の課題: 二値の調査変数の場合への拡張(数値例では提案手法のままでもうまくいっているけれど)。ベイズ信用区間を出す場合の事前分布。母集団についての補助情報が\(X\)の母平均ではなくて確率標本だった場合にどうするか。
———-
　なるほどね。
　要するにこういうことだ。提案指標は、補助変数の母平均がわかっているというカリブレーション的な状況下で、非確率標本にどのくらいの選択バイアスがあるかを、補助変数で100%補正できる場合から全然補正できない場合までの幅という形で示す指標である。標本選択メカニズムについての明示的なモデルもないのにバイアスのサイズがわかるなんて奇妙な感じがするが、かわりに、補助変数と調査変数が標本と非標本でそれぞれ二変量正規分布するという、それはそれで強力なモデルに依存している。おかげで、補助変数と調査変数の相関が弱いときにはお手上げになる。