Chrostowski, L., Chlebicki, P., & Beresewicz, M. (2025) nonprobsvy – An R Package for modern methods for non-probability surveys. arXIv.

　非確率標本に基づく母集団特性の推定のためのRパッケージとして、CRAN Task Viewの “Official Statistics & Survey Statistics“には、NonProbEstパッケージとnonprobsvyパッケージというのが載っている。前者についてはR Journalに記事があった。後者のvignetteにあたるらしきプレプリントを見つけたので読んでみた。こっちのほうが新しい機能を積んでいるようだ。

1. イントロダクション
　[…大幅中略…]
　そんなわけで、非確率標本のためのパッケージとしては、RだとNonProbEst, Pythonだとbalanceとinpsがある。他にRでは、MRPをrstanarmでやったり、一般化標本選択モデルのGJRMパッケージというのがあったりする。
　NonProbEstはちょっと古い。GJRMの提供する標本選択モデルは観察研究の選択バイアスを修正するために広く使われているけれど、母集団特定の推定の際にどう使えるかについては理論がない。[GJRMって初めて聞いた… NMAR下でのアウトカムと欠損の同時モデルみたいなもの？]
　というわけで、nonprobsvyパッケージを作りました。

2. 非確率標本のための手法
2.1 基本的セットアップ
　母集団\(U = \{1, \ldots, N\}\), 補助変数\(\mathbf{x}_i\), 調査変数\(y_i\)。非確率標本\(S_A\), 確率標本\(S_B\)、デザインウェイト\(d^B_i = 1/\pi^B_i\)。標本包含インジケータ\(R^A_i, R^B_i\)、未知の傾向スコア\(\pi^A_i=P(R^A_i = 1 | \mathbf{x}_i, y_i) = P(R^A_i = 1 | \mathbf{x}_i)\)。推定対象 \(\mu_y = \frac{1}{N} y_i\)。[最近このセッティングの論文ばかり読んでいるのでメモもぶっきらぼうになる]
　以下に登場する多くの手法は以下の仮定に基づく。

• A1. \(R^A_i\)と\(y_i\)は\(\mathbf{x}_i\)のもとで独立 (つまりMARメカニズム)。
• A2. 母集団の全個体で \(\pi^A_i \gt 0\)。
• A3. 補助変数の下で包含インジケータは独立。

　\(S_A\)と\(S_B\)のあいだに重複なし、\(y_i\)に測定誤差なし。

2.2 予測ベース手法

　まずは予測推定量について。
　以下のセミパラメトリックモデルを仮定する。$$ E_\xi(y_i | \mathbf{x}_i) = m(\mathbf{x}_i, \beta) $$ $$ V_\xi(y_i | \mathbf{x}_i) = v(\mathbf{x}_i) \sigma^2$$ モデルは線形かもしれないし、一般化線形モデルのように非線形かもしれない。パラメータは疑似最尤推定法で推定できる。
　Wu(2022)が述べているように、よく使われる予測推定量の形式がふたつある。\(m_i = (\mathbf{x}_i, \beta_0), \hat{m}_i = m(\mathbf{x}_i, \hat{\beta})\)として、$$ \hat{\mu}_{y, PR1} = \frac{1}{N} \sum_{i=1}^N \hat{m}_i $$ $$ \hat{\mu}_{y,PR2} = \frac{1}{N} \left( \sum_{i \in S_A} y_i – \sum_{i \in S_A} \hat{m}_i + \sum_{i=1}^N \hat{m}_i \right) $$ \(m(\mathbf{x}_i, \beta) = \mathbf{x}_i^\top \beta\)なら, \(\mu_\mathbf{x} = \frac{1}{N} \sum_{i=1}^N \mathbf{x}_i, \bar{x}_A = \frac{1}{n_A} \sum_{i \in S_A} \mathbf{x}_i\)として、$$ \hat{\mu}_{y, PR1} = \mu^\top_\mathbf{x} \hat{\beta}$$ $$ \hat{\mu}_y = \frac{n_A}{N} (\bar{y}_A – \bar{x}_A^\top \hat{\beta}) + \mu^\top_\mathbf{x} \hat{\beta} $$ となる。つまり、母集団合計ないしその推定値だけがわかっておればよいことになる。

　マス代入推定量について。\(S_B\)のすべての\(y_i\)が欠損と考えて代入する。$$ \hat{\mu}_{y, MI} = \frac{1}{\hat{N}_B} \sum_{i in S_B} d^B_i y^*_i $$ マス代入推定量は、もし回帰代入が決定論的なのであれば予測推定量と同じである。
　このパッケージでは以下をご用意している。

• GLMに基づく手法(MI-GLM)。その性質についてはKim, Park, Chen, Wu(2021 JRSS:A)を参照。このパッケージではgaussian, binomial, poissonをご用意している。
• カーネル平滑化による手法(MI-NPAR)。詳しくはChen, Yang, & Kim(2022 J.SurveyStat.Methodol.)をみてね。このパッケージではstats::loessを使った方法をご用意している。いずれはnpかKernSmoothパッケージを使った手法もご提供するよ。[npとはノンパラメトリックカーネル法のパッケージらしい]
• 最近隣法を用いた手法(MI-NN)。もともとRivers(2017 Conf.)が標本マッチングと呼んでいた手法。詳しくはYang, Kim, & Hwang (2021 SurveyMethodol.)をみよ。共変量の数が多いときにバイアスが生じる。
• 予測平均マッチング(MI-PMM)。詳しくはChlebicki, Chrostowski & Beresewicz(2024 arXiv)をみよ。簡単にいえばこういう手法である。まず\(S_A\)にモデル\(m(\mathbf{x}_i, \beta)\)をあてはめる。次に、\(S_A, S_B\)の両方に予測値\(\hat{m}_i\)を振る。で、\(S_B\)の全個体を、\(\mathbf{x}_i\)じゃなくて\(\hat{m}_i\)で\(S_A\)の個体に最近隣マッチングさせる(ないし、\(S_B\)側は\(\hat{m}_i\), \(S_A\)側は\(y_i\)で最近隣マッチングさせる)。こうすると共変量の数が多くてもバイアスが生じない。[へえええ？？？ なんだかナイーブな手法のようにみえるけど…]

　分散推定について。[パスするけど、MI-GLM, MI-NPARについては分析的に出すのとブートストラップで出すのをご用意、MI-NN, MI-PMMについて後者のみ、とのこと]

2.3 逆確率ウェイティング
　IPW推定量にはHT形式とHajek形式がある。$$ \hat{\mu}_{y,IPW1} = \frac{1}{N} \sum_{i \in S_A} \frac{y_i}{\hat{\pi}^A_i} $$ $$ \hat{\mu}_{y,IPW2} = \frac{1}{\hat{N}^A} \sum_{i \in S_A} \frac{y_i}{\hat{\pi}^A_i} $$ 詳しくは、Lee(2006 J.OfficialStats.), Biffignandi & Bethlehem(2021 Book), Elliott & Vallient (2017), Chen et al.(2020), Wu(2022)をみよ。
　\(\hat{\pi}^A_i\)は擬似対数尤度$$ l^*(\gamma) = \sum_{i \in S_A} \log \frac{\pi(\mathbf{x}, \gamma)}{1-\pi(\mathbf{x}, \gamma)} + \sum_{i \in S_B} d^B_i \log(1-\pi(\mathbf{x}_i, \gamma)) $$ で得る。ロジット関数を仮定すると擬似スコア関数はこうなる[式省略]。\(qp\)同時ランダム化の観点からの期待値をとると、パラメータの真値を\(\gamma_0\)として\(E_{qp}(U(\gamma_0)) = \mathbf{0} \)である。このことは、\(\hat{\alpha}\)が\(\alpha_0\)の一致推定量であることを意味する[いきなり\(\alpha\)ってのが出てきたけど\(\gamma\)の誤植ではないかしらん]。
　一般推定方程式の観点から見ると… [カリブレーションの式が導出できるという話。云々。このへん、Wu(2022)の内容と同一と思うのでメモ省略]

　分散推定は…[パス]

2.4 二重頑健推定
　DR推定量の一般的な式は$$ \tilde{\mu}_{y, DR} = \frac{1}{N} \sum_{i \in S_A} \frac{y_i – m_i}{\pi^A_i} + \frac{1}{N} \sum_{i=1}^N m_i $$ なんだけど、推定にあたっては2つのアプローチがある。

　その1、Chen et al.(2020)のアプローチ。
　\(d^A_i = \pi(\mathbf{x}_i, \gamma)^{-1}\)とする[あれ？ \(\pi(\mathbf{x}_i, \hat{\gamma})^{-1}\)じゃないの？]。$$ \hat{\mu}_{y, DR1} = \frac{1}{N} \sum_{i \in S_A} d^A_i \left( y_i – m(\mathbf{x}_i, \hat{\beta}) \right) + \frac{1}{N} \sum_{i \in S_B} d^B_i m(\mathbf{x}_i, \hat{\beta}) $$ $$ \hat{\mu}_{y, DR2} = \frac{1}{\hat{N}_A} \sum_{i \in S_A} d^A_i \left( y_i – m(\mathbf{x}_i, \hat{\beta}) \right) + \frac{1}{\hat{N}_B} \sum_{i \in S_B} d^B_i m(\mathbf{x}_i, \hat{\beta}) $$ バイアスとMSEの観点からは後者のほうがよい[え？ 分散は小さいだろうけれど、バイアスの観点からもそうなの？]
　2つのモデルをバラバラに推定する。使う変数を変えてもよい。

　その2、Yang, Kim, & Song (2020 JRSS:B)のアプローチ。
　[数式が長いし難しいのでメモは省略するけど、\(\hat{\mu}_{DR}\)のバイアスを最小化するような\(\beta, \gamma\)を一気に求めるのだそうだ。ただし、2つのモデルの次元が同じでなければならないので、最初に変数選択をするとのこと。へー。これはたぶん元論文を読んだ方がわかりやすいだろう]
　
　分散推定。
　\(\hat{\mu}_{y, DR1}\) [HT型のほうね] の分散推定量は閉形式で出せる。ところが\(\hat{\mu}_{y, DR2}\) [Hajek型のほう] は出せない。仕方がないのでこのパッケージでは\(\hat{N}_A, \hat{N}_B\)を\(N\)に置き換えて出しちゃっている。それよかブートストラップ法のほうがお勧め。

2.5 変数選択アルゴリズム
　結局のところ変数選択がカギだ。モデルの安定性と計算のしやすさにも影響するし、いらん変数をいれていると分散が増すかもしれんし。

　一番ポピュラーなのは罰則法だ。LASSOとか、smoothly clipped absolute devivation(SCAD)とか、minimax concave penalty(MCP)とか。損失関数が適切なら効かない変数の係数を退化させることができる。
　非確率標本の場合、損失関数が外的データソースを説明するように修正する[?]。たとえばYang et al.(2020)はこう提案している。アウトカムモデルを\(m(\mathbf{x}_i, \beta(\lambda_\beta))\)として(\(\lambda_\beta\)はチューニングパラメータ)、$$ Loss(\lambda_\beta) = \sum_{i=1}^N R^A_i [y_i – m(\mathbf{x}_i, \beta(\lambda_\beta))]^2 $$ とする[標本におけるアウトカムモデルの誤差の二乗和を最小化するようにチューニングするわけね]。いっぽう傾向スコアモデルについては、IPW推定量の場合の損失関数は$$ Loss(\lambda_\gamma) = \sum_{j=1}^p \left( \sum_{i=1}^N \left[ R^A_i – \frac{R^B_i \pi(\mathbf{x}_i, \gamma(\lambda_\gamma))}{\pi^B_i} \right] x_{ij} \right)^2 $$ とする[えーっと、ある共変量について、\(S_A\)に入ってる人の値の合計と、\(S_B\)側から求まる「傾向スコアと値の積」の母合計のHT推定量との差を二乗和を出す。これを共変量を通じて合計している。つまり、\(S_A\)での共変量の合計を\(S_B\)で当てようとしたときの推定誤差を損失関数にしているわけだ。なるほど]。
　カリブレーション推定量の場合の損失関数は…[略]

　Yang et al.(2020)は、罰則がSCADでIPW推定量の場合についてのみ検討している。本パッケージはその他いろいろ提供します。
　[申し訳ないけど話のポイントがわかってない。アウトカムモデルと傾向スコアモデルの説明変数セットは同じ、どちらも正則化項を入れる、チューニングパラメータはそれぞれについて決める、ってこと？ もしそうなら、説明変数セットがモデル間で違っててもよくない？]

3. 主要な関数とパッケージの機能
[ここからはパッケージの機能説明。見出しのみメモする]

3.1 nonprob関数
これが主力関数。引数data, selection(傾向スコアモデルのformula), outcome(アウトカムモデルのformula), target (targetにする変数のformula)を指定する。参照標本のデザインについてsurveyパッケージのsvydesign2クラスのオブジェクトを渡せる。[…]

3.2 推定量のタイプの制御
3.3 分散推定の制御

4. データ分析例
[これは必要になったら読もう…]

5. クラスとS3Methods
[パス]

6. 要約と今後の課題
[…]
　今後の課題: モデル・ベースのカリブレーション[あれ？ カリブレーション推定もカバーしているんでしょ？ あ、そうじゃなくて、アウトカムモデル側でカリブレーションするのはカバーしてないってこと？]。確率標本と非確率標本に重複がある場合への対応。無視可能でない標本選択への対応。カリブレーション・ウェイトとの一致性の維持の理論[?]。複数の非確率標本の統合。二重頑健推定の経験尤度アプローチ。debiased/double機械学習のような手法の統合。ビッグデータ変数の測定誤差の扱い。ブートストラップの拡張。
————–
　よくわからん点もあったけど、勉強になりましたですー。
　理論のほうは基本的にWu(2022)に沿った内容だと思うんだけど、二重頑健推定のパラメータ推定を両方のモデルで一気にやるという提案があるとは知らなんだ。へえええ。