Wu, C. (2022) Statistical Inference with Non-probability survey samples. Survey Methodology, 48(2), 283-311.
仕事の都合で読んだ奴。非確率標本からの母集団特性推定についての研究レビュー。掲載誌はカナダ統計局の発行。
本文は35ページ、その後5人の論者による質疑応答がついて、全部で100ページ近い。いや、読みますよ、読みますけどね、いろいろ辛いなあ。
1. イントロダクション
[…] というわけで、非確率標本の分析についてレビューし議論します。
2. 仮定と推論フレームワーク
目標母集団を\(U = \{1, 2, \ldots, N\}\)とする。調査変数を\(y_i\)、補助変数を\(\mathbf{x}_i\)とする。母平均\(\mu_y = \frac{1}{N} \sum_{i=1}^N y_i\)を推定したい。
非確率調査標本のデータセットを\(\{(y_i, \mathbf{x}_i), i\in S_A\}\)とする。そのサイズを\(n_A\)とする。
2.1 仮定
標本包含インジケータを\(R_i = I(i \in S_A)\)とする。$$ \pi^A_i = P(i \in S_A | \mathbf{x}_i, y_i) = P(R_i = 1 | \mathbf{x}, y_i) $$とする。欠損データ分析の用語に従いこれを傾向スコアと呼ぶが、参加確率という人もいる。
Chen, Li & Wu (2020 JASA) に従い、以下を仮定する。
- A1. \(R_i\)と\(y_i\)は\(\mathbf{x}_i\)の下で独立。すなわち \( (R_i \perp y_i) | \mathbf{x}_i\)。
- A2. 目標母集団の全ユニットが非ゼロの傾向スコアを持つ。すなわち \(\pi^A_i \gt 0, i = 1, \ldots, N\)。
- A3. \(R_1, R_2, \ldots, R_N\)は\(\mathbf{x}_1, \mathbf{x}_2, \ldots, \mathbf{x}_N\)の下で独立。
A1は欠損データ分析でいうところのMARに近い。A2は揉める話で、7章で検討する。A3はクラスタ抽出の時には破られる。
さて、理想的には母集団について補助変数が既知であるとよいのだが、もうちょっと現実的に次のように仮定する。
- A4. サイズ\(n_B\)の確率標本\(S_B\)が存在し、デザイン・ウェイトを\(d^B_i\)としてデータセット\(\{(\mathbf{x}_i, d^B_i), i \in S_B\}\)が手に入っている(\(y\)については未知でよい)。
以下では\(S_B\)を参照確率標本と呼ぶ。補助変数が\(S_A\)と\(S_B\)の両方で手に入っているという点がポイント。
2.2 推論フレームワーク
変動の源は3つある。
- \(S_A\)における包含・参加の傾向スコアのモデル \(q\)。
- アウトカムの回帰 \(y | \mathbf{x}\) ないし代入のモデル \(\xi\)。
- \(S_B\)の抽出デザイン\(p\)。
推論フレームワークは3つある。
- モデル・ベースの予測アプローチ。\(\xi p\)フレームワークとも呼べる。\(\xi\)と\(p\)の同時ランダム化の下でのフレームワーク。
- 推定された傾向スコアを用いた逆確率ウェイティング。\(q p\)フレームワークとも呼べる。\(q\)と\(p\)との同時ランダム化の下でのフレームワーク。
- 二重頑健推論。\(q p\)フレームワークないし\(\xi p\)フレームワーク(どちらとも決められない)。
[同時ランダム化と云っている意味がよくわかんないんだけど… たぶん、たとえばモデルベース推論ならば推定量の不確実性はモデル\(\xi\)の不確実性と\(S_B\)の抽出デザイン\(p\)に起因する不確実性を含む、という意味だろうと思う]
3. モデル・ベース予測アプローチ
このアプローチは、補助情報の量とモデルの信頼性がカギになる。補助情報がない場合は、モデルとしては$$ E_\xi(y_i) = \mu_0, \ \ V_\xi(y_i) = \sigma^2 $$ がリーズナブルだが、その下でのモデルベース予測推定量は$$ \hat{\mu}_y = \bar{y}_A = \frac{1}{n_A} \sum_{i \in S_A} y_i $$ なわけで、これは確かにモデルの下で不偏なんだけど、受け入れがたいだろう。
3.1 セミパラメトリック・アウトカム回帰モデル
以下では\(\mathbf{x}\)の先頭を1とする。
こういうセミパラメトリック・モデルが考えられる。$$ 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$$ 平均関数\(m(\cdot|\cdot)\), 分散関数\(v(\cdot)\)の形式は既知とする。
仮定A1.の下で、$$ E_\xi(y_i | \mathbf{x}_i, R_i = 1), \ \ V_\xi(y_i | \mathbf{x}_i, R_i = 1) = V_\xi(y_i | \mathbf{x}_i)$$ つまり、モデルが母集団で正しいなら\(S_A\)でも正しい。
\(\beta\)の真値を\(\beta_0\)として、その擬似最尤推定量は$$ S(\beta) = \sum_{i \in S_A} \frac{\partial m(\mathbf{x}_i, \beta)}{\partial \beta} \{ v(\mathbf{x}_i) \}^{-1} \{y_i – m(\mathbf{x}_i, \beta) \} = 0 $$ の解として与えられる。[個々の個体の平均関数からの偏差に重み付けて足し上げたのが0になるような\(\beta\)を求めればいいよって話ね。で、重みとは平均関数を\(\beta\)で微分してから分散関数で割ったものだ。うん、えーと… なんでこうなるんだっけ？ みたことあるような式なんだけど思い出せない]
3.2 予測推定量のふたつの一般形式
完全な補助情報\(\{\mathbf{x}_1, \ldots, \mathbf{x}_N\}\)があるとしよう。このときの推定量の形式は2つある。以下では\(m_i = m(\mathbf{x}_i, \beta_0), E_\xi(\mu_y) = \frac{1}{N} \sum_{i=1}^N m_i \)と略記する。また\(\hat{m}_i = m(\mathbf{x}, \hat{\beta})\)とする。
- 一般形式1。$$ \hat{\mu}_{y1} = \frac{1}{N} \sum_{i=1}^N \hat{m}_i $$ とする[つまり、予測モデルを使って母集団の全個体の\(y\)を予測して単純平均すりゃいいんじゃんってわけね]。もし線形回帰モデル \(m(\mathbf{x}, \beta) = \mathbf{x}^\top \beta\)だったら、\(\mu_x = \frac{1}{N} \sum_{i=1}^N \mathbf{x}_i\)として $$ \hat{\mu}_{yi} = \mu_x^\top \hat{\beta} $$ と書ける。
- 一般形式2。\(\mu_y = \frac{1}{N} \left(\sum_{i \in S_A} y_i + \sum_{i \notin S_A} y_i \right) \)であり、わかんないのは後半だけなんだから、それを\(\hat{m}_i\)で予測して$$ \hat{\mu}_{y2} = \frac{1}{N} \left( \sum_{i \in S_A} y_u – \sum_{i \in S_A} \hat{m}_i + \sum_{i=1}^N \hat{m}_i \right) $$ とする。線形回帰モデルだったら、標本平均\(\bar{y}_A, \bar{\mathbf{x}}_A\)を使って$$ \hat{\mu}_{y2} = \frac{n_A}{N} (\bar{y}_A – \bar{\mathbf{x}}^\top_A \hat{\beta}) \mu^\top_x \hat{\beta} $$ と書ける。
\(\hat{\beta}\)がOLS推定量だったら、残差の合計は0だから最後の式のカッコの中が消え、結局両者は同じとなる。[逆にいうと、線形モデル・OLS推定量じゃなかったら違ってくるわけだ]
完全な補助情報がなくても、参照確率標本を使って\(\sum_{i=1}^N \hat{m}_i\)を\(\sum_{i \in S_B} d^B_i \hat{m}_i\)で置き換え、\(\mu_x\)を\(\hat{\mu}_x = \frac{1}{\hat{N}_B} \sum_{i \in S_B} d^B_i \mathbf{x}_i, \ \hat{N}_B = \sum_{i \in S_B} d^B_i \)で置き換えればよい。
3.3 マス代入
伝統的には、非確率標本と参照確率標本を使った\(\mu_y\)のモデルベース予測推定量はマス代入推定量mass imputation estimatorと呼ばれてきた[なんて訳せばいいんだろうか…]。\(S_B\)の個体において\(y_i\)が全部欠損だったので代入して平均しました、と考えるわけだ。代入した値を\(y^*_i\)とすると、推定量は$$ \hat{y}_{y MI} = \frac{1}{\hat{N}_B} \sum_{i \in S_B} d^B_i y^*_i $$ とかける。MIってのは多重代入の略じゃないから気を付けてね。もし\(y_i^* = \mathbf{x}^\top \hat{\beta}\)だったら前の節と同じことになる。
なんでこんな話をしてるかというと、ひとたび欠損値代入だと思えばいろんな可能性が開けてくるからだ。たとえば、参照確率標本の個体を非確率標本の個体とマッチングしマッチした奴から\(y_i\)をもらってきちゃうとか(最近隣代入推定量という)。\(\xi\)がノンパラメトリック・モデルになっているといってもよい。
決定論的回帰によるマス代入に、本章冒頭のセミパラメトリック・モデルを使ってもよい。つまり、まず\(S_A\)から\(\hat{\beta}\)を得て(それはA1の下で一致推定量である)、次に\(S_B\)について\(y^*_i = m(\mathbf{x}, \hat{\beta})\)を得て、それらを重み付き集計して\(\hat{\mu}_{y MI}\)を得るわけだ。これはモデルベース予測推定量のふたつの一般形式のいっぽうと同一である。
[えー、同一？ 一般形式1は、母集団の全個体について予測して単純平均、一般形式2は、母集団の全個体について予測して単純合計し、標本部分の合計について実際の合計と予測された合計の差を求めてそれを足して調整し、母集団サイズで割る、だ。いっぽうここで説明しているのは、参照確率標本の全個体について予測し、その予測値を使って母平均をHajek推定するといってるんだから、一般形式1と漸近的に一致する、ってことなんじゃない？]
こうしてみると、仮定A1は、いわゆる「モデル移入可能性」を含意している。非確率標本で作ったモデルを参照確率標本での予測に使える、っていうことである。[なるほど]
マス代入アプローチでは、\(S_A\)の\(y_i\)は\(\beta\)の推定にしか用いられない。それでいいの？ という疑問が湧くだろう。Wu & Thompson (2020, 書籍)には「逆標本マッチング」推定量というのが出てくる。これは非確率標本を参照確率標本と最近隣マッチングして、マッチした相手\(j\)が持っている\(d^B_j\)を\(d^*_i\)として \(\hat{\mu}_{y A} = \frac{1}{\hat{N}^*} \sum_{i \in S_A} d^*_i y_i\) を求めるというものだ。その理論的特性についてはまだよくわかっていない。
[はあああ？ あなたなにいってんの？ 参照確率標本が単純無作為標本だったらどうすんのよ。\(d^*_i\)はすべて同一になり、どうマッチングしようが母平均の推定量は非確率標本の単純平均になるじゃない？ なにか読み間違えているのだろうか]
逆標本マッチングに際して、\(\hat{p}_i, \hat{p}_j\)のカーネル距離を\(K_{ij}\)として \(d^*_i \propto \sum_{j \in S_B} K_{ij} d_j\)とする、というアイデアもある[\(\hat{p}_i, \hat{p}_j\)ってなんのことだか説明がないけど、非確率標本の個体を彦星、参照確率標本の個体を織姫として、参照確率標本の抽出デザインにおける彦星の包含確率と、数多くの織姫たちの包含確率のカーネル距離を足し上げたのを重みにする、ということかな。うん、まあそれならわかる。メモは省略するけど、いろいろ研究がある由]
[ともあれ、マス代入アプローチについてはいろいろ理論的課題があっていろいろ研究されているよ、という話が1パラグラフ。省略]
[あれれ？ モデルベース推論の話ってこれで終わり？ MRPは出てこないのか。マルチレベル回帰モデルだっていろいろありうる\( m(\mathbf{x}_i, \beta) \) のひとつに過ぎないんだから、とりたてて言及する必要なし、ってことかなあ]
4. 傾向スコア・ベースのアプローチ
今度は非確率標本の傾向スコア\(\pi^A_i = P(R_i = 1|\mathbf{x}, y_i) \)に基づくアプローチ。理屈の上からは、母集団の要素はみな\(\pi^A_i\)を持っているわけだ。それを非確率標本の全個体について推定したい。
4.1 傾向スコアの推定
A1の下では\(\pi^A_i = \pi(\mathbf{x}_i)\)と書ける。しかし関数形はわかんない。仕方がないので、次の3つのパラメトリックな形式をよく使う。
- 逆ロジット関数 \(\pi^A_i = 1-\frac{1}{1+\exp(\mathbf{x}^\top_i \mathbf{\alpha})} \)
- 逆プロビット関数 \(\pi^A_i = \Phi (\mathbf{x}^\top_i \mathbf{\alpha}) \)
- 逆補対数対数関数 \(\pi^A_i = 1 – \exp(-\exp(\mathbf{x}^\top_i \mathbf{\alpha})) \)
ほかに、関数形を考えずにノンパラメトリックな方法を使うという手もある。
では、傾向スコアをどうやって推定するか。
- まずは擬似最尤法から。
母集団について補助変数が既知で、A3を仮定できるなら、傾向スコアモデルのパラメータ\(\mathbf{\alpha}\)の完全対数尤度関数は$$ l(\mathbf{\alpha}) = \log \left( \prod_{i=1}^N (\pi_i^A)^{R_i} (1-\pi^A_i)^{1-R_i} \right) = \sum_{i \in S_A} \log \left( \frac{\pi^A_i}{1-\pi^A_i} \right) + \sum_{i=1}^N \log(1-\pi^A_i) $$ となるわけで、これを最大化すればよい。補助情報が参照確率標本でしかわかってないなら、第二項を差し替えて$$ l^*(\mathbf{\alpha}) = \sum_{i \in S_A} \log \left( \frac{\pi^A_i}{1-\pi^A_i} \right) + \sum_{i \in S_B} d^B_i \log(1-\pi^A_i)$$これを最大化する\(\mathbf{\alpha}\)が疑似最尤推定量となる。
これを微分した\( \mathbf{U}(\alpha) = \partial l^*(\alpha) / \partial \alpha\)を擬似スコア関数という。一般に、擬似スコア関数 \(\mathbf{U}(\alpha)\)は真値\(\alpha_0\)において(qpランダム化の下での)期待値がゼロ、すなわち\(E_{qp}(\mathbf{U}(\alpha_0)) = \mathbf{0}\)である。これを、\(\hat{\alpha}\)が\(\alpha_0\)のqp一致推定量であるという。
というわけで、\(\mathbf{U}(\alpha) = 0\)を解けばよろしい。もし\(\pi^A_i\)の形式が逆ロジット関数ならば、$$ \mathbf{U}(\alpha) = \sum_{i \in S_A} \mathbf{x}_i – \sum_{i \in S_B} d^B \pi(\mathbf{x}_i, \alpha) \mathbf{x}_i $$ となる。 - [ここはちょっと話がそれている箇所だと思う]
Valliant & Dever (2011 Sociol.Methods.Res.)は、\(S_A\)と\(S_B\)をプールして傾向スコアを推定するという方法を提案した。重複をいとわずプールした標本を\(S_{AB}\)とする。\(R^*_i\)は\(S_A\)出身なら1, \(S_B\)出身なら0とする。ウェイト\(d_i\)として、\(S_A\)出身者には1、\(S_B\)出身者には\(d^B_i \left( 1-\frac{n_A}{\hat{N}_B} \right)\)を与える (\(S_{AB}\)を通じた合計を1にしている)。で、\(R^*_i\)のロジスティック回帰モデルをウェイトつきで推定する。
実はこの方法だと、\(S_A\)が単純無作為標本でないかぎり、傾向スコアモデルのパラメータの一致推定量にはならないことが示されている。Chen et al.(2020 JASA)をみよ。[へー]
この話の教訓は、標本をプールする路線は基本的に難しいということである。Valliantらは\(S_A\)出身者に\(d_i = 1\)を与えたが、これは要するに非確率標本の個体について交換可能性を仮定してしまっているわけだ。
Valliantらはくじけずに調整ロジスティック傾向スコア(ALP)というのを提案してきた(Wang, Vallient, & Li, 2021)。まずValliant & Dever と似た方法で、\(S_{AB}\)について調査ウェイト付きのロジスティックモデルをつくる。ここで\(d_i\)は\(S_A\)出身なら1, \(S_B\)出身なら\(d^B_i\)とする。次に、\(\hat{\pi}^A_i = \frac{\hat{p}_i}{1-\hat{p}_i} \)とする。[…中略…] しかしこれも概念的におかしくて […]。 - 推定方程式ベースの方法。
さきほどの\(\mathbf{U}(\alpha)\)を、一般推定方程式のシステムで置き換えることを考えよう。[さあ面妖な話になってきました…!!! 深呼吸！]
ユーザが指定した関数の、\(\alpha\)と同じ長さのベクトルを\(\mathbf{h}(\mathbf{x}, \alpha)\)とする。$$ \mathbf{G}(\alpha) = \sum_{i \in S_A} \mathbf{h}(\mathbf{x}_i, \alpha) – \sum_{i \in S_B} d^B_i \pi(\mathbf{x}_i, \alpha) \mathbf{h}(\mathbf{x}_i, \alpha)$$ とする。\(\mathbf{h}(\mathbf{x}, \alpha)\)がなんであれ、\(E_{qp}(\mathbf{G}(\alpha_0)) = \mathbf{0}\)である。
[なにをいっているのかというと、「それが0になるという方程式を解けばその解が推定量になる」ような関数を作りたいのである。この\(\mathbf{G}(\alpha)\)がそれだ。話を簡単にするために参照確率標本が単純無作為標本である場合について考えると、第2項は\(\sum_{i \in S_B} \pi^A \mathbf{h}(\mathbf{x}_i, \alpha)\)、つまり単純無作為標本から非確率標本と同じような抽出デザイン(?)で抽出しなおして合計すると云っていることになり、期待値は第1項の期待値と等しくなる。なるほど]
というわけで、\(\pi^A(\mathbf{x}, \alpha)\)の選択したパラメトリック形式とユーザがなんらか決めた\(\mathbf{h}(\mathbf{x}, \alpha)\)の下で\(\mathbf{G}(\alpha) = \mathbf{0}\)を解けば\(\alpha\)の推定量が得られる。
最尤推定量の最適性により、このやりかたはふつう効率性がよりひくい。しかしある種の場面では役に立つ。
たとえば、\(\mathbf{h}(\mathbf{x}, \alpha) = \frac{\mathbf{x}}{\pi(\mathbf{x}, a)}\)としよう。すると$$ \mathbf{G}(\alpha) = \sum_{i \in S_A} \frac{\mathbf{x}_i}{\pi(\mathbf{x}_i, \alpha)} – \sum_{i \in S_B} d^B_i \mathbf{x}_i $$ ところで、さっきの擬似最尤法の話で、\(\pi^A_i\)の関数形がもし逆ロジット関数だったら、擬似スコア関数は $$ \mathbf{U}(\alpha) = \sum_{i \in S_A} \mathbf{x}_i – \sum_{i \in S_B} d^B \pi(\mathbf{x}_i, \alpha) \mathbf{x}_i $$ であった。つまり、\(\mathbf{G}(\alpha)\)は傾向スコアモデルがロジスティック回帰であるときの擬似スコア関数を「歪めた」形になっている。さらに第二項をよくみると、これは\(\mathbf{x}\)の母合計の推定量であって、それがわかってんなら、参照確率標本がなくても傾向スコアを推定できるのである。
[なるほどーーーーー。たしかに、直観として、非確率標本の個体の包含確率推定値はどんなものであってほしいかというと、それを使ってなにかの変数の母合計のHT推定量を求めたとき既知の「正解」にぴたっと合うようなものであってほしいな、と思うもんね。この話はそれを正当化しているのか] - ノンパラメトリック手法。
仮に、有限母集団について\(\{(R_i, \mathbf{x}_i), i = 1, \ldots, N\}\)が手に入っているとしよう [台帳の全ユニットについて補助変数がわかっているとしよう]。\(\pi^A_i = \pi(\mathbf{x}_i)\)のNadaraya-Watsonカーネル回帰推定量というのを降臨させる。$$ \tilde{\pi} = \frac{ \sum_{j=1}^N K_h(\mathbf{x} – \mathbf{x}_j) R_j}{\sum_{j=1}^N K_h(\mathbf{x} – \mathbf{x}_j)} $$ ただし\(K_h(t) = K(t/h)\)はバンド幅\(h\)のカーネルである。
[はい深呼吸。心を広く持ち、学者どものわけのわからん発言を許容しましょう。なにをゆうとるのかというと、補助変数の2つのデータ点が似ているかどうかを返してくれる奴がいるとして、あるデータ点について、その点と標本包含されたすべてのデータ点との非類似性の和を求め、その最大が1となるように調整し、これを標本包含確率の推定値とする、という話である。標本包含されたデータ点との非類似性が高かったら、おまえなんか標本包含されねえよばーか、と推測するわけである]
実際には母集団の全個体について\(\mathbf{x}_i\)がわかるわけではない。でも、この式の分子は非確率標本からわかるし、分母は参照確率標本から推定できる。こうして、傾向スコア推定量 $$ \hat{\pi}^A_i = \hat{\pi}(\mathbf{x}_i) = \frac{\sum_{j \in S_A} K_h(\mathbf{x}_i – \mathbf{x}_j) }{\sum_{j \in S_B} d^B_j K_h(\mathbf{x}_i – \mathbf{x}_j) } $$ が得られる。同時\(qp\)フレームワークの下で一致推定量であることを示せる。シミュレーションでも、正規カーネル、ふつうのバンド幅でうまくいくことがわかっている。[へええええええ。Yuan, Ki, & Wu(2022 WorkingPaper)というのが挙げられている。これ、だれかパッケージ作ってないかしらん。自力でも意外に簡単に書けそう] - 回帰木ベースの手法。
Chu & Beaumont (2019 Conf.) のTrIPW法というのがある。非確率標本と参照確率標本をCARTみたいなアルゴリズムに食わせて学習させる。なお、プールした標本\(S_{AB}\)を単純に使うナイーブな方法はうまくいかないから注意するように。
4.2 逆確率ウェイティング
おめでとう、あなたは\(\hat{\pi}^A_i\)を手に入れました。では、\(\mu_y\)をどうやって推定するか。
逆確率ウェイテッド(IPW)推定量には2種類ある。$$ \hat{\mu}_{IPW1} = \frac{1}{N} \sum_{i \in S_A} \frac{y_i}{\hat{\pi}^A_i} $$ $$ \hat{\mu}_{IPW2} = \frac{1}{\hat{N}^A} \sum_{i \in S_A} \frac{y_i}{\hat{\pi}^A_i}, \ \ \hat{N}^A = \sum_{i \in S_A} \frac{1}{\hat{\pi}^A_i}$$ デザインベース推定理論では、1本目をHT推定量、2本目をHajek推定量という。後者のほうを使うべしといわれている(たとえ\(N\)が既知であっても)。
仮定A1, A2, そしてモデル\(\pi^A_i = \pi(\mathbf{x}_i, \alpha_0)\)が正しければ、\(\hat{\mu}_{IPW1}\)は一致推定量である。なぜならば… [メモ省略。これは非確率標本に限った話じゃないんじゃないかな]
4.3 二重頑健推定
IPW推定量の弱みは傾向スコアモデルの妥当性である。ところが、傾向スコアモデルの誤指定にも頑健な二重頑健(DR)推定量というのが提案され、成功を収めている。
傾向スコアモデル\(q\)の下で非確率標本について\(\pi^A_i\)があり、アウトカム回帰モデル\(\xi\)のもとで母集団について\(m_i = E_\xi(y_i |\mathbf{x}_i)\)があるとする。DR推定量は$$ \tilde{\mu}_{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$$ 第二項のモデルベース予測を第一項で調整しているわけである。DR推定量は\(q\)と\(\xi\)のどっちかの指定が正しければ不偏推定量である。
現実的には2つの方法がある。$$ \hat{\mu}_{DR1} = \frac{1}{N} \sum_{i \in S_A} \frac{y_i – \hat{m}_i}{\hat{\pi}^A_i} + \frac{1}{N} \sum_{i \in S_B} d^B_i \hat{m}_i$$ $$ \hat{\mu}_{DR2} = \frac{1}{\hat{N}^A} \sum_{i \in S_A} \frac{y_i – \hat{m}_i}{\hat{\pi}^A_i} + \frac{1}{\hat{N}^B} \sum_{i \in S_B} d^B_i \hat{m}_i$$ [あ、HT推定量とHajek推定量だね…] 後者を使うべし。
なお、もし\(S_A\)が単純無作為標本だったとしても、DR推定量はモデルベース推定量の一般形式2 $$ \hat{\mu}_{y2} = \frac{1}{N} \left( \sum_{i \in S_A} y_u – \sum_{i \in S_A} \hat{m}_i + \sum_{i=1}^N \hat{m}_i \right) $$ に帰着しないことに注意。[なるほど]
4.4 擬似経験尤度アプローチ
擬似経験尤度(PEL)法というのを紹介しよう。この20年ほどで研究が進んだ方法である。
[… なんかしらんが、尤度関数じゃなくてPEL関数というのを最大化するんだそうな。得られる推定量は\(\hat{\mu}_{IPW2}\)と同一になるらしい。へー、Hajek推定が正当化されるのか。どういうご利益があるのかよくわかんなかったんだけど、信頼区間とか仮説検定とかの際に助かるとか、そんな感じの話だと思う]
5. quota surveyと事後層別
[quota surveyってなんて訳せばいいんだろうか？ 定訳は割当法調査だと思うけど、そうは訳したくないんですよね、超わかりにくいと思う]
quotaサーベイというのは、あらかじめ\(n_A\)とその下位母集団別のサイズを決めておいて埋まるまで集めますという方法である。通常、クオータ以外には抽出のコントロールがなされない。
以下では、\(S_A\)はクオータサーベイ標本、\(\mathbf{x}\)はクオータの設定に使ったカテゴリカル変数の集合とする。クオータの数を\(K\)とする。クオータ\(k\)の標本を\(S_{Ak}\), 標本サイズと母集団サイズを\(n_k, N_k\)とする。
A1の下では\(\pi^A_i = n_k / N_k\)なので、$$ \hat{\mu}_{IPW2} = \frac{1}{\hat{N}^A} \sum_{k=1}^K \sum_{i \in S_{Ak}} \frac{y_i}{\hat{\pi}^A_i} = \sum_{k=1}^K \hat{W}_k \hat{y}_k $$ と書ける。結局、標準的な事後層別推定量と同じである。
この推定量が妥当な母集団推定値を提供する条件をまとめておこう。
- \(\mathbf{x}\)が任意調査への個体のの参加行動を特徴づけていること。
- それぞれの下位母集団では個体の包含が比較的にランダムであり、特定の集団が除外されていたりしないこと。
- 層ウェイト\(\hat{W}_k = \hat{N}_k / \hat{N}^A\)が手に入ること。
- ハードコア無回答者、つまり任意調査に絶対参加してくれないひとたちが、調査参加者と似た\(y\)を持っていること。
IPW推定量は小さな傾向スコアに弱いので、あえて非確率標本を傾向スコア順に並べて\(K\)層に等分して層別推定量を使うという手もある。参照確率標本の側ではどうやるかというと… [はっはっは。面白いけどメモ省略]。\(K\)はどのくらいがいいかというと、層の標本サイズが30以上になるくらいがいいんじゃないでしょうか。年配の方なら覚えているでしょうけど、かつての古き良き日々、\(n \geq 30\)だったら「標本サイズが大きい」と称したものでした。[ほんとにこう書いてある。ははは]
6. 分散推定
[残念ながら時間切れにつきまるごとパス。見出しをみると、マス代入推定量の分散推定、IPW推定量の分散推定、二重頑健推定量の分散推定、の順に話が進むらしい]
7. 仮定についてもう一度考える
7.1 仮定A1
A1、すなわち\(\pi^A_i = P(R_i = 1|\mathbf{x}_i, y_i) = P(R_i | \mathbf{x}_i) \)はもっとも致命的な仮定である。データ収集の前に、標本包含に関連するであろう潜在的要因について検討しておく必要がある。収集後は、\(S_A\)の経験分布関数(ないし積率)を\(S_B\)からの重み付きのそれらと比べるべし。\(y\)と似ている\(z\)を選び、共変量\(\mathbf{u}\)を選び、\(S_A\)から条件つきモデル\(z | \mathbf{u}\)をつくり\(S_B\)からもウェイトを使って\(z|\mathbf{u}\)を作り、両者が似ていることを確かめるとよい。もし違ってたら\(z\)が\(y\)の補足変数になるか、あるいは仮定が正しくない。
7.2 仮定A2
A2は確率標本で行くところのアンダーカバレッジの問題に似ている。この仮定も外見よりずっと厳しい。非確率標本の場合、そもそも包含インジケータ\(R\)が確率変数だというのがひとつの仮定にすぎないのである。
A2が破られる典型的なシナリオがふたつある。
- 確率論的アンダーカバレッジ。\(S_A\)は\(U_0 = \{i | i \in U and \pi^A_i \gt 0\}\)から選ばれており、\(U_0\)は目標母集団の\(U\)の確率標本だという場合。A2は破られない。
- 決定論的アンダーカバレッジ。ある種の特徴を持つ個体が標本に決して含まれない。この場合は簡単な修正方法はない。Chen(2020, 博論)をみよ。[おいおいおい！ えらいぶん投げ方だな]
7.3 仮定A3
A3はあんまし重要でない。これが破られてても擬似スコア関数\(\mathbf{U}(\alpha)\)は不偏である。ちょっと効率性が下がるだけ。
7.4 仮定A4
目標母集団についての確率標本を探してくるのは難しくないとしても、ほしい補助変数が入っているやつを探してくるのは難しかろう。これはよほどのお金持ちの悩みだけど、確率標本が複数あったらどれをどう結合するかという問題もある。実務的には以下を薦めたい。非確率標本における参加行動を特徴づける、ないし調査変数への予測力を持つ、重要な補助変数が手に入るかどうかをチェックすること。そのうえでの考慮事項は、まず非確率標本と共通する変数がたくさんあること。次に大きさ。最後に、データ収集モードが非確率標本と同じであること。
8. 結び
[2013年のAAPORタスクフォースの報告書の話…]
[研究小史…]
非確率標本による統計的推測は、「データ統合」、すなわち複数のデータソースから得たデータを結合するというより一般的なトピックの一部分である。またマルチフレーム調査については研究が進んでいて…[略。いくつか文献がreferされている。よさそうなのをメモしておくと、Wu(2004, CanadianJ.Stat.), Kim & Rao (2012 Biometrika), Lohr & Raghunathan(2017 Stat.Sci.), Thompson(2019 Int.Stat.Rev.), Kim & Tam (2021 Int.Stat.Rev.), Yang & Kimg (2020 JapaneseJ.Stat.DataSci.)]
本論文の本質的メッセージのひとつは、非確率標本の分析における妥当性と効率性の概念である。妥当性とは点推定量との一致性を指す。効率性とは点推定量の漸近分散によって測られる。妥当性と効率性についての議論は適切な推論フレームワークと厳密な手続きを必要とする。確率標本についての伝統的なデザインベースフレームワークもモデルベースフレームワークも非確率標本にはフィットしない。
云々。
——————-
いやー、枠組みが勉強になった。デザインベース推論とモデルベース推論という風に考えるより、(だいたいそれらに対応するけど)傾向スコア推定とアウトカム予測、そして二重頑健推定、という風に考えればいいのか。なるほど。
傾向スコア推定のための擬似最尤法と推定方程式ベース手法についての解説も大変勉強になった。要するに確率標本と非確率標本をマージして個体がどっちの出身かを当てるロジスティック回帰モデルを組めばいいんでしょ、という風にざっくり考えていたが、もっときちんと考えないといけないんだなあ。
実務的な教訓としては、傾向スコア推定もアウトカム予測もいろいろ方法があるけれど、結局のところ質の良い共変量がカギ、黄金の弾丸はない、ということなんでしょうね。
ところで、Rには関連しそうなパッケージとしてNonProbEstとnonprobsvyがある(どっちもCRAN Task Viewに載っている)。ほかにもあるかもしれない。ちょっと調べておこう。
質疑応答編は別のエントリでメモする。引き続き頑張るぞ！たぶん！