Godambe, V.P, & Thompson, M.E. (1986) Parameters of Superpopulation and Survey Population: Their Relationships and Estimation. International Statistical Review, 53(2), 127-138.
調査データの分析について調べていると、母平均のHajek推定量 (標本の個々の値をその個体の標本包含確率の逆数で重みづけて平均した量)は母平均のデザイン不偏推定量ではないけれど、実はある種のモデルのもとでモデル不偏推定量なのだ… という話がでてきて、そこでこの論文がよく引用されている。
第一著者のGodambeさんという人は推定方程式アプローチの理論で有名な人らしい。第二著者はHorvitz-Thompson推定量のThompsonとは別人の模様。
こういうのは、もう、一行づつ読んでいくしかない。あきらかに私の理解能力を超える論文だが、人生においては、敗れるとわかっている戦いに挑まなければならないこともあるのだ…
1. イントロダクション
本論文の概要
標本調査で、超母集団モデルのパラメータなり、そのモデルに基づく母集団特性関数なりを推定したい場面を考える。たとえば、超母集団モデルが$$ y_i = \beta x_i + \epsilon_i$$ で、\(\epsilon_i\)は平均\(E(\epsilon_i) = 0\), 分散\(V(\epsilon_i) = \sigma^2\)に独立に従うとする[原文では、Eではなくて花文字のE, Vではなくてvが使われているが、わかりにくいので勝手に変える]。目的は標本から\(\beta\)を推定することだったり、有限母集団からの\(\beta\)のOLS推定値\(\beta_N\)を推定することだったりする。
\(\beta\)みたいなパラメータを推定する最適な手続きを考える際には、標本に基づく推定量\(\hat{\beta}_s\)が、所与の抽出デザインのもとで\(\beta_N\)の(ほぼ)不偏な推定量になっているか、が問題になる。抽出デザインとモデルの下で、デザイン不偏でありかつ同時MSEが最小である推定量のことを最適と呼ぶ。
本論文では推定関数の理論を用いて、最適性基準の応用範囲を拡張する。
超母集団モデルと推定関数
長さ\(N\)の母集団ベクトルを\(\mathbf{y}\)とする(要素\(y_i\)はベクトルかもしれない)。\(\mathbf{y}\)は分布\(\xi\)から生成されていて、\(\xi\)はクラス\(C = \{\xi\}\)のメンバーであることがわかっているとする。この\(C\)のことを超母集団モデルという。
\(\theta\)が超母集団パラメータだとしよう。推定関数の理論では、有限母集団ベクトル\(\mathbf{y}\)に基づく関数\(g(\mathbf{y}, \theta)\)が以下を満たすとき、それを不偏推定関数という。すべての\(\xi \in C\)について$$ E_\xi[g(\mathbf{y}, \theta(\xi))] = 0 $$ ただし\(E_\xi\)は\(\xi\)のもとでの期待値を表す。
[はい深呼吸！ そんなに難しいことはいってない。推定関数とは、それが0になるという方程式を解けばその解が推定量になるような関数のことだ。最尤法では密度関数\(p(\mathbf{y}, \theta)\)を最大化する\(\theta\)を求めるけど、実際には\(\frac{d}{d\theta} \log p(\mathbf{y}; \theta) = 0\)を解くじゃないですか、あの左辺を推定関数と呼んでいるわけね。ここでいっているのは、モデルはなんだかわからんけれど、ありうるすべてのモデルを通じて期待値をとったとき、推定関数に真の\(\theta\)を放り込んで得られる量の期待値は0であってほしい、ということだ。うん、そりゃまあそうですね]
推定関数の最適性
不偏推定関数\(g^*(\mathbf{y}, \theta)\)は、以下の条件を満たすときに最適であるという。その条件とは、すべての\(\xi \in C\)について適切な正則性をみたしつつ、他のすべての不偏推定関数のなかでもっとも $$ \frac{E_\xi[g^2]}{ ( E_\xi[\partial g / \partial \theta]_{\theta = \theta(\xi)} )^2 } $$ が小さいこと、である。そのとき方程式\(g^*(\mathbf{y}, \theta) = 0\)を最適推定方程式といい、その解を最適推定値という。
[そうそう！ ここが全然わかんないんだよな。私の理解を超える。たぶん、推定量の分散が小さいことに相当する基準だと思うんだけど…]
推定方程式の解の意義
\(\mathbf{y}= (y_1, \ldots, y_N)\)が調査母集団ベクトルなら、\(g^*(\mathbf{y}, \theta_N)=0\)の解\(\theta_N(\mathbf{y})\)は以下の2つの特徴を同時に持つ。
- \(\mathbf{y}\)のすべての要素が既知のとき、\(\theta_N(\mathbf{y})\)は超母集団パラメータ\(\theta\)の推定値である。
- \(\mathbf{y}\)が未知のとき、\(\theta_N(\mathbf{y})\)は調査母集団のパラメータの定義となる。
[前者はいいとして、後者の言葉遣いがちょっとわかりにくい。たとえば回帰モデルで、「もし母集団が利用可能だったら得られていたはずの\(\beta\)の推定値」のことをセンサス・パラメータっていうことがあるじゃないですか。そういう意味でパラメータと言っているのだと思う。たとえば超母集団モデルが\(y_i = \beta + \epsilon_i\)だとして、なんらかの推定方程式アプローチによる\(\beta\)の推定量\(\hat\beta\)を構成したとして、それは\(\beta\)の推定量であると同時に、未知の有限母集団平均そのものの定義でもあるよね… ということかな。
うーん、やっぱりちょっと腑に落ちない。デザイン・ベース・アプローチなら、有限母集団平均は定数であってモデルで定義されるものではないじゃないですか。モデル・ベース・アプローチでも、有限母集団の母平均とはデータ生成プロセスを\(N\)回動かして得た\(y_i\)の実現値たちの平均だ、と考えるのが普通ではないかと思う(だから「有限母集団平均を予測する」という言い方になる)。私の誤解かもしれないが、この論文って意外に変わったことを云ってないだろうか？]
2. 最適推定関数と線形最適推定関数
\(C = \{\xi\}\)が超母集団モデルで、そのもとで\(y_1, \ldots, y_N\)が独立だとしよう。いくつかの重要な事例では、最適推定関数\(g^*\)が以下の形をとる。$$ g^* = \sum_{i = 1}^N \phi_i(y_i, \theta)$$ ただし\(\phi_i\)は\( E_\xi[\phi_i] = 0\)であるようななんらかの実関数で、\(C\)によって一意に定まる。
[書き方が小難しくて困るが、最適推定関数が個体レベルのなにかの関数の合計になることがあるよ、ということであろう]
最適推定関数が存在しないこともある。そうした場合には、線形最適推定関数というより狭い概念を使うことがある。推定関数 \(g(\mathbf{y}, \theta)\)が線形であるとは、なんらかの実関数\(a_i(\theta)\)を使って$$ g(\mathbf{y}, \theta) = \sum_{i=1}^N a_i(\theta) \phi_i(y_i, \theta) $$ と書けるときである。\( E_\xi[\phi_i] = 0\)だから自動的に不偏推定関数となる。
[推定関数が、\(y_i\)を関数にいれて(返し値の期待値はゼロ)、パラメータで決まる個人別定数で重みづけて足し上げるという関数であったら、それを線形不偏推定関数というよ、ということであろう]
ここからは、次のようなクラス\(C\)に注目しよう。
- すべての\(\xi \in C\)に関して\(y_1, \ldots, y_N\)が独立。[iidとはいってない点に注意]
- すべての\(\xi \in C\)に関して\( E_\xi[\phi_i(y_i, \theta(\xi)] = 0\) [個人レベルで不偏推定関数があるってことね]
- \(g^* = \sum_{i=1}^N \phi_i(y_i, \theta)\)が\(\theta\)の最適推定関数ないし線形最適推定関数 [最適推定関数が個人レベル推定関数の積み上げの形になっているってことね]
推定関数が線形最適であることの十分条件
推定関数\( g^* = \sum_{i=1}^N \phi_i(y_i, \theta) \)が線形最適である十分条件は以下を満たすことである。すべての\(i\)について$$ \frac{E_\xi[\partial \phi_i / \partial \theta]_{\theta = \theta(\xi)}}{E_\xi[\phi^2_i]} = \mathrm{constant}(\theta(\xi)) $$ [証明が書いてある。略]
さらに、\(\xi\)が上に述べたクラスに属し、\(y_1, \ldots, y_N\)がiidで、すべての\(i\)について\(\phi_i = \phi\)なら、\(g^* = \sum_i \phi(y_i, \theta)\)が\(\theta\)の線形最適推定関数であることを示せる。Godambe & Thompson (1978, 1984)をみよ。
[なんかわからんけど要約すると、(1)\(y_i\)がiidで、(2)\(E_\xi[\phi(y_i, \theta)] = 0\)であるような個人レベル関数\(\phi(y_i, \theta)\)があって、(3)その関数を\(\theta\)で偏微分して得た関数の真の\(\theta\)における期待値が、その関数の二乗の期待値と比例していたら、\(g^* = \sum_i \phi(y_i, \theta)\)の解を\(\theta\)の推定量にするといいよ、それは不偏推定量だし、(最初に決めたなんだかよくわからん基準において)最適だよ、ということであろう]
線形最適推定方程式の解の意義
\(g^*\)が最適ないし線形最適なら、$$ \sum_{i=1}^N \phi_i(y_i, \theta_N) = 0 $$ の解として定義される\(\theta_N\)は、もし\(y_1, \ldots, y_N\)が利用可能ならば超母集団パラメータ\(\theta\)の推定値だし、そうでなければ調査母集団ベクトル\(\mathbf{y}\)のパラメータである。
事例
- \(C\)が以下を満たすすべての\(\xi\)からなるとしよう。
- \(y_1, \ldots, y_N\)は独立に分布する [iidとはいってない]
- すべての\(i\)について\( E_\xi[y_i] = \theta(\xi)\)
- \( E_\xi[y_i – \theta(\xi)]^2\)が\(i\)と独立
このとき、\(g^* = \sum_i \phi_i, \ \phi_i = y_i-\theta\)は線形最適である。さらに、\(y_1, \ldots, y_N\)がiidなら、この\(g^*\)は\(\theta\)の最適推定関数である。[よくわからん。\(y_1, \ldots, y_N\)が独立だけどiidじゃなかったら、\(g^* = \sum_i (y_i – \theta)\)は線形最適な関数だけど最適推定関数ではないってこと？]
\(\mathbf{y}\)に基づく\(\theta(\xi)\)の推定値\(\theta_N\)とは、調査母集団平均 $$ \theta_N = \frac{1}{N} \sum_{i=1}^N y_i $$ となる。\(\theta_N\)はそれ自体が調査母集団ベクトル\(y\)のパラメータでもある。 - 冒頭のモデル \( y_i = \beta x_i + \epsilon_i\)の\(\beta\)を\(\theta\)にしよう。\(g^* = \sum_i \phi_i(y_i, \theta), \ \phi_i(y_i, \theta) = x_i(y_i – \theta x_i)\)は線形最適であり、$$ \theta_N = \frac{\sum_{i=1}^N x_i y_i}{\sum_{i=1}^N x^2} $$ は\(\mathbf{x}, \mathbf{y}\)に基づく\(\theta\)の最適推定値であり、かつ調査母集団のパラメータである。
3. 超母集団パラメータと調査母集団パラメータの同時推定
設定
有限母集団\(P = \{1, \ldots, i, \ldots, N\}\)の部分集合である標本を\(s\)とし、\(S = \{s\}\)の確率分布を抽出デザイン\(p\)とし、$$ p(s) = Prob(\mathrm{sample \ is \ } s) $$ とする[原文ではP, Sは花文字]。標本\(s\)について\(y\)の値を観察して得たデータを $$ \chi_s = \{(i, y_i): i \in s\} $$ とする。
超母集団モデル\(C = \{ \xi \}\)が2節冒頭の3条件を満たすとしよう[\(y_i\)が独立で、個人レベルの不偏推定関数\(\phi_i(y_i, \theta)\)があって、\(\theta\)の(線形)最適推定関数\(g^*\)がその積み上げになっている]。\(g^* = 0\)の解を\(\theta_N\)とする。
デザイン不偏な推定関数
\(\theta_N\)の[標本からの]推定が、\(h(\chi_s, \theta) = 0\)の解として得られるとしよう。以下を要請するのが自然である。$$ E_p[h(\chi_s, \theta)] = \sum_{i=1}^N \phi_i(y_i, \theta) $$ ここでの期待値は抽出デザインの下での期待値である[原文は単に\(E\)だが、頭が混乱するので\(E_p\)とメモする]。いま述べたことは、1節冒頭で述べたデザイン不偏性基準に対応している。特に、抽出デザイン\(p\)での包含確率は $$ \pi_i = \sum_{s: i \in s} p(s) \gt 0 $$ であるから、関数 $$ h(\chi_s, \theta) = \sum_{i \in s} \frac{\phi_i(y_i, \theta)}{\pi_i} $$ はこの条件を満たす。
デザイン不偏な最適推定関数
では、この条件を満たす\(h\)のうち最適なのはどれか。それは以下の最小化として定義できる。$$ \frac{E_\xi[h^2(\chi_s, \theta(\xi))]}{ \left( E_\xi E_p \left[ \frac{\partial h}{\partial \theta} \right]_{\theta = \theta(\xi)} \right)^2} $$ [うううう… ここが理解できていない箇所である。信じて写経した]
上の条件とあわせるとこうなる。$$ E_\xi E_p \left[ \left( \frac{\partial h}{\partial \theta} \right)_{\theta = \theta(\xi)} \right] = E_\xi \left[ \frac{\partial}{\partial \theta} \sum_{i=1}^N \phi_i \right]_{\theta = \theta(\xi)} $$ [最小化する式の分母がどうなるかという話をしているんだけど思うけど、添え字\(\theta = \theta(\xi)\)は期待値記号のなかにつけても外につけても同じことなの？ ううううう…]
以下の点に注意されたい。
- 3条件における \(\sum_i \phi_i\)はある定数倍という形でしか決まらないが、その定数が決まれば、\(h\)の定数倍が決まる。
- 上の式の右辺は\(h\)とは独立であるから、結局、\(E_\xi[h^2(\chi_s, \theta(\xi))]\)を最小化すればよいことになる。
- \(E_p \left[ \left( h-\sum_{i=1}^N \phi_i \right)^2 \right] = E[h^2] – \left( \sum_{i=1}^N \phi_i \right)^2 \)だから、結局、$$ E_\xi E_p \left[ \left( h(\chi_s, \theta(\xi)) – \sum_{i=1}^N \phi_i(y_i, \theta(\xi)) \right)^2 \right]$$ を最小化すればよいことになる。これは結局、通常の抽出理論における最適性基準と整合している。
[わ、わからん… ついに完全に置いて行かれた…]
\(\theta_N\)の最適推定の定理
方程式 $$ \sum_{i \in s} \frac{\phi_i(y_i, \theta)}{\pi_i} = 0 $$ の\(\theta\)の解\(\hat{\theta}_s\)によって\(\theta_N\)を最適に推定できるのはいつか。以下の定理がその条件を示している。
定理1. 超母集団モデル\(C = \{\xi\}\)の下で\(\phi_i(y_i, \theta)\)が独立であり、それぞれの\(i\)について\(E_\xi[\phi_i(y, \theta(\xi))] = 0\)だとしよう。すべての\(\mathbf{y}, \theta\)について\(E_p[ h(\chi_s, \theta) ] = \sum_{i=1}^N \phi_i(y_i, \theta) \)を満たす推定関数 \(h(\chi_2, \theta)\)を考える。$$ \frac{E_\xi[h^2(\chi_s, \theta(\xi))]}{ \left( E_\xi E_p \left[ \frac{\partial h}{\partial \theta} \right]_{\theta = \theta(\xi)} \right)^2} $$ を最小化する\(h\)は下式である。$$ h^*(\chi_s, \theta) = \sum_{i \in s} \frac{\phi_i(y_i, \theta)}{\pi_i} $$ 証明: [略]
定理2. 定理1の下で、標本サイズを\(n\)に固定した時に$$ E_\xi E_p \left[ \left( \sum_{i \in s} \frac{\phi_i(y_i, \theta(\xi)}{\pi_i} \right)^2 \right] = \sum_{i = 1}^N \frac{E_\xi[\phi^2]}{\pi_i} $$ を最小化するデザインをの最適抽出デザインと呼ぶと、それは以下の包含確率を持つデザインである。$$ \pi_i \propto \sqrt{ E_\xi[\phi^2_i(y_i, \theta(\xi))] } $$
[…ちょっと中略…]
4. 例
例1.
超母集団モデル\(C = \{ \xi\} \)が、個々の\(\xi\)の下で\(y_1, \ldots, y_N\)が独立で、\(E_\xi[y_i] = \theta(\xi)\)で、\(y_i\)の分散が\(\sigma^2_i\)に比例し、\(\sigma^2_i\)が既知だとしよう。このとき、線形最適推定関数は$$ \phi_i(y_i, \theta) = \frac{y_i – \theta}{\sigma^2} $$ 調査母集団パラメータは $$ \theta_N = \sum_{i=1}^N \frac{y_i}{\sigma^2_i} / \sum_{i=1}^N \frac{1}{\sigma^2_i} $$ となる。最適推定方程式の解として得られる最適推定量は$$ \hat{\theta}_s = \sum_{i \in s} \frac{y_i}{\pi_i \sigma^2_i} / \sum_{i \in s} \frac{1}{\pi_i \sigma^2_i} $$となる。最適デザインは $$ \pi_i \propto \sqrt{ E_\xi \left[ \left( \frac{y_i – \theta(\xi)}{\sigma_i^2} \right)^2 \right] } = \frac{1}{\sigma_i} $$ となる。
例2.
例1.とほかは同じで、\(E_\xi(y_i) = \theta(\xi) x_i\)であるとする。\(x_i\)は母集団について既知とする。[おっとー。頭が大混乱したが、もはや母平均の推定ではない。\(y_i = \theta x_i + \epsilon_i\)というモデルを考えていて、この\(\theta\)を推定したいという話になっているのだ]
線形最適推定関数は $$ \phi_i(y_i, \theta) = \frac{y-\theta x_i) x_i}{\sigma^2} $$ 調査母集団パラメータは $$ \theta_N = \sum_{i=1}^N \frac{y_i x_i}{\sigma^2_i} / \sum_{i=1}^N \frac{x^2_i}{\sigma^2_i} $$ 最適推定量は$$ \hat{\theta}_s = \sum_{i \in s} \frac{y_i x_i}{\pi_i \sigma^2_i} / \sum_{i \in s} \frac{x_i^2}{\pi_i \sigma^2_s} $$ 最適デザインは $$\pi_i \propto \frac{x_i}{\sigma_i} $$ となる。
特殊ケースをふたつ考えよう。
- \(\sigma^2_i \propto x_i\)の場合は、$$ \theta_N = \sum_{i=1}^N y_i / \sum_{i=1}^N x_i, \ \hat{\theta}_s = \sum_{i \in s} \frac{y_i}{\pi_i} / \sum_{i \in s} \frac{x_i}{\pi_i}$$
- \(\sigma^2_i \propto x_i^2\)の場合は、$$ \theta_N = \sum_{i=1}^N y_i \sum_{i=1}^N \frac{y_i}{x_i}, \ \hat{\theta}_s = \sum_{i \in s} \frac{y_i}{\pi_i} / \sum_{i \in s} \frac{x_i}{\pi_i}$$
[ああ、なるほどね。超母集団モデルの構造の部分が一緒でも、攪乱項の分散についての仮定が変わると\(\theta_N\)は変わるのだ。だって\(\theta_N\)はモデルのパラメータではなくてセンサス・パラメータだから]
最後に、最適デザインについて。分析的調査(超母集団パラメータの推定が主目的である調査)では、より安定している個人の包含確率を高くしたほうがいいから、\(\sigma_i\)が小さい人の\(\pi_i\)を大きくしたほうがいい。これが上の例1, 例2である。いっぽう、累積的調査(調査母集団平均の推定が主目的である調査)では、より変動が大きい人の\(\pi_i\)を大きくしたほうがよい。Godambe(1982)をみよ。[うおおお、なにかすごく大事な話をしておられる… しかしよく理解できない。後者の場合はどういう推定方程式が最適なの？]
5. パラメータが複数ある場合への拡張
[回帰モデルのパラメータを推定したいというような場合への拡張。さすがに力尽きたのでパス]
6. 考察
[ここから逐語訳する。もはや気分は修行僧]
3節では、2節冒頭の3条件で示した超母集団モデルから導出される調査母集団パラメータ\(\theta_N\)の最適推定を示した。しかし、3節の定理1のなかの\(\phi_i\)を\(a \phi_i\)に代えても、やはり定理1は成り立つ。たとえば4節の例1についてみてみよう。もし\(a_i = \sigma_i^2\)なら、\(\sum_i a_i \phi_i = \sum_i(y_i – \theta)\)である。そのときの調査母集団平均\(\theta’_N\)は方程式\(\sum_{i=1}^N (y_i – \theta’) = 0\)で定義される。最適推定関数は$$ h^*(\chi_s, \theta) = \sum_{i \in s} \frac{\sigma^2 \phi_i}{\pi_i} = \sum_{i \in s}\frac{y_i – \theta}{\pi_i} $$ となり、調査母集団平均の推定値は$$ \hat{\theta}_s = \sum_{i \in s} \frac{y_i}{\pi_i} / \sum_{i \in s} \frac{1}{\pi_i} $$ となる[Hajek推定量だあああ]。3条件のモデルの動機づけられたパラメータは\(\theta_N\)であって\(\theta’_N\)ではないが、定理1を広く解釈すれば、それは\(\theta’_N\)に関しても上記の推定関数を与えてくれている。
伝統的な不偏最小期待分散推定量と、本論文の推定量を比べてみよう。
調査母集団平均の不偏最小期待分散推定値は、特定の\(\theta = \theta(\xi)\)を持つ\(\xi\)に限定していえば、$$ e = \frac{1}{N} \left( \sum_{i \in s} \frac{y_i-\theta(\xi)}{\pi_i} + \theta(\xi)N \right) $$ である。\(\theta(\xi)\)は未知だから、\(e\)がいくら最適でもそれ自体は無意味だが、伝統的理論における最適性の結果の特殊ケースが、以下の事実から得られる。\(\pi_i \gt 0\)が指定され、抽出デザインが\( Prob(s: \sum_{i \in s} \frac{1}{\pi_i} = N \} = 1\)であるとき[抽出確率の逆数の標本合計が母集団サイズになるとき、つまりHajek推定量とHorvits-Thompson推定量が等しくなるときね]、推定値\(e\)は\(\theta\)から独立になり、\(\hat{\theta}_s = e\)となる。4節の例2で議論したモデルからも類似した結果が得られる。上記の条件を満たす抽出デザイン(補足変数\(x\)があるなら、それを考慮したうえで上記の条件を満たす抽出デザイン)のなかには、層別抽出デザインのさまざまな種類が含まれる。このこと、そして本節冒頭の\(h^*(\chi_s, \theta)\)の最適性がモデル\(C\)全体をおいて維持されるという事実から、次のことがわかる。本論文で構築した理論は、伝統的理論の一般化なのである。
[まだまだ難しい話が続くけど、パス]
—————–
そうじゃないかとは思っていたが、やはり難しくてよくわからんかった。ともあれ、ある種のセッティングとある種の最適性概念のもとでHajek推定量が最適推定量だといえる、ということがわかりました。そういうのを「わかった」といってよいならの話だが。
素朴な疑問過ぎて恥ずかしいが、こんなに大がかりな話が必要なの？ Hajek推定量をモデル・ベースで正当化するためには、推定方程式アプローチのような難しい話がどうしても必要なのだろうか？ ごくふつうに、\(Y_i\)についてのこれこれのモデルの下で、Hajek推定量がモデル・パラメータの最良線形不偏推定量になりますとか、そのモデルから生成される\(Y_i\)の\(N\)個の実現値の平均の最良線形不偏予測量になりますとか、そういう風な、俺にわかるようなかんたんな正当化はできないものなの？