読了: Delovoye & Savje (2020) Horvitz-Thompson推定量の一致性

 標本抽出について調べていると、値を抽出確率の逆数で重みづけて合計して母集団サイズで割ると(Horvitz-Thompson推定量)、それは母平均の不偏推定値だよ、いっぽう母集団サイズじゃなくて重みの和で割ると(Hajek統計量)、それは不偏じゃない、でも一致推定量だよ、なんていう話が出てくる。で、不偏性の証明は出てくるけど、一致性の証明は出てこない。おまえらユーザには理解できんだろうから省略するよ、というニュアンスがひしひしと伝わってくる。なんだかイライラする。いま文化大革命とか起きたら、私は紅衛兵となってキャンパスに突入し統計家を群衆のもとに引き出して自己批判を迫るかもしれない。(すいません冗談です)
 まあご配慮は正しいんだけどさ。推定量の漸近的挙動なんて途方に暮れるじゃないですか。聞かされても絶対理解できるわけないじゃないですか。でもちょっと覗き見したいんですよね。背伸びして大人の世界に触れてみたい、というか。(いいおっさんが… 我ながらキモイ)

Delovoye, A., Savje, F. (2020) Consistency of the Horvitz-Thompson estimator under general sampling and experimental design. Journal of Statistical Planning and Inference, 207, 190-197.

 そんなこんなで、本文が短いので読んでみた。HT推定量の一致性の一般的条件を示すのだそうです。Google様いわく、被引用件数11件。
 さあ、怖いものみたさでゴー!

1. イントロダクション
 サイズ\(N\)の有限母集団における$$ \mu = \frac{1}{N} \sum_{i = 1}^N y_i$$ を、その下位集合\(\mathbf{S}\)の\(y_i\)を使って推定する、という問題を考える。\(\mathbf{S}\)がランダムならいわゆるHorvitz-Thompson推定量が使える。デザインのもとで、$$ \hat{\mu} = \frac{1}{N} \sum_{i \in \mathbf{S}} \frac{y_i}{\pi_i} $$ ただし\(\pi_i = Pr(i \in \mathbf{S})\)である。これ、開発したのはHorvitz & Thompson(1952)とNarain(1951)だけどね。[へー。後者の掲載誌はJ. Indian Soc. Agricultural Statist. だそうだ。それで発見が遅れたんだろうね…]

 HT推定量は標本抽出でも使うし、因果推論でも使う。後者の場合、\(y_i\)は潜在アウトカムである。たとえば、\(y_i(1)\)を処理下のアウトカム、\(\mu_1\)をその母平均、\(y_i(0)\)を統制下のアウトカム、\(\mu_0\)をその母平均として、ATEは\(\mu_1 – \mu_0\)である。\(\mu_1, \mu_0\)をHT推定量で推定するわけだ。

 本稿の目的はHT推定量の漸近特性の検討である。

2. HT推定量の一致性
2.1 抽出推論
 漸近レジームととして\(N \rightarrow \infty\)を考える。この考え方の欠点は、収束速度が標本サイズに関してではなくて母集団サイズに関して表現されるという点なんだけど、平均包含確率を\( \bar{\pi} = \frac{1}{N} \sum_{i=1}^N \pi_i\)として\(E[|\mathbf{S}|] = \bar{\pi} N\)だから両者はつながっている。\(\bar{\pi} \gt 0\)とする (\(\bar{\pi} \rightarrow 0\)は許す)。
 HT推定量は以下のように書き換えられる。包含インジケータを\(S_i = \mathbf{1}[i \in \mathbf{S}]\) とい、\(\pi_i = 0\)のとき\(w_i = 0\)、でなければ\(w_i = 1/\pi_i\)として[\(\pi_i = 0\)は許すのか。非確率標本も視野にあるわけね]、$$ \hat{\mu} = \frac{1}{N} \sum_{i=1}^N S_i w_i y_i$$ 話を手っ取り早くするために、\( \tilde{\pi}_i = \pi_i / \bar{\pi}, \tilde{w}_i = \bar{\pi} w_i\)と基準化する。

[さあ、いよいよやってきますよ、難しい話が…!]

条件1. 次のような\(p, q\)が存在する: \(p \gt 2, q \gt 1, pq \geq p + 2q\)であり、関心ある変数の\(p\)次積率と抽出ウェイトの\(q\)次積率が有界となる、すなわち$$ \left[ \frac{1}{N} \sum_{i=1}^N |y_i|^p \right]^{1/p} \leq k_y$$ $$ \left[ \frac{1}{N} \sum_{i=1}^N |\tilde{w}_i|^q \right]^{1/q} \leq k_\pi$$ となる。定数\(p, q, k_y, k_\pi\)は漸近系列を通じて固定とする。

 この積率条件は、関心ある変数とデザインについての制御を与えている。
 最初のパートは、すごく大きな\(y_i\)を持つユニットが母集団ではレアであることを保証している。一部のユニットで変数が無限に大きくなってもかまわないが、そういうユニットは母集団における消えゆく部分でなければならない。この条件が破られるのは、たとえば、母集団のある部分が、消失は十分に遅く、かつ\(y_i \rightarrow \infty\)となるのは十分に速いような場合である。
 二つ目のパートは、デザインがなんらかのユニットを抽出する頻度が非比例的に低いことがないことを保証している。つまり、包含確率がその平均から極端に逸脱していることがないことを要求している。この条件が破られるのは、\(\bar{\pi} \rightarrow 0\)の速度よりも十分に速く、母集団の十分に大きな部分が\(\pi_i \rightarrow 0\)であるような場合である。
 ふたつの積率条件は\(pq \geq p + 2q\)でつながっている。これは、関心ある変数についての制御がきつかったら、ウェイトについての制御は緩くてよい、ということを意味している。制御のきつさをどう配分するかは適用場面によるけれど、関心ある変数の4次積率が有界だと想定するのが一般的だし、直感的にわかりやすい。\(p=4\)なら\(q \geq 2\)なわけで、つまり\(\tilde{w}_i\)の2次積率が有界なわけだ。そういうデザインの例を挙げると、\(N^{0.5}\)以下の最大の整数個のユニットで\(\pi_i = 1/4N^{0.25}\), 残りすべてのユニットで\(\pi_i = 1/4\)、というデザインがあてはまる。

定義1. ユニット\(j\)が抽出された下でのユニット\(i\)の包含確率を基準化した \(\tilde{\pi}_{i|j} = Pr(i \in \mathbf{S} | j \in \mathbf{S}) / \bar{\pi} \) について考える。もし\(\pi_j = 0\)なら\(\tilde{\pi}_{i|j} = \tilde{\pi}_i\)とする。
 平均デザイン依存性を以下のように定義する。$$ D(r) = \left[ \frac{1}{N^2} \sum_{i = 1}^N \sum_{j \neq i} |\tilde{\pi}_{i|j} – \tilde{\pi}_i|^{1/r} \right]^r $$

補題1. 条件1の下で、下式が成り立つ。$$ Var(\hat{\mu}) \leq \frac{k^2_y k_\pi}{\bar{\pi}N} + k^2_y k_\pi D \left(1- \frac{1}{p} – \frac{1}{q} \right) $$
証明 [たぶん読んでもわからんのでパス]

 [補題1が示しているのはこういうことだと思う。確率標本か非確率標本かを問わず、\(y_i\)と一次包含確率が、条件1に示した(たぶんマイルドな)制約のもとにある、とだけ仮定しよう。そのときHT推定量の分散は上界を持つ。上界は、もし二次包含確率が一次包含確率の積なら(つまりユニットの抽出が互いに独立なら)、期待標本サイズの逆数に、謎の定数を掛けたものになる。つまり分散の上界は期待標本サイズに反比例して小さくなる。ユニットの抽出に依存性がある場合、その依存性のぶんだけ、上界が少し上にずれる。うーん、ここまでややこしい話になっているのは\(\pi_i = 0\)を許しているからだろうな。もし許さないならばもっと簡単に書けるはずだ]

 補題1の上界の第1項は、期待標本サイズ \(E[|\mathbf{S}|] = \bar{\pi} N\)との関係を表している。第2項を無視すれば、分散は標本サイズにおいてconventional linear rateで消失する[どういう意味だかわからんが、まあ雰囲気はわかるよ]。この項をみると、\(\bar{\pi}N \rightarrow \infty\)となるようなデザインでないかぎりこの推定量はconcentrateしない、ということがわかる。
 第2項は、ユニットのペアの間の抽出依存性をとらえている。\(Pr(i \in \mathbf{S}|j \in \mathbf{S}) – Pr(i \in \mathbf{S})\)は、ユニット\(j\)が抽出されたということを知った時、ユニット\(i\)が抽出される確率がどれだけ高くなるかという、我々の知識の利得を測っている。これを包含確率と同じように基準化すると\(\tilde{\pi}_{i|j} – \tilde{\pi}_i\)となり、これを全ペアで平均すると\(D(r)\)になる。\(r=1\)なら通常の算術平均だし、\(r \lt 1\)なら大きな依存性を強調する形で平均することになる。要するに、デザインによって導入された全体的な依存性を測っているわけだ。
 平均的デザイン依存性は低いほうが望ましい。なぜなら、分散の制御がより強くなるからだ。特に、後述する条件2で示すようにこの量が消失するとき、依存性は十分に弱くなり、有効標本サイズが名目標本サイズとともに成長することが保証される。
 慣用的なデザインでこの依存性がどのような速度で消失するか、ここで考えておこう。まずベルヌーイ抽出やポアソン抽出では、ユニットは独立に抽出されるから\(D(r) = 0\)である。SRSの場合は[…中略…] \(D(r) = O(1/\bar{\pi}N)\)である。クラスタ抽出の場合、[…中略…]。クラスタ数が固定されたクラスタ抽出のようなデザインでは、平均デザイン依存性が固定されるので、有効標本サイズは\(N\)とともに成長しない。

条件2.(デザイン依存性が弱いこと) \(D \left(1-\frac{1}{p}-\frac{1}{q} \right) \rightarrow 0\)

 補題1と条件2を合わせると、推定量の標本分布が集中することが保証される。ただし、どこに集中するかは制御されない。問題は、一部のユニットが決して観察されないかもしれないという点である。それらのユニットは母平均には影響するが推定量には影響しない。最悪の場合、すべてのユニットで\(\pi_i = 0\)なら\(Var(\hat{\mu}) = 0\)であり、\(\mu\)がなんであれ推定量は0である[ははは]。
 なにかの構造的仮定をおいて補外するのでもないかぎり、収束点を制御するためには、デザインによって母集団から除外されていしまう部分が\(\mu\)に影響を与えない程度に小さい、ということを保証するよりほかにない。

定義2. \(e_i = \mathbf{1}[\pi_i = 0]\)とする。除外ユニットが母集団に占める割合を\(\bar{e} = \frac{1}{N} \sum_{i=1}^N e_i\)とする。

補題2. 条件1の下で、以下が成り立つ。$$ |\mu – E[\hat{\mu}| \leq 2k_y \bar{e}^{1-1/p} $$
証明: [読んでもわからんだろうからパス]

 [補題2が示しているのはこういうことだと思う。引き続き、確率標本か非確率標本かは問わず、\(y_i\)と一次包含確率が条件1に示した(たぶんマイルドな)制約のもとにある、とだけ仮定する。そのときHT推定量のバイアスは上界を持つ。上界は、謎の定数 \(k_y\)と\(p\)、そして「一次包含確率が0の人の割合」\(\bar{e}\)で決まる。なお、確率抽出がそうであるように\(\pi_i = 0\)を許さないなら、\(\bar{e} = 0\)となるから、それだけでHT推定量は不偏推定量になる]

 ふたつの補題は、標本分布の位置と幅についての制御を提供している。これらを合わせると、母平均からの偏差についての制御が得られる。本論文の主たる結果に到達した。

命題1. 条件1のもとで、MSEの上界は以下となる。$$ E[(\hat{\mu} – \mu)^2] \leq 4k^2_y \bar{e}^{2-2/p} + \frac{k^2_y k_\pi}{\bar{\pi} N} + k^2_y k_\pi D \left( 1-\frac{1}{p} – \frac{1}{q} \right) $$
証明 MSEを二乗バイアスと分散に分解し、補題1と補題2の上界を適用する。

系1. 条件1と条件2に加えて \(\bar{e} \rightarrow 0\) かつ \(\bar{\pi} N \rightarrow \infty\)ならば、二次平均において[MSEでみたときに、ということだと思う]、HT推定量は母平均に収束する。

2.2 因果推論
 因果推論の文脈ではHT推定量をIPW推定量ということがある。この文脈では、標本が無限の超母集団からiid標本であるといえるような漸近レジームに焦点が当てられてきた。このレジームでは、デザインベースの観点が持つニュアンスの一部が失われてしまう[標本抽出の文脈では抽出がiidでないほうがむしろ普通だしね]。そのかわり、包含確率と推定量の効率性に焦点が当てられてきた。
 しかし例外もある。Aronow & Middleton (2013)は、因果的な量のHT推定量の有限標本特性についてデザインベース・フレームワークで検討している。以下の議論は彼らの大標本での結果を補足するものである。
 [以下、メモは省略するけど、ふたつのHT推定量の差についてそのMSEの上界を示している。MSEが悪いほうの群のデザインのMSEに支配されるのだそうだ]

3. 結語
3.1 先行する結果
 HT推定量の収束についての研究としては Robinson(1982)があって、[… 中略…] 本研究はその一般化と捉えらえる。別のアプローチとしてIsaki & Fuller(1982)というのもあるけど [… 中略…] デザインと関心ある関数の両方にある合成条件をかけていて推論しにくい[???]。我々の貢献は、これらの異なるアプローチのいいとこどりだといえる。

3.2 実務的含意
 本結果は主に、ある種の悪条件なセッティングに関連している。たとえば関心ある変数の分布の裾が重かったり、デザインが歪んでいたりする場合である。実務家はできるかぎりそういう状況を避けるべきだ。推定量は(たとえ依然として一致性があるとしても)効率が悪くなるだろう。実務家にとって本結果は、悪条件が避けられないとき、とくにデザインをあまり制御できない場合に有益であろう。

  • そういう場面のひとつとして、統計的考慮と実務的考慮の間にトレードオフがあるデザインが挙げられる。USセンサス局のAmerican Community Surveyは、年次のセンサスを月次で補う調査だが、アラスカの一部地域はアクセスが難しいので台帳から除外されている。
  • ロシアのサンクトペテルブルグでの研究で、注射ドラッグ・ユーザのHIV発生率を推定しているのがあるけれど、母集団のかなりの部分がホームレスなので接触できず、包含確率を自由に決められない。仕方がないので対象者駆動型の抽出デザインを使っている[スノーボール抽出ってことね]。他にドラッグユーザの知り合いがいないユーザは最初のほうにしか出てこないわけだ。
  • 因果推論でもこういうことはある。[…略…]
  • こうしたデザインのゆがみはフィールド実験ではよくある話で[…略…]
  • 自然実験ならなおさらそうで[…略。いちいち実例を示していて面白いんだけど、メモを取るのがめんどくさい]

 本稿の結果は、このようなタイプの研究において関心が持たれるだろう。
 特に、条件1と条件2で捉えられている、一致性の主な条件のあいだの密接なつながりに注目されたい。このつながりが示しているのは、関心ある変数の挙動が良いことがわかっていたら(たぶん一様な上界を持っているなという場合なら)、調査なり実験ありのデザインには余裕が生まれるので、実務家はほかの大事な事柄にもっと目を向けてよい、ということだ。この結果は、実務家がデザインを制御できない場面で、自然に生じたデザインがいったいどんなデザインだったら推論ができるのかを示している。

付録. ほかのデザインベース推定量の一致性
 Hajek推定量 $$ \hat{\mu}_{HA} = \frac{\sum_{i \in \mathbf{S}} \frac{y_i}{\pi_i} }{\sum_{i \in \mathbf{S}} \frac{1}{\pi_i} } $$ について。
命題A1. 条件1と条件2が成り立ち、\(\bar{e} \rightarrow 0\)かつ\(\bar{\pi}N \rightarrow \infty\) のとき、Hajek推定量は母平均の一致推定量である。
証明 \(h = \frac{1}{N} \sum_{i \in \mathbf{S}} \frac{1}{\pi_i} \)とすれば\( \hat{\mu}_{HA} = \hat{\mu} / h \)である。
 系1により、上述の条件下で\(\hat{\mu}\)は\(\mu\)に収束する。\(h\)は、すべてのユニットで\(y_i = 1\)であるような母集団に関するHT推定量であるから[おお、なるほどね]、系1により、1に収束する。
 連続写像の定理により証明を完了できる。[なんだかしらんがまあ雰囲気はわかるよ]

 差分推定量について。[略]
—————–
 あらためてメモしておくと、こういうことだ。
 目的変数の値を\(y_i\)、一次包含確率を基準化してつくったウェイトを\(\tilde{w_i}\)とする。ある\(p\)と\(q\)について、\( p \gt 2, q \gt 1, pq \geq p + 2q, \left[ \frac{1}{N} \sum_{i=1}^N |y_i|^p \right]^{1/p} \leq k_y, \left[ \frac{1}{N} \sum_{i=1}^N |\tilde{w}_i|^q \right]^{1/q} \leq k_\pi \)が成り立つとする。ユニット抽出のユニット間依存性を測る関数\(D(r)\)があり(\(r\)は大きな依存性を強調する程度を表すパラメータ)、母集団において一次包含確率が0であるユニットの割合を\(\bar{e}\)とする。抽出デザインを固定し、標本サイズと母集団サイズを大きくしていくにつれ、\(D \left(1-\frac{1}{p}-\frac{1}{q} \right)\)が0に近づき\(\bar{e}\)が0に近づくとしよう。このとき、HT推定量もHajek推定量も母平均に近づく。

 なるほどね。おおざっぱにいうと、確率抽出でないどころか、たとえ「絶対に標本にならない奴ら」が母集団にいるような抽出デザインであっても、目的変数があまり極端に動かない、ウェイトがあまり極端に動かない、ある人が標本が選ばれるかどうかと別の人が選ばれるかどうかとのあいだにあまり強い依存性がない、という条件の下では、HT推定量やHajek推定量は母平均の一致推定量となるわけだ。ここで一致推定量といっているのは、抽出デザインは標本割合を含めて一定で、標本サイズがぐんぐん大きくなり(ということは母集団サイズもぐんぐん大きくなり)、「絶対に標本にならない奴ら」の割合がぐんぐん0に近づくときに、推定量は母平均にぐんぐん近づいていくよ、夕陽に向かって走っていくよ、さあ僕と一緒にどこまでもいこうウラジミール、という意味である。(すいません冗談です)
 うん。そりゃまあそうだよな。結論だけみれば、いっちゃなんだがあたりまえな話のように聞こえるんだけど、でもこういうことをきちんと証明できるというのは凄いことだ。