Yang, S., Kim, J.K., Song, R. (2020) Doubly robust inference when combining probability and non-probability samples with high dimensional data. Journal of the Royal Statistical Society Series B: Statistical Methodology, 82(2), 445–465.
非確率標本と確率標本があるとき、前者をうまいことウェイティングしましょうとか、前者でモデルを組んで後者にあてはめて母集団特性を予測しましょうといった方法があるけれど、最近ではその延長線上に、両方やって二重頑健推定しましょうという話もある。ひゅー、かっこいい、今流行りの因果推論みたいだ。私も市場調査みたいな地味な仕事じゃなくて、web広告の最適化とかでぶいぶいいわせられるかもしれない。よーし転職してタワマンの上の方に住んで港区女子と不倫するぞ！ (←貧困なイメージ)
まあとにかくそういうとき、かつ、共変量がたくさんありすぎて困っちゃうとき(どんどんSFっぽくなっていくね)、罰則付き回帰で変数選択するぞ、という論文。
1. イントロダクション
[…]
非確率標本と確率標本を統合する方法は大きく3つにわけられる。
- 傾向スコア調整。Lee & Valliant(2009 Sco.Meth.Res.), Valliant & Dever(2011), Elliott & Valliant(2017), Chen, Li, & Wu(2018), Stuart et al.(2011 JRSS:A, 2015 Prev.Sci.)をみよ。Buchanan et al.(2018 JRSS:A)はRCTの結果の一般化に傾向スコア調整を使っている。O’Muircheartaigh & Hedges(2014 Appl.Statist.)は非ランダム化実験の分析に傾向スコア層別を使っている。傾向スコア調整の弱点は傾向スコアモデルの誤指定である。
- カリブレーション・ウェイティング。Deville & Sarndal(1992 JASA), Kott(2006 Surv.Methodol.), Chen, Valliant, & Elliott(2018 Surv.Methodol.). Chen et al.(2019 Appl.Statist.)をみよ。
- mass imputation. Breidt et al.(1996 J.Ind.Soc.Agri.Statist.), Kim & Rao(2012 Biometrika), Chipperfield et al.(2012 Aust.NZ.J.Statist.), Yang & Kim(2018 arXiv)をみよ。
非確率調査Aで(X,Y)を、確率調査BでXを観察しているとする。分析者はYの母平均を推定するため、抽出とアウトカムの両方の予測子である共変量をコントロールしようとする。X変数がいっぱいあるときは変数選択が重要になる。予測の変数選択と違って、データ統合のための変数選択の研究は少ない。
- Gao & Carroll(2017): 擬似尤度アプローチ。すべての尤度を正しく指定する必要があり、モデル誤指定に弱い。
- Chen, Valliant, & Elliott (2018): lasso回帰を使ったモデルベース・カリブレーション。アウトカムモデルの誤指定に弱い。
本論文では二重頑健な変数選択法を提案します。[手法の概要紹介、2パラグラフ。メモ省略]
2. 基本セットアップ
2.1 表記: 2標本の場合
有限母集団がある。サイズ\(N\)は既知、インデクス集合を\(\mathcal{U} = \{1, \ldots, N\}\)とし、\(\mathcal{F}_N = \{(X_i, Y_i):i\in \mathcal{U}\)とする。母平均\(\mu = \frac{1}{N} \sum_{i=1}^N Y_i\)に関心がある。
確率標本A(サイズ\(n_A\)), 非確率標本B(サイズ\(n_B\))がある。Aでは\(\mathcal{O}_A = \{(d_{A,i} = \pi^{-1}_{A,i}:i \in \mathcal{A}\}\)が観察されている(ただし\(\pi_{A,i} = P(i \in \mathcal{A})\))。Bでは\(\mathcal{O}_B = \{(X,Y):i \in \mathcal{B}\}\)が観察されている。選択インジケータを\(I_{A,i}, I_{B,i}\)とする。
2.2 識別のための仮定
有限母集団を生成する超母集団モデル\(\zeta\)を考え、\(Y\)条件付き確率分布を\(f(Y|X)\)とする。
仮定1.
(a) \(P(I_B = 1 | X, Y) = P(I_B | X) \)。これを抽出スコア\(\pi_B(X)\)と呼ぶ。
(b) すべての\(X\)について\(\pi_B(X) \gt N^{\gamma-1} \delta_B \gt 0\), ただし\(\gamma \in (\frac{2}{3}, 1]\)。
(a)は\(m(X) = E(Y|X) = E(Y|X, I_B = 1)\)を含意する。
(b)は\(\pi_B(X)\)の下界を特定している。先行研究での標準的条件は、\(\pi_B(X) \gt \delta_B \gt 0\)という意味での強い正値性を指定していた。しかしそれは\(n^{-1}_B = O(N^{-1})\)を意味し、調査の実践においては制約的かもしれない。ここでは我々はこの条件を緩和し、1未満でありうる\(\gamma\)について\(n^{-1}_B = O(N^{-\gamma})\)を許容する。[なんで包含確率の下界を決めないといけないのかよくわかんないんだけど、おそらく推定量の漸近特性について考える際に必要になるのだろう]
仮定1の下で、\(E(\mu\)\)は\(E[I_A m(X)]\)ないし\(E[I_B Y / \pi_B(X)]\)に基づいて識別可能となる。しかし仮定1はデータからは検証できない。リサーチャーは\(I_B, Y\)のたくさんの予測子について検討する羽目になることが多い。
2.3 既存の推定量
\(m(X), \pi_B(X)\)のモデルを\(m(X^\top \beta), \pi_B(X^\top \alpha)\)とする[線形モデルに決めちゃうんだ…]。これまでに提案されてきた\(\mu\)の推定量を紹介しよう。
- 抽出確率の逆数のスコアでウェイティング。$$ \hat{\mu}_{IPW}(\hat{\alpha}) = \frac{1}{N} \sum_{i=1}^N \frac{I_{B,i}}{\pi_B(X^\top_i \hat{\alpha})} Y_i $$ \(\hat{\alpha}\)を得る方法はいろいろある。
- Valliant & Dever(2011): 結合データ(ウェイトは\(d_{A,i}\)と1)で抽出スコアモデルを推定する。\(n_B\)が小さいときに妥当。
- Elliott & Valliant(2017): 結合データにおける標本選択のオッズ\(O_B(X) = P(I_B = 1|X, \mathcal{O}_A \cup \mathcal{O}_B) / P(I_B = 0 |X, \mathcal{O}_A \cup \mathcal{O}_B)\)を求め、\(\pi_B(X) \propto P(I_A=1 | X) O_B(X)\)とする。もし\(X\)が標本Aのデザイン変数に対応していなかった場合は、\(P(I_A = 1 | X)\)の追加モデルを仮定する必要がある。さらに、このアプローチにおいては\(X\)が高次元の時の変数選択が難しい。
- \(\hat{\alpha}\)を得る方法として、我々は\(\alpha\)の以下の推定方程式を用いる: $$ \sum_{i=1}^N \left\{ \frac{I_{B,i}}{\pi(X^\top_i \alpha)} – \frac{I_{A,i}}{\pi_{A,i}} \right\} h(X_i; \alpha) = 0$$ Kott & Liao(2017 J.Surv.Statis.Methodol.) は\(h(X;\alpha)=X\), Chen, Li, & Wu (2018)は\(h(X;\alpha) = \pi(X^\top; \alpha) X\)としている。
[やめて、推定方程式で書かないで… 問題を難しくしないでほしい…
Chen, Li, & Wu で与えられている推定方程式をこの論文の表記に直して書くと、$$ \sum_{i=1}^N \left\{ I_{B,i} – I_{A,i} d_{A,i} \pi(X_i^\top; \alpha) \right\} X_i = 0 $$ だと思うんだけど、上の式に一致する？ → あ、一致してるわ]
- 標本Aに基づくアウトカム回帰。$$ \hat{\mu}_{reg}(\hat{\beta}) = \frac{1}{N} \sum_{i=1}^M I_{A,i} d_{A,i} m(X^\top_i \hat{\beta}) $$ \(\hat{\beta}\)は\(\mathcal{O}_B\)に基づいて得る。
- カリブレーション・ウェイティング。$$ \hat{\mu}_{cal} = \frac{1}{N} \sum_{i=1}^N w_i I_{B,i} Y_i $$ ただし、\(\{w_i: i \in \mathcal{B}\}\)は以下のどちらかの制約を満たす。McCounville et al.(2017), Chen, Valliant, Elliott (2018), Chen et al. (2019)をみよ。
- \(\sum_{i \in S_B} w_i X_i = \sum_{i \in S_A} d_{A,i} X_i\)
[日本語訳: 非確率標本における\(X_i\)の\(w_i\)で重みづけた合計が、確率標本における\(X_i\)の抽出ウェイトで重みづけた合計と一致するようにする。下添字に\(S_B, S_A\)というのが出てくるけど、ここで初出だと思う。\(\mathcal{B}, \mathcal{A}\)の誤植かな？]
この場合、アウトカムモデルの線形性(つまりなんらかの\(\beta^*\)で\(m(X) = X^\top \beta^*\)が成り立つこと)、ないし抽出ウェイトの逆数のウェイトの線形性(つまりなんらかの\(\alpha^*\)で\(\pi_B(X)^{-1} = X^\top \alpha^*\)が成り立つこと)、に依存する。非連続変数の場合この線形条件はおそらく維持されない。これらのケースでは\(\hat{\mu}_{cal}\)は歪む。
[ここ、全然わかんない… カリブレーション・ウェイトでしょ？ アウトカムモデルとは関係なくない？] - \(\sum_{i \in S_B} w_i m(X_i; \hat{\beta}) = \sum_{i \in A} d_{A, i} m(X_i; \hat{\beta})\)
[日本語訳: 非確率標本におけるアウトカムモデルに基づく\(Y\)の期待値を\(w_i\)で重みづけた合計が、確率標本における同じアウトカムモデルに基づく\(Y\)の期待値を抽出ウェイトで重みづけた合計と一致するようにする]
この場合、\(m(X;\beta)\)の正しい指定に依存する。
- \(\sum_{i \in S_B} w_i X_i = \sum_{i \in S_A} d_{A,i} X_i\)
- 二重頑健推定。$$ \hat{\mu}_{dr}(\hat{\alpha}, \hat{\beta}) = \frac{1}{N} \sum_{i=1}^N \left[ \frac{I_{B,i}}{\hat{\pi}_B(X^\top_i \hat{\alpha})} \{ y_i – m(X^\top_i \hat{\beta})\} + I_{A,i} d_{A,i} m(X^\top_i \hat{\beta}) \right] $$
3. 高次元データにおける方法論
不要な共変量をモデルに含めてしまうと計算が不安定になったり推定誤差が大きくなったりするので、高次元データの場合は変数選択が必要になる。
任意の \( \mathcal{J} \subseteq \{1, \ldots, p\}\)について、\(\alpha\)の要素のうちそのインデクスが\(\mathcal{J}\)に入っている要素で作ったサブベクトルを\(\alpha_\mathcal{J}\)としよう。\(\mathcal{J}\)の補集合を\(\mathcal{J}^c\)とする。
任意の\(\mathcal{J}_1, \mathcal{J}_2 \subseteq \{1, \ldots, p\}\)と行列 \(\Sigma \in \mathbb{R}^{p \times p}\)について、\(\mathcal{J}_1\)に入っている行, \(\mathcal{J}_2\)に入っている列で作った部分行列を\(\Sigma_{\mathcal{J}_1, \mathcal{J}_2}\)とする。
以下では共変量の分散はほぼ1に標準化されているものとする。
仮定2. 標本Bの抽出メカニズム\(\pi_B(X)\)はロジスティック回帰モデル \(logit\{\pi_B(X^\top \alpha)\} = X^\top \alpha\)に従う。
仮定3. アウトカム平均関数 \(m(X)\)はリンク関数\(m(\cdot)\)の一般化線形回帰モデル \( m(X) = m(X^\top \beta)\)に基づく。
[表記が超・わかりにくい。非確率標本のいわば真の標本包含確率のようなものを \(\pi_B(X)\), それを予測するロジスティック回帰モデルを\(\pi_B(X^\top \alpha)\), アウトカムの期待値を\(m(X)\), それを予測するGLMモデルを\(m(X^\top \beta)\)と書いているのだと思う。同じ関数名を違う意味で使うのやめてくれませんか]
以下では、\(\alpha^*\)をKLダイバージェンスを最小化するパラメータとする。すなわち$$ \alpha^* = arg \min_{\alpha \in \mathbb{R}^p} E \left[ \pi_B(X) \log \left\{ \frac{\pi_B(X)}{\pi_B(X^T \alpha)} \right\} + \{1 – \pi_B (X)\} \log \left\{ \frac{1-\pi_B(X)}{1-\pi_B(X^\top \alpha)} \right\} \right] $$ また以下とする。$$ \beta^* = arg \min_{\beta \in \mathbb{R}^p} E[\{ Y – m(X^\top \beta)\}^2] $$ \(X_B(X^\top \alpha), m(X^\top \beta)\)は作業モデルで、誤っているかもしれない。正しく指定されれば、\(\pi_B(X) = \pi_B(X^\top \alpha^*), m(X) = m(X^\top \beta^*)\)となる。
[パラメータの真値のようなものを\(\alpha^*, \beta^*\)と呼んでいるんだと思う。KLダイバージェンスとMSE最小化の観点から定義しているけれど、これは実際に計算できるって話じゃないよね？]
提案手法。2ステップからなる。
ステップ1, 抽出スコアモデルとアウトカムモデルにおける重要な変数を選ぶ。
\(\alpha, \beta\)の推定関数を以下とする: $$ U_1(\alpha) = \frac{1}{N} \sum_{i=1}^N \left\{ \frac{I_{B,i}}{\pi_B(X^\top_i \alpha)} – \frac{I_{A,i}}{\pi_{A,i}} \right\} X_i $$ $$ U_2(\beta) = \frac{1}{N} \sum_{i=1}^B I_{B,i} \{ Y_i – m(X_i^\top \beta) \} X_i$$ これをくっつけて \(\theta = (\alpha^\top, \beta^\top)^\top, U(\theta) = (U_1(\alpha)^\top, U_2(\beta)^\top)^\top\)とする。[一本目はChen,Li,&Wuと同じね。2本目はアウトカムモデルを線形モデルとみなした場合のOLS推定量に対応する推定関数だ。だって残差二乗和 \( \sum_{i \in B} (Y_i – X_i^\top)\beta^2 \) を\(\beta\)で微分すると \( -2 \sum_{i \in B} (Y_i – X_i^\top) X_i \) だから。アウトカムモデルはほんとは一般化線形モデルなんだけど、変数選択の際には線形とみなしているわけだ]
\(p\)が大きいとしよう。罰則付き推定関数を考える。$$ U^p(\alpha, \beta) = U(\alpha, \beta) – \left( \begin{array}{c} q_{\lambda_\alpha}(|\alpha|) sgn(\alpha) \\ q_{\lambda_\beta}(|\beta|) sgn(\beta) \end{array} \right) $$ ここで \(q_{\lambda_\alpha}(|\alpha|)\)というのは、\(\alpha\)の各要素の絶対値をなんらかの連続関数\(q_{\lambda_\alpha}(\cdot)\)で変換して積み上げたベクトル。\(q_{\lambda_\alpha}(|\alpha|)sgn(\alpha)\)というのは要素ごとの積である。で、\(\alpha, \beta\)に共通して\(q_\lambda(x) = d p_\lambda (x) / dx\)で、\(p_\lambda(x)\)はfolded concave smoothly clipped absolute deviation penalty functionとする[わからん！]。
結局こういう式になる: $$ q_\lambda(|\theta|) = \lambda \left\{ I(|\theta| \lt \lambda) + \frac{(a \lambda – |\theta|)_+}{(a-1)\lambda)} I (|\theta| \geq \lambda) \right\} $$ ただし\(a \gt 0\)とする。\( (\cdot)_{+} \)は打ち切り線形関数で、\(x \lt 0\)のとき\((x)_+ = 0\), でなければ\((x)_{+} = x\)である。Fan \amp; Li (2001 JASA)に従い \(a=3.7\)とする。こうして選ばれた変数の和集合を\(X_\mathcal{C}\)とする。
[日本語でいうと、話を簡単にするために\(\alpha_1\)が正だとして、もし閾値の3.7倍を上回るようだったらなにもしない。閾値の3.7倍は下回るが閾値を超えているならば、推定関数の出力のうち\(\alpha_1\)のところから、下回った分/2.7 だけ引く。閾値を下回っていたら閾値を引く。あっれー… パラメータのある要素の推定値をゼロに縮小させたいとき、最小化したい関数にその分の罰則を加えるよね？ 推定関数の場合は罰則項を(足すんじゃなくて)引くの？ 頭が混乱してきたよ…]
ステップ2. \(X_\mathcal{C}\)だけで、2節で示した二重頑健推定量を得る。これを\(\hat{\mu}_{dr}(\hat{\alpha}, \hat{\beta})\)とする。
\(X_\mathcal{C}\)には、抽出スコアモデルで重要な変数もアウトカムモデルで重要な変数もふくまれているはずである。だからどっちかのモデルが正しく特徴づけられていれば漸近不偏である。そうでない場合はこれこれのアイアスが生じる[式省略]。そこで、これを最小化するためにこういう推定関数を使う[式省略]。
ステップ1と2では推定関数が異なる。ステップ1では\(\alpha\)と\(\beta\)の選択が分離されているので、どっちかのモデルに誤指定があっても安定する。
ステップ2では両方のモデルが誤指定されている可能性を考慮して漸近バイアスを最小化するための同時推定関数を使っているが、だからとして、両方誤指定されていても不偏推定量が得られるわけではないことに注意。
なお、他の戦略として以下がありうる。
- それぞれのモデルで選択された変数の和集合を使うのではなくて、それぞれの変数セットを使ってフィッティングする。そうしなかった理由は、ステップ2の同時方程式が同じ次元数でないといけないから。さらに、傾向スコアモデルにアウトカムの予測子を含めるとバイアスは増えずATE推定の精度が上がるということがわかっているので、変数選択の際には両方を考えたほうがいいだろう。
- 抽出スコアモデルには関連しているがアウトカムモデルには関連していない予測子は取り除く。たしかに、もし両方のモデルが正しく指定されているならばそっちのほうがいいだろう。でも、アウトカムモデルに誤指定があるかもしれない。
4. 計算
[さっすがにそこまでは付き合いきれない… まるごとパス]
5. 変数選択と推定の漸近的結果
[さすがは一流誌、そこまでやるのね… もちろんパス]
6. シミュレーション研究
[すいませんパスします]
7. 応用例
[Pew Research CenterのデータをBRFSSという確率標本でなんとかするという例。パス]
————-
全21ページ中9ページまでしか読んでないけど、読み終えたということにしておこう。おつかれさまでした！
素朴すぎる感想だけど、「いま手元に確率標本と非確率標本があって共通の共変量を測定しているの、でも共変量が多すぎるの、困る～」という場面がどのくらいあるのかがわからない。「共変量が超しょぼすぎて、こんなんで調整して意味あるのかどうかわかんない、泣きそう～」という場面がほとんどじゃないですかね。まあ、それは私がそういう場面ばかり経験しているということだろうけれども。
幸いなことに、提案手法はRのnonprobsvyパッケージで実装されている。今度試してみよう。