Rueda, M., Ferri-García, R., & Castro L. (2020) The R package NonProbEst for estimation in non-probability surveys. The R Journal, 12(1), 406-418.
　ここんとこ非確率標本に基づく推定の話を調べていて、含蓄の深い話が続いて心底くたびれた。リハビリのため、あまり深い話が出てこないであろう資料を手に取った。
　非確率標本のためのRパッケージNonProbEstの紹介記事。要はソフトの機能紹介で、後半はコード例である。癒されるねえ。

1. イントロダクション
　オンラインの非確率調査における選択バイアスを修正する方法について、詳しくは Elliott & Valliant(2017) をみてください。重要なアプローチは以下の3つです。

• 擬似デザインベース推論(擬似ランダム化)。ウェイトをつくって選択バイアスを補正する。まず、反応確率を推定してHT推定量なりHajek推定量なりを使うという方法がある。反応確率推定の方法としては傾向スコア調整がある。
• 標本マッチング。
• 予測モデル。目標変数を従属変数として非確率標本でモデルを作り、これを確率標本にあてはめて予測する。

NonProbEstはこれらの方法のなかのいくつかを実装しております。

2. 統計的手法
　有限母集団を\(\mathcal{U} = \{1, \ldots, k, \ldots, N\}\)とする。そのサブセットとしてオンライン母集団\(\mathcal{U}_V\)があり、そこから自己選択された任意型非確率標本\(s_V\)(サイズ\(n_V\))があるとする。
　関心ある変数を\(y\)、その母合計を\(Y\)とする。補助情報がなかったらHT推定量 $$ \hat{Y}_{HT} = \sum_{k \in S_V} w_{vk} y_k $$ が使われる。もっともシンプルには\(w_{vk} = N/n_V\)。これは\(s_V\)を\(\mathcal{U}\)のSRS標本だと捉えているわけだ。当然バイアスを受ける。主要なのは選択バイアスとカバレッジバイアス(\(\mathcal{U}\)と\(\mathcal{U}_V\)のずれ)である。
　ウェイティングの成功のカギは強力な補助情報である。このパッケージで補助情報の形式を3つに分けている。

• InfoTP。補助情報の母合計がわかっている。
• InfoES。確率標本の全個体で補助情報の値がわかっている。
• InfoEP: 母集団の全要素で補助情報の値がわかっている。

[3節以降ではInfoTP, InfoSP, InfoUPという名前にしれっと変わっている。ちゃんと校正してほしいなあ。略語についての説明はないが、TはTotal, Sは標本, Uは母集団を表しているのだろう]
以下では補助情報の形式別に、バイアスを扱う主な方法を説明する。

3. InfoTP [補助情報の母合計がわかっているときね]
3.1 カリブレーション
　個体\(k\)の補助変数ベクトルを\(\mathbf{x}_k\)とし、その母合計\(\mathbf{X}\)が既知とする。\(Y\)のカリブレーション推定量とは、オリジナルの標本ウェイト\(w_{vk}\)を$$ \sum_{k \in s_V} w_k \mathbf{x}_k = \mathbf{X} $$ を満たすよう最小限に修正したウェイト\(w_k\)を使うというもの。具体的には、擬似距離\(G\)を定義して、$$ min_{w_k} \{ \sum_{k \in s_V} G(w_k, w_{vk})\} $$ を最小化する。擬似距離としては線形距離が用いられることが多い。もしカイ二乗距離なら、GREG推定量 $$ Y_{reg} = \sum_{s_V} w_k y_k = \sum_{s_V} d_k y_k + (\mathbf{X} – \sum_{s_V} w_{vk} \mathbf{x}_k)^\top \hat{B}_{s_V} $$ $$ \hat{B}_{s_V} = \frac{ \sum_{s_V} w_{vk} \mathbf{x}_k y_k }{\sum_{s_V} w_{vk} \mathbf{x}_k \mathbf{x}^\top_k } $$ となる。[\(d_k\)ってなによ…?]
　この方法でバイアスが取れるのはMARのときだけで、NMARではうまくいかないといわれている。

4. InfoSP [確率標本の全個体で補助情報の値がわかっているときね]
4.1 傾向スコア調整
　傾向スコア調整自体はRosenbaum & Rubin(1983)にはじまっており、標本抽出の文脈ではRubin(1986)にさかのぼるが、オンライン調査への適用は2000年代初頭から。Taylor et al.(2001 IJMR)をみよ。

　[以下、ちょっと納得いかないので細かくメモする]
　非確率調査と参照調査の対象者を結合し、ある対象者がどちらの調査に参加したかを表す二値変数\(I_{s_V}\)についての予測モデル(ロジスティック回帰とすることが多い)を構築して、そこから、ある個人が非確率調査に参加する傾向を得る。[んんん？ この言い方で合っているの？ 推定対象は母集団における個体が非確率調査に参加する確率ですよね? 手順の説明なのであれば、これはWuさんが批判しているValliant & Deverの方法ではないだろうか]
　対象者が非確率標本に参加している確率を $$ \pi(\mathbf{x}) = \frac{1}{\exp(-\gamma^\top \mathbf{x}) + 1} $$ とする。で、推定された反応傾向\(\hat{\pi}(\mathbf{x}_k)\)を得る。推定量の作り方にはいろいろあって、

• Valliant(2019 J.SurveyStat.Methodol.): $$ \hat{Y}_{PSA1} = \sum_{k \in s_V} y_k w_k^{PSA1}$$ $$ w_k^{PSA1} = \frac{w_{V_k}}{\hat{\pi}(\mathbf{x}_k)} $$　[\(w_{V_k}\)ってなんなのさ。オリジナルの標本ウェイト\(w_{v_k}\)のことか？]
• Schonlau & Couper (2017 Stat.Sci.) : $$ \hat{Y}_{PSA2} = \sum_{k \in s_V} y_k w_k^{PSA2}$$ $$ w_k^{PSA2} = \frac{1-\hat{\pi}(\mathbf{x}_k) }{\hat{\pi}(\mathbf{x}_k) } $$
• Valliant & Dever (2011 Soc.Methods&Res.): \(\hat{\pi}(\mathbf{x}_k)\)で標本を並べて\(g\)クラスに等分し、平均\(\bar{\pi}_g\)を求めて、$$ \hat{Y}_{PSA3} = \sum_g \sum_{k \in s_{V_g}} y_k w_k^{PSA3}$$ $$ w_k^{PSA3} = w_{V_k} / \bar{\pi}_g $$
• Lee & Valliant (2009 Soc.Methods&Res.): \(g\)クラスにわけるところまで同じ。参照標本側(R)のウェイトをつかう。$$ \hat{Y}_{PSA4} = \sum_g \sum_{k \in s_{V_g}} y_k w_k^{PSA4}$$ $$ w_k^{PSA4} = w_{V_k} f_g$$ $$ f_g = \frac{ \frac{\sum_{k \in s_{R_g}} w_{R_k}}{\sum_{k \in s_R} w_{R_k}} }{ \frac{\sum_{k \in s_{V_g}} w_{V_k}}{\sum_{k \in s_V} w_{V_k}} } $$

PSAはバイアスをうまく取り除けることもあるけれど、分散の増大も招く。カリブレーション調整と組み合わせるという提案もある。上記のほかにFerri-Garcia & Rueda (2018 SORT)をみよ。
　PSAにおける分散推定としては…[大事な話だろうけどそこまで頭が回らない。パス]

4.2 統計的マッチング
　この路線はRivers(2007)にはじまる。$$ \hat{Y}_{SM} = \sum_{s_R} \hat{y}_k w_{R_k} $$ とする。では\(\hat{y}_k\)はどうするか。もっとも多いのは$$ \hat{y}_k = \mathbf{x}^\top_k \hat{\beta}$$ なんだけど、他の方法も提案されている。
　この手法の欠点は、非確率標本の精度が落ちること。

5. InfoUP [母集団の全個体で補助情報の値がわかっている]
母集団\(\mathbf{y} = (y_1, \ldots, y_)^\top\)が超母集団からえた実現値だと考えて、超母集団モデル $$ Y_k = m(\mathbf{x}_k) + e_k $$ を考える。で、どうにかして $$ \hat{y}_k = E_m(y_k | \mathbf{x}_k) $$ を得る。添え字の\(m\)は具体的なモデルのもとでのって意味ね。
　推定量の作り方はいろいろある。

• モデル・ベース推定量: $$ \hat{Y}_m = \sum_{k \in s_V} y_k + \sum_{k \in \bar{s_V}} \hat{y}_k $$
• モデル・アシステッド推定量: $$ \hat{Y}_{ma} = \sum_{k \in \mathcal{U}} \hat{y}_k + \sum_{k \in s_V} (y_k – \hat{y}_k) w_{V_k} $$
• モデル・カリブレット推定量: $$ \hat{Y}_{mcal} = \sum_{k \in s_V} y_k w_k^{CAL} $$ \(w_k^{CAL}\)は、\(\sum_{k \in s_V} w_k^{CAL} \hat{y}_k = \sum_{k \in \mathcal{U}} \hat{y}_k \) となるようにカリブレートしたウェイト。

6. 非確率標本における機械学習アルゴリズムの使用
　いろんな提案がある。近年のモデル・アシステッド推定量については Breidt et al.(2017 Stat.Sci.)とかをみよ。いくつか挙げると:

• Phipps et al(2012 Annals.App.Stat.): 無回答傾向のPSAへの分類木・回帰木の活用
• Buskirk & Kolenikov (2015 SurveyMethods): 無回答傾向のPSAへのランダム・フォレストの活用
• Chen et al.(2019 JRSS:C): LASSOによる調整
• Buelens et al.(2018 Int.Stat.Rev.): レビュー

7. RパッケージNonProbEst
[関数の紹介。ちゃんと読んでないけど、提供されているのはたぶん上に出てきた奴らだけで、二重頑健推定はないと思う。機械学習系の手法ってどうなってんの？ という興味があったのだが、マニュアルをみると、どうやら具体的なモデルの指定についてはcaretパッケージにお任せしているようだ]

8. NonProbEstを使った非確率標本による推定
[コード例。パス]
——————–
最近、なんか疲れるのが早くなってきてしまった。私ももう30を超えているからなあ… (嘘ではないがまぎらわしい言明)