Seber, G.A.F., Salehi, M.M. (2013) “Adaptive Sampling Design: Inference for Sparse and Clustered Populations.” Springer. Chapter 5. Inverse Sampling Methods.
ちょっと関心があって、適応的抽出デザインについてのモノグラフ(全6章)を読んできた(1章, 2章, 3章)。いよいよ、この本を読み始めた目的である第5章、逆抽出についての解説である。
なお、4章 “Primary and Secondary Units”はどうやら適応的クラスタ抽出が多段とか層別とかになっている場合の話らしい。適応的クラスタ抽出には差し当たって関心ないのでスキップ。6章 “Adaptive Allocation” は、おそらく層別抽出の標本サイズ配分を適応的にするという話だと思う。興味はあるけれど、差し迫った用事はないので今回はパスとする。
というわけで、この章を読み終えたら解散である。がんばるぞ！
1. イントロダクション
逆抽出を考えたのはHaldane(1945)だといわれている。一般的な定義としては「観察がある条件を満たすまで観察を求め続ける方法」だが、この定義だとたとえば固定コスト系列抽出(Pathak, 1976 Ann.Stat.)とか制約付き適応的クラスタ抽出なんかも逆抽出だということになってしまう[←どういう手法なのだろうか]。ここでは伝統的な逆抽出に焦点を当てる。すなわち、事前に決めた数の個人が観察されるまで抽出を続ける方法である。
2. Murthy推定量の新たな証明
1章でMurthy推定量というのを導入した。ここではSelehi & Seber (2001, Aust.NZ.J.Stat.)による新しい証明を示す。
標本の重複のないユニットの順序なしの集合を\(s_R = \{i_1, i_2, \ldots, i_\nu\}\)とする。[どういう話かというと、母集団のユニットに1から順に連番が振ってあり、これを背番号と呼ぶとして、そこから標本を抽出し、(復元抽出の場合は)重複を取り除き、含まれる背番号を順不同で集めたのが\(s_R\)である。\(\nu\)は必ずしも標本サイズでないことに注意]
最初のユニットとして\(i\)番目のユニットが選ばれたら1となるインジケータ変数を\(J_i\)とする。[母集団の全ユニットが持つ変数ね。\(i\)は背番号]
\(E[J_i] = p_i\)とする。[包含確率じゃなくて抽出確率だ。これも注意しないとな]
次の式は、trivialだが、\(\mu\)の不偏推定量である。[\(\mu\)ってのは\(y\)値の母平均でしょうね。\(N\)は母集団サイズ] $$ \hat{\mu} = \frac{1}{N} \sum_{i=1}^N \frac{y_i}{p_i} J_i$$
3章の定理1より、\(d_R = \{(i, y_i) : i \in s_R\}\)を値として持つ確率変数\(D_R\)は、\(\theta = (y_1, y_2, \ldots, y_N)^\top\)の十分統計量である。ということは、Rao-Blackwellの定理より、次の不偏推定量がつくれる。$$ \hat{\mu}_{RB} = E[\hat{\mu} | D_R] $$ $$ = \frac{1}{N} \sum_{i=1}^N \frac{y_i}{p_i} E[J_i | D_R] $$ $$ = \frac{1}{N} \sum_{i=1}^N \frac{y_i}{p_i} Pr(J_i = 1 | D_R) $$ [\(Pr(J_i = 1 | D_R) = \frac{Pr(J_1 = 1, D_R = s_R)}{Pr(D_R = s_R)} \)より] $$ = \frac{1}{N} \sum \frac{y_i}{p_i} \frac{Pr(J_i = 1, s_R)}{P(s_R)} $$ [\(\frac{Pr(J_i = 1, s_R)}{P(J_i = 1)} = Pr(s_R | J_i = 1)\)より $$ = \frac{1}{N} \sum_{i=1}^N \frac{P(s_R | J_i = 1)}{P(s_R)} y_i $$ だけど、もし\(i\)が\(s_R\)に含まれてなかったら分子は0だから] $$ = \frac{1}{N} \sum_{i=1}^\nu \frac{P(s_R | i)}{P(s_R)} y_i $$ これは1章で示したMurthy推定量である。
\(\frac{P(s_R | i)}{P(s_R)}\)は、\(s_R\)に含まれているユニットについて確率がわかっていればわかるのであって、\(\nu\)が確率変数かどうかは関係ない、というのが、この導出のミソである。[ああ、なるほどね。推定量の不偏性を導出する標本サイズを定数とみなしてしまうと、標本サイズが固定されているようなデザインでしか使えなくなるでしょ、という話だ]
この定式化でも、1章で示した分散公式とその不偏推定量が使える。Salehi & Seber(2001)をみよ。
3. 逆抽出デザイン
以下ではSelehi & Seber (2004 Aust.NZ.J.Stat.)に従い、母集団を\(C\)を満たすユニットの集合\(\mathcal{P}_M\)(サイズ\(M\), 平均\(\mu_M\), 分散\(\sigma^2_M\))と、満たさないユニットの集合\(\mathcal{P}_{N-M}\)に分ける。 [あああ、そうか、伝統的な逆抽出に焦点を当てるというのはこういうことか。私の抱えている問題の場合、よく考えると二段抽出になっていて、PSUの逐次抽出がSSUの適応的抽出の結果によって停止する。なのでこういう定式化はできない。がっくり…]
母集団から非復元のSRSをやって、\(\mathcal{P}_M\)から\(k\)個のユニットを得たら停止ということにする。最終標本サイズ\(n_1\)は確率変数になる。標本のユニットのインデクスの集合も分割して\(\mathcal{S}_M, \mathcal{S}_{N-M}\)とする。
\(\hat{\mu}_{RB}\)の\(\frac{P(s_R | i)}{P(s_R)}\)の部分を評価したい。どうするか。
[読んでて混乱したのだが、ここからの話は、むしろ$$ \hat{\mu}_{RB} = \frac{1}{N} \sum_{i=1}^N \frac{y_i}{p_i} \frac{Pr(J_i = 1, s_R)}{P(s_R)} $$ に基づいている。\(i\)が\(s_R\)に入ってなかったら分母は0なので、$$ \hat{\mu}_{RB} = \frac{1}{N} \sum_{i \in s_R} \frac{y_i}{p_i} \frac{Pr(J_i = 1, s_R)}{P(s_R)} $$ といってもよい。\(p_i = 1/N\)なので $$ \hat{\mu}_{RB} = \sum_{i \in s_R} y_i \frac{Pr(J_i = 1, s_R)}{P(s_R)} $$ といってもよい]
話を簡単にするために、\(C\)を満たす標本ユニットの集合を\(s_C\)として、インデクス\(i=1, \ldots, k\)を振りなおし、残りの標本ユニット集合を\(s_{\bar{C}}\)として、インデクス\(i=k+1, \ldots, n_1\)を振りなおす。
\(s_R\)は何通りあるか。最後に抽出したユニットは必ず\(C\)を満たしているから\(k\)通り。残りのユニットの順列は\((n_1-1)!\)通りある。ゆえに\(s_R\)は\(k \times (n_1 – 1)!\)通りある。
事象\(\{J_i = 1, s_R\}\)は何通りあるか。\(i\)が\(s_C\)に入っている場合、最後に抽出したユニットも\(C\)を満たしているから\(k-1\)通り。残りのユニットの順列は\((n_1-2)!\)通り。ゆえに\(\{J_i = 1, s_R\}\)は\((k-1)\times (n_1-2)!\)通りある。入っていない場合、最後に抽出したユニットは\(C\)を満たしているから\(k\)通り、残りのユニットの順列は\((n_1-2)!\)通り。ゆえに\(\{J_i = 1, s_R\}\)は\(k \times (n_1-2)!\)通りある。
[話を戻して、\(\hat{\mu}_{RB}\)の式の分子と分母をそれぞれ「通り」数で置き換えれば] $$ \hat{\mu}_{RB} = \sum_{i = 1}^k y_i \frac{ (k-1)\times (n_1 – 2)!}{k \times (n_1 – 1)!} + \sum_{i = k + 1}^{n_1} y_i \frac{ k \times (n_1 – 2)!}{k \times (n_1 – 1)!} $$ [分数のところを整理して$$ = \sum_{i = 1}^k y_i \frac{k-1}{k(n_1-1)} + \sum_{i = k + 1}^{n_1} y_i \frac{1}{n_1 – 1} $$ となり、] \(\bar{y}_M = \frac{1}{k} \sum_{i \in \mathcal{S}_M} y_i, \bar{y}_{N-M} = \frac{1}{n_1 – k} \sum_{i \in \mathcal{S}_{N-M}} y_i \)として $$ = \bar{y}_M \frac{k-1}{n_1-1} + \bar{y}_{N-M} \frac{n_1-k}{n_1-1} $$ \(\hat{P} = \frac{k-1}{n_1-1}\)として $$ = \hat{P} \bar{y}_M + (1-\hat{P}) \bar{y}_{N-M} $$ なお\(\hat{P}\)は\(M/N\)の不偏推定量である。
[さんざん苦労したけど、ある意味で当たり前のことがわかったわけだ。逆抽出標本に基づく母平均の不偏推定量は、条件\(C\)に合致したユニットの標本平均と、しなかったユニットの標本平均を、合致率の不偏推定量で重みづけた値になる。合致率の不偏推定量が\(k/n_1\)ではなくて\((k-1)/(n_1-1)\)だというところが注意点である]
分散の不偏推定量は…[1頁にわたる説明。スキップしたけどそれほど難しそうではない]。なお、この式はGreco & Naddeo (2007)が超幾何分布を使って導出した式と同じである。
\(\hat{P}\)の分散推定量は…[パス]
[…中略…]
母集団の下位部分についてサイズがわかっているという場合もある。そういう場合の母平均のHT推定量や、事後層別した層別母平均の不偏推定量も提案されている。[…後略…]
4. 一般的逆抽出デザイン
停止ルールにもいろいろある。次の3通りについて考えよう。
- 1ユニットづつ抽出していき、標本のなかのレア・ユニットの数が\(k\)になったら停止。いわゆるふつうの逆抽出である。
- まずサイズ\(n_0\)のSRS標本をがばっとととる。そのなかにレア・ユニットが\(k\)個あったらそこで停止。足りなかったら追加していき、\(k\)個になったら停止。
- まずサイズ\(n_0\)のSRS標本をがばっとととる。そのなかにレア・ユニットがひとつでもあったらそこで停止。足りなかったら追加していき、\(k\)個になったら停止。[ああ、なるほどね。レアユニットの出現率にしか関心がないんならそれでいいような気がするね]
(1)(2)の不偏推定量と分散はChristman & Lan (2001, Biometrics), (3)の不偏推定量と分散はSelehi & Seber(2004)を参照。
さらにSelehi & Seber(2004)では、(2)を変形して、金か時間が尽きたときも停めるというデザインについて考え、これを一般逆抽出(GIS)と呼んでいる。(1)(2)はGISの特殊ケースなわけである。
GISデザインのもとでの母平均の推定量について考えよう。最終標本サイズを\(n_1\), 上限を\(n_2\)とする。Murthy推定量は…[めんどくさいのでメモは省くけど、要は、\(n_0\)で止まった、\(n_1\)で止まった、\(n_2\)でようやく止まった、の3つのシナリオに場合分けするだけである。なあんだ]。分散は…[パス]
Moradi et al.(2007 J.Prob.Stat.Sci.)は、分母がゼロになりうるような推定量で比を推定するという問題について検討している。逆抽出デザインが自然ですよね。彼らはテイラー線形化を使って、GISの下での比の漸近不偏推定量とその分散を求めている。イラン政府がクルド人地区のはちみつ生産量を推定しようとして…[面白いエピソードだけど、力尽きつつあるのでパス]。そこで彼らはGISの下での回帰推定量をつくっていて…[略]
5. 多重逆抽出
母集団をサイズ既知の下位母集団に分割できるんだけど、しかしユニットがどの下位母集団に属しているかは抽出してみないとわかんない、というような場面は結構ある。事後層別だと事後層のなかにユニットが少なすぎるのが出てくるかもしれない。そこでChang et al.(1998 J.Stat.Planning&Inference)は多重逆抽出(MIS)というデザインを考えている。Salehi & Chang (2005 J.Stat.Planning&Inference) は母集団全体ならびに層ごとの\(y\)の合計のMurthy推定量を提案している。
[何度か読み返してようやく意味がわかった。なるほどね、これは興味深い話だ。GISもそうだけど、実務でいうところのブースト標本の話である]
5.1 クオータ抽出
クオータ抽出は実務家の誤用で悪名高いが、有限母集団からのMISの一種とみなせる。ただし、クオータ抽出は主観的に誤用されてきたが、MISには適切な停止ルールと推定量がある。
まずSRS標本を得て、いくつかの下位母集団が少なすぎたら、そいつらが目標に達するまでひとつづつ取り続ける、としよう。実務家たちはこれをSRS標本とみてしまうんだけど、これは間違い。
Sarraf-Zadegan et al.(2003)を心血管疾患の研究を例にあげよう。彼らはイランの都市イスファハンを93個のPSUに分割し(各PSUにはおよそ1000世帯はいっている)、25個のPSUをSRSでとった。で、各PSUから5-10%の世帯をとり、各世帯から大人をひとりランダムにとり、年代別に標本サイズを調べ、事前に設定した目標に足りなかったら世帯を増やした。彼はこれをクオータ抽出と云っているけれど、これはMISである。[なるほどね]
推定量についてはSalehi & Chang (2005)をみてほしいが、以下で理論的な説明をお届けしよう。
5.2 トランケートのある多重逆抽出
母集団\(\mathcal{P}\)が\(L\)個の下位母集団\(P_h\)にわかれている。非復元標本をとり、すべての\(h\)で標本サイズ\(n_h\)を目標\(m_h\)より大きくしたい。そこで、まず\(n_0\)のSRS標本をとって、一個ずつ増やしていき、すべての\(h\)で目標に到達するか、標本全体のサイズが\(n_T\)になるかまで続ける。
各\(h\)の標本のうちサイズが\(m_h\)個である層[目標に到達した層]をまとめて\(\mathcal{S}_M\)、残りを\(\mathcal{S}_\bar{M}\)とする。\(\mathcal{S}_M\)のcardinalityを\(k\)とし\(\hat{p} = (k-1)/(n-1)\)とする。[cardinality??? なぜこんな難しい言葉を使うのかよくわからない。目標に到達した層の標本ユニット数のことではないかと思うのだが…]
以下とする。
- \( [Q_1] = \{ n = n_0\} \)
- \( [Q_2] = \{ n_0 \lt n \lt n_T\} \ \mathrm{or} \ \{n = n_T \ \mathrm{and} \ n_h \geq m_h \ \mathrm{for \ all} \ h \} \)
- \( [Q_3] = \{ n = n_T \ \mathrm{and} \ m_h \lt m_h \ \mathrm{for \ at \ least \ one} \ h \} \)
[初期SRS標本で目標に到達した[Q1], 追加標本で目標に到達したか、上限に到達した瞬間に目標に到達していたら[Q2], 目標に到達してなかったら[Q3]だ]
すると推定量は$$ \hat{\mu} = \begin{cases} \sum_{h=1}^L \frac{n_h}{n_0} \bar{y}_h & \mathrm{if} \ [Q_1] \\ \sum_{h \in \mathcal{S_M}} \frac{n_h \hat{p}}{k} \bar{y}_h + \sum_{h \in \mathcal{S}_{\bar{M}}} \frac{n_h(1-\hat{p})}{n-k} \bar{y}_h & \mathrm{if} [Q_2] \\ \sum_{h=1}^L \frac{n_h}{n_T} \bar{y}_h & \mathrm{if} \ [Q_3] \end{cases} $$
[どんなにいかつい話かと思ったが、意外にふつうの話であった。はい深呼吸。[Q1][Q3]の場合は話は簡単で、層別の標本平均を、層の割合\(n_h/n\)で重みづけする。問題は[Q2]のケースで、この場合は、目標に到達した層と到達しなかった層を分ける必要がある。目標に到達した層については、普通に考えたらその層の標本平均を\(n_h/n\)で割ったのが重みだが、そうではなくて、目標に到達した層たちの標本ユニット数\(k\)で割り(つまり目標に到達した層たちの中での割合にし)、さらに\(\hat{p} = (k-1)/(n-1)\)、つまり目標に到達した層たちが母集団に占める割合の推定値を掛ける。目標に到達しなかった層についても同様で、目標に到達しなかった層たちの標本ユニット数\(n-k\)で割り、目標に到達しなかった層たちが母集団に占める割合の推定値\(1-\hat{p}\)を掛ける。ポイントはちゃんと場合分けすること、そして逆抽出の停止ルールで停止した時は到達した層としなかった層をわけることだな]
分散推定量は…[略]
下位母集団\(h\)の母平均の推定量は、層\(h\)の標本集合を\(S_h\)として、$$ \hat{\mu}_h = \begin{cases} \frac{n_h}{n_0} \bar{y}_h & \mathrm{if} \ [Q_1] \\ \frac{n_h \hat{p}}{k} \bar{y}_h & \mathrm{if} \ [Q_2] \ \mathrm{and} \ S_h \in \mathcal{S}_M \\ \frac{n_h (1-\hat{p})}{n-k} \bar{y}_h & \mathrm{if} \ [Q_2] \ \mathrm{and} \ S_h \in \mathcal{S}_\bar{M} \\ \frac{n_h}{n_T} \bar{y}_h & \mathrm{if} \ [Q_3] \end{cases} $$ [あっれーーー？ \(\bar{y}_h\)は層\(h\)の標本平均、\(\mu_h\)は層\(h\)の母平均だと思うんだけど？？ [Q1]のケースでなぜ、 \(\bar{y}_h\)に層\(h\)が標本において占める割合\(n_h / n_0\)を掛けているの？ なにか勘違いしているのだろうか…
元論文のSelehi & Chang (2005)を見ると、層別の母合計の推定の話をしていて、[Q1]では\( N \frac{n_h}{n_0} \bar{y}_h \)となっている。\(N_h\)の推定量は\(N \frac{n_h}{n_0}\)だからこれは理解できる。でも、層別の母平均の推定量は\(\bar{y}_h\)そのものではなかろうか。もとの式の\(N\)をうっかり機械的に取っ払っちゃっただけなんじゃないだろうか。
いやいや、ちょっとまて。そもそも、どんな停止ルールで抽出を止めたとしても、調査変数の値で止めたわけではないのだから、各層の母平均の不偏推定量は\(\bar{y}_h\)でよくない？ 私、なんか勘違いしているのだろうか…]
分散推定量は…[略]
—————-
やれやれ、読み終えたぞ。
クオータ抽出のようなヤクザな抽出でも、ちゃんとした定式化に従って逆抽出していれば不偏推定できる、という点が勉強になった。しかも、推定量は単に抽出停止シナリオ別に場合分けしているだけで、個々の式はそれほど複雑でない。逆抽出といってもそんなに怯えることはないんですね。勉強になりましたです。