Fattorini, L. (2006) Applying the Horvitz-Thompson criterion in complex designs: A computer-intensive perspective for estimating inclusion probabilities. Biometrika, 93(2), 269-278.

　仕事の関連で調べものをしていて、適切なキーワードがわからず迷走していたんだけど、この論文のイントロ部分にあまり期待せず目を通し、探していたタイプの研究がついに目の前に現れたことに気が付いた。長かった。Google様いわく引用回数93。

　知らない言葉がたくさん出てくるので、イントロ部分はほぼ全訳する。
　なお、いきなり「1次・2次包含確率」という言葉が出てくるけど、包含はinclusionの訳。包含確率とは抽出確率とか選択確率というのと同じ意味だろう。本文を読み進めて気が付いたのだが、ここで1次・2次といっているのは多段抽出の1段目・2段目のことじゃじゃなくて、ある抽出単位の抽出確率が1次、ある抽出単位ペアの同時抽出確率が2次、ということらしい。そう思って読むと以下の記述も理解できる。わけわからんと思っても我慢して読み進めてみるものね…

　母合計のHorvitz-Thompson[以下HTと略記]推定量は不偏推定量で、その標本抽出分散は1次・2次包含確率に依存する。さらに、2次包含確率が正であることを保証するデザインであれば、標本抽出分散はHT分散推定量で不偏推定される。固定サイズデザインの場合はそのかわりにSen-Yates-Grundy統計量が用いられることもある。いくつかの2次包含確率がゼロである場合、ないし2次包含確率の計算ができない場合については、Berger(2004)がいくつかの選択肢を提案しているが、Hansen-Hurvitz[以下HHと略記]分散推定量がもっともよく用いられている。Wolter(1985)はHH推定量がほとんどの固定サイズデザインに対して保守的であることを示している。

　HT推定量の主な利点は、固定サイズデザインにおいて調査変数とほぼ比例的であるような補足サイズ変数に1次包含確率が比例しているときに有効性があるという点である。いわゆるπPSデザイン[自治体の抽出確率を自治体の人口に比例させるようなデザインのことだろう]のようないくつかの選択メカニズムでは、1次包含確率がサイズ変数に比例する。しかし、πPSデザインを計画する際には標本サイズはあらかじめ固定されていないといけない。いっぽう場合によっては、リサーチャーは調査時点でいくつ抽出できるのかわからないこともある。そのため事例によってはPPS系列抽出が好まれる。この手法では、それぞれの抽出において、これまで抽出されていない単位の選択確率が、抽出されていない単位の合計サイズに対するその単位のサイズの比率によって与えられる。このスキーマでは、単位は系列的に抽出され、標本サイズは調査期間が終わった時点で選ばれた単位数となる。

　残念なことに、たいていの系列抽出スキーマでは、1次・2次包含確率を陽に導出できない[あ、やっぱしそうなのか…]。そのためHT推定量は使えない。さらに、Raj(1956)推定量も非許容的であり[それよかましな推定量が存在するってことかな]、単位が標本に入ってくる順序に依存する。Rao-Blackwellのバージョンというのもあるけれど、さほど大標本でなくても計算的に無理である。

　しかし、母集団がwith-frameであり[有限サイズの抽出台帳が存在するってこと？]、抽出スキーマが未知の母集団特性に依存しないならば、標本抽出を独立に十分な回数反復し、単位ないし単位のペアが標本に入ってきた回数の割合に基づいて1次・2次包含確率を推定できるだろう。以下ではこれらの推定値を経験包含確率と呼ぶ。真の包含確率のかわりにこれらの推定値を使うことがHT統計量を求めることができる。
　本論文はHT推定量ならびにHT分散推定量・HH分散推定量における経験包含確率の使用について扱う。

　というわけで本題。
　母集団サイズを\( N \)とする。各単位が持つ調査変数(非負)の値を\(y_1, \ldots, y_N\)とする。母合計\(T\)に関心があるとする。非復元抽出スキーマを考える。1次包含確率を\( \pi_j > 0 (j=1,\ldots,N) \)、2次包含確率を\( \pi_{jh} (h > j=1,\ldots,N) \)とする。標本\(S\)は母集団インデクスのサブセットと見ることができる。
　1次包含確率の陽な導出ができないとき、HT推定量 $$ \hat{T} = \sum_{j \in S} \frac{y_j}{\pi_j} $$ は使えない。しかし抽出スキーマが\( y_1, \ldots, y_N \)とか他の道の母集団特性に依存していないなら、1次包含確率をモンテカルロ推定できるだろう。\(S\)を選んだ時と同じルールで台帳から選択するのを独立に繰り返して\(M\)個の標本\(S_1, \ldots, S_M\)を得る。単位\(j\)が\(M\)個の標本に入った回数を\(X_j\)として、\(\pi\)の常に正な推定量は $$ p_j = \frac{X_j + 1}{M +1} $$である。\(p_j\)は\(p_j\)の一致推定量である。HT推定量を $$ \hat{T}_M = \sum_{j \in S} \frac{y_j}{p_j} $$と修正しよう。\(M\)が増大につれ\( \hat{T}_M \)が \( \hat{T} \)にほとんど確実に収束することは自明である。実際、\( \hat{T}_M \)は\(M \rightarrow \infty\)に伴い漸近的に\( \hat{T} \)と等価であり、漸近的に不偏であり、そのMSEは\(\hat{T}\)の分散に収束するということを証明できる。付録を読め。

　\( \hat{T}_M \)の分散の推定について。2次包含確率が常に正かどうかはふつう抽出スキーマをみればわかる。
　常に正な場合、その推定量は$$ \pi_{jh} = \frac{X_{jh}}{M + 1} $$ で、$$ v^2_M = \sum_{j \in S} y^2_j \left( \frac{1}{p^2_j} – \frac{1}{p} \right) + 2 \sum_{j \in S} \sum_{h > j} y_j y_h \left( \frac{1}{p_j p_h} – \frac{1}{p_{jh}} \right) $$ が\(var(\hat{T}_M\))の漸近不偏推定量になる[うぉぉぉ…]。付録を読め。
　常に正ではない場合は、保守的な分散推定として…[ \(s^2_M\)というのが導入される。メモ省略] …となる。付録を読め。

　前述のように、\( \hat{T}_M\)は\(\hat{T}\)と漸近的に等価である。では、どれだけの\(M\)があれば、\( \hat{T}_M\)が\(\hat{T}\)とほぼ同程度の有効性を持ち、\( \hat{T}_M\)のバイアスが無視できるようになり、\( var(\hat{T}_M \)のふたつの推定量のバイアスが\(var(\hat{T}\)のバイアスと同程度になるのか。
　有効性の損失の指標を $$ L(M) = \frac{|MSE(\hat{T}_M) – var(\hat{T})|}{var(\hat{T}} $$としよう。すると下式が示せる。 \(\pi_0 = \min \pi_j, CV(\hat{T}) = \frac{\sqrt{var(\hat{T})}}{T}\)として$$ L(M) = \leq \frac{9}{(M+2) \pi_0} \left( 1+\frac{1}{(CV(\hat{T}))^2} \right) $$ [おおおおすげえええ。えーと、抽出確率が均一でない標本から母合計を推定する際、ウェイトを求める際にほんとの抽出確率ではなくてモンテカルロ推定した抽出確率を使うと推定量の分散がちょっと上がっちゃうんだけど、その程度は、モンテカルロ反復数が大きいときに小さくなり、真の抽出確率の最小値が大きいときに小さくなり、ほんとの抽出確率を使った推定量の分散が大きいときに小さくなる、ってことね。へえええ]
　さらに、\(\hat{T}_M\)の絶対相対バイアスは $$ ARB(\hat{T}_M) = \frac{|E(\hat{T}_M) – T|}{T} \leq \frac{1}{(M+2)\pi_0} $$ となる。
　\( v^2_M \) のARBは…[メモ省略]
　\( s^2_M \)の ARBは…[メモ省略]

　以上から次のことがわかる。有効性の損失を減らすためには計算的努力が必要。\(\hat{T}_M\)のバイアスは\(\hat{T}\)の有効性が高いときに大きくなり、包含確率が小さい単位が存在するときに大きくなる。場合によっては現実的に無理なモンテカルロ反復が必要になってしまう。

　\(\hat{T}\)のかわりに\(\hat{T}_M\)を使ってもいいよということを確認するための小規模なシミュレーションを行おう。
　(0,1)の一様分布からN=400の母集団を生成します。抽出単位のサイズを\(x_j = y_j + u_j, \ u_j \sim Un(-y_j/2, y_j/2) \)ってことにします(相関は0.86になる)。単位をサイズ順に並べます。で、n=20ないしn=40の標本を1000回抽出します。抽出方法は、非復元単純無作為、系統抽出、サイズで層を定義した層別抽出、Sampfordのrejective抽出[なにそれ…]。各標本についてのモンテカルロ反復は\(M = 10^6\)とします。
　結果。有効性とバイアスの損失は\(M=10^6\)もあれば無視できる。ただしn=20の系統抽出では有効性の損失が4%に及ぶ。[…中略…]

　空間抽出への適用。
　空間にちらばっている単位を抽出するとき、非復元単純無作為抽出をやるんだけど、空間自己相関が怖いので、ひとつ抽出するたびに残りの抽出単位のうちすでに抽出済みの奴に隣接している奴の抽出確率を下げる(上げる)というやり方を採ることがある[へええええ]。またサイズに合わせて抽出確率を決めることもある。
　そこで… [面白そうなんだけどいまちょっと時間がないのでパス]

　事例。
　2001年夏、北イタリアはトレンティーノの森で、Centro di Ecologica Alpinaが炭素インベントリを行った[よくわからんが、地中の炭素の蓄積を調べることらしい]。地域を1キロ四方の6200マスに分け、航空写真で森が3472マスあると突き止めた。知りたいのはこの地域の炭素の総量である。
　あるマスの炭素量を調べるだけで大変なので、とても全部は調べられない。そこで標本調査をやった。別の調査で、各マスの樹木量のおおまかな推定値が得られている。樹木量は炭素量と強く関連するのでこれを補足変数とする。夏の間に何マス調べられるかわからないので系列的に抽出する。2001年5月下旬から9月下旬にかけ、抽出確率を樹木量に比例させたPPS抽出で各地を回りました。結局調べたのは150マス。
　150マス調べたというのを所与として各マスの1次包含確率を求めようとすると、各マスについてそれが抽出されるC(3471, 149)パターンと150!の順序の数え上げが必要にある。これは無理なので、かわりに同じPPSスキーマでサイズ150を抽出するのを\(M=10^6\)回反復した。計算には140時間かかった。[ええええええ？そんなにかかったの？ 3741件から150件をサイズ比例抽出するのを100万回繰り返すだけでしょ？]

　…いやー、複雑な方法で抽出した標本にウェイティングをかけたいんだけど抽出確率がわかんないのでシミュレーションで殴ろうという発想は俺だけじゃないよな、絶対あるにちがいないと思って、あれこれ探してたんだけど、ようやく見つけることができた。Horvitz-Thompson estimatorというキーワードでヒットするとは。
　モンテカルロ反復数と推定量の性質との関係を、まさか解析的に導出しやがるとは思わなかった。さすがは統計学者だ。
　とはいえ、私の抱えている問題はちょっと違う問題なので (台帳から均一確率・非復元で系列抽出するんだけど、実験条件の割付が共変量で複雑に決まってて割付確率の算出が難しく、抽出しなかった行の共変量はわからない、というタイプの問題なので)、この論文の方法は直接には適用できない。誰かプロの統計学者に研究してもらえないだろうか… あいにくお金はないもんで、報酬は熱い抱擁とキスとかで…