Chen, Y., Li, P., & Wu, C. (2019) Doubly Robust Inference With Nonprobability Survey Samples. Journal of the American Statistical Association.
相変わらず非確率標本からの母集団特性推定について調べているのだけれど、このたびは、先日読んだ Wu(2022)でフィーチャーされていた、Wuさんチームの論文を読んでみた。筆頭著者は博士課程の院生さん。
1. イントロダクション
[…]
非確率調査の標本と、参照確率標本から得られる補助情報とを用いて行う統計的推測についての一般的フレームワークを提出する。
2. 傾向スコアと標本マッチング
2.1 基本的セッティングと傾向スコア
母集団を\(\mathcal{U} = \{1, 2, \ldots, N\}\)、個体\(i\)の補助変数の値を\(\mathbf{x}_i\)、反応変数の値を\(y_i\)とする。デザインベース・フレームワークに立ち、\(\mathcal{F}_N = \{ (\mathbf{x}_i, y_i), i \in \mathcal{U} \}\)を定数とみる。\(\mu_y = \frac{1}{N} \sum_{i=1}^N y_i \)とする。
サイズ\(n_A\)の非確率標本\(\mathcal{S}_A\)があって\( \mathbf{x}_i, y_i\)がわかっている。包含インジケータを\(R_i\)とする。傾向スコアを$$ \pi^A_i = E_q(R_i | \mathbf{x}_i, y_i) = P_q(R_i = 1 | \mathbf{x}_i, y_i)$$ とする。添え字の\(q\)は傾向スコアモデル(非確率標本選択のモデル)を指す。
以下のように仮定する。
- A1. \(R_i\)と\(y_i\)は\(\mathbf{x}_i\)の下で独立。
- A2. すべての\(i\)において\(\pi^A_i gt 0\)。
- A3. \(\mathbf{x}_i, \mathbf{x}_j\) (\(i \neq j\)) の下で\(R_i\)と\(R_j\)は独立。
A1とA2を合わせて強い無視可能性条件という。もっとも”ignorable”という表現はよろしくない。分析者が無視していいわけじゃないから[なるほど？]。ちなみにMARっていうのもよくない、randomly missingではないから[それはほんとにそうだよね]。
参照確率標本\(\mathcal{S}_B\)があって\(\mathbf{x}_i\)が利用できるとしますね。もちろん\(y_i\)はわかんない。\(\mathcal{S}_A\)と\(\mathcal{S}_B\)とは独立とする。包含確率\(\pi^B_i\)、その逆数であるウェイト\(d^B_i\)は既知とする。
Lee(2006 J.OfficialStat.), Lee & Valliant(2009 Soc.Methods&Res.)は、\(\mathcal{S}_A, \mathcal{S}_B\)をプールして、\(i \in \mathcal{S}_A\)のときに1, \(i \in \mathcal{S}_B\)のときに0になるインジケータ\(\tilde{R}_i\)をつくり、\(\tilde{R}^A_i = P(\tilde{R}_i = 1 | \mathbf{x}_i)\)を推定しようとした。しかしこれは\(\pi^A_i\)の正しい推定ではない。そのことはValliant & Dever(2011 Soc.Methods&Res.)も指摘している。
[今回調べていて、非確率標本の包含確率の推定をこの\(\tilde{R}^A_i\)でやっているのが結構あって、びびっていたのであった。普通に考えるとそうはならないと思う(それってケース・コントロール研究でオッズ比じゃなくてリスク比をみるようなもんじゃないですか？)。プロの研究者がそう考えたんだから、それはそれでなんか理屈があるのだろうかと考え込んでいたのだが、やっぱり単なる勘違いだったのかなあ…]
2.2 標本マッチングと予測アプローチ
すでにモデルベースの予測アプローチについての検討が進んでいる。母集団がモデル$$ y_i = m(\mathbf{x}_i) + \epsilon_i $$からの無作為標本で、\(m(\mathbf{x}_i) = E_\xi(y_i | \mathbf{x}_i)\)だとする (添字\(\xi\)は予測モデルのもとでの期待値であることを指す)。誤差項は\(\epsilon_i\)は期待値0, 分散\( v(\mathbf{x}_i) \sigma^2 \)に独立に従うとする。
A1の下で\(E_\xi(y_i | \mathbf{x}_i, R_i = 1) = E_\xi(y | \mathbf{x}_i) = m(\mathbf{x}_i)\)となる。これは非確率標本で推定できる。\(m(\mathbf{x}_i) = \mathbf{x}_i^\top \beta, v(\mathbf{x}_i) = 1\)ならば\(\beta\)はこうなって[式省略]、参照確率標本でも\(\hat{y}_i\)が得られるので\(\mu_{REG}\)が得られる[式は省略するけどHajek推定]。これはモデルと抽出デザインが正しければほぼ不偏である。[\(N\)が既知ならHT推定できてそれは厳密に不偏だけど、でもやっぱしこっちのほうがいいよという話。メモ省略]
回帰予測推定量はいわゆる「マス代入」の特殊ケースである。代入値の作り方としては最近隣法もあるよね。ノンパラな手法である。\(y_i\)をベクトルにすることもできる[あー、なるほどね！ モデルベースなら目的変数ごとにモデリングしなければならないと思い込んでいたが、そうでもないね！]
ともあれ、\(\mathcal{S}_B\)の全個体になんらかの代入値\(y^*_i\)を代入してHajek推定するのを\(\hat{\mu}_{SM}\)としよう(SMはsurvey matchingの略)。
3．二重頑健推定の手続き
3.1 傾向スコア推定
[この節は勝手に小見出しをつける]
[提案手法]
仮にですが、\(\mathcal{U}\)で\(\mathbf{x}_i\)がわかっているとしますね。で、\(\pi^A_i\)がパラメトリックに\(\pi^A_i = \pi(\mathbf{x}_i, \theta_0\)と書けるとします。
\(\theta_0\)の最尤推定量\(\theta\)を作りましょう。対数尤度関数はどうなるでしょうか？ はいそこの君、黒板に書きなさい。[そうは書いてないけどさ] $$ l(\theta) = \sum_{i=1}^N \left( R_i \log \pi^A_i + (1-R_i) \log(1-\pi^A_i) \right) $$ $$ = \sum_{i \in S_A} \log \left( \frac{\pi(\mathbf{x}_i, \theta)}{1-\pi(\mathbf{x}_i, \theta)} \right) + \sum_{i=1}^N \log (1-\pi(\mathbf{x}_i, \theta)) $$ はい、よくできました。[なぜ非確率標本での対数尤度と母集団での対数尤度に無理やりわけたのか、この次の段落でわかる]
現実には、\(\mathcal{U}\)では\(\mathbf{x}_i\)はわからず、\(\mathcal{S}_B\)でしかわかんないだろう。そこで、次の擬似対数尤度関数を最大化する。$$ l^*(\theta) = \sum_{i \in S_A} \log \left( \frac{\pi(\mathbf{x}_i, \theta)}{1-\pi(\mathbf{x}_i, \theta)} \right) + \sum_{i \in S_B} d^B_i \log (1-\pi(\mathbf{x}_i, \theta)) $$ つまり、母合計の項を\(S_B\)によるHT推定量に差し替えたわけである。[なるほどね]
\(\pi^A_i\)のパラメトリック・モデルがロジスティック回帰モデル$$ \pi(\mathbf{x}_i, \theta_0) = \frac{\exp(\mathbf{x}_i^\top \theta_0)}{1+\exp(\mathbf{x}_i^\top \theta_0)} $$ ならば、擬似対数尤度関数はこうなる: $$ l^*(\theta) = \sum_{i \in S_A} \mathbf{x}_i^\top \theta – \sum_{i \in S_B} d^B_i \log(1+\exp(\mathbf{x}_i^\top \theta)) $$ これを最大化すればいいわけだ。
というわけで、$$ U(\theta) = \frac{\partial}{\partial \theta} l^*(\theta) = \sum_{i \in S_A} \mathbf{x}_i – \sum_{i \in S_B} d^B_i \pi(\mathbf{x}_i, \theta) \mathbf{x}_i $$ として、スコア方程式\(U(\theta) = 0\)を解けばよい。Newton-Raphson法で解ける。[…]
[Valliant & Dever批判]
Valliant & Dever(2011) は \(\mathcal{S}_{AB} = \mathbf{S}_A \cup \mathbf{S}_B\)を使った重み付きロジスティック回帰を考えた。重みとして、\(i \in S_A\)には\(d_i = 1\), \(i \in S_B\)には\(d_i = d^B_i \frac{1-n_A}{\hat{N}_B}\)を振る。
ところが、この方法の対数尤度関数は $$ l^*(\theta) + \sum_{i \in S_A} \log(1-\pi_i^A) – \frac{n_A}{\hat{N}_B} \sum_{i \in S_B} \log(1-\pi^A_i) $$ となってしまう。この方法では一致推定量にならない(\(S_A\)がSRSなら別だけど)。
[推定方程式から見ると]
次のクラスの推定方程式を考える。$$ \sum_{i=1}^N R_i \mathbf{h}(\mathbf{x}_i, \theta) – \sum_{i=1}^N \pi_i(\theta) \mathbf{h}(\mathbf{x}_i, \theta) = 0$$ ここで\(\pi_i(\theta) = \pi(\mathbf{x}_i, \theta)\) である。この第2項に\(d^B_i\)を追加し、さらに\( \mathbf{h}(\mathbf{x}_i, \theta) = \mathbf{x}_i \) とすれば、$$ \sum_{i \in S_A} \mathbf{x}_i – \sum_{i \in S_B} d^B_i \pi(\mathbf{x}_i, \theta) \mathbf{x}_i = 0$$ となり、提案手法に一致する。
では、\( \mathbf{h}(\mathbf{x}_i, \theta) = \pi(\mathbf{x}_i, \theta)^{-1} \mathbf{x}_i \) とすればどうなるか？ $$ \sum_{i \in S_A} \frac{\mathbf{x}_i}{\pi_i(\theta)} = \sum_{i=1}^N \mathbf{x}_i $$ つまりカリブレーションになる。
Kim & Riddles(2012 SurveyMethodol.) いわく、線形回帰モデルを維持するならば、このクラスの推定方程式のなかで\( \mathbf{h}(\mathbf{x}_i, \theta) = \pi(\mathbf{x}_i, \theta)^{-1} \mathbf{x}_i \) とするのが最適な推定になる。Wu & Sitter(2001 JASA), Wu(2003 Biometrika)をみよ。
[うーん。この話、Wu(2002)でも読んだんだけど、いまいちピンとこない。推定方程式の理論についてよく理解できていないからだろう]
[データ結合問題として捉えると]
[略]
3.2 逆確率ウェイテッド推定量
ご存じHT推定は、ウェイトが抽出デザインで決まる場合であった。いっぽうIPW推定量は傾向スコアモデルに基づく。サーヴェイの文脈では擬似ランダム化アプローチともいう。
傾向スコアの推定値\(\hat{\pi}^A_i\)が得られたら、\(\mu\)のふたつの推定量をつくれる。$$ \hat{\mu}_{IPW1} = \frac{1}{N} \sum_{i \in A} \frac{y_i}{\hat{\pi}^A_i} $$ $$ \hat{\mu}_{IPW2} = \frac{1}{\sum_{i \in S_A} (\hat{\pi}^A_i)^{-1}} \sum_{i \in A} \frac{y_i}{\hat{\pi}^A_i} $$
有限母集団の系列\(\mathcal{U}_v\)があるとしよう。サイズを\(N_v\)とする。それぞれの有限母集団に対応する非確率標本\(S_{A,v}\)(サイズ\(n_{A,v}\))と確率標本\(S_{B,v}\)(サイズ\(n_{B,v}\))がある。\(v \rightarrow \infty\)と共に三つとも\(\infty\)に近づくとする。以下では\(v\)のことは書かず、\(N \rightarrow \infty\)と略記する。
定理1. 仮定A1-A3, Supplimentary Materialに書いた6個の正則性条件、傾向スコアのロジスティック回帰モデル、を仮定すると、以下がいえる。[ふたつの推定量の誤差と分散が\(O\)-記法(っていうんですか)で示されている。メモ省略。こういう頭のいい人だけがかかるひどい水虫とかがあるといいな]
3.3 二重頑健推定量
IPW推定量はモデルの誤指定に弱い。Tan(2007)をみよ。
[…] 欠損値データ分析の最近の研究でさかんに検討されている二重頑健推定量を導入しよう。
パラメトリックモデル \( E_\xi(y | \mathbf{x}) = m(\mathbf{x}, \beta) \)について考える(添え字\(\xi\)はアウトカム回帰のモデルを指す)。二重頑健(DR)推定量は$$ \hat{\mu}_{DR} = \frac{1}{N} \sum_{i=1}^N \frac{R_i(y_i – m(\mathbf{x}_i, \beta))}{\pi(\mathbf{x}_i, \theta)} + \frac{1}{N} \sum_{i=1}^N m(\mathbf{x}_i, \beta) $$ 仮定A1の下で、\(\beta\)の推定量は\(S_A\)に基づく最小二乗推定量なり最尤推定量なりでよい。
Wu & Sitter(2001)は、母集団の\(\mathbf{x}_i\)がわかっている状況での「一般化差分推定量」というモデル・アシステッド推定量を提案していたんだけど、上の推定量はそれと同じである。
では、現在のセッティングにあわせると、\(d^A_i = 1/\pi(\mathbf{x}_i, \theta)\)として、$$ \hat{\mu}_{DR1} = \frac{1}{N} \sum_{i \in S_A} d^A_i (y_i – m(\mathbf{x}_i, \theta)) + \frac{1}{N} \sum_{i \in S_B} d^B_i m(\mathbf{x}_i, \beta) $$ $$ \hat{\mu}_{DR2} = \frac{1}{\hat{N}^A} \sum_{i \in S_A} d^A_i (y_i – m(\mathbf{x}_i, \theta)) + \frac{1}{\hat{N}^B} \sum_{i \in S_B} d^B_i m(\mathbf{x}_i, \beta) $$ この理論的特性はどうなっておるんかというと… [要するに、傾向スコアモデルとアウトカム回帰モデルのどっちかが正しく指定されていれば一致推定量である。再び\(O\)-記法が乱れ飛ぶ世界の話になるのでスキップ。なおIPW推定量とDR推定量の効率性の比較は難しい問題なのだそうだが、回帰モデルがフィットしていればDR推定量のほうが分散が小さいとのこと]
4. 分散推定
[この辺の話は！プロの方に全部お任せしたく！何卒よしなによろしくお願いいたします！]
5. シミュレーション研究
サイズ20000の有限母集団があって、$$ z_{1i} \sim Bernoulli(0.5)$$ $$z_{2i} \sim Uniform(0,2)$$ $$ z_{3i} \sim Exponential(1) $$ $$ z_{4i} \sim \chi^2(4)$$ $$ x_{1i} = z_{1i} $$ $$ x_{2i} = z_{2i} + 0.3 x_{1i}$$ $$ x_{3i} = z_{3i} + 0.2(x_{1i}+x_{2i}) $$ $$ x_{4i} = z_{4i} + 0.1(x_{1i} + x_{2i} + x_{3i}) $$ $$ y_i = 2 + x_{1i} + x_{2i} + x_{3i} + x_{4i} + \sigma \epsilon_i, \ \epsilon_i \sim_{iid} N(0,1) $$ だとする。\(\sigma\)は\(y\)と\(\mathbf{x}^\top \beta\)の相関が0.3, 0.5, 0.8になるように選ぶ。
ここから\(S_A\)を抽出する。真の傾向スコアモデルは $$ \log \left( \frac{\pi^A_i}{1-\pi^A_o} \right) = \theta_0 + 0.1 x_{1i} + 0.2 x_{2i} + 0.1 x_{3i} + 0.2 x_{4i} $$ \(theta_0\)は標本サイズが\(n_A\)になるようにうまく決める。この\(\pi^A_i\)でポアソン抽出する。[ポアソン抽出で標本サイズを狙い撃ちするのは結構骨だと思うんだけど、なんか頑張ったんでしょうね。まあそこは本題じゃないしね]
いっぽう、\(S_B\)(サイズ\(n_B\))は\(\pi^B_i \propto z_i = c + x_{3i}\)でPPS抽出する。\(c\)は\(z_i\)の最大値と最小値の比が50になるようにうまく選ぶ(ウェイトの分散をコントロールしている)。
以下の4つのシナリオを考える。
- TT: アウトカム回帰モデル\(\xi\)も傾向スコアモデル\(q\)も正しい。
- FT: \(\xi\)をうっかり\(m(\mathbf{x}_i, \beta) = \beta_0 + \beta_1 x_{1i} + \beta_2 x_{2i} + \beta_3 x_{3i} \)と指定(\(x_{4i}\)を落としている)。\(q\)は正しい。
- TF: \(\xi\)は正しい。\(q\)からうっかり\(x_{4i}\)を落としている。
- FF: 両方のモデルからうっかり\(x_{4i}\)を落としている。
標本サイズ\((n_A, n_B\)はそれぞれ500か1000、計4通り。
実験1。推定量の相対バイアスとMSEを調べる。5000試行のシミュレーション。選手は以下の通り: \(\hat{\mu}_{IPW1}, \hat{\mu}_{IPW2}, \hat{\mu}_{REG}, \hat{\mu}_{DR1}, \hat{\mu}_{DR2}\)。ゲスト選手として、標本平均\(\hat{\mu}_A = \bar{y}\)、 \(\mathcal{S}_A\)と\(\mathcal{S}_B\)の積上データを使ってウェイトなして\(\tilde{R}_i\)の確率を推定してIPW推定した\(\hat{\mu}_{C1}, \hat{\mu}_{C2}\)にもご参加いただきます。
結果。[面倒なので表は見ずに本文中のコメントをメモする]
- IPW推定量の性能は、TT, FTなら(つまりモデル\(q\)が正しければ)良好、そうでないと悪い。
- 予測推定量\(\hat{\mu}_{REG}\)の性能は、TT, TFなら(つまりモデル\(\xi\)が正しければ)とても良好、そうでないと悪い。
- 二重頑健推定量はTT,FT,TFで良好。
- FFではどの推定量も悪い。
- ゲスト選手の皆さんの性能はどのシナリオでも悪い。[はっはっは、性格悪いなあ… \(\mathcal{S}_B\)はよりによって共変量をサイズ変数にしたPPS抽出という一番やばいやつなのに、それを単純無作為抽出だとみなしているんだから、そりゃあ\(\hat{\mu}_{C1}, \hat{\mu}_{C2}\)はぼろ負けするだろうさ]
実験2。分散推定量の挙動を調べる。[すいませんがスキップ]
6. Pewのサーベイデータへの適用
[非確率標本は、Pew Research Centreの調査データ、参照標本はBRFSSとCPSという公的データ。すいませんけどスキップ]
7. Additional Remarks
本論文の基盤になるのは仮定A1-A3である。しかし実務的には、補助変数\(\mathbf{x}\)に抽出メカニズムを特徴づける要素がすべて入っているのかを判断するのは難しい。一般的原則としては、性年代とかSESとか、調査参加傾向と反応変数の両方を説明しそうな変数を入れておくこと。
非確率的な手法で収集されたデータがどんどん増えている昨今、多くの統計エージェンシーが採用しているサーベイ中心的なアプローチは、新しいデータの時代に適応して進化しなければならない。[そう、それは痛感しております]
母集団に傾向スコアがゼロの個体がいるというシナリオの場合(確率標本でいうアンダーカバレッジ問題)、標本が表現している母集団はなにかという点について注意深い評価が必要である。問題の深刻さは、カバーされていない個体の割合、そして母集団の2つのパートの反応変数の差によって決まる。本論文のフレームワークをこの問題に拡張するのは今後の課題である。
—————–
難しいところを何か所かスキップしておいてなんだが、大変わかりやすく書かれた論文であったと思う。先にWuさんのレビューを読んでいたからかもしれないけれど、楽しく読み終えることができた。JASAなどという統計学のトップジャーナルであっても、あんまり怖れずに気楽に読んじゃえばいいのである。いや、てんで歯が立たない奴とか晦渋すぎる奴とかもあるけどさ。
実務側としては、推定量と分散推定については先生方にお任せしますんで、共変量の選び方、ならびに「共変量はこれで足りておるのだろうか」問題に対するなんらかの示唆を頂きたいところなんだけど、それって難しい問題なんでしょうね、きっと。
ちょっと疑問に思っているのは、この論文でいう仮定A1とA2、すなわち、母集団における共変量の下での標本包含と目的変数の条件付き独立性という仮定と、母集団における標本包含確率の正値性という仮定を、なぜ別建てにしているのだろうか、という点である。A1が成立していさえすれば(つまり、\(\mathbf{x}_i\)の下で\(R_i\)と\(y_i\)が条件付き独立ならば)、たとえA2が破られていても(母集団に包含確率がゼロの人がいたとしても)、IPW推定量なり予測推定量なり二重頑健推定量なりはモデルが正しければ不偏性を持つのではなかろうか。いや、A1とA2では実務的な意味合いが異なる(A1は無回答のMAR仮定、A2は台帳のアンダーカバレッジに相当する)というのはよくわかるんだけど。