Dorfman, A.H., Valliant, R. (1997) The Hajek Estimator Revisited. Proceedings of the Section on Survey Methods Research, American Statistical Association.
標本ウェイティングについて調べ物をしていると、googleの検索結果にはいつもこれが上位に挙がってくる。内容は難しそうだし、そもそも論文ではないし、毎回「これはちがうだろう」と避けていたのだが、あまりによく見かけるので、これもなにかの勉強だろうと思い、メモをとりながら読んでみた。正直、かなり早い段階で後悔しはじめたのだが…
著者はこの段階ではUS労働統計局の人。第二著者は昨年読んだ Elliot & Valliant (2017) のValliantさんで、たまたま見つけた略歴によれば、役所をやめたあとしばらくWestatに勤め、その後大学に転じて現在はミシガン大の先生らしい。
1. Hajek推定量とは
要素数\(N\)の母集団\(U\)の平均\(\bar{Y} = \frac{1}{N} \sum_i^N Y_i\)について推定したい。サイズ\(n\)に固定した標本を得る。ところが\(U\)を通じて既知な正の値\(x_i\)がある。そこで、\(\bar{x} = \frac{1}{N} \sum_i^N x_i\)として、標本包含確率\(\pi_i=\frac{nx_i}{N\bar{x}}\)でpps抽出した。このときあなたは、Hajek推定量 $$ \hat{\bar{Y}}_{Haj} = \frac{\sum_s \frac{Y_i}{\pi_i}}{\sum_s \frac{1}{\pi_i}} $$ を使いますね？ 使いますよね？ [畳みかけるスタイル]
Sarndal-Swenson-Wretman本いわく、pps不偏なHT推定量$$ \hat{\bar{Y}}_{HT} = \frac{1}{N} \sum_s \frac{Y_i}{\pi_i}$$ なんかより、Hajek推定量のほうが気が利いている。Hajek推定量は最適推定方程式の理論からも導出できる。つまり、\(\bar{Y}\)を超母集団モデル\(Y_i \sim_{iid} (\mu, \sigma^2)\)の下で生み出された有限母集団パラメータだと考えると、Hajek推定量が出てくるのよ。Godambe & Thompson (1986) の例1をみろ。
本稿では、調査における抽出についてのモデルベースというか予測アプローチの観点から、Hajek推定量とHT推定量について検討する。
2. 背景: 単純釣り合い標本
(モデル)
\(Y_i, x_i\)がこうなってるとすんじゃんか。\(\delta_j\)を二値変数として $$ Y_i = \sum_{j = 0}^{J} \delta_j \beta_j x_i^j + \epsilon_i v_i^{1/2}, \ \ \epsilon \sim_{iid} (0, \sigma^2) $$ これを\(M(\delta_0, \delta_1, \ldots, \delta_J; v)\)って書くね。そのBLUPを\(\hat{\bar{Y}}(\delta_0, \delta_1, \ldots, \delta_J; v)\)って書くね。
たとえば、\(M(0, 1:x)\)ってのは$$ Y_i = \beta_i x_i + \epsilon_i x_i^{1/2} $$ のこと。そのBLUP は比推定量\(\hat{\bar{Y}}_R = \bar{Y}_s \frac{\bar{x}}{\bar{x}_s}\)。
(頑健性の調べ方)
model-based opusにおける主要な関心事のひとつは頑健性である[opusってなんのことだろう?]。その検討のためには、単純な作業モデルともっと複雑な真モデルを考える。
真モデルが\(M(\delta_0, \delta_1, \ldots, \delta_j: v)\)だったとき、比推定量\( \hat{\bar{Y}}(0,1:x) \) はどうなるか。こうなります。\(\bar{x}^{(j)}_s = \frac{1}{n}\sum_s x^j_i, \bar{x}^{(j)}_s = \frac{1}{n} \sum_i^N x^j_i\) と書く。モデルベースだということを表すために期待値に添字\(M\)をつける。$$ E_M[\hat{\bar{Y}}(0,1:x) -\bar{Y}] = \bar{x} \sum_j^J \delta_j \beta_j \left[ \frac{\bar{x}^{(j)}_s}{\bar{x}_s} – \frac{\bar{x}^{(j)}}{\bar{x}} \right] $$ 当然ながら\(j=1\)の項は効かない。
[はい深呼吸！ 真モデルが\(x\)の\(J\)次多項式の範囲にあるときに\(x\)でうっかり比推定しちゃったらどんなバイアスが起きますか、という話をしているのである。たとえば真モデルが$$Y_i = \beta_1 x_i + \beta_2 x_i^2 + \epsilon_i v_i^{1/2} $$ であるときにうっかり\(\hat{\bar{Y}}_R = \bar{Y}_s \frac{\bar{x}}{\bar{x}_s}\)を使っちゃったら、バイアスは$$ \bar{x} \beta_2 \left[ \frac{\bar{x}_s^{(2)}}{\bar{x}_s} – \frac{\bar{x}^{(2)}}{\bar{x}} \right] $$ になる]
(単純釣り合い標本)
もしすべての\(j\)について\( \bar{x}_s^{(j)} = \bar{x}^{(j)} \) であるような標本があったら、比推定量は不偏になる。そういう標本を釣り合い標本という。
[比推定量に対応するモデルが正しいとき、ないし、実はモデル\(x\)の\(J\)次多項式なんだけどすべての\(j\)について\(x_i^j\)の標本平均と一致しているとき(比推定量はやっぱり不偏)、その標本を釣り合い標本という、というわけね。要は、「目的変数に効いてる共変量について、母集団からの共変量シフトがない」ということであろう]
以下ではこういうのを単純釣り合い標本と呼ぶ。本項でいいたいのは、この条件を満たすよう標本を注意深く抽出すれば、ある種のモデルエラーを防げる、ということである。
(Hajek推定量と比推定量の頑健性)
抽出スキーマがSRSなら、Hajek, HT, 標本平均はみな同じですよね。もし標本が釣り合いだったら、比推定量も標本平均になる。さらに、SRSでは、標本は釣り合い標本だと期待される。つまり、いま確率抽出スキーマ下での期待値を\(E_\pi[\cdot]\)と書くとして、$$ E_\pi[\bar{x}_s^{(j)}] = \frac{1}{n} \sum_i^N \pi x_i^j = \bar{x}^{(j)} $$ である。このように、SRSとHajek/HTの組み合わせは、モデル・ベース理論からも追加の支持を受けているように思われるわけだ。
しかし、これは事実だろうか？ 推定量の利点を評価するには、釣り合いが維持されなかったときに不偏性が維持されるか、を調べなければならないだろう。
比推定量はその点でHajek推定量より敏感性が低い。たとえば、真のモデルが\(M(1,1,1,:v)\)だとしよう[つまり\( Y_i = \beta_0 + \beta_1 x_i + \beta_2 x^2_i + \epsilon_i v_i^{1/2}\)]。\(\bar{x}_s = (1+e_1)\bar{x}, \bar{x}_s^{(2)} = (1+e_2)\bar{x}^{(2)} \)とする[つまり、補助変数はあるけど標本平均と二乗の標本平均が真値からずれている]。Hajek/HT推定量のモデルバイアスは[補助変数を無視するから] $$ E_M[\bar{Y}_s -\bar{Y}] = \beta_1 \bar{x} e_1 + \beta_2 \bar{x}^{(2)} e_2 $$ だ。いっぽう比推定量のモデルバイアスは $$ E_M[\hat{\bar{Y}}_R – \bar{Y}] = -\beta_0 e_1 + \beta_2 \bar{x}^{(2)} (e_2 – e_1)$$ である。ふつう、Hajek推定量のバイアスの第二項は比推定量のバイアスの第二項より大きな絶対値を持つ。というのは、多くの場合\(e_1\)と\(e_2\)は同じ符号なので、比推定量のバイアスの第二項はキャンセルされるだろうからだ。第一項もたぶんHajekのほうが大きい。比推定量を使いたい時というのは\(Y\)が\(x\)と連動しているように思われるときだからだ。確かに非釣り合いは怖いけど、Hajek推定量のほうがもっと敏感なのである。
というわけで、もしHajek推定量を使うなら、SRSに頼らず、注意深く釣り合いをとって抽出すべきだ。「モデル・エラーを防ぐ最善の方法はランダム化だ」というよくある主張は、別のモデルで何が起きるのかという視点から見ると必ずしも正しくないわけである。
3. 重み付き釣り合い標本
モデル\(M(0, 1:x)\)のもとでは、釣り合い標本＋比推定という抽出戦略はもっとも効率的な戦略ではない。\(x_i\)が大きいときに\(Y_i\)の分散が大きくなるわけだから、\(x_i\)が大きい要素を強めに抽出すべきだからだ。
[はいはい、なるほどね。話の筋がみえてきたような気がする。層別抽出でいうところの最適配分みたいな話をしたいわけね]
(モデルとBLUP)
以下では、長さ\(n\)の1のベクトルを\(\mathbf{1}_s\)、長さ\(N-n\)の1のベクトルを\(\mathbf{1}_r\)、長さ\(N\)の1のベクトルを\(\mathbf{1}_N\)と書く。
以下のモデルを考える。$$ E_M(\mathbf{Y}) = \mathbf{X} \beta$$ $$ var_M(\mathbf{Y}) = \mathbf{V} \sigma^2$$ \(\mathbf{V}\)は対角行列とする。これを\(M(\mathbf{X}:\mathbf{V})\)と書く。補足変数の行列を標本と非標本に分けて、\(\mathbf{X} = (\mathbf{X}_s^\top, \mathbf{X}^\top_r)^\top\)とする。
上のモデルの下でのBLUPはどうなるかというと、$$ \hat{\bar{Y}}(\mathbf{X}:\mathbf{V}) = \frac{1}{N} ( \mathbf{1}^\top_s \mathbf{Y}_s + \mathbf{1}^\top_r \mathbf{X}_r \hat{\beta} ) $$ $$ \hat{\beta} = \mathbf{A}_s^{-1} \mathbf{X}^\top_s \mathbf{V}^{-1}_{ss} \mathbf{Y}_s$$ $$ \mathbf{A}_s = \mathbf{X}^\top_s \mathbf{V}_{ss}^{-1} \mathbf{X}_s$$ ただし\(\mathbf{V}_{ss}\)は標本の共分散行列で、\(n \times n\)の対角行列である。[ここは著者を信じて写経した]
(重み付き釣り合い標本)
さて、以下の条件を満たす標本を\(root(\mathbf{W})\)釣り合い標本 \(B(\mathbf{B}:\mathbf{W})\)と呼ぼう。\(\mathbf{W}\)を\(N\times N\)行列、そこから標本を取り出したのを\(\mathbf{W}_s\)として$$ \frac{1}{n} \mathbf{1}^\top_s \mathbf{W}^{-1/2}_s \mathbf{X}_s = \frac{\mathbf{1}^\top_N \mathbf{X}}{\mathbf{1}^\top_N \mathbf{W}^{1/2} \mathbf{1}_N} $$ モデルが \(\mathbf{V} = \mathbf{W} \sigma^2\)のとき、\(root(\mathbf{W})\)釣り合い標本が最適となる。
[細かいところはわかんないけど雰囲気はわかる。モデル上で攪乱項の分散が個体ごとに決まっているとき、それを考慮して「この標本なら共変量シフトが起きない」といえるような標本を\(root(\mathbf{W})\)釣り合い標本と呼んでるわけだ]
もし\(\mathbf{W} = \mathbf{I}\)なら、\(B(\mathbf{B}:\mathbf{I})\)は列に\(\mathbf{X}\)について釣り合いが取れた標本の集合、つまり\( \mathbf{1}^\top_s \mathbf{X}_s / n = \mathbf{1}^\top_N \mathbf{X}_N / N \)であるような標本の集合となる。もし\(Y\)のモデルが\(x\)の多項式なら、\(\bar{x}^{(j)}_s = \bar{x}^{(j)}\)を満たす標本集合、つまり2節でいうところの釣り合い標本となる。
(BLUPの分散の下界)
Royall(1992)は以下を示している。\(\mathbf{X}\)の列によって生成される線形多様体、つまり、\(\mathbf{X}\)の列のすべての線形結合で張ったベクトル空間を\(m(\mathbf{X})\)とする。
定理1. \(M(\mathbf{X}: \mathbf{V})\)の下で、\(\mathbf{V1}_N \in m(\mathbf{X})\), \(\mathbf{V}^{1/2} \mathbf{1}_N \in m(\mathbf{X})\) がともに満たされるとき、$$ var_M \left[\hat{\bar{Y}}(\mathbf{X}: \mathbf{V}) – \bar{Y} \right] \geq N^{-2} \left[ n^{-1} (\mathbf{1}^\top_N \mathbf{V}^{1/2} \mathbf{1}_N)^2 – \mathbf{1}^\top_N \mathbf{V} \mathbf{1}_N \right] \sigma^2$$ 等号成立の必要十分条件は\(s \in B(\mathbf{X}: \mathbf{V})\)で、このとき\(\bar{v}^{(1/2)} = \bar{v}^{(1/2)} = N^{-1} \sum_{i=1}^{N} v_i^{1/2} \)として$$ \hat{\bar{Y}}(\mathbf{X}:\mathbf{V}) = \frac{\bar{v}^{(1/2)}}{n} \sum_s \frac{Y_i}{v_i^{1/2}} $$となる。
[常人にわかる言葉に翻訳しよう。目的変数が説明変数\(\mathbf{V}\)の線形モデルでできていると考えた分析者が、そのモデルの下で母平均をBLUに推定したとする(パラメータのBLUEではなくて有限母集団平均のBLUPということであろう)。上の定理は、そのときの推定量の分散の下界を与えている。下界が達成されるのは、標本が攪乱項の真の共分散行列\(\mathbf{V}\)のもとでの\(root(\mathbf{V})\)釣り合い標本であるときで、BLUな推定量は、標本の\(Y_i\)を各個体の真のSDで割って合計し、それに真のSDの母平均を掛けて\(n\)で割った量だ。
あれれ？ きちんと式で書くと $$ \hat{\bar{Y}} = \left( \frac{1}{N} \sum_U v^{1/2}_i \right) \times \frac{1}{n} \sum_s \frac{Y_i}{v^{1/2}_i} $$ だよね？ デザインベース推定量と比べてみよう。\(\pi’_i = c \pi_i\)とすると、その母平均は\(c\frac{n}{N}\)だから、HT推定量は$$ \hat{\bar{Y}}_{HT} = \frac{1}{N} \sum_s \frac{Y_i}{\pi_i} = \frac{c}{N} \sum_s \frac{Y_i}{\pi’_i} = \left(\frac{1}{N} \sum_U \pi’_i \right) \times \frac{1}{n} \sum_s \frac{Y_i}{\pi’_i}$$ となる。このモデルからみてBLUなのはHT推定量なの？ 実はHajek推定量だという落ちかと思ってた。
あっれーー？ だんだん頭が混乱してきた。この論文は有限母集団の推定について考えているんだよね？ モデルベースではあるけれど、推定対象はあくまで\(\bar{Y} = \frac{1}{N} \sum_U Y_i\)で、モデルのパラメータじゃないよね？ デザインベースなら推定量の分散には \(fpc= 1-\frac{n}{N}\) みたいな項が必ず入るじゃないですか。モデルベースでも母平均のBLUPの分散にはfpcが入るはずじゃない？ でも定理で示された下界って、これfpcっぽい項が入ってなくない？ 俺にわからないだけで実は入ってんの？]
次の点に注目してほしい。重み付き釣りあい条件の下では、推定量それじたいも、その分散も、\(\mathbf{X}\)に明示的には依存しない。\(\mathbf{X}\)の列が同じなら、\(\mathbf{X}\)自体がすり替わっても推定量は依然として不偏だし効率性は下がらない。
(推定量の比較)
- 例1. 作業モデルが\(M(1,1,1:x^2)\)だとしよう[\(Y_i = \beta_0 + \beta_1 x_i + \beta_2 x^2_i + x_i \epsilon_i \)ね]。\(x, x^2\)がモデルに入っているから定理の条件は満たされている。分散の下界は$$ \frac{1}{N^2} \left[ \frac{(N \bar{x})^2}{n} – \sum_i^N x^2 \right] \sigma^2 $$ である[あ、\(x_i = 1\)なら\(\frac{1}{N^2} (\frac{N^2}{n} – N)\sigma^2 = (1-\frac{n}{N}) \frac{\sigma^2}{n}\)だ。fpcはいってた。失礼しました]
下界に到達するのは$$ \bar{x}^{(j-1)}_s = \frac{\bar{x}^{(j)}}{\bar{x}} $$ が\(j=0,1,2\)で成り立つことである。\(j=1\)では自明。\(j=0\)は、標本の\(x\)の調和平均が、母集団の算術平均に等しいことを指す。[おおお、そうなのか…]
以下、めんどくさいので、これを単に(\(J\)次で)x-釣り合いだと呼ぼう。先行研究では「\(\pi\)-釣り合い」と呼ばれている。
BLUPは\(\hat{\bar{Y}} = \bar{x} \sum_s \frac{y_i}{nx_i} \)であり(「平均の比」推定量)、\(\pi_i \propto x_i\)で\(n\)が固定ならHT推定量である。\(j=0\)で釣り合ってるんならつまりHajek推定量でもある。x-釣り合いの下では、モデルのBLUPと、pps抽出下のHT推定量・Hajek推定量が等しくなるわけだ。 - 例2. 作業モデルが\(M(0, 1:x)\)だとしよう。[…中略…]
- より一般的に、\(n\)固定, サイズ変数\(v^{1/2}_i\)のpps抽出、つまり包含確率が$$ \pi_i = \frac{n v^{1/2}_i}{\mathbf{1}^\top_N \mathbf{V}^{1/2} \mathbf{1}_N} $$ である抽出を考えよう。これをpp(\(v^{1/2})\)抽出と略記する。この抽出は以下の性質を持つ。
- \(X\)がどうであれ、デザイン期待値において重み付き釣り合いが成立する。
- \(root(v)\)釣り合い標本のBLUPはHT推定量になる。
- \(\frac{n}{\mathbf{1}^\top_s \mathbf{V}^{-1/2} \mathbf{1}_s} = \frac{1}{N} \mathbf{1}^\top_N \mathbf{V}^{1/2} \mathbf{1}_N \)ならばBLUPはHajek推定量になる。
釣り合い標本を選ぶ方法のひとつは、上記のpps抽出計画を用い、かつ特定の積率においてきちんと釣り合いが取れていない標本を避けることである。図1をみてほしい。[いや、図なんて載ってないよ？！ まったくもう… 後略]
4. 非釣り合い標本における推定量の比較
例1の状況に注目する。以下、所与の共分散行列\(\mathbf{V}\)について \(M(\mathbf{X}_V; \mathbf{V}), \mathbf{X}_V = (\mathbf{V}^{1/2} \mathbf{1}_N, \mathbf{V1}_N) \)のことを最小モデルという。
分散が\(x^2\)に比例するとき、最小モデルは\(M(0,1,1:x^2)\)である。このモデルに対応するBLUPは、釣り合いが取れていようがいなかろうが$$ \hat{\bar{Y}}_{BLU} = \frac{1}{n(\bar{x}^{(2)}_s – (\bar{x})^2)} \left( (\bar{x} \bar{x}^{(2)}_s – \bar{x}^{(2)} \bar{x}_s) \sum_s \frac{Y_i}{x_i} + (\bar{x}^{(2)} – \bar{x} \bar{x}_s) \sum_s Y_i \right)$$ である。
さて、実は$$ \bar{x}^{j-1}_s = \frac{\bar{x}^{(j)}}{\bar{x}} (1+e_j)$$ だとしよう。この例では\(e_1 = 0\)である。\(J\)次多項式モデルの下でのモデルバイアス\(E_M(\hat{\bar{Y}} – \bar{Y})\)を、HT推定量、Hajek推定量、最小モデルBLUPを求めてみると以下となる。$$ Bias_{Haj} \sim \sum_{j=0}^J \beta_j \bar{x}^{(j)} (e_j – e_0) $$ $$ Bias_{HT} = \sum_{j=0}^{J} \beta_j \bar{x}^{(j)} e_j $$ [BLUPのバイアスの式はめちゃくちゃややこしいのでメモ省略]
もし切片以外のすべての係数が0なら、釣り合いからの逸脱に強いのはHajek推定量であることがわかる。これは冒頭で紹介したSSWのHajek推しに対応している。でも、\(Y\)が\(x\)に定数でない連続的な依存性を持っているなら、むしろHTのほうがバイアスが減る。たぶん\(e_0\)とそれ以外の\(e_0\)の符号が違うからである。そしてBLUPの式をよく見ると、これはHTよりバイアスが小さい。[？？？ よくわからん…]
5. シミュレーション研究
[力尽きた。まるごとパス]
6. 結論
Hajek推定量は(そしてHajek推定量ほどではないけどHT推定量も)、所与の抽出計画の下で、重みつき釣り合いからの逸脱に対して敏感である。抽出計画が期待値において重み付き釣り合いを達成していたとしても、やはり敏感なのである。例外は補足変数と目的変数が無関連だった時だけである。
抽出計画に対応する最小推定量を用いるのが望ましい。標本を重み付き釣り合い標本に制約するのがベストである。不釣り合い標本の場合、最小モデルに対応し切片を修正したBLU推定量[?]が、補足変数に関連する目的変数についても、また共通の平均を持つ複数の目的変数についても、効率的になる。
———-
途中から「なぜ俺はこれを読んでいるんだ…」という実存的な疑念にとらわれたが、意地になって読み終えた。
単純釣り合い標本を重み付き釣り合い標本へと拡張したあたりから話についていけなくなってしまったので、内容をほとんど理解できていないのだけれど、標本において運悪く共変量シフトが起きたときに推定量のバイアスがどうなるか、という話をしているのだと思う。場合によってはHajek推定量がいいとはいえなくなり、共変量を組み込んだモデルに基づいて推定量を考えたほうがいい、という結論なのだろうな(事後層別みたいなものであろう)。しらんけど。