Gelman, A., Little, T.C. (1997) Poststratification Into Many Categories Using Hierarchical Logistic Regression, Survey Methodology, 23(2), 127-135.
00年代の選挙予測でブイブイいわせ、いまでは調査データ分析の標準的ツールにまで成り上がってしまった、ご存じミスターP(MRP, マルチレベル回帰・層化)の元論文としてよく引用されるもの。掲載誌はカナダ統計局の地味ーな雑誌である。
MRPの解説はいまや巷に溢れているから、別にいまこれを読むことはないんだろうけど、調べ物のついでに読んだ次第である。温故知新！
1. イントロダクション
世論調査で事後層別のウェイトつけるのってふつうですよね。でもあんまし細かく層別しすぎると回答者数が小さくなっちゃいますよね。そこで反応を事後層別変数で条件付けたモデルを作ったりしますよね。レイキングってありますけど、あれはモデルでいえば高次交互作用をゼロにしていることになるわけです。事後層別カテゴリが階層構造になってたら階層モデルを組むという手もあります。Longford(1996)は小地域推定に階層線形モデルを使うのを提案してます。
本論文では二値変数の事後層別推定値を得るための階層ロジスティック回帰モデルを提案する。普通の事後層別と違って、詳細な母集団情報を使って細かーく層別きる。小地域推定にも便利です。
2. モデル
2.1 標本抽出と事後層別情報
\(R\)個のカテゴリカル変数があって水準数は\(J_1, \ldots, J_R\)だとする。セル数を\(J = \prod_{r=1}^R J_r\)とする。セル\(j\)の母集団サイズを\(N_j\)とする。
関心ある二値変数を\(y\)とし、そのセル\(j\)での母平均を\(\pi_j\)とする。$$ Y = \frac{\sum_j N_j \pi_j}{\sum_j N_j} $$ とする。母集団は十分に大きいとして有限母集団修正は無視する。
調査をやりました。セル\(j\)の個体数は\(n_j\)でした。無回答は無視可能でした。
2.2 事後層別の文脈での回帰モデリング
”yes”確率\(\pi_j\)について$$ logit(\pi_j) = X_j \beta $$ と考える。
ここで\(\beta\)に一様事前分布を与えたならば、\(X\)に応じて、これは古典的なウェイティングと密接に関連する(正確に言うと非線形変換が入っているところがちがうけど)。
- \(X\)が\(J \times J\)単位行列なら、セル\(j\)にウェイト\(N_j / n_j\)をつけているのと同じ、つまり単純な事後層別と同じである。
[なにをいうておるのかというと…式にlogit()が入っているので話がややこしいが、これを無視して\(\pi_j = X_j \beta\) だとしよう。\(X_j\)が単位行列なら\(\pi_j = \beta_j\)である。その推定量はyesと答えた人数を\(y_j\)として\(\hat{\beta}_j = y_j/n_j\)だ。で、\(N=\sum_j N_j\)として\(\hat{Y} = \frac{1}{N} \sum_j \hat{\beta}_j N_j\)を求めればよい。書き換えると \(\hat{Y} = \frac{1}{N} \sum_j y_j (N_j/n_j)\)だ。という話をしているのであろう] - \(X\)が\(J \times (\sum_{r=1}^R J_r)\)のインジケータ行列だったら、one-wayのマージンにレイキングしたのと同じである。
- \(X\)にいろいろ交互作用を入れた構造を作れば、そういう構造でプーリングした事後層別ウェイティングになる。
- \(X\)が単なる1のベクトルだったらそれは標本平均である。
2.3 部分プーリングのための階層回帰モデリング
セルを通じた部分的なプーリングのため、混合効果モデルを組んでみよう。上の式の\(\beta\)を$$ \beta = (\alpha, \gamma_1, \ldots, \gamma_L)^\top $$ とする。\(\alpha\)はプールしない係数の下位ベクトル、\(\gamma_l\)は残りは階層モデル $$ \gamma_{kl} \sim_{ind} N(0, \tau^2_l), \ \ k=1, \ldots, K_l$$ の\( \gamma_{kl}\)のサブベクトルである。[難しいことを云うね… 3節の事例にあてはめれば、要するに、\(\beta\)とは全国における係数と地域別切片を縦に積んだベクトルである。\(l=1,2,3,4\)というのが地域(Northeast, South, North Central, West)であり、\(k=1,…,K_1\)というのがNortheastの州(12州)である]
反応をセルに分ける\(n \times J\)行列を\(C\)とする。\(C_{ij} = 1\)は回答者\(i\)がセル\(j\)に属することを表す。\(Z = CX\)とする。するとこう書ける。$$ y_i \sim Bernoulli(p_i) $$ $$ logit(p_i) = Z \beta $$ $$ \beta \sim N(0, \Sigma_\beta) $$ ここで\(\Sigma^{-1}_\beta\)は対角行列で、\(a\)の箇所は0、\(\gamma_l\)の箇所は\(\tau^{-2}_l\)である。
2.4 モデルの下での推論
経験ベイズ作戦でよろしい。まず、データのもとでハイパーパラメータ\(\tau_l\)を推定する。次に、データと\(\tau_l\)の推定値の下で回帰係数\(\beta\)をベイズ推定する。で、\(\pi = logit^{-1}(X \beta)\)を求める。で、最後に\(N_j \pi_j\)を合計して母集団特性の推定値を出す。
ま、そうじゃなくてフル・ベイジアンでやってもいいけどな。この違いは技術的な話なのでここでは触れない。この論文のポイントはそこじゃない。
各セルの\(\beta\)の推定値はシュリンクする。その程度は\(n_j\)が小さいときに大きく、ロジスティック回帰モデルによる\(\bar{y}_j\)の予測が実際のそれと離れていた時に大きくなる。もちろん、\(\tau_l\)が小さいときにシュリンケージは大きくなる。\(\tau_l\)に予測力がなく0に近いときも大きくなる。
3. 適用例
3.1 サーベイ・データ
1988年US大統領選(Bush vs. Dukakis)の直前2週間にCBSがやった7個の全国世論調査(48州)を使う。Bush支持を\(y_i = 1\), Dukakis支持を\(y_i = 0\)とし、意見なし(15%)は捨てる。
CBSはウェイトをレイキングで作っていた。変数は{地域(5), 性別(2), エスニシティ(2), 年代(4), 学歴(4)}で、レイキングに当たっては、5つの主効果と、性xエスニシティ、年代x学歴を含めていた[つまり、周辺分布5つと同時分布2つをターゲットしてIPWを推定したわけね]。
さて、我々はこの5変数を固定効果としてモデルに入れ、さらに48州のインジケータも入れる。
3.2 事後層別のための人口データ
人口データとして、性(2)xエスニシティ(2)x年代(4)x州(48)をセルにした\(N_j\)を用意する[おっとお？ 学歴がないぞ]。Public Use Micro Survey(PUMS)データを使う。
3.3 結果
7つの調査に対して次の4つの手法を適用した。
- 古典的推定。CBSとほぼ同様で、地域, 性, エスニシティ, 年代, 学歴, 性xエスニシティ, 年代x学歴についてレイキングしたウェイトを使う。
- 回帰推定。上記7つの変数と州インジケータを使う。
- 回帰推定。上記7変数を使う。
- 階層モデルを使う。上記7変数と州の効果を入れる。2節の表現では、\(L=4\)で、\(k\)が州となる。[ちょっと待って、7変数の効果は固定なの？ まじ？ いまどきMRPっていったら個人特性の係数もランダム傾きにするよね]
モデルを推定し、それぞれの係数についての事後シミュレーション・ドローを結果をPUMSに基づいて重みづけし、各州のBush投票率の事後層別推定値を得た。
結果を実際の投票結果と比べる。4つの方法による推定は、国レベルではほぼ同じ。州レベルでは階層モデルの圧勝。
[ここからちょっと面白いけど、混乱してきたので逐語訳]
興味深いことに、階層モデルはデータを国レベル平均へと十分にシュリンクさせていないようにみえる。実際の選挙結果は、予測値が低い箇所では予測値より高く、予測値が高い箇所では予測値より低い。推定されたパラメータ\(\hat{\tau}_l\)が真の値よりおそらく高くなっていることをアンダーシュリンケージという。アンダーシュリンケージが生じるひとつの理由は、無視可能でない無回答のパターンが週によって異なるせいで、週別の比率の観察された変動が、実際の平均的世論の変動だけでなく無回答パターンの違いによっても引き起こされていた、というものである。Little & Gelman (1996 tech.rep.), Krieger & Pfeffermann (1992 SurveyMethodology)をみよ。アンダーシュリンケージは推定された値を最適なシュリンケージ水準と比べることで定量化できるが、その比較のためには真値が必要になる。
[調査別の分析など、細かそうな議論。めんどくさいのでパス]
4. ディスカッション
事後層別は選択確率の不均一性と無回答を補正するための標準的手法である。モデリングの観点からは、共変量集合への事後層別は共変量で条件づけた回帰モデルと密接に関連している。しかし多すぎる共変量へのレイキングは無理がある。本論文では階層モデルを適用した解決を示した。
関連する課題: 異なる組織による調査の結合。一般化線形モデルに拡張して連続的反応を扱う。
[…]
本論文で述べた手法はレイキング型の事後層別調整の改善で合って、無視可能でない無回答の修正を意図しているわけではない[あっそうなんだ]。しかし、この手法ではたくさんの変数への調整ができるから、無視可能性の想定がより合理的であるようなモデルが可能になる。面白いことに、この手法では古典的ウェイティングと違って、たくさんのセルを作ったほうがより良いモデルになる。
—–
わかりやすい論文であった。私のなかでGelman先生とLittle先生はどちらもわかりやすい論文を書くことで有名なのだが、その二大巨頭の共著だもんな。
この論文、MRPの始まりとしてよく引用されるのだが、実のところかなりマユツバだと思っていた。後になって「そういえば俺たちの書いたあの論文はいうなればMRPだったよな」って気が付いたんじゃないの？ という疑念があったのである。このたびはじめて目を通してきて、こりゃあ確かにMRPだと思った。疑ってスイマセンでした。
とはいえ、いまでいうMRPとはかなり違いがある。地域レベル共変量を入れていないし、個人レベル共変量の効果はおそらく固定で、地域レベル切片だけをランダムにしているようにみえる。ついでにいうと、フルベイズで推定するんじゃなくて経験ベイズでいいよといっていたり、付録でわざわざEMアルゴリズムによる推定方法を示している点も面白い。このあたりは時代の違いですね。WikipediaによればStanの初リリースは2012年だそうだ。
MRPの適用に対する個人的な疑問のひとつとして、センサス側で層別変数になってない個人レベル共変量をわざわざモデルに入れていることがあるけど、あれって意味あるの？ という点がある。いや、その係数に関心があるんなら別ですけど、母集団特性の推定という観点からは意味ないんじゃない？ という疑問である。この論文の例でいうと学歴がそうである。いま、高学歴な人ほどBush支持の確率が高いとする(知らんけど)。アーカンソー州の標本は母集団に比べて運悪く高学歴者に偏ってしまっているが我々はそれを知らないとしよう。モデルに学歴を含めなかった場合、アーカンソー州のランダム切片が高めに歪む。よって、事後ドローにおいてBush支持率は高めに歪む。いっぽう、モデルに学歴を含めた場合、アーカンソー州のランダム切片は歪まずに済むが、州レベル支持率を推定するための事後ドローにおいては、アーカンソー州についてのドローは学歴が高い人に偏るので、支持率の推定値はやはり高めに偏る。結局どちらも偏る。だったらモデルはシンプルにしておいた方がいいんじゃないだろうか？
というわけで、層別変数でない個人レベル共変量をモデルに入れるのはどうなんだろうかと疑っているのだが、実はこの論文でもそうなっていたことがわかった。すいませんでした。うーん、でも納得いかないなあ。