仕事の都合で、調査データのウェイティングという複雑怪奇な問題について考える羽目になることがあるのだけれど、なにか役にたつ資料はないものかと探していると、Statistical Scienceの2017年の調査特集号に突き当たることが多い。去年目を通した Elliot & Valliant(2017)、先日読んだBreidt & Opsomer(2017)もこの号に掲載された論文だった。
四の五の言わず、この号に載った9本の論文(残り7本)に片っ端から目を通しちゃえばいいのではないか、と思って、とりあえず(おそらくは編者によるのであろう)イントロをパラパラめくり始めたら、これが意外に面白く、かつ難しく… 結局、最初からメモを取りながら読み直すことになった。
Skinner, C., Wakefield, J. (2017) Introduction to the Design and Analysis of Complex Survey Data. Statistical Science, 32(2), 165-175.
以下ではcomplex survey, complex sampling designといった用語を「複雑調査」「複雑抽出デザイン」という風に直訳する。ここでいるcomplexとは、なんて訳せばいいのかわからずいつも困ってしまうのだけれど、要するに、層別抽出とかクラスタ抽出とかのことです。
1. イントロダクション
[…中略…]
歴史的には、調査での標本抽出は統計学の他のトピックとは別のものだとみなされてきました。調査デザインと調査データ分析をやっていたのはほかの分野とは違う人々でしたし、応用統計学の主流であったモデル・ベース推論とは異なり、調査デザインと結びついた推定手法は多くの場合デザイン・ベース・アプローチ(ランダム化アプローチ)を中心としてきました。デザイン・ベース・アプローチには明確な論拠があるのですが、統計学における慣習的なモデル・ベースの訓練を受けてきた人にとってはわかりにくいことがあります。この分離に拍車をかけてきたのが、ほとんどの統計ソフトに複雑調査手法が含まれていなかったという点です。
しかし、近年では調査での標本抽出と応用統計学的分析手法との間に多大な交配が起きてきたように思います。本論文、そして本特集号の目的は、調査での標本抽出における近年の発展をより広い読者に向けて発信することで、こうした考え方の共有を支援することです。
統計的推論について検討する前に、複雑デザインのいくつかの特徴について概観しておきましょう。
- 我々は確率抽出についてのみ考えます。確率抽出とは、収集されうる潜在的標本を通じたある確率分布を通じて特徴づけられるデザインです。関心ある母集団のそれぞれの単位が非ゼロの被選択確率を持ちます。
個々のデザインを、単純無作為抽出(SRS)からの逸脱という観点からとらえることもできるでしょう。逸脱したくなる動機として、効率的推定、台帳の性質による実務的制約、データ収集コスト、などが挙げられます。実務におけるさまざまなデザインについてはValliant, Dever, & Kreuter(2013 書籍)をみてください。 - […層別単純無作為抽出、一段クラスタ抽出、二段クラスタ抽出、層別多段クラスタ抽出についてごく簡単に説明…] 確率抽出デザインについて、詳しくはTille & Wilhelm (2017, 本特集号)をみてください。
2. 調査データの分析のための推論アプローチ
きちんと定義されたサイズ\(N\)の有限母集団\(U\)があり、単位\(k \in U\)が持っている関心ある調査変数の値が\(y_k\)だとします。標本\(S \subset U\)が選択確率\(p(s)\) で選択されます(\(\sum_s p(s) = 1\))。\(k \in S\)である\(y_k\)しか観察できません。単位\(k\)の抽出確率を\(\pi_k = \sum_{s:k\in S} p(s)\)とし、いわゆるデザイン・ウェイトを\(d_k = \pi_k^{-1}\)とします。おおざっぱにいうと、ウェイトはその単位が「表現している」母集団単位の数です。
たいていの文献では、最初は有限母集団合計についての推論を扱うんですが、ここでは有限母集団平均\(\bar{y}_U = \frac{1}{N} \sum_{k \in U} y_k\)についての推論に焦点をあてましょう。
調査データに基づく推論には大きくふたつのパラダイムがあります。デザイン・ベース・アプローチ(ランダム化ベース・アプローチ)とモデル・ベース・アプローチです。前者はランダム性の源を確率抽出に伴うランダム化と捉え、後者は母集団の値\(y_k\)を生成するモデルだと捉えます。もっとも後者でもランダム化実験についてはランダム化に基づく推論をすることがありますが。
前者についての人気ある教科書として Lohr (2010, “Sampling: Design and Analysis”, 2nd ed.)があります。
2.1 デザイン・ベース・アプローチとモデル・ベース・アプローチ
2.1.1 デザイン・ベース・アプローチ
このアプローチでは、母集団の値\(y_1, \ldots, y_N\)を固定された定数と捉えます。
- 母平均の標準的な推定量はHT推定量です。\(N\)を既知として $$ \bar{y}_{HT} = \frac{\sum_{k \in S} d_k y_k}{N} $$ です。デザイン・ウェイティングをかけているわけですが、その主要な動機はバイアスの除去です。Haziza & Beaumont (2017, 本特集号)をみてください。
- もうひとつの推定量はHajek推定量です。\(N\)が既知かどうかは問いません。$$ \bar{y}_{HJ} = \frac{\sum_{k \in S} d_k y_k}{\hat{N}} $$ ただし\(\hat{N} = \sum_{k \in S} d_k\)です。
バイアスやそのほかの積率は、\(U\)からの繰り返し抽出という観点から評価します。期待値を\(E_s[\bar{y}_{HT}]\)という風に、また分散を\(var_s(\bar{y}_{HT})\)という風に書きましょう。添え字\(s\)はありうる下位集合を通じた平均であることを強調するためにつけています。
インフォーマルな定義ですが、推定量についての次の2つの基準について定義しておきましょう。
- デザイン不偏: 期待値が真値と等しいこと
- デザイン一致: 標本サイズが大きくなるにつれて、デザインのバイアスと分散がゼロに近づくこと
Hajek推定量はバイアスがありますがデザイン一致です。後述する理由で、\(N\)が既知であってもHT推定量よりHajek推定量が好まれることがあります。Hajek推定量をみると、伝統的な調査標本抽出の驚くべき側面に気づきます。比の推定へのこだわりです。
デザイン・ベース推定量の性質を導出する際、カギとなるのは選択の二値インジケータ\(I_k\)を定義することです。\(E_s[I_k] = \pi_k, var_s (I_k) = \pi_k(1-\pi_k)\)となります。\(E_s[I_k I_l] = \pi_{kl}, cov_s(I_k, I_l) = \pi_{kl} – \pi_k \pi_l = \Delta_{kl}\)とします。
多くの場合、抽出は非復元で行われ、\(\pi_{kl} \neq \pi_k \pi_l\)です。この点はモデル・ベース・アプローチとは対極的です(モデル・ベース・アプローチでは値は仮説的な無限母集団から抽出されるので、ふつう独立だと仮定されます)。
では、HT推定量のバイアスについてみてみましょう。$$ E_s[\bar{y}_{HT}] = \frac{1}{N} E_s \left[ \sum_{k \in S} d_k y_k \right] = \frac{1}{N} E_s \left[ \sum_{k \in U} d_k I_k y_k \right] $$ $$ = \frac{1}{N} \sum_{k \in U} \pi^{-1}_k E_s [I_k] y_k = \bar{y}_U $$ このようにデザイン不偏です。\(I_k\)を導入して\(U\)を通じて総和をとる、というトリックに注目してください。
このデザイン不偏性は逆確率ウェイティングというテクニックに由来しており、\(\pi_k \gt 0\)が必要となります。直感的にもわかりやすいですね。絶対抽出されない単位があったら不偏推定量は手に入らないでしょう。
HT推定量の分散は$$ var_s(\bar{y}_{HT}) = \frac{1}{N^2} var_s \left( \sum_{k \in S} d_k y_k \right) = \frac{1}{N^2} var_s \left( \sum_{k \in U} d_k I_k y_k \right) $$ $$ = \frac{1}{N^2} \sum_{k \in U}\sum_{l \in U} d_k d_l cov_s (I_k, I_l) y_k y_l $$ $$ = \frac{1}{N^2} \sum_{k \in U} \sum_{l \in U} \Delta_{kl} \frac{y_k}{\pi_k} \frac{y_l}{\pi_l} $$ [あっれえ？ 恥ずかしながら、2行目に行くところがわからない…]
分散の不偏推定量はこうなります: $$ \hat{var}_s (\bar{y}_{HT}) = \frac{1}{N^2} \sum_{k \in S} \sum_{l \in S} \frac{\Delta_{kl}}{\pi_{kl}} \frac{y_k}{\pi_k} \frac{y_l}{\pi_l} $$ ここでカギになるのは\(\pi_{kl} \gt 0\)であるということです。
上に示した分散とその不偏推定量は、モデル・ベース陣営の人にはミステリアスに見えるかもしれません。\(y_k\)の分散に依存していないからです。
上に示した分散の不偏推定量は(そして関連するSen-Yates-Grundy推定量も)いくつかの欠点を持っています。負になりえますし、\(\pi_{kl}\)は未知であることが多いです。そこで、デザインを「復元あり」抽出だと見立てて推定したり(Loarの教科書を見てください)、ジャックナイフ法やブートストラップ法のようなリサンプリング法で推定したりします。
上記の表現のうちいくつかについて例示するために、次のようなSRSを考えてみます。\(s\)が\(n\)個の単位を持っていたら\(p(s) = C(N, n)^{-1}\), そうでなかったら\(p(s)=0\)とします。\(\pi_k = n/N, d_k = N/n, \pi_{kl} = \frac{n}{N}\frac{n-1}{N-1}\)です。このとき、以下が示せます。$$ \hat{N} = \sum_{k \in S} d_k = N $$ $$ \bar{y}_{HT} = \bar{y}_{HJ} = \sum_{k \in s} \frac{y_k}{n} $$ $$ var_s(\bar{y}_{HT}) = \left( 1-\frac{n}{N} \right) \frac{S^2_y}{n} $$ ただし\(S^2_y = \frac{1}{N-1} \sum_{k \in U}(y_k – \bar{y}_U)^2 \)です。\(1-\frac{n}{N}\)を有限母集団修正といいます。\(S^2_y\)のデザイン不偏推定量は$$ s^2_y = \frac{1}{n-1} \sum_{k \in S}(y_k – \bar{y}_{ht})^2 $$ です。このように、推定量の分散の式は(有限母集団修正を別にすれば)おなじみの形になるわけです。
[この段落の文意がつかめない。なぜここでこんな話をするんだろう？ 有限母集団であってもSRSだったら推定量の分散の不偏推定量が\(y_k\)の分散に依存する形になるよ、ということを云っているの？ ]
2.1.2 モデル・ベース推論
このアプローチでは、\(y_k\)が確率変数\(Y_k\)の実現値だと考えます。\(Y_k\)はなんらかのモデルに従っていて、母集団は仮説的な無限の超母集団からの抽出だと考えるのです。
伝統的な統計的モデリングでは、関心があるのは有限母集団の特性というよりモデル・パラメータですが、有限母集団の特性についての推論もモデリングの枠組みで考えることができます。関心の対象を\(\bar{Y}_U = \frac{1}{N} \sum_{k \in U} Y_k\)としましょう。これは確率変数なので、これから考える\(\hat{\bar{Y}}\)のことを、\(\bar{Y}_U\)の推定量ではなくて予測量と呼びましょう。
\(\hat{\bar{Y}}\)の評価の際の主な基準として次の2つがあります。(1)バイアス \( E_M[\hat{\bar{Y}} – \bar{Y}_U] \)。(2)モデルに関する分散\(E_M[ (\hat{\bar{Y}} – \bar{Y}_U)^2] \)。なお、ここでは頻度主義の枠組みで考えています。ベイジアンの枠組みなら、\( \bar{Y}_U \)の事前分布に基づく推論を行うところですが、ここでは触れません。Gelman(2007), Little(2013, J.Off.Stat.)をみてください。
予測アプローチの例として、以下のモデルについて考えましょう。$$ \mu = E_M(Y_k), \ var_M(Y_k) = \sigma^2 $$ 予測量は標本平均 $$ \hat{\bar{Y}}_n = \frac{1}{n} \sum_{k=1}^n Y_k $$となります。これは不偏予測量です[数式略]。予測の分散は、有限母集団の分散を\(\sigma^2\)として $$ E_M [ ( \hat{\bar{Y}}_n – \bar{Y}_U )^2] = E_M \left[ \left( \frac{1}{n} \sum_{k=1}^n Y_k – \frac{1}{N} \sum_{k=1}^N Y_k \right)^2 \right] $$ $$ = \left( 1 – \frac{n}{N} \right) \frac{\sigma^2}{n} $$ \(\sigma^2\)に\(S^2_y\)を代入すれば、SRSのデザイン・ベース推論の時と同じ分散になります。
モデル・ベース・アプローチでは、複雑抽出スキーマがなんの役割も果たしていないように見えるかもしれません。しかし、それは次の役割を持ちます。
- 複雑抽出スキーマは母集団の構造(たとえば層別やクラスタリング)に依存するでしょう。モデル・ベース推論ではこの構造をモデルで捉えることが本質的に重要です。
- ふつうモデル・ベース推論では母集団レベルのモデルが観察にも適用されると仮定します。これは暗黙的ですが強い仮定です。モデルベース・アプローチの下での予測は\(I_k\)に条件づけられていますから、この仮定は、\(I_k =1\)のもとでの\(Y_k\)の条件付き分布が\(I_k = 0\)の下での\(Y_k\)の条件付き分布と同じだということを意味します。このとき、抽出は非情報的であるといいます。この仮定が成り立たないときに選択バイアスが生じえます。確率抽出の主たる利点は\(I_k\)と\(Y_k\)の独立性の保証によって選択バイアスを防ぐことにあります。その保証がなされていない標本から推定を行うのはリスキーです。
[うーん、なにをいってんだかわかりにくい説明だなあ。要するに、\(I_k\)が\(Y_k\)と独立でないときは、\(I_k = 1\)のときしか\(Y_k\)を観察できないということを\(Y_k\)のモデルに組み込まないといけないよ、ということでしょうか]
2.2 パラダイムのスイッチング
[この節は非常に勉強になった。その反面、私にとっては難解なところが多いので、丁寧にメモする。勝手に小見出しをつける]
2つのアプローチを比較するために、デザイン・ベース推定量をモデル・ベース基準で検討したり、その逆を検討したりしてみましょう。
(モデル・ベースの観点から見たデザイン・ベース推定量のバイアス)
まず、HT推定量とHajek推定量の、モデル\(\mu = E_M[Y_k]\)の下でのバイアスについて検討します。推定量の式に入っている\(y_k\)を\(Y_k\)に置き換えます。$$ E_M[\bar{y}_{HT} – \bar{Y}_U] = \left( \sum_{k=1}^n d_k – N \right) \frac{\mu}{N} $$ $$ E_M[\bar{y}_{HJ} – \bar{Y}_U] = 0$$ [なに！？ なにが起きたの！？
深呼吸して、モデルベース陣営の発想に切り替えよう。\(Y_1, \ldots, Y_n\)は観察されていようがいまいが確率変数である。愚かなデザインベース陣営の奴らは\( \frac{1}{N} \sum_{k=1}^n d_k Y_k \)をHT推定量\(\bar{y}_{HT}\)と称している。推定対象\(\bar{Y}_U\)に対するそのバイアスについて考えると、$$ E_M[\bar{y}_{HT} – \bar{Y}_U] = E_M \left[\frac{1}{N} \sum_{k=1}^n d_k Y_k – \frac{1}{N} \sum_{k \in U} Y_k \right] $$ 確率変数は\(Y_k\)だけなので$$ \frac{1}{N} \sum_{k=1}^n d_k E_M[Y_k] – \frac{1}{N} \sum_{k \in U} E_M[Y_k] $$ \(E_M[Y_k] = \mu\)より $$ = \frac{1}{N} \mu \sum_{k=1}^n d_k – \frac{1}{N} N \mu $$ $$ = \left(\sum_{k = 1}^n d_k – N \right) \frac{\mu}{N} $$ 次に、あいつらは\( \frac{\sum_{k=1}^n d_k Y_k}{\sum_{k=1} d_k} \)をHajek推定量\(\bar{y}_{HJ}\)と称している。そのバイアスについて考えると $$ E_M[\bar{y}_{HT} – \bar{Y}_U] = E_M \left[ \frac{\sum_{k=1}^n d_k Y_k}{\sum_{k=1} d_k} – \frac{1}{N} \sum_{k \in U} Y_k \right] $$ \(E_M[Y_k] = \mu\)より $$ = \frac{\sum_{k=1}^n d_k}{\sum_{k=1}^n d_k} \mu – \frac{1}{N} N \mu = 0$$ となる… ってことでしょうか]
このように、このモデルのもとで、HT推定量は\(\sum_{k=1}^n d_k = N\)のときしか不偏でありませんが (たとえばSRSのときがそうです)、Hajek推定量は常に不偏です。
(デザイン・ベースの観点から見たモデル・ベース推定量のバイアス)
今度は、モデル不偏推定量\(\hat{\bar{Y}}_n\)のデザインバイアスについて考えましょう。式の中の\(Y_k\)を\(y_k\)に置き換えます。
[再び深呼吸して、デザインベース陣営の発想に切り替える。母集団の\(y_1, \ldots, y_k\)は定数だ。モデルベース陣営の馬鹿どもは\(\frac{1}{n} \sum_{k \in S} y_k \)をモデル不偏推定量\(\hat{\bar{Y}}_n \)などと呼んでいる。その期待値について考えよう] $$ E_s[\hat{\bar{Y}}_n] = E_s \left[\frac{1}{n} \sum_{k=1}^N I_k y_k \right] = \frac{1}{n} \sum_{k=1}^N \pi_k y_k $$ これは一般的にはデザイン不偏ではありません。Haziza & Beaumont (2017, 本特集号)が示しているように、これがデザイン不偏になるのは、\(\pi_k\)と\(y_k\)が無相関であるときのみです。モデルベースの言葉でいえば、抽出が無情報であるときです。このことは、モデルの誤指定によってもたらされる(デザイン)バイアスについて示しています。ウェイトなしの推定量\(\hat{\bar{Y}}_n\)は、モデル\(\mu = E_M[Y_k]\)の下で(そして無情報抽出という仮定のもとで)、\(\bar{Y}_U\)の不偏推定量です。しかし、これらの仮定が誤っていた時、バイアスから守られていません。これがデザイン・ウェイトつきの推定量とのちがいです。
(モデル・ベースの観点から見たデザイン・ベース推定量の分散)
今度は分散について考えてみましょう。モデルを\(var_M(Y_k) = \sigma^2_k\)とし、\(Y_k\)と\(Y_l\)はモデルの下で独立だとします。HT推定量について考えると$$ E_M[ (\bar{y}_{HT} – \bar{Y}_U)^2 ] = \sum_{k=1}^n (N^{-1} d_k – N^{-1})^2 \sigma^2_k + \sum_{k = n+1}^N N^{-2} \sigma^2_k $$ となります[ここはわからんので写経した]。あきらかに、デザインベースの観点からみた\( var_s(\bar{y}_{HT}) \)とは異なります。なお、Hajek推定量の場合は$$ E_M[ (\bar{y}_{HT} – \bar{Y}_U)^2 ] = \sum_{k=1}^n (\hat{N}^{-1} d_k – N^{-1})^2 \sigma^2_k + \sum_{k = n+1}^N N^{-2} \sigma^2_k $$ となります。
デザイン・ベース推定量のモデル分散について検討することは、その効率性を評価する助けになります。たとえば次のモデルについて考えてみましょう:$$ Y_k = \theta \pi_k + \epsilon_k $$ 誤差項\(\epsilon_k\)は\(E_M[\epsilon_k] = 0, var_M(\epsilon_k) = \pi^2_k \sigma^2\)に独立に従うとします。WLS推定量\(\hat{\theta}\)は、$$ \sum_{k=1}^n \frac{(y_k – \theta \pi_k)^2}{\pi^2_k \sigma^2} $$を最小化する推定量ですが、これはHT推定量\(\bar{y}_{HT}\)に対応します。上記のモデルは正しくなくてもHT推定量の特性は成り立ちますが、このことは、推定量がうまくいくのはどんな状況かということを教えてくれます[\(Y_k\)と\(\pi_k\)が関連しているときにHT推定量はWLS推定量になるよということね。]。これをモデル・アシステッド・アプローチとみなすこともできます。[??? どういうことだろうか]
(Basuの象)
\(y_k\)が\(\pi_k\)と関連していないときにHT推定量がうまくいかなくなるという極端な例として、有名なBasuの象の例があります。サーカスのオーナーが50頭の象の重さを、一頭だけ体重を測って推定しようとします。3年前の記録ではSamboの体重が平均的体重だったので、Samboの体重を測ってそれを50倍しようとします。サーカスの統計家はデザイン不偏推定量にこだわり、いや、Samboを\( \pi_k = 99/100\)で、残りの49頭を\(\pi_k = (1/100) \times (1/49) \)で抽出するという計画を提案します。選ばれた象の体重を\(y\)として、HT推定量はもしSamboが選ばれたら\( (100/99) y\), そうでなかったら\( (100 \times 49) y\)だということになります。これは明らかに変です。\(y_k\)が\(\pi_k\)に比例していないというのがこの話のポイントです。
(Hajek推定量が好まれる理由)
モデリングの観点からみると、\(E_M[Y_k] = \theta, var_M(y_k) = \pi^2_k \sigma^2\)というモデルからHajek推定量が得られます。ここから、反応が(抽出確率と比例するのではなくて)ほぼ一定のとき、HT推定量よりHajek推定量のほうがよいことが期待されます。これが、多くの事例でHajek推定量が勧められる理由です。
(まとめ)
まとめましょう。モデル・ベース推定量をデザイン・ベースの観点から見たり、その逆を考えたりすることは有益です。推定量がうまくいくのはどんな状況かについての洞察が得られるからです。本節で得られた結論をいくつかおさらいしましょう: HT推定量はデザイン不偏だが特定のモデルの下ではバイアスを持ちます。モデル分散はデザイン分散と対応しないことがあります。
3. 補足変数とモデル・アシステッド推定
モデル・ベースの観点からは、モデルの共変量に補足変数(母集団単位についての補足情報)を含めるのは自然だし当然です。デザインで使われた補足変数で条件づけることで、抽出を非情報的にすること(つまり共変量のもとで\(Y_k\)と\(I_k\)を独立にすること)ができます。
デザインベースの観点からも、\(Y_k\)の回帰モデルの共変量としてなにが適切かを考えることで、推論を「アシスト」することができます。こういうモデル・アシステッド推論の教科書の決定版としてSarndal, Swensson, & Wretman (1992) があります。
話を単純にするため、\(y_k\)と関連すると我々が信じている単一の変数\(x_k\)があり、その平均\(\bar{x}_U\)を我々は知っているとしましょう。「作業モデル」として下式を考えます。$$ Y_k = B_0 + B_1 x_k + \epsilon_k$$ $$ B_1 = \frac{ \sum_{k=1}^N (x_k – \bar{x}_U)(y_k – \bar{y}_U) }{ \sum_{k=1}^N (x_k – \bar{x}_U)^2 } = \frac{R S^2_y}{S^2_X} $$ $$ B_0 = \bar{y}_U – B_1 \bar{x}_U$$ \(\epsilon_k\)は独立に\(E_M[\epsilon_k] = 0, var_M(\epsilon_k) = \sigma^2 \)に従うとします。
作業モデルの\(B_0+B_1 x_k\)の部分は\(\bar{y}_U\)の推定量のつもりです。\(B_0, B_1\)のデザインベース推定量は$$ \hat{B}_1 = \frac{ \sum_{k \in S} (x_k – \bar{x}_S)(y_k – \bar{y}_S) }{ \sum_{k \in S} (x_k – \bar{x}_S)^2 } = \frac{r s^2_y}{s^2_X}$$ $$ \hat{B}_0 = \bar{y}_S – \hat{B}_1 \bar{x}_S$$ 比がはいっているのでデザイン不偏ではありませんがデザイン一致ではあります。結局、\(\bar{y}_U\)の推定量はこうなります: $$ \hat{B}_0 + \hat{B}_1 \bar{x}_U = \bar{y}_S + \hat{B}_1(\bar{x}_u + \bar{x}_s) $$ これが伝統的な「回帰推定量」です。SRSの下ではデザイン一致です。単に\(\bar{y}_S\)を使うよりも精度が上がります。その分散は[…中略…] \(y_k\)と\(x_k\)の相関の二乗に応じて減少します。
一般的な複雑デザインでは、\(B_0, B_1\)の推定にデザイン・ウェイトを使います。Breidt & Opsomer (2017, 本特集号)が論じているように、これはGREGモデルにアシストされた推定量$$ \bar{y}_{GREG} = \hat{B}_0 + \hat{B}_1 \bar{x}_U = \bar{y}_{HT} + \hat{B}_1 (\bar{x}_U – \bar{x}_{HT}) $$ となります。HT推定量が、\(x_k\)の母平均とその標本推定量との差によって調整されているわけです。
GREG推定量はデザイン一致です。バイアスはありますが、大標本ではHT推定量より精度が高くなります。\(B_0\)をゼロにすれば比推定量になります。
モデル・アシステッド・アプローチは、サンドイッチ分散推定量を使ったロバスト推定と似ているところがあります。作業モデルを指定しますが、一致性は非常に弱い想定の下で保障されておりモデルに依存しません。
GREG推定量や関連する推定量は重みつき推定量として表現できます。いろいろな調整が可能ですし、ウェイトの構築は複雑になることがあります。Haziza & Beaumont (2017 本特集号), Chen et al.(2017 本特集号)をみてください。
4. モデル・パラメータ
本節では、有限母集団の平均の推定の延長ではとらえられないような、モデル・パラメータの推定に焦点をあてます。詳細はLumley & Scott (2017 本特集号)をみてください。
- モデル・ベースの観点からは、調査変数のモデルが関心あるパラメータ\(\theta\)だけでなく、デザインに用いられた補足変数も含んでいるとき、\(\theta\)の推論においては抽出スキーマは無視可能になり、標準的なウェイトのないアプローチを使用できます。この場合、抽出デザインも有限母集団も無視できて、単に標本データを生成すると思われるモデルについて考えるだけでよくなります。その一方で、デザイン変数で条件づけることでモデルは関心の範囲から逸れかねません。Skinner, Hold, & Smith (1989 書籍)はデザイン変数に条件づける分析を「非累積」分析、条件づけない分析を「累積」分析と呼び、この2つが全く異なる分析上の目的を持つことを示しています。[へー。デザインを無視して分析するのにもそれなりの目的がありうるってことかしらん]
- デザイン・ベースの観点からは、そもそも関心あるパラメータをどうやって定義すんの、というのが問題になります。一般的なアプローチは、\(\theta\)を含む(超母集団)モデルを指定することなのですが、モチベーションとしてはそれは単なる「作業モデル」であって推論に使うものではありません。この目的を心にとめつつ、「センサス・パラメータ」\(\theta_U\)なるものを定義します。それは(\theta\)のなんらかの推定量で、もしも有限母集団全体が仮想的なセンサスによって観察されたら得られるはずのものです。
たとえば、失業のモデリングに関心があり、\(\theta\)はある時点にある労働者が失業している確率を表すパラメータだとします。\(\theta_U\)はこの母集団における実際の失業率で、モデリング上の仮定の下で、大標本ならば\(\theta\)を近似します。\(\theta\)ではなく\(\theta_U\)が関心の対象だ、と考えるわけです。[ああ、なるほどね、言っている意味がわかったよ]
この考え方の魅力は、\(\theta_U\)が有限母集団の量であり、デザインベース推論の直接的な関心の対象であるという点です。いっぽう、具体的な推定アプローチの観点から見ると、センサス・パラメータの定義はいささか恣意的です。やはり多くの場合は、\(\theta\)を目標とみるのが好ましいでしょう。その場合でも、\(\theta_U\)として使うつもりだったのと同じ推定量を\(\hat{\theta}\)として使うのは合理的ですが、分散推定やそのほかの推定手続きにおいては、デザイン・ベース推論とモデル・ベース推論を組み合わせた修正が必要になります。
たとえば次の場合について考えましょう。超母集団の平均\(\mu = E_M[Y]\)についての推論を行いたいです。有限母集団平均\(\bar{y}_U\)をセンサス・パラメータとみなし、デザイン・ベースの観点から、HT推定量\(\bar{y}_{HT}\)を\(\bar{y}_U\)と\(\mu\)の両方の点推定量として用います。これは\(\bar{y}_U\)に対してデザイン不偏であるだけでなく、モデルベースとデザインベースを組み合わせた観点からみて、\(\mu\)に対して不偏でもあります。つまり、(モデルベース的に\(\bar{y}_U\)を\(\bar{Y}_U\)と書き換えて) $$ E[\bar{y}_{HT}] = E_M[E_s[\bar{y}_{HT}]] = E_M[\bar{y}_U] = \mu $$ です。さて、分散について考えましょう。不確実性は、サイズ\(n\)の標本の選択だけでなく、サイズ\(N\)の母集団の選択においても生じます。話を単純にするためにSRSだとして、$$ var(\bar{y}_{HT}) = E_M[var_s(\bar{y}_{HT})] + var_M(E_s[\bar{y}_{HT}]) $$ […中略…] $$ = \frac{\sigma^2}{n} $$ この例では、分散の式はあたかも超母集団からの無作為標本のようになります。さらに、\(N\)が\(n\)より十分に大きければ第二項は無視できるので、デザイン・ベースの分散推定でもよいことになります。[めんどくさい議論のわりに結論はあたりまえだ…]
モデル・パラメータの推定についてのより一般的な議論として、Graubard & Korn (2002 Stat.Sci.)をみてください。
5. 非確率抽出と無回答
ここまで無視してきましたが、たいていの調査では無回答が生じます。さらにいえば、最近は非確率抽出も増えてます。Schonlau & Couper (2017 本特集号)をみてください。
ここでは(項目無回答ではなくて)単位無回答について考えます。
無回答と非確率抽出はどちらも、デザイン上コントロールされていない標本選択を含んでおり、モデリング上の想定が必要になります。2つのアプローチを紹介します。詳しくはElliott & Valliant(2017 本特集号)をみてください。
- 準ランダム化。あたかも確率抽出で得られたかのように標本を表現することを目指します。無回答の場合でいうと、無回答を抽出の第二段階とみなし、デザイン・ベースのウェイトを作ったりします。他の確率標本を使って疑似ウェイトを作ることもあります。
- 超母集団モデル。それについて条件づければ標本選択が非情報的になるような補足情報を探します。
6. 分散推定
複雑調査データによる区間推定にはふつう漸近不偏性を用います。[…中略…] Breidt & Opsomer (2017 本特集号)をみてください。そこでカギになるのは、適切な点推定量とその分散推定量を同定することです。
多くの場合、ウェイト構築は調査主体が行い、それを推定にどう組み込むかは分析者が考えます。分散推定についてもこのような分離ができないでしょうか。データにすべての\(\pi_{kl}\)をつければよいのですが、ふつうは無理です。
- ひとつのアプローチとして、現実のデザインを、復元抽出を行う類似のデザインで近似するというやりかたがあります。層別多段抽出で一般的なアプローチです。層、PSU、ウェイトがあれば分散推定できます。Valliant, Dever, & Lreuter (2013 書籍)をみてください。
- 複雑な非線形推定量に関しては線形化という手法が用いられています(主流の統計学でいうところのデルタ法です)。ほかにも近似方法があります。Berger & Tille (2009 “Handbook of Statistics”)をみてください。
- もう一つのアプローチはreplicationによる分散推定です。ブートストラップとジャックナイフが有名です。[… めんどくさいのでメモ省略…]
- 最後に、モデルベースのアプローチもあります。たとえばクラスタを混合モデルで表現し、分散はサンドイッチ推定します。[…]
7. 結び
触れられなかった話題として異なるデータソースの結合があります。Lohr & Raghunathan (2017 本特集号)をみてください。
[…]
—————
細部の説明がちょっと雑でイラっとしたところもあるんだけど、特集号のイントロだから仕方がないだろう。
2節後半の、デザインベース推定量(Horvitz-Thompson推定量とHajek推定量)をモデルベースの観点から評価する、ないしその逆、というところがとても勉強になった。デザインベースでみると不偏推定なのにモデルベースでみるとそうでない、なんてことがあるのね。ここの話をもう少し詳しく勉強したいんだけど、なにをみればいいのだろうか。