Tille, Y. (2016) Unequal probability inverse sampling. Survey Methodology, 42(2), 283-295.
3年位前からずっとあれこれ思い悩んでいた問題があって、このたびその問題について自分なりにまとめた説明を作っていて、たまたま見つけた関連論文をなにかの足しになるかなと思って眺めていたところ、不意に、これじゃん！ と気が付いた。一見そうは見えないけれど、まさにこれこそが自分が抱えている問題ではないか。ぞわぞわっと鳥肌が立ち、慌ててきちんと読み始めた。
マーケティング・リサーチにおけるささやかな問題について悩んでいるはずだったのに、なんということでしょう、まさかの標本抽出法の論文である。できれば関わり合いを持ちたくない分野だ。
二段目がポアソン抽出の逆抽出であるような標本抽出法についての解説。カナダにJob Vacancy and Wage Surveyというのがあり、そこでの標本抽出を例に挙げている。著者は先日読んだ標本抽出についての論文の著者である。
1. 問題
[略]
2. 問題の定式化
\(N\)個の企業からなる母集団を\(U\)とする。ここから企業\(i\)をなんらかの抽出デザインで抽出する。項目数\(M\)の職種リスト\(L\)があり、そのうちその企業にある職種のリストを\(F_i\)(サイズ\(Mp_i\))、ない職種のリストを\(D_i\)とする。なんらかの抽出デザインで、存在する職種を\(r\)個までを抽出する。抽出に失敗した数を\(X_i\)とする[ここではわかりにくいが、3章で意味がわかる]。
企業\(i\), 職種\(k\)の平均賃金を\(y_{ik}\), 従業員数を\(z_{ij}\)とする。目的は母集団における職種別平均賃金$$ \bar{Y}_k = \frac{ \sum_{i \in U | k \in F_i} z_{ij} y_{ik}}{ \sum_{i \in U | k \in F_i} z_{ik}} $$ である。
\(U\)から包含確率\(\pi_{1i}\)で抽出された企業の標本を\(S_1\)とする。[面食らったんだけど、\(U\)のメンバー\(i\)に包含確率\(\pi_{1i}\)を与えるような抽出デザインによって得たなんらかのサイズの標本、ってことですね。具体的な抽出デザインにはコミットしていない。この書き方だと非確率標本の可能性さえあるわけだが、もちろんそういう話ではない]
企業\(i\)から包含確率\(\pi_{k|i}\)で抽出された職種の標本を\(S_i\)とする。[なるほどね。私はこの問題を、「個々の対象者に割付可能な実験条件の同定し、割付可能な条件のなかから確率的に選択する」という2段階の割付として捉えていたのだが、この人の言い方では、\(F_i\)から\(r\)個を確率抽出するという問題なのだ]
\(\bar{Y}_k\)の推定量として、次の「比」タイプの推定量を用いることができる。[式はややこしいが、よく見ると単なるHajek推定量なのでびびらなくてよろしい] $$ \hat{\bar{Y}} = \left( \sum_{i \in S_1 | k \in (S_i \cap F_i)} \frac{z_{ik} y_{ik}}{\pi_{1i} \pi_{k|i}} \right) / \left( \sum_{i \in S_1 | k \in (S_i \cap F_i)} \frac{z_{ik}}{\pi_{1i} \pi_{k|i}} \right) $$
問題は、逆抽出型のデザインだと包含確率がわからんということだ。
3. 復元SRS
[頭が大変こんがらがったのだが、この節はさっき定式化した問題ではなく、それを簡略化した問題について考えている]
企業\(i\)が\(L\)のうち持っている職種の割合を\(p_i\)とする。\(L\)からの復元SRSを、\(r\)個の職業が同定されるまで続けるとしよう。このとき\(X_i\)は負の二項分布に従う。すなわち\(X_i \sim NB(r, p_i)\)である。確率分布は $$ Pr(X_i = x_i) = C(r+x_i – 1, x_i) p^r_i (1-p_i)^{x_i} $$ さらに $$ E(X_i) = \frac{r(1-p_i)}{p_i} $$ $$ var(X_i) = \frac{r(1-p_i)}{p^2_i} $$ [原文では組み合わせの数を\(nCr\)という風に書いているけれど、読みにくいので\(C(n,r)\)と書いた。えーと、なにを云っているのかというと、\(X_i\)は「その企業にある職種を\(r\)回ドローするまでに、その企業にない職種をドローしてしまう回数」だから負の二項分布に従いますよね、という話である。しかし成功した\(r\)回のドローには重複があるわけで、つまりこの節で考えているのは本当の問題ではない]
リスト\(L\)の項目\(k\)が抽出された回数を\(A_{ik}\)とする。これは多項分布に従い、抽出回数を\(n\)として、$$ Pr(A_{ik} = a_{ik} ) = \frac{n!}{M^n} \prod_{k \in L} \frac{1}{a_{ik}!} $$ $$ \sum_{k \in L} a_{ik} = n$$ となる。[ああそうか… 確かに、職種ドローの結果ってカテゴリ数\(M\)の多項変数のベクトルだもんね。\(n\)回の試行で、\(A_{i1} = x_1, A_{i2} = x_2, …, A_{iM} = x_M\) となる同時確率は、$$ Pr(A_{ik} = x_k, k \in L) = \frac{n!}{x_1! \cdots x_k!} p_1^{x_1} \cdots p_M^{x_M} $$ どの\(p_k\)も\(1/M\)だから $$ = \frac{n!}{x_1! \cdots x_k!} \left( \frac{1}{M} \right)^{ x_1+\cdots+x_M } = \frac{n!}{M^n} \prod_{k=1}^M \frac{1}{x_k} $$ なるほどね、あなたのいうとおりだ]
この多項ベクトルが、母集団のある所与の部分においてある固定的なサイズに条件付けられているならば[\(L\)からドローし続けたんだけど、結果的に\(F_i\)から\(r\)個がドローされたならば、ってことね。\(r\)種類とはいっていない]、$$ Pr \left( A_{ik} = a_{ik}, k \in F_i | \sum_{k \in F_i} A_{ik} = r \right) $$ $$ = Pr \left( A_{ik} = a_{ik}, k \in F_i \ \mathrm{and} \ \sum_{k \in F_i} A_{ik}= r \right) / Pr \left(\sum_{k \in F_i} A_{ik} = r \right)$$ $$ = \left( \frac{n! (1-p_i)^{(n-r)}}{(n-r)! M^r} \prod_{k \in F_i} \frac{1}{a_{ik}!} \right) / \left( \frac{n!p^r_i(1-p_i)^{n-r}}{r!(n-r)!} \right) $$ [分子は\(F_i\)の各要素のドロー回数の同時確率なんだけど、\(n\)個のドローのうち\(n-r\)回は失敗しなければならないので、そこに\(\frac{(1-p_i)^{(n-r)}}{(n-r)!}\)が掛けられている。分母は、\(n\)個のドローのうち\(r\)個が\(F_i\)側に落ちる確率だから二項分布である。うんうん。で、分子と分母にある\( \frac{n!(1-p_i)^{(n-r)}}{(n-r)!} \)が消えるから…] $$ = r! \left( \frac{1}{Mp_i} \right)^r \prod_{k \in F_i} \frac{1}{a_{ik}!} $$ ただし $$ \sum_{k \in F_i} a_{ik} = r$$ ここからわかるように、母集団の一部について\(A_{ik}\)の合計が条件づけられているとしても、分布は多項分布のままであり、復元SRSデザインのままである。
[さて、我々は企業\(i\)がドローされた下での職種\(k\)の包含確率を知りたい。でもここでは復元SRSを考えておるので…]
\(r\)個の職種を得るまで復元SRSを続けた場合、$$ E(A_{ik} | X_i) = \begin{cases} \frac{r}{Mp_i} & \mathrm{if} \ k \in F_i \\ \frac{X_i}{M-Mp_i} & \mathrm{if} \ k \in D_i \end{cases} $$である。
[ここで躓いて一旦投げ出したのだけれど、一晩寝て読み直したらあっさりわかった。職種リストのうちその会社に存在する職種を青、存在しない職種を赤とする。職種数は\(M\)個、うち青職種は\(Mp_i\)個ある。ここから青職種が\(r\)個ドローされるまで復元SRSしたんだから、つまり各職種に等確率にチップを配り青職種に配られたチップ枚数が\(r\)になるまで続けたといっているのである。青職種に撒かれた\(r\)は定数で、撒いたチップ枚数\(n\)と赤職種に撒かれた枚数\(X_i\)は確率変数である。
で、話を簡単にするために、赤職種に撒かれた枚数\(X_i\)で条件づける。すると、ある青職種\(k\)に撒かれたチップ枚数の期待値は、当然ながら\(r/Mp_i\)である。ある赤職種\(k\)に撒かれたチップ枚数の期待値は、当然ながら\(X_i / (M – Mp_i) \)である。なあんだ、あたりまえじゃんか。\(X_i\)で条件づけているから逆抽出であることは忘れてよいというのがポイントである]
復元抽出の場合、求めるべき\(\pi_{k|i}\)は正確には包含確率ではなく、\(A_{ik}\)の期待値である。\(k \in L\)について$$ \pi_{k|i} = EE(A_{ik}|X_i) = \frac{r}{Mp_i}$$ [ここもわかんなくって困惑した。なんで期待値を二回取っているの？ \(k\)が赤なのか青なのかが確率的に決まるってこと？と思ったんだけど、そうすると確率\(p_i\)で赤、\(1-p_i\)で青なんだから、\(A_{ik}\)の期待値は\((r+X_i)/M\)である。あっれー？
と困惑していたのだが、4章を読んで、\(E[E(A_{ik}|X_i)] \)という意味だと分かった。\(k \in F_i\)ならば\(X_i\)がどうだろうが関係ない。\(k \in D_i\)ならば、$$ \sum_{x_i \in \{0, 1, \ldots\}} \frac{x_i}{M-Mp_i} Pr(X_i=x_i) = \sum_{x_i \in \{0, 1, \ldots\}} \frac{x_i}{M-Mp_i} \frac{(r + x_i – 1)!}{x_i!} p_i^r (1-p_i)^{x_i} $$ で、よくわかんないけど、これが結局\(\frac{r}{Mp_i}\)になるということなのではないかと思う。
よくわかってないのに文句をつけるのもなんですが、\(k \in D_i\) の場合なんてどうでもよくないですか？ いま\(\pi_{k|i}\)について考えているのはHajek推定量の式のなかに出てくるからであって、その際には\(k \in F_i\)な\(k\)だけが問題になるわけだしさ]
さて、我々は\(M, r, X_i\)を知っているが\(p_i\)を知らない。[ここで慌てて冒頭部を読み直した。これ、そういう話か！ てっきり企業\(i\)にまず\(L\)をみせて\(F_i\)と\(D_i\)に分けてもらうのかと思っていた。全部見せるのではなく一部分を見せるという話なのだ]
\(E(X_i)=X_i\)を解いて積率法で\(p_i\)を推定すると[つまり\(X_i\)の実現値が\(E(X_i)\)と等しいとみなせば、\( E(X_i) = \frac{r(1-p_i)}{p_i} \)より] $$ X_i = \frac{r(1-\hat{p}_i)}{\hat{p}_i} $$ したがって、ひとつめの推定量として $$ \hat{p}_{i1} = \frac{r}{X_i+r}$$ がある。最尤推定でも同じ推定量が得られる。しかしこの推定量にはバイアスがある。[よくわからんが分母に確率変数\(X_i\)が入っているからだろう。はなから逆数を推定すればいいのに]。\(r \geq 2\)ならば、\(p_i\)の不偏最小分散推定量は $$ \hat{p}_{i2} = \frac{r-1}{X_i + r -1} $$ である。しかし\(1/\hat{p}_{i1}\)は\(1/p_i\)の不偏推定量である。
というわけで、はなっから逆数を推定すると $$ \hat{(1/\pi)}_{k|i} = \begin{cases} \frac{M \hat{p}_{i2}}{r} & = \frac{M(r-1)}{r(X_i + r – 1)} & \mathrm{if} \ k \in F_i \\ \frac{M(1-\hat{p}_{i2})}{X_i} & = \frac{M}{X_i + r – 1} & \mathrm{if} \ k \in D_i \end{cases} $$ [式自体には納得だけど、疑問点が2つ。\(p_{i2}\)を使ってるのはなぜだろうか？ なぜ\(k \in D_i\)のときについても推定しているの、使い道なくない？]
残念ながら、復元SRSでは同じ職種が複数回選ばれてしまうし、\(Mp_i\)が小さいときには抽出が長くなるかもしれない。
4. 非復元SRS
というわけで、こんどは職種の非復元SRSを考える。歯を食いしばれ。
まず変わってくるのは、\(X_i\)がNBDではなくて、負の超幾何分布に従うという点である。\(NH \sim (M, r, Mp_i)\)と書く。
確率分布はこうなる。$$ Pr(X_i = x) = p(xl M, r, Mp_i) = \frac{C(x+r-1, x) C(M-x-r, Mp_i-r)}{C(M, Mp_i)}$$ $$ E(X_i) = \frac{Mr(1-p_i)}{Mp_i+1}$$ $$ var(X_i) = \frac{rM(1-p_i)(M+1)(Mp_i-r+1)}{(Mp_i+1)^2 (Mp_1+2)} $$ [なにも考えずに写経した。まあそういう日もあるよね]
今度は、\(A_{ik}\)は0か1であって、$$ Pr(A_{ik} = a_{ik}, k \in L) = C(M,n)^{-1}$$ ただし $$ \sum_{k \in L} a_{ik} = n$$ である。[そらそうだ。ドローされた職種がある組み合わせになる同時確率は、\(M\)個から\(n\)を取り出す組み合わせの数の逆数である]
\(F_i\)における個数で条件づけると、$$ Pr(A_{ik}, k \in F_i | \sum_{k \in F_i} A_{ik} = r) = C(Mp_i, r)^{-1}$$ ただし $$ \sum_{k \in F_i} a_{ik} = r$$ となる。分布は依然として非復元SRSのままである。[導出のプロセスはメモを省略した]
[途中の進め方は復元SRSの場合とな時なので省略して…]
というわけで、包含確率の逆数の推定量は以下となる。 [なあんだ。最終的には復元SRSと同じなのだ] $$ \hat{(1/\pi)}_{k|i} = \begin{cases} \frac{M \hat{p}_{i2}}{r} & = \frac{M(r-1)}{r(X_i + r – 1)} & \mathrm{if} \ k \in F_i \\ \frac{M(1-\hat{p}_{i2})}{X_i} & = \frac{M}{X_i + r – 1} & \mathrm{if} \ k \in D_i \end{cases} $$ このウェイトはMurthy(1957)がすでに示している。Selehi & Seber (2001, Australian & NZ J. Stat.)をみよ。
5. 不均一確率による復元抽出
これを不均一確率抽出に一般化する。復元抽出であればそんなに難しくない。
職種をドローする際の抽出確率\(p_{ik}\)を所与とする。\(\sum_{k \in L} p_{ik} = 1\)である。\(P_{i} = \sum_{k \in F_i} p_{ik}\)とする。[なるほど。SRSのときの\(p_i\)に相当するものだね]
\(X_i\)の分布は復元SRSと同じ。
\(A_{ik}\)の同時分布と、\( \sum_{k \in F_i} A_{ik} = r \)と条件づけたときの同時分布は、復元SRSのときとは違うけど、まあ似たような式になる。
これまでと同様、\(X_i\)で条件づけた\(A_{ik}\)の期待値を求めて、\(P_i\)を推定して… 結局、以下となる。$$ \hat{(1/\pi)}_{k|i} = \begin{cases} \frac{\hat{P}_{i2}}{r p_{ik}} & = \frac{r-1}{(X_i + r – 1)r p_{ik}} & \mathrm{if} \ k \in F_i \\ \frac{1-\hat{P}_{i2}}{X_i p_{ik}} & = \frac{1}{(X_i + r – 1)p_{ik}} & \mathrm{if} \ k \in D_i \end{cases} $$
6. 不均一確率による非復元抽出
いよいよ本題である。不均一な確率の下での非復元抽出。[よくわかんないけど、\(\pi_{k|i}\)が所与、ないし補足変数\(b_k\)が所与なのだろうと思う]
6.1 非復元の系列抽出
逆デザインを決める前に、まずはデザインを決めないといけない。以下のデザインが考えられる。
- Ohlsson(1995) いうところの系列ポアソン抽出。次の手順で抽出する。残念ながら固定された包含確率を厳密に満たすわけではないのだけれど、とても正確な近似になっている。
- 各職種について、\([0,1]\)の一様確率変数\(u_{ik}\)を生成する。
- \(u_{ik} / \pi_{k|i}\)が小さいものから\(n\)個を選ぶ。[各対象者について\(n\)職種をドローするってこと? それが\(F_i\)かどうかは問わないわけ?]
- Sampford(1962), Pathak(1964)の方法。ここでは、包含確率が厳密に満たされる解を紹介する。
- 厳密に正な補足変数\(b_k\)を用意し、それに比例していて、かつ和が\(n\)になるような包含確率\(\pi_{k|i}(n)\)を求める。私によるものを含めいくつかアルゴリズムがある。Rのsamplingパッケージのinclusionprobabilities()を使うとよい。
- この包含確率に従って\(n\)個を系列抽出する。どうやるかというと、職種の完全なリストからはじめて、ひとつづつ消していく。ステップ\(j (=1,\ldots,N)\)で職種\(k\)を削る確率を、\( 1- \frac{\pi_{k|i} (N-j)}{\pi_{k|i}(N-j+1)}\)とする。こうやると、逆向きにみたとき、最初の\(n\)個の要素は確率\(\pi_{i|k}(n)\)で標本に包含されていることになる。[へえー。どうでもいいけどそうなんですか]
6.2 確率不均等な逆デザイン
では逆デザインの定義である。[明確に書いてないけど、たぶん系列ポアソン抽出でなく、\(\pi_{k|i}(n)\)を求めるほうの方法を採用することが前提になっているのだと思う]
選ばれた職種が\(r\)個になるまで削除をつづけるとなると、\(X_i\)の確率分布も、条件付き包含確率\(E(A_{ik} | X_i)\)も推定できなくなる。しかし心配はいらない。5章の\(\hat{(1/\pi_{k|i}}\))の式のなかの\(p_{ik}\)を\( \frac{\pi_{k|i}(r + X_i)}{r + X_i} \)におきかえるだけでよい。
[つまり、ドロー回数\(n\)を横軸、包含確率\(\pi_{k|i}(n)\)を縦軸に取った曲線の、実際のドロー回数\(r + X_i\)における傾きを抽出確率とした復元抽出だと思ってよい、ということだ。なぜだろう？ わからない…]
7. 考察
かくして、復元・非復元の均等・不均等確率の抽出について問題が解けた。
企業の包含確率を決める際には、存在する職種が多そうな企業を高めにするのがよかろう。そうしないと、職種の包含確率が過度にばらついてしまう。
——————–
いやあ、勉強になった… なにより、問題の定式化が勉強になった。
整理しておこう。この論文は、二段抽出で、二次抽出単位が層別されていて、一次抽出単位ごとの二段目の抽出要素数に制約があり、かつ関心の対象が層別母平均であるときに、そのHorvitz-Thompson推定量なりHajek推定量なりに用いる重みをいかにして求めるか、という話である。
調査対象者に対する刺激割付の問題として考えると、刺激の数が\(M\)、うちその人に割付可能な刺激の数が\(Mp_{i}\)で、割付可能刺激のうち\(r\)個までの刺激を、なんらかの確率的な方法で割り付けた時(この論文の場合は、要はPPS抽出みたいな割付方法)、刺激\(k\)の割付可能対象者が\(k\)に割り付られる確率の逆数\(1/\pi_{ik}\)をどうやって推定するか、という話なのであった。
いやね、それね、私がこの3年ずーっと悩み続けている問題ですよ、先生。