前回に引き続き、
Breidt, F.J., Opsomer, J.D. (2017) Model-assisted survey estimation with modern prediction techniques. Statistical Science, 32(2), 190-205.
のメモ、その第二回。いよいよ、著者の云うところのモデル・アシステッド推定とはなんなのかが明らかになる。話の先取りになるけど、要するにGREG推定量のことなのだ。ナアンダというなかれ、俺はこういう風に考えたことなかったよ。
4. モデル・アシステッド推定
差分推定量は標本とは独立な手法\(m(\cdot)\)を必要としますが、実務的にはなかなか困難です。モデル・アシステッド推定アプローチでは、予測のための「作業モデル」を導入します。多くの場合、次の形で書けるモデルです。$$ y_k = \mu(\mathbf{x}_k) + \epsilon_k $$ \(\{\epsilon\}\)は平均ゼロの確率変数です。ここでは母集団の\(\{y_k\}_{k \in U}\)を、超母集団の確率的モデルの実現値としてモデル化します。
ここで重要なのは、このモデルは別に正しくなくてもいい、ということです。ある程度の予測力があればそれでいいのです。
推定・推論のための一般的なレシピを紹介しましょう。
- もし \(\{\mathbf{x}_k, y_k\}_{k \in U}\)が母集団全体について観察可能なら、\(\mu(\cdot)\)を推定する標準的な統計手法が\(m_N(\cdot)\)となります。それは母集団に依存しますが標本に依存しません。[想像しづらい状況だが、論理的な可能性として述べているのだろう]
- 標本しか観察できない場合、\(m_N(\cdot)\)を\(\hat{m}(\cdot)\)によって推定します。
- 差分推定量の式に\(\hat{m}(\cdot)\)をプラグインし、モデル・アシステッド推定量を得ます。$$ DIFF(y, \hat{m}) = \sum_{k \in U} \hat{m}(\mathbf{x}_k) + \sum_{k \in s} \frac{y_k – \hat{m} (\mathbf{x}_k)}{\pi_k} $$ [あれ、まだ母集団の全個体の\(\mathbf{x}_k\)が利用可能だと仮定しているのか… いずれは集計値しかわからない場合にもっていかないといけないはずだが]
- 推論のためには、大標本において推定量が漸近的正規性を持つと仮定します。分散の推定は$$ \hat{V}(DIFF(y, \hat{m})) = \sum_{k,l \in U} \Delta_{kl} \frac{y_k – \hat{m}(\mathbf{x}_k)}{\pi_k} \frac{y_l – \hat{m}(\mathbf{x}_l)}{\pi_l} \frac{I_k I_l}{\pi_{kl}} $$
このレシピは驚くべき柔軟性と応用可能性を備えています。次の節では、モデル・アシステッド推定量の漸近特性を得る方法を一般的に述べましょう。
5. サーベイの漸近解析 I
モデル・アシステッド・推定量の式を書き換えます。$$ DIFF(y, \hat{m}) = \sum_{k \in U} \hat{m}(\mathbf{x}_k) + \sum_{k \in s} \frac{y_k – \hat{m}(\mathbf{x}_k)}{\pi_k} $$ $$ = \sum_{k \in U} m_N (\mathbf{x}_k) + \sum_{k \in U} \frac{(y_k – m_N(\mathbf{x}_k)) I_k}{\pi_k} + \sum_{k \in U} (\hat{m}(\mathbf{x}_k) – m_N(\mathbf{x}_k))\left( 1 – \frac{I_k}{\pi_l} \right) $$ $$ = DIFF(y, m_N) + (残り)$$ つまり、モデル・アシステッド推定量とは、\(m_N(\cdot)\)に基づく差分推定量(それは手に入りませんが不偏です)と残りの部分の和です。残りの部分が無視できる程度に小さいことを示す、というのが課題です。それをどうやるかは推定方法によりますが、それがうまくいけば、モデル・アシステッド・推定量は直ちに差分推定量の漸近特性を受け継ぎます。
モデル・アシステッド推定量は漸近不偏です。なぜなら差分推定量が(作業モデル\(\mu(\cdot)\)がどうであれ)不偏だからです。モデル・アシステッド推定量は、HT推定量がデザイン一致・漸近正規である条件と似た条件を満たせば、デザイン一致で漸近正規です。分散は差分推定量の分散の式に\(\hat{m}(\cdot)\)をプラグインした式で一致推定できます。そして、ありがたいことに、\(HT(y)\)より分散が小さいのです。
[夢のある話をありがとう。そろそろ現実的な話をしようぜ]
次の節からは、このモデル・アシステッド推定量の特殊ケースについて検討します。
6. 一般化回帰推定
モデル・アシステッド推定量のもっともよく知られているクラスとして、一般化回帰推定量、別の名をGREG推定量があります。
作業モデルを分散不均一な重回帰モデルにします。$$ y_k = \mathbf{x}^\top_k \beta + \epsilon_k$$ $$ \{\epsilon_k\} \ \mathrm{uncorrelated} \ (0, \sigma^2_k) $$
もし母集団全体が観察可能なら、\(\beta\)はWLS推定すればいいでしょう。つまり $$ m_N (\mathbf{x}_k) = \mathbf{x}^\top \mathbf{B}_N = \mathbf{x}^\top_k \left( \sum_{j \in U} \frac{\mathbf{x}_j \mathbf{x}^\top_j}{\sigma^2_j} \right)^{-1} \sum_{j \in U} \frac{\mathbf{x}_j y_j}{\sigma^2_j} $$ です。しかし標本しか観察できないので、$$ \hat{m}(\mathbf{x}_k) = \mathbf{x}^\top_k \hat{\mathbf{B}} = \mathbf{x}^\top_k \left( \sum_{j \in s} \frac{\mathbf{x}_j \mathbf{x}^\top_j}{\sigma^2_j} \right)^{-1} \sum_{j \in s} \frac{\mathbf{x}_j y_j}{\sigma^2_j} $$ をプラグインして $$ DIFF(y, \mathbf{x}^\top_k \hat{\mathbf{B}}) = \sum_{k \in U} \mathbf{x}^\top_k \hat{\mathbf{B}} + \sum_{k \in s} \frac{y_k – \mathbf{x}^\top_k \hat{\mathbf{B}}}{\pi_k} $$ これがGREG推定量です。
よく知られている例として、事後層別推定量(\(\mathbf{x}_k\)はカテゴリカル共変量の水準のインジケータのベクトル)、古典的サーベイ比推定量(\(x_k\)はスカラーでモデルは原点を通る分散不均一回帰)、古典的サーベイ回帰推定量(\(x_k\)はスカラーでモデルは分散均一な回帰)、が挙げられます。[おおおおおお… 古色蒼然とした奴らが勢ぞろいだ]
GREG推定量の特性について考えてみましょう。$$ DIFF(y, \mathbf{x}^\top_k \hat{\mathbf{B}}) = \sum_{k \in U} \mathbf{x}^\top_k \hat{\mathbf{B}} + \sum_{k \in s} \frac{y_k – \mathbf{x}^\top_k \hat{\mathbf{B}}}{\pi_k} $$ $$ = \sum_{k \in U} \mathbf{x}^\top_k \mathbf{B}_N + \sum_{k \in U} \frac{(y_k – \mathbf{x}^\top_k \mathbf{B}_N) I_k}{\pi_k} + (\hat{\mathbf{B}}_N – \mathbf{B}_N)^\top (t_x – HT(\mathbf{x})) $$ [第1項、元の式の\(\hat{\mathbf{B}}\)を\(\mathbf{B}_N\)にすり替えたので、\(\sum_{k \in U} (\mathbf{x}^\top_k (\hat{\mathbf{B}} – \mathbf{B}_N))\)を足してつじつまを合わせないといけない。第2項、総和をとる範囲を標本から母集団に代えるとともに同じすり替えをしているので、\(\sum_{k \in s} (\mathbf{x}^\top_k (\hat{\mathbf{B}} – \mathbf{B}_N)) / \pi_k \)を引かないといけない。これをまとめたのが第3項になっている] $$ = DIFF(y, \mathbf{x}^\top_k \mathbf{B}_N) + (残り)$$ 残りの部分は、\( \sum_{i=1}^p(\hat{B}_i – B_{N,i})(t_{x_i} – HT(x_i)) \)と書き換えられます。回帰係数の差の部分は \(O_p((N \pi^*_N)^{-1/2})\)であり、\(x_i\)とその全体との差の部分は\(O_p(N(N \pi^*_N)^{-1/2})\)です。これらの速度がすべての\(i\)に等しくかかるとすれば、残りの項は\(N \pi^*_N \rightarrow \infty\)とともに \( o_p(N(N \pi^*_N)^{-1/2}) \)です。[このくだり、学力不足で理解できない…]
したがって、GREGの漸近挙動は差分推定量と似ています。分散不均一回帰モデルの質がどうであろうが、それは漸近不偏で、平均平方一致です。その分散は[…式省略…]と漸近的に等しく、有限母集団の回帰の残差の分散が\(\{y_k\}\)の分散より小さい限り、その分散は\(HT(y)\)の分散より小さいです。
6.1 GREGウェイティングとカリブレーション
GREGの便利な特徴として次の点が挙げられます。$$ DIFF(y, \mathbf{x}^\top_k \hat{B}) = \sum_{k \in s} \frac{y_k – \mathbf{x}^\top_k \hat{B}}{\pi_k} + \sum_{k \in U} \mathbf{x}^\top_k \hat{B} $$ $$ = \sum_{k \in s} \left( \frac{1}{\pi_k} + (t_x – HT(\mathbf{x}))^\top \left(\sum_{k \in s} \frac{\mathbf{x}_k \mathbf{x}^\top_k}{\pi_k} \right)^{-1} \frac{\mathbf{x_k}}{\pi_k} \right) y_k$$ $$ = \sum_{k \in s} \omega_{ks} y_k $$ \( \{ \omega_{ks} \} \)をGREGウェイトといいます。これは元のHTウェイト \( \{ \pi_k^{-1} \} \)を、\(\mathbf{x}\)の情報を使って修正したウェイトです。GREGウェイトは、(\(t_x\)を別にすれば) \( \{\mathbf{x}_k \}_{k \in s}\)にしか依存しません。\(y\)もみていないので多目的サーベイの際に便利です。
[思い出した。このくだり、ウェイティング関連の論文で何度か読んだんだけど、そのたびにキツネにつままれたような気分になるのだ。そして今回もそうだ… 進歩がない…
式の変形をきちんと追いかけていないんだけど、大きな視点でいうと、\(DIFF(y, \mathbf{x}^\top_k \hat{B})\)が\(\sum_{k \in s} w_{ks} y_k\)という形に書き換わるよね、\(w_{ks}\)を定義する式には\(y_k\)は出てこないし、(\(t_x\)をつかってもいいのなら)標本の外側は見に行かなくていいよね、ということさえ理解できればいいのだと思う。
第1項の\(\sum_{k \in s} \frac{y_k}{\pi_k}\)の部分はもちろんこの形になっている。
第1項の\(\sum_{k \in s} \frac{\mathbf{x}^\top_k \hat{\mathbf{B}}}{\pi_k} \)の部分は、よくわかんないけど重回帰モデルのOLS推定量というのは昔から\( \frac{\mathbf{X}^\top}{\mathbf{X}^\top \mathbf{X}} \mathbf{Y}\)だと相場が決まっていて、つまり標本の\(y_k\)のなんらかの重みづけ和なんだから、たとえそれが分散不均一回帰モデルであってもはやりそうなのだろう、つまり結局は\(\sum_{k \in s} w_{ks} y_k\)という形になるのだろう。
第2項が一番不思議で、母集団の総和記号があるんだけど、さっきGREG推定量の変形をしたときと同様に\(\sum_{x \in U} x^k = t_x\)に帰着してしまい、やはり\(\sum_{k \in s} w_{ks} y_k\)という形になるのだろう。
と、むりやり納得して先に進もう…]
GREGウェイトを標本の\(\mathbf{x}\)に適用してみましょう。$$ DIFF(\mathbf{x}^\top, \mathbf{x}^\top_k \hat{B}) $$ $$ = \sum_{k \in s} \omega_{ks} \mathbf{x}^\top_k $$ $$ = \sum_{k \in s} \left( \frac{1}{\pi_k} + (t_x – HT(\mathbf{x}))^\top \left( \sum_{k \in s} \frac{\mathbf{x}_k \mathbf{x}^\top_k}{\pi_k} \right)^{-1} \frac{\mathbf{x}_k}{\pi_k} \right) \mathbf{x}^\top_k$$ $$ = HT(\mathbf{x}^\top) + (t_x – HT(\mathbf{x}))^\top \left( \sum_{k \in s} \frac{\mathbf{x}_k \mathbf{x}^\top_k}{\pi_k} \right)^{-1} \sum_{k \in s} \frac{\mathbf{x}_k \mathbf{x}^\top_k}{\pi_k} $$ $$ = t^\top_x$$ [おおおお…. いや、そりゃ当然そうならないと困るけどさ。でもかっこいいっすね]
つまりGREGウェイトは、ウェイティングされた標本推定値が既知の母集団合計に一致するようにカリブレートされているわけです。
6.2 GREGの近年の拡張
- Cardot & Josserand (2011)は、関数データ(つまりスカラー\(y_k\)のかわりに関数\(y_k(t)\)があるようなデータ)のHT推定について議論しています… [ひええええ。怖ろしい。私はたぶん使わないのでパス]
- Beaumont, Haziza & Ruiz-Gazen (2013)は、影響あるデータ点の効果を取り除いた修正版HT推定量を提案しています。[なるほどね、そういうニーズはあるよね… でもパス]
6.3 GREGの弱点
GREGウェイトはすごく大きくなったり負になったりすることがあります。この弱点を乗り越えることがモチベーションになっている推定手法がたくさんあります。本論文では、特に「回帰っぽい」手法について強調したいと思います。それらの多くは、ウェイトをトリミングしたり平滑化したりします。また、分散不均一回帰でないモデルを使ったりします。たとえば線形混合モデルとか、ノンパラ回帰とか、セミノンパラ回帰とかです。ないし、もっと幅広い統計的手法や機械学習手法を使ったりします。
7. モデル・アシステッド推定のための線形混合モデル
これまで長さ\(p\)の既知の共変量ベクトル\(\mathbf{x}_k\)について考えてきましたが、これに加えて、長さ\(K\)の既知の共変量ベクトル\(z_k\)もあるとしましょう。母集団の行列を\(\mathbf{X}_U = [\mathbf{x}^\top]_{k \in U}, \mathbf{Z}_U = [\mathbf{z}^\top]_{k \in U}\)とします。横につなげて\(\mathbf{C}_U = [\mathbf{X}_U, \mathbf{Z}_U]\)とします。
次の線形混合モデルを作業モデルとしましょう。$$ \mathbf{y}_U + \mathbf{X}_U \beta + \mathbf{Z}_U \mathbf{b} + \epsilon_U$$ $$ E \left( \left[ \begin{array}{c} \mathbf{b} \\ \epsilon_U \end{array} \right] \right) = \mathbf{0} $$ $$ Var \left( \left[ \begin{array}{c} \mathbf{b} \\ \epsilon_U \end{array} \right] \right) = \sigma^2 \left[ \begin{array}{cc} \lambda^{-2} \mathbf{Q} & \mathbf{0} \\ \mathbf{0} & \mathbf{R} \end{array} \right] $$ \(\mathbf{Q}\)は既知で正則、\(\lambda\)は事前に固定とします[なにいい？？？ いったい何を企んでいるの？]。\(\sigma^2\)は未知です。
これはもちろんとても広いモデルクラスであって、重回帰やANOVAはもちろん、多くの縦断モデル、罰則付きスプライン、係数変動モデル、セミパラ加法モデル、低ランク・クリギング、などなどを含みます。
複雑な調査、特に小地域推定においては、線形混合モデルが広く使われています。ウェイトの安定化のためにもよく使われます。ここではGREGの拡張としての線形混合モデルに焦点を当てましょう。
仮に有限母集団全体が観察されていたならば、\(\beta\)のBLUEと\(\mathbf{b}\)のBLUP (最良線形不偏予測子)は以下となります。$$ \mathbf{B}_N = (\mathbf{C}^\top_U \mathbf{C}_U + \Lambda)^{-1} C^\top_U \mathbf{y}_U$$ $$ \Lambda = \mathrm{blockdiag}(\mathbf{0}, \lambda^2 \mathbf{Q}^{-1}) $$ これに対応する\(y_k\)のBLUPは、\( \mathbf{c}^\top_k = (\mathbf{x}^\top_k, \mathbf{z}^\top_k) \)として$$ m_N(\mathbf{c}_k) = \mathbf{c}^\top_k \mathbf{B}_N $$ となります。対応するサーベイ重みづけバージョンは$$ \hat{m}(\mathbf{c}_k) = \mathbf{c}^\top_k \hat{\mathbf{B}} $$ $$ \hat{\mathbf{B}} = \left( \sum_{k \in s} \frac{\mathbf{c} \mathbf{c}^\top_k}{\pi_k} + \Lambda \right)^{-1} \sum_{k \in s} \frac{\mathbf{c}_k y_k}{\pi_k} $$ この\( \hat{m}(\mathbf{c}_k)\)を差分推定量にプラグインすればいいわけです。
得られる推定量は、GREG同様、\(\sum_{k \in s} w_{ks} y_k\)という形式で書けます。ウェイトは\(\mathbf{x}\)の母集団合計に対してカリブレートされますが\(\mathbf{z}\)の母集団合計に対してはカリブレートされません。[なるほど… ウェイティングとしてみると、カリブレーションを緩和したことになるわけか]
——-
ここまでの話をかんたんにまとめておこう。
抽出確率が不均一な調査データに基づき、なんらかの目的変数の、母集団合計を推定する場面について考える(母平均の推定はそのすぐ先にあるわけで、実はすごく応用範囲の広い話をしている)。一番簡単なのはHT推定、つまり、標本の値を抽出確率の逆数で重みづけて合計することだ。でもいかんせん推定量の分散が大きい。なんとかならんのか。
仮に母集団の全メンバーについて、目的変数と関連するなにかの変数がわかっているとする(年齢とか)。そしたら、どうにかして目的変数を年齢で予測するモデルを用意する。別にデータ生成プロセスを反映していなくても、ある程度当たるモデルならなんでもよい。で、母集団の全メンバーについて年齢をつかって目的変数を予測して合計する。さらに、標本の各個体について、年齢による予測と観察された値とのずれを求め、これを抽出確率の逆数で重みづけて合計する。この二つを足す。これが差分推定量だ。
完全なデザイン・ベース推定(デザイン情報しか使わない)でもないし、完全なモデル・ベース推定(目的変数と共変量の関係をガチでモデル化する)でもない。いわばモデルのアシストを受けてデザインベースで推定するアプローチだ。頭いいでしょ？
問題は、目的変数を年齢で予測するモデルをどうやってつくるんだい、という話だ。だって母集団のすべてについて目的変数がわかっているわけじゃないから(もしわかっていたら母集団合計を推定したいと思わない)、そんなモデルを母集団で作れるわけがない。せっかくうまいこと思いついたのにね。残念だったね。
そこで登場するのが、よし、標本でモデルを作っちゃえ、というアイデアである。その代表例が、調査データ解析の本によく出てくるGREG推定量である(モデルは分散不均一回帰モデル)。見ようによっては、HT推定で使うウェイトを、補助変数の標本でのウェイト付き合計が母集団合計(既知)とぴったり一致するようにカリブレートする手法であるともいえる。
しかし、古典的なGREGだけが能じゃない。要は、目的変数を補助変数でもって多少なりともあてられさえすればよいのだ。他にもいろんな可能性があるじゃないの。さあ旅立とう、複雑怪奇なモデリングの世界へと。
以下次号！