Little, R.J.A. (1993) Post-Stratification: A Modeler’s Perspective. Journal of the American Statistical Association, 88(423), 1001-1012.
統計学者Little先生による、調査における事後層別(post-stratification)についての解説。伝統的なデザイン・ベースの説明ではなく、モデル・ベースの観点から説明するというのがミソである。
良く引用される論文で、前から気になっていたのだがずるずると後回しになっていた。このたびちょっと調べ物をしていて、その流れで思い出し、試しに読んでみた。別にいま読まなくてもいいっちゃいいんだけど、この一週間体調を崩して寝込んでいたので、そのリハビリを兼ねている。
途中で気が付いたけど、これ、招待講演を元にした論文なんですね。講義を思わせるちょっとカジュアルな書きぶりである。
1. イントロダクション
たとえば人口調査はふつう年齢で層別抽出できない(年齢はインタビューしてみないとわかんない)。でもセンサスデータから年齢の集計値をもってきて事後層別することはできる。標本サイズを\(r\), 事後層\(h\)の標本サイズを\(r_h\), センサスから得た層の人口割合を\(P_h\)として、ウェイトは\(w_h = r P_h / r_h\)となる(ウェイトの合計が\(r\)になるようにスケーリングしている)。なお、臨床試験では層別なしのランダム化デザインで得たデータを層別して分析することを事後層別というけど、この論文で扱うのはサーベイの話ね。
事後層別によってバイアスを取り除いたり制度を増大させたりできる。非常に広く用いられており、多くの公的調査において重要な役割を果たしている。しかし、事後層別の統計的な諸特性にはあまり注意が向けられていない。
事後層別についての従来の研究はランダム化の観点からなされていた。つまり、母集団の値は固定で、抽出分布を通じた推定量のバイアスとMSEを問題にしてきた[デザインベースってことね]。この観点から見れば事後層別というのはウェイティング調整のひとつである。
本論文では事後層別を予測モデリングの観点から捉える。つまり、母集団の値はモデルのもとでの確率変数で、有限母集団特性についての推論をモデルのもとでの予測分布に基づいて行うという見方である。事後層別とは、事後層別変数\(Z\)が母集団の全個体について既知で、そのうち\(n\)個だけについて調査変数を測定できるとき、欠損をどう埋めるか、という問題である。
2. 単一の周辺への事後層別
2.1 事後層別平均
[ここは、まあ、おさらいのパートである。気が楽でよい]
調査変数を\(Y\), 事後層別変数を\(Z\)、事後層\(h\)の有限母集団平均を\(\bar{Y}_h\)として、\(\bar{Y} = \sum_h P_h \bar{Y}_h\)の推論について考えよう。サイズ\(n\)の単純無作為標本から回答者\(r\)人を得た。\(\mathbf{P} = \{ P_1, \ldots, P_H\}, \mathbf{n} = \{n_1, \ldots, n_h\}, \mathbf{r} = \{r_1, \ldots, r_H\}\)とし、\(\mathbf{P}, \mathbf{r}\)を既知とする。なお当面は\(\mathbf{n}\)はどうでもよい。
無回答は無視可能、つまり、事後層\(h\)の回答者は標本の無作為下位標本とみなせるとしよう。形式的に書くと $$ r_h | \mathbf{n} \sim bin(n_h, \varphi_h) $$ である。もしすべての\(h\)について\(\varphi_h = \varphi\)なら、これはMCARであって、無回答バイアスの除去は不要となる(分散を下げるためになお事後層別してもいいけど)。MCARでない場合、事後層別によってバイアスを取り除ける。
\(\bar{Y}\)の推定量としてはふつう以下が用いられる: $$ \bar{y}_{ps} = \sum_h^H P_h \bar{y}_h = \frac{1}{r} \sum_i^r w_i y_i$$ ただし\(w_i\)は、\(i\)が\(h\)に属しているとして\(w_i = r P_h / r_h\)ですね。
2.2 事後層別平均のランダム化分散
さて、\(\bar{y}_{ps}\)のランダム化分散として適切なのはなにか[いわゆるデザイン分散のことね]。これは論争的話題である。
無条件な標本抽出分散は$$ var(\bar{y}_{ps}) = E[var(\bar{y}_{ps} | \mathbf{r})] + var(E[\bar{y}_{ps} | \mathbf{r}]) $$ 2項目の\(var()\)の内側、\(E[\bar{y}_{ps} | \mathbf{r}]\)は定数\(\bar{Y}\)であるから、第2項はゼロである。というわけで、第1項をどう近似するかを考える。Cochranの教科書とかはそうなっている。
いっぽう、Holt & Smith(1979 JASA)いわく、もし\( \{ \mathbf{r} \}\)が期待値からすごく離れたらどうするんだ[えーと、ある層の無回答が極端に多いといった悲しい事態までを通じた分散を求めるの？ それって悲観的すぎるでしょ、ってことかな]。むしろ、条件付き分散を考えたほうが良い。事後層における\(Y\)の分散を\(S^2_h\)として $$ var(\bar{y}_{ps} | \mathbf{r}) = \sum \frac{P^2_h(1-f_h) S^2_h}{r_h} $$ ただし\(1-f_h\)は有限母集団修正。この式は層別平均の分散の通常の表現である。
両者の違いは大標本ならたいしたことがないんだけど、事後層が多くなると無視できない。
このことから次の点がわかる。ランダム化アプローチというのは、関心ある量についてのの直接的情報はもたないがその精度を示しているような補足変数[ここでは\(\{\mathbf{r}\}\)ね]を条件づけるべきかという点であいまいである。実のところ、\(\{\mathbf{r}\}\)は直感的には補足変数だが、調査のランダム化理論における補足変数についてのフォーマルな理論はあるのか。私は寡聞にして聞いたことがない。
[Little先生、デザイン・ベース陣営をディスるの巻。要するに、\(\bar{y}_{ps}\)の分散として、無条件な分散と\(\mathbf{r}\)で条件づけた分散を考えることができる、ってことね。これがデザイン・ベース・アプローチの欠点だというのは反論がありうるんじゃないかと思うぞ]
2.3 事後層があるときの\(\bar{Y}\)についてのベイジアン推論
\(\bar{Y} = \sum P_h \bar{Y}_h\)についてのベイジアン推論のためには、\(Z\)のもとでの\(Y\)のモデルが必要である。\(\bar{Y}\)の事後分布は事後平均$$ E(\bar{Y}|data) = \sum P_h E(\bar{Y}_h | data) $$ で要約される。これはランダム化推論における推定量と似ている。精度評価においては事後分散 $$ var(\bar{Y}|data) = \sum P^2_h var(\bar{Y}_h | data) + \sum \sum_{h \neq k} P_h P_k cov(\bar{Y}_h, \bar{Y}_k | data) $$ が重要である。
基本的正規事後層別モデル(BNPM)$$ (y_i | z_i = h, \mu_h, \sigma^2_h) \sim_{ind} G(\mu_h, \sigma^2_h) $$ $$ p(\mu_h, \log \sigma_h) = const$$ のもとでの推論について考えよう。なお、\(z_i\)とは\(i\)が属する事後層、\(G(a,b)\)は平均\(a\), 分散\(b\)の正規分布、\(p(\mu_h, \log \sigma_h)\)は無情報ジェフリーズ事前分布である。正規性を仮定したけどそれはたいした問題じゃない。ポイントは、それぞれの事後層が別の平均と分散を持っていて、層の中ではiidだということである。
このモデルの下で、\(\bar{Y}\)の事後分布はt分布の混合分布になる。平均と分散は$$ E(\bar{Y} | \mathbf{Z}, \mathbf{Y}_s) = \bar{y}_{ps} $$ $$ var(\bar{Y} | \mathbf{Z}, \mathbf{Y}_s) = \upsilon_{ps} = \sum \frac{P^2_h (1-f_h) \delta_h s^2_h}{r_h} $$ ただし、\(s^2_h\)は\(Y\)の標本分散。\(\delta_h = \frac{r_h-1}{r_h-3}\)は小標本修正。\(\mathbf{Z}\)は\(Z\)についての標本とセンサスのデータ。\(\mathbf{Y}_s\)は標本の\(y_i\)のベクトルで長さ\(r\)である。
Holt & Smith (1979)がいうように、このモデルによる分析は条件づけられた分散を支持する。式を見比べると、確かに条件づけられた分散の形になっている。
2.4 ウェイトなし平均との比較
事後層別平均は大標本では効率性が高いけど、小標本では不安定だし、事後層に回答者がいなかったら定義できなくなる。そこで、\(w_i = r P_h/r_h\)をトリミングなり平滑化なりしようということになる。極論をいえば、すべて1に平滑化すればウェイトなし平均\(\bar{y}\)になる。
\(\bar{y}_{ps}\)と\(\bar{y}\)の比較に際しても、\( \{ \mathbf{r} \} \)で条件づけるかどうかという疑問が生じる。$$ var(\bar{y}_{ps} | \mathbf{r}) = \sum P^2_h(1-f_h) S^2_h / r_h$$ と、母分散を\(S^2\)とし\(f = r/N\)として $$ var(\bar{y}) = (1-f) S^2/r$$ を比べるのは変である。だって後者では\( \{ \mathbf{r} \}\)は確率変数だから。[ああ、そりゃそうだな]
モデリングの観点からいうと、ウェイトは小標本での事後層別という問題にとっての原因ではなくて症状である。真のポイントは、先ほど定式化した基本的正規事後層別モデルが、\(r_h\)が小さい事後層で\(Y\)の分布をうまく予測できないという点にある。\(\bar{Y}\)の事後分布についての推論の仕方を修正しないといけないのである。
正攻法は、\(Z\)のもとでの\(Y\)の分布について修正したモデルを作り、そのもとでの事後平均と事後分散を求めることである。別の方法は、\(Z\)マージンについての情報の一部ないしすべてを捨ててモデル化することである[よくわかんないけど、層の併合のこと？]。
まず、事後層を通じて\(Y\)の分布が同じだというナルモデルを考えよう。$$ (y_i | z_i = h, \mu, \sigma^2) \sim_{ind} G(\mu, \sigma^2)$$ $$ p(\mu, \log \sigma) = const$$ ここから以下が得られる: $$ E(\bar{Y} | \mathbf{Z}, \mathbf{Y}_s) = \bar{y}$$ $$ var(\bar{\mathbf{Y}} | \mathbf{Z}, \mathbf{Y}_s) = \upsilon = \delta(1-f)s^2/r$$ ただし\(s^2\)は標本分散で\(\delta = \frac{r-1}{r-3}\)。
今度は、\(Z\)についてのデータを無視すると $$ (y_i | \mu, \sigma^2) \sim_{ind} G(\mu, \sigma^2)$$ $$p(\mu, \log \sigma) = const$$ ここからも \( E(\bar{Y} | \mathbf{Y}_s) = \bar{y}, var(\bar{\mathbf{Y}} | \mathbf{Y}_s) = \upsilon \)が得られる。前者は後者を含意する(逆はいえない)。MCARであれば両者は同じことだが、MCARでない場合に話がちがってくる。
このように、基本的正規事後層別モデルのもとでは事後平均と分散は\(\bar{y}_{ps}, \upsilon_{ps}\)になるし、上の2つのモデルのもとでは\(\bar{y}, \upsilon\)になる。これは、ランダム化フレームワークで\(\bar{y}_{ps}\)の条件付き分散と\(\bar{y}\)の無条件分散が得られたのと同じである[あっ！ そういう話の展開だったのね！]。
このように、ベイジアンアプローチでは、精度については2つの自然な指標 \(\upsilon_{ps}, \upsilon\)が得られる。
2.5 事後層の併合アルゴリズム
事後層カウントを使う、無視する、の中間に、小さな事後層は併合する、というのがある。そのアルゴリズムについて考える。
まずは、MCARを仮定したアルゴリズムについて。基本的正規事後層別モデルの下で事後分散\(\upsilon_ps\)の期待値を小さくするのが目標である。
層\(i, j\)を併合したときの事後平均は、\( \bar{y}_{(i+j)} = \frac{r_i \bar{y}_i + r_j \bar{y}_j}{r_i + r_j}\) として$$ \bar{y}^{(ij)}_{ps} = \sum_{h \neq i,j} P_h \bar{y}_h + (P_i + P_j) \bar{y}_{(i+j)} $$ 事後分散はこうなる。\(\delta_{ij} = \frac{r_i + r_j – 1}{r_i + r_j – 3}, f_{ij} = \frac{r_i + r_j}{N_i + N_j}\), 併合した層の標本分散を\(s^2_{ij}\)として$$ \upsilon^{(ij)}_{ps} = \upsilon_{ps} – \Delta \upsilon_{ij}$$ $$ \Delta \upsilon_{ij} = \frac{P^2_i \delta_i(1-f_i)s^2_i}{r_i} + \frac{P^2_j \delta_j(1-f_j)s^2_j}{r_j} – \frac{(P_i + P_j)^2 \delta_{ij} (1-f_{ij} s^2_{ij}}{r_i + r_j} $$ [落ち着け、そんなに難しいことはいっていない]
期待値をとろう。修正項\(\delta_i, \delta_{ij}\)は無視する。$$ E[\Delta \upsilon_{ij}] = \frac{P^2_i (1-f_i) \sigma^2_i}{r_i} + \frac{P^2_j (1-f_j) \sigma^2_j}{r_j}- \frac{(P_i + P_j)^2 (1-f_{ij} \sigma^2_{ij}}{r_i + r_j} $$ \(Z\)と\(Y\)が独立なら、簡単にはこうなる。$$ E\left[ \frac{\Delta \upsilon_{ij}}{\sigma^2} \right] \approx \frac{r_i r_j}{r^2(r_i + r_j)} (w_i -w_j)^2 $$ つまり、事後層のウェイトが異なるのならば分散は下がる。いっぽう\(Z\)と\(Y\)が独立でないなら、かわりに\(\sigma^2_{ij}\)が大きくなる。
というわけで、お勧めの手続きはこうだ。
- 近隣が等質になるように事後層を並べる[??? 等質というのはたぶん、\(\mu_i, \sigma^2_i\)が似てそうだということだと思う]。事後層がもともと順序つきであればこのステップはいらない。
- 隣接層ペアのうち、\( E[\Delta \upsilon_{ij}] \)が最大のペアを選んで併合する。
- \( E[\Delta \upsilon_{ij}] \)がいずれも小さくなるか、層の数がreasonableになるまで続ける。
[… この後も解説が続くが、あんまし関心なくなってきたのでパス]
[事例。パス]
今度は、バイアスが存在するときの併合アルゴリズムについて。
[面白そうな話ではあるが、いまちょっと関心ないのでパス。事例もあるけどパス]
2.6 その他の併合戦略
事前分布をいじるという手もある。層の平均が交換可能だとみなせるなら、それらは共通の分布からのiidだとモデル化できる。つまりランダム効果モデルである。
たとえば、すべての\(h\)で\(\sigma_h = \sigma\)として $$ p(\mu_h | \sigma) \sim_{ind} G(\mu, \tau^2)$$ $$ p(\mu, \log \tau^2) = const$$ とすれば、\(\bar{y}_{ps}\)は\(\bar{y}\)に平滑化される。そうでなくて$$ p(\mu_h | \sigma_h) \sim_{ind} G(\mu, \tau^2)$$ $$ p(\mu, \log \tau^2, \log \sigma^2_h) = const$$ とすれば、\(\bar{y}_{ps}\)は\(s^{-2}_h\)で重みづけた平均へと平滑化される。ただし、交換可能性を仮定するというコストを支払っていることに注意。
ウェイトを直接に平滑化するという手もある。上限(たとえば3)を決めてトリミングするとか。
回答率をモデル化するという手もある。標本の回答率を\(\varphi_h\)とし、選択かつ回答率を\(\pi_h = \delta \varphi_h\)としよう。\(\pi_h\)がベータ分布に従うと仮定して…[パス]
こうしたアプローチは予測モデリングの形式ではないことに注意。モデル不確実性は無回答じゃなくて、\(Z\)のもとでの\(Y\)の分布にあるからだ。MCARの場合を考えるとわかりやすい。\(\pi_h\)の真のモデルは\(\pi_h = \pi\)だから、\(\bar{Y}\)の推定量は\(\bar{y}\)になっちゃうでしょ？ でも事後層がアウトカムの予測に役立つならばこれは効率的推定量でない。
3. 事後層別子が複数あるとき
[集中力が途切れてきたので細かーくメモする]
3.1 イントロダクション
事後層別子が複数ある場合、事後層が空だったりスパースだったりすることが多くなるので、基本的正規事後層別モデルを修正する必要性が高まる。また、可能な修正の幅も広くなる。
以下では事後層別子が2つ(\(Z_1, Z_2\))の場合について考える。\(Z_1 = h, Z_2 = k\)の母集団割合を\(P_{hk}\)とする。$$ \bar{Y} = \sum_h \sum_k P_{hk} \bar{Y}_{hk} $$ となる。基本的正規二元事後層別モデルは$$ (y_i | z_{1i} = h, z_{2i} = k, \mu_{hk}, \sigma^2_{hk}) \sim_{ind} G(\mu_{hk}, \sigma^2_{hk})$$ $$ p(\mu_{hk}, \log \sigma_{hk}) = const$$ となる。
3.2 基本的な正規二元事後層別モデルの下での推論
まず、\(\bar{Y}\)の事後平均について。モデルより、$$ E[\bar{Y} | \mathbf{Z}, \mathbf{Y}_s] = \bar{y}_{ps} = \sum_h \sum_k P_{hk} \bar{y}_{hk}$$ である。状況によっては、周辺分布 \( \{ P_{h+} \}, \{ P_{+k}\} \)だけがわかっているということもある。
- もし無回答についても\(Z_1, Z_2\)がわかっているのなら、つまり\(\mathbf{n} = \{ n_{hk} \}\)が観察されているのならば、無作為抽出のもとでは$$ \{ n_{hk} \} | n \sim MNOM[\{P_{hk}\}, n]$$ と考えられる。\(\{ P_{hk}\}\)の事前分布をジェフリーズとすると、\(\{n_{hk}\}\)のもとでの事後分布はパラメータ\(\{n_{hk} + 1/2\}\)のディリクレ分布である。さらに\( \{ P_{h\cdot} \}, \{ P_{\cdot k} \}\) が所与だとすると、事後分布はもうややこしくて書けないが、事後平均は近似的に、\( \{ n_{hk} \}\)を\( \{ P_{h\cdot} \}, \{ P_{\cdot k} \}\)にレイキングした結果になる: $$ \{ \hat{P}^{(1)}_{hk} \} = Rake[ \{ n_{hk} \}; \{P_{h\cdot}\}, \{P_{\cdot k}\} ] $$
- 無回答について\(Z_1, Z_2\)がわかんない、つまり\(\mathbf{n}\)が未知だとしよう。こんどは\(\{r_{hk}\}\)でどうにかしないといけない。こういうモデルを考える。$$ \{ r_{hk} \} | r \sim MNOM[\{P_{hk} \varphi_{hk} / \bar{\varphi} \}, r]$$ ただし \(\bar{\varphi} = \sum_h \sum_k P_{hk} \varphi_{hk} \)。このモデルは識別できないけど、\(\{r_{hk}\}\)を周辺にレイキングして$$ \{ \hat{P}^{(2)}_{hk} \} = Rake[ \{ r_{hk} \}; \{P_{h\cdot}\}, \{P_{\cdot k}\} ] $$ とすれば、これは$$ \phi_{hk} = \alpha_h \beta_k$$ という仮定の下で最尤である。Little & Wu (1991 JASA)をみよ。
まあとにかく、どうにかして\( \{ \hat{P}_{hk} \} \)を得れば、あとは$$ \hat{y} \sum_h \sum_k \hat{P}_{hk} \bar{y}_{hk} $$ とすればよい。
次に、\(\bar{Y}\)の事後分散について。
[…力尽きた。パス]
3.3 小さなセルの扱い
[ほんとに集中力がなくなってきたので、この節のみほぼ逐語訳でゴー]
本節では、事後層別子が2つのときのBNPM[基本的正規事後層別モデル] のもとでの推論を修正する方法について概略を示す。
ひとつのprincipledなベイジアン・アプローチは、3.1節末尾に示したBNPMモデルにおける平らな事前分布を変え、事後層間で借りてくるようにすることである。修正の仕方は状況による。ひとつのアプローチは、\(Z\)における\(Y\)の高次交互作用をランダム効果とみなすことである。たとえば、$$ \mu_{hk} = \mu + \alpha_h + \beta_k + \gamma_{hk}$$ とし、\( \{\mu, \alpha_h, \beta_k\} \)の事前分布は平らにするが、\(\gamma_{hk}\)の事前分布は\(N(0, \sigma^2_\gamma)\)とする、というもである。その結果は、主効果が固定で交互作用効果がランダムな混合ANOVAモデルになる。\(\bar{Y}_{hk}\)の予測は、有限母集団修正を省けば $$ \tilde{y}_{hk} = u_{hk} \bar{y}_{hk} + (1-u_{hk}) \hat{y}_{hk} $$ $$ \hat{y}_{hk} = \hat{\mu} + \hat{\alpha}_h + \hat{\beta}_{k} $$ $$ u_{hk} = \frac{n_{hk} \hat{\sigma}^2_\gamma}{ n_{hk} \hat{\sigma}^2_\gamma + \hat{\sigma}^2_{hk}} $$ という形になる。\(u_{hk} \)はシュリンケージ・ファクターである。\(\hat{\sigma}^2_\gamma, \hat{\sigma}^2_{hk}\)は最尤法なりモーメント法なりで求める。
BNPMモデルの事前分布の修正は、それぞれのアウトカムに対してチューニングする必要がある。だから、たくさんの変数についての大きな調査では現実的でないだろう。もっと単純な戦略として、事後層の部分的情報を無視するという方法がある。よく用いられるのは、標本ないし回答者カウントをマージンにレイキングしてしまい、\(Z_1, Z_2\)の同時分布の情報を無視してしまうというものである。回答者カウントより標本をレイキングするほうがよい(回答確率についてモデルがないから)。
レイキングを正当化できるのは、結果変数の平均が近似的に加法的構造を持っているときである。$$ \bar{Y}_{hk} = \mu + \alpha_h + \beta_k$$ ならば、$$ E[\bar{Y} | data) = \hat{\mu} + \sigma_h P_{h\cdot} \hat{\alpha}_h + \sigma_k P_{\cdot k} \hat{\beta}_k$$ であり、ここには\(Z_1, Z_2\)の周辺分布しかないから、レイキングしても情報は失われない。[…中略…]
交互作用がありそうなときはどうするか。その場合は併合戦略が良いだろう。\(Z_1\)は主要な層別子でアウトカムともっとも密接に関連し、\(Z_2\)は二次的な層別子だとしよう。この場合は、\(Z_1\)の各値ごとに\(Z_2\)のクラスの併合を行うとか、\(Z_2\)を速く動かすserpentine orderでセルを並べて併合するとか[\(Z_1\)が変わったら折り返すような並びをつくるってことだろう]。
無回答の場合のように、MCAR仮定が破られる場合、併合戦略は無回答バイアスを制約することに焦点を当てるべきである。無回答率を事後層別変数に回帰して、無回答率の推定値について等質であるような事後層を併合するという方法があるだろう。事後層別子の数が多いときに便利な戦略である。
4. 結論
本論文で触れられなかったこと: 他の量(母集団における傾きとか)についての推論。もっと複雑な抽出デザインの話。事後層別をかける先のマージンが怪しいとき(センサスと調査で収入の定義がちょっと違うとか)。標本カウントじゃなくて補足変数の合計にレイキングする話。
調査についてのモデリング的視点にご同意いただけない統計家のみなさんも、こうした場合のモデルについての研究から利益を得ることができるだろう。モデラーのみなさんにはこういういいたい。事後層別、さらにいえばサーヴェイ手法というのは、豊かで手付かずな領域であり、多くの解かれるべき問題が残っている、と。[ヒューヒュー]
————-
いやー、面白かった。いろんな話題が登場して飽きさせない。大変勉強になりました。
途中で疑問に思ったことをメモしておく。
非確率標本に基づいて母集団特性を推定したい時、母集団と標本における共変量の同時分布を使って、ある個体が標本に属するかどうかを予測する傾向スコアモデルを作り、傾向スコアの逆数をウェイトにして標本を集計分析することがあるでしょ？ バリバリにモデリングしているけど、関心ある変数についてのモデリングではない。よって、あれはモデルベース・アプローチ(Little先生の言い方では「予測モデリング」)ではない… ということであっているだろうか？
たぶんそうなんだろうな。でも、デザインベース・アプローチ(Little先生のいう「ランダム化」)と呼ぶのも違和感がある(標本抽出デザインの情報を使っていない)。デザイン・ベースとモデル・ベースというのは母集団推定アプローチの理念型に過ぎず、個別具体的な手法については必ずしもクリアに分類できない、ということだろうか。