Seber, G.A.F., Salehi, M.M. (2013) “Adaptive Sampling Design: Inference for Sparse and Clustered Populations.” Springer. Chapter 2. Adaptive Cluster Sampling.
適応的抽出についてのモノグラフ、全6章のうちの第2章。適応的クラスタ抽出についての章だが、この抽出デザインが主題というより、この抽出デザインを題材として適応的抽出の際の主要な推定量を導出する、という主旨だと思う。
苦手分野だけど、たったの12ページだ。がんばろう！ と気合をいれて…
本章からunitのことをユニットと訳す。単位だとわかりにくいし個体だと混乱するので(区画がunitで、変数はその区画にいる動物かなんかの個体数だという場合があるので)。
1章で説明されてたんだけど、適応的クラスタ抽出とは、母集団のユニット間で「近隣」がどうかが定義されていて(全ユニット間に無向グラフみたいなのがあって)、ユニットのなかにレアなやつがあるときに、レアなユニットを多く抽出する方法。まず一定数のノードをSRSで抽出し(初期標本)、各ノードについて、もしそれがレアなユニットだったらその近隣すべてを抽出する。これを繰り返すと、抽出されたノードがどんどん広がっていくが(マインスイーパーみたいなイメージですね)、あるところで止まる。このクラスタというかネットワークの一番外側にあるユニットはレアなユニットではない。なお、初期標本のユニットがたまたま近くにあったら結果的に同じクラスタに属してしまうこともある。
1. 不偏推定
1.1 表記
母平均と母分散を定義しときますね。$$ \mu = \frac{1}{N}\sum_{i=1}^N y_i $$ $$ \sigma^2 = \frac{1}{N-1} \sum_{i=1}^N (y-\mu)^2 $$ サイズ\(n_1\)の単純無作為標本を手に入れ、これに観察を適応的に付け加えて\(n\)までもっていくものとする。\(n\)は確率変数である。
1.2 インジケータ変数
ユニット\(j = 1, 2, \ldots, r\)のインジケータ変数を\(I_j\)とし、\(\pi_j =E[I_j]\)とする。\(var[I_j] = \pi_j(1-\pi_j)\)である。\(\pi_{jk} = P(I_j = 1 \cap I_k = 1)\)とする。\(cov[I_j, I_k] = \pi_{jk} – \pi_j \pi_k\) である。
定数\(z_j\)について$$ Z = \sum_{j=1}^r z_j I_j $$ という形をした\(Z_j\)に関心があるとき[母合計とかね]、$$ E[Z] = \sum_{j=1}^r z_j \pi_j $$ $$ var[Z] = \sum_{j=1}^r \sum_{k=1}^t z_j z_k (\pi_{jk} – \pi_j \pi_k) $$ $$ \hat{var}[Z] = \sum_{j=1}^r \sum_{k=1}^r z_j z_k I_j I_k \left( \frac{\pi_{jk} – \pi_j \pi_k}{\pi_{jk}} \right) $$ である。[導出は省略したが、まあここまでは大丈夫だな]
1.3 修正HT推定量
[たぶんここが山場である。深呼吸…]
適応的クラスタ抽出でユニット\(i\)を抽出したとき、そいつが属するネットワークから境界のユニットを抜いたネットワークを\(A_i\)、それに属する個体数を\(m_i\)とする。
[細かいことだけど… 適応的クラスタ抽出では抽出したユニットがネットワーク上でクラスタを形成するよね。あるクラスタから境界クラスタを取り除くという意図はわかる(逆抽出における最後のユニットが定義上必ずonであるようなもので、境界のユニットは定義上必ずoffである)。では、境界ユニットは「抽出した」ことになるのだろうか？ なるよね、実際のところ抽出したんだから。ってことは、境界ユニット\(i\)は\(A_i\)に属さないことになるわけ？ たぶんなるんだろうな、きっと自分自身しか属していないクラスタという扱いになるのだろう]
HT推定量を使うためには、最終的な標本にはいった各ユニットの選択確率を知る必要がある。残念ながら、初期標本にはいっていたユニットを除き、選択確率はわからない。境界ユニットの場合、他の非選択クラスタに属している可能性もある[あっそうか… 別の始点からめくっていってもめくられていた可能性があるわけか]。というわけで、初期標本に入っていなかった境界ユニットは無視する。[まじか。めくらなかったことにするわけ？ つまり「抽出しなかった」という扱い？ 母集団が小さくても無視する気だろうか。強気だなあ]
母集団に存在するネットワークの数を\(K\)とする[おそらく境界ユニットはひとりクラスタとみなすのだろう]。、\(k\)番目のネットワークが持つユニット数を\(x_k\)とする(ネットワーク\(k\)に属する個体の\(m_i\)は\(x_k\)となる)。ネットワーク\(k\)の選択確率を\(\alpha_k\)とする(属する個体の選択確率も\(\alpha_k\)である)。
では、Thompson(1990 JASA)が提案した推定量を紹介しよう。ネットワーク\(k\)の\(y\)の合計をy^*_k とする。標本におけるネットワーク数を\(\kappa\)とする。ネットワークに振るインジケータ変数を\(J_k\)とする。$$ \hat{\mu}_{HT} = \frac{1}{N} \sum_{i=1}^K \frac{y^*_k J_k}{\alpha_k} = \frac{1}{N} \sum_{i=1}^\kappa \frac{y^*_k}{\alpha_k} $$ [単なるHT推定量だ…]
ネットワーク\(k\)が選択されない確率を\(p_k\)とすると、$$ p_k = \frac{C(N-x_k, n_1)}{C(N, n_1)} $$ である。[ええと、すべてのありうる初期標本に占める、「そのネットワークに含まれないユニットから選んだ初期標本」の割合だ]
したがって、$$ E[\hat{\mu}_{HT}] = \frac{1}{N} \sum_{k=1}^K y^*_k = \frac{1}{N} \sum_{i=1}^N y_i = \mu$$ であり、\(\hat{\mu}_{HT}\)は不偏推定量である。[これは上の\(p_k\)の式とは無関係に、単に\(E[J_k] = \alpha_k\)から示しているのだと思う]
分散はどうなるか。[導出。 \(var[\hat{\mu}_{HT}], \hat{var}[\hat{\mu}_{HT}] \)が示されている。めんどくさいのでメモ省略。導出のプロセスは追わなかったが、そんなに難しくはなさそう]
[以上は非復元抽出を想定していたが、復元抽出でもこの推定量は不偏だという説明。略]
上記の話はすべて\(n_1\)に依存している。逆にいうと、\(n_1\)はパイロット調査で決定できない。
1.4 修正HH推定量
[普段使わないので忘れがちだが、Hansen-Hurwitz推定量というのはサイズ\(n\)の復元抽出標本からの推定量のことで、包含確率じゃなくて抽出確率\(p_i\)を使う。母合計の推定量は\(\frac{1}{n} \sum_{i=1}^n \frac{y_i}{p_i} \)となる]
もうひとつの推定量として、HH推定量を一般化した次の推定量がある。最終標本のユニット\(i\)がめくられた回数を\(f_i\)として、$$ \hat{\mu}_{HH} = \frac{1}{N} \sum_{i=1}^N y_i \frac{f_i}{E[f_i]} $$ とする。
[しばらく途方に暮れたけど不意に腑に落ちた。めくられた回数ね！ どういう話かというと、ネットワークのノードをまず\(n_1\)個めくり、赤かったノードの隣接ノードをめくり、赤かったノードの隣接ノードをめくり…というのを繰り返していくんだけど、もうめくっているノードをまためくりたくなることがあるわけで、それをカウントしておくのだ。いやまて、そんなら隣接に赤ノードを持つ赤ノードは必ず二回めくられることにならない？ いやそうじゃないな、エッジを逆戻りするのはたぶんナシだ。同じクラスタに2個初期ノードが落ちてた時には2回めくりたくなるよね、ということだろう]
\(f_i\)は超幾何分布\((N, m_i, n_1)\)に従うので\(E[f_i] = n_1 m_i / N\)である。
[復習しよう。超幾何分布(N,K,n)とは、K 個の成功状態をもつサイズNの母集団からn個の要素を非復元抽出したときに含まれている成功状態の数Xがkである確率の分布である。期待値はnk/Nである。さて\(f_i\)とは、サイズ\(N\)のネットワークに、俺ノード\(i\)が属するクラスタに属する仲間ノードが\(m_i\)個あり、さてそこから初期ノードを\(n_1\)個抽出した時、仲間ノードが抽出される個数だ(仲間が\(f_i\)個抽出されたら俺は\(f_i\)回めくられるから)。なるほどね、超幾何分布\((N, m_i, n_1)\)に従っている]
よって、$$ \hat{\mu}_{HH} = \frac{1}{n_1} \sum_{i=1}^N \frac{y_i f_i}{m_i} $$ $$ = \frac{1}{n_1} \sum_{i=1}^{n_1} \frac{1}{m_i} \sum_{j \in A_i} y_j $$ [待って待って、2本目の式では\(i\)の意味が変わっている。最初の総和記号の\(i\)は標本ユニットでなくて初期標本ユニットのインデクスだ。初期標本のユニットごとに、そいつの所属クラスタにおける\(y\)の平均をもとめ、それを単純平均しているわけだ。当然ながらダブルカウントもある]
\(A_i\)における\(m_i\)個の観察における平均を\(w_i\)、その平均を\(\bar{w}\)として$$ = \frac{1}{n_1} \sum_{i=1}^{n_1} w_i = \bar{w} $$ この分散とその不偏推定量は… [メモ省略]
\(\hat{\mu}_{HH}\)の別形式を与えよう。ネットワーク\(k\)において全ユニットの\(w_i\)は等しいので\(\bar{v}_k\)と書こう。ネットワークが出てくる回数[初期標本ユニットが落ちた回数ってことかな]を\(b_k\)として、$$ \hat{\mu}_{HH} = \frac{1}{n_1} \sum_{k=1}^\kappa b_k \bar{v}_k $$ […中略…] ネットワーク\(k\)における\(y\)の合計を\(y^*_k\)として $$ = \frac{1}{N} \sum_{k=1}^K y^*_k \frac{b_k}{E[b_k]} $$ ただし、\(b_k\)は超幾何分布\((N, x_k, n_1)\)に従うので\(E[b_k] = n_1 x_k / N\)である。
[最初の式と同じ形式に書き換えたわけだ。母集団ユニットでみてもいいしネットワークでみてもいいわけね]
[以上は非復元抽出を想定していたが、復元抽出でもこの推定量は不偏だという説明。略]
1.5 HT推定量とHH推定量の比較
どっちの推定量の分散もネットワーク内分散の影響を受けないが、SRSの標本平均は影響を受けるわけで、つまり、ネットワーク内の分散が高くなると、これらの推定量の効率性はSRSに比べて相対的に高くなる。[これは推定量というよりデザインの話だね]
効率性の点からはHTのほうが、計算しやすさの点ではHHのほうがよい。Salehi(2003 Environ.Eco.Stat.)を見よ。[ここでいう推定量の効率性って何のことだろうか]
1.6 慣用的手法と比べた効率性ゲイン
適応的クラスタ抽出の効率性はふつう高いけど、それは母集団の空間分布によるわけで、ゲインが保証されているわけではない。効率的なのは、ネットワーク内の分散が母分散に近いとき[?]、最終標本割合が初期標本割合に近いときなんだけど、後者が成り立つときはふつうネットワーク内分散がちいさいわけで、この2つの要因は対立しうる。
適応的クラスタ抽出には二つの問題がある。(1)効率性が事前にわからない。(2)最終標本サイズが事前にわからない。後者について、予算が吹っ飛ぶのを避けるには停止ルールを設けるという手がある(制約付き適応的クラスタ抽出という)。推定はブートストラッピングでやるという提案がある。Hanselman, et al.(2005 Fisheries Bulletin)をみよ。[シミュレーション研究。中略]
ほかに、ネットワークのサイズがすごく違うときにどうするかとか、最終標本サイズを制御するためにはどうするかといった提案がある。
1.7 いくつかの適用例
[生態学っぽい論文が列挙されている。メモ省略]
1.8 不完全検出可能性
適応的クラスタ抽出でレアな種を抽出する場合、検出可能性は不完全である。[初期標本でヒットしないかもという話だろうか]
この問題を解決するために…[面倒なので省略]
2. Hajek推定量
\(\mu\)の別の推定量を紹介しよう。HT推定量の\(N\)を、\(N\)を知っている知らないは別にして、\(N\)の推定量に置き換える。ネットワーク\(k\)のユニット数を\(x_k\)として$$ \hat{N} = \sum_{k=1}^\kappa \frac{x_k}{\alpha_k}$$ $$ \hat{\mu}_{HJ} = \frac{1}{\hat{N}} \sum_{i=1}^N \frac{y^*k J_k}{\alpha_k} = \frac{1}{\hat{N}} \sum_{k=1}^\kappa \frac{y^*_k}{\alpha_k} $$ 分散推定は、\(z_k = y^*_k – x_k \mu\)としてテイラー線形化を使って $$ var[\hat{\mu}_{HJ}] \approx \frac{1}{N^2} \sum_{r=1}^K \sum_{s=1}^K (\alpha_{rs} – \alpha_r \alpha_s) \frac{z_r}{\alpha_r} \frac{z_s}{\alpha_s} \equiv V_{HJ}$$ その近似的推定量は…[略]
この\(\hat{\mu}_{HJ}\)はHajek(1971)にはじまるといわれることが多い。不偏性は近似的にしか成り立たない。標本サイズ固定で\(\alpha_k\)が\(y\)に比例する場合はHT推定量のほうが効率的だが、標本サイズがランダムな時にはHJ推定量のほうが効率的かもしれない。Sarndal et al.(1992)をみよ。
[HT vs HJについての詳しい話、特にHJのバイアスの話が読みたかったんだけど、やっぱりSarndal本をreferして終わりか… 残念。]
3. 信頼区間
\(\hat{\mu}_{HT}, \hat{\mu}_{HH}\)の漸近挙動はどうなっているのか。漸近正規性の理論的条件は複雑なんだけど、要するに、HTは小さなクラスタがたくさんあるとき漸近正規であり、HHはもっと漸近正規になりやすく、特にクラスタ平均がわりかし同じ時に有利である。[…中略…]
3.1 ブートストラップ信頼区間
[略]
3.2 経験尤度信頼区間
ノンパラメトリック尤度比関数を使った経験尤度法で信頼区間を構成するという手もある。どういう話かというと、観察データを使って分布についての考察を制約し、既知のパラメトリックな形式を指定することなく、推定量の母集団分布についての限定的な仮定を置くのである。このアイデアはHartley & Rao (1968)がはじまりで、彼らはこれをscale-loadアプローチと呼んだ。経験尤度と呼んだのはOwen(1988)である。
Salehi et al.(2010 Environ.Eco.Stat.)では適応的クラスタ抽出について、HH推定量とHJ推定量それぞれについて擬似経験尤度関数を考えた。信頼区間が出せるよ。HHではブートストラップ信頼区間と同程度だった。HJではHTのブートストラップ信頼区間よりよかったんだけど、HJのブートストラップ信頼区間との比較はこれからだ。
4. 非復元選択されたネットワーク
[疲れてきたのでパス]
5. さらなる拡張
- 初期標本がSRSじゃなくて、たとえばPPS抽出だったらどうなるか。[そう！それ読んで思った！ すでに研究があるそうな。メモ省略]
- 初期標本で測定した値が小さい順にユニットをとっていったらどうなるか。[…中略…]
- インターネットみたいな隠れた母集団からの適応的抽出とか、スノーボール・サンプリングとか、ランダム・ウォーク法とか、ネットワーク・サンプリングとか。モデルベースの手法やベイジアンモデルが提案されている。[いくつか論文がreferされているが、やはりこの分野もS.K.Thompsonという人が有名らしい]
- 単変量の値じゃなくてベクトルが観察される場合。[…略]
—————-
いやー、修正Hansen-Hurwitz推定量の話がものすごく面白かった。日曜日に最適な大人の頭の体操である。期待していた話(Hajek推定量のバイアスの話)はなかったんだけど、面白かったから良しとしよう。
素朴な疑問なんだけど、適応的クラスタ抽出って、母集団全体のネットワークやサイズがわからん場合にもいけるんだろうか。最後の節に書いてあったインターネットとかの話がそれなのかな。