読了:Gao et al.(2019) MRP(マルチレベル回帰・層化)に構造化事前分布をいれる

Gao, Y., Kennedy, L., Simpson, D., Gelman, A. (2019) Improving multilevel regression and poststratification with structured priors. arXiv:1908.06716v2. 30 Sep 2019.
 しばらく前に読んだ奴。たしか勉強のつもりで読んだのだと思う。
 最近の選挙予測でブイブイいわせているらしき、Mr.P こと Multilevel Regression and Poststratification (日本語ではなんていうんだろう? マルチレベル回帰・層化?) に、構造を持つ事前分布をいれるという論文。Mr.Pの生みの親 Andrew Gelman さんも著者に入っている。たぶん未公刊。

 いまみたらarXivに改訂版があがっていた。なんか内容が大幅改善されてるっぽい… INLAとかも使ってるっぽい.. . くそう…

1. イントロダクション
 MRPは非代表標本を大母集団へと調整する方法として人気上昇中である。小地域推定や恣意的抽出のような、伝統的なデザインベースのアプローチが苦戦していた分野でもうまくいく。
 伝統的なデザインベースの事後層化ウェイトとどうちがうかというと、partial poolingを用いるところである。単純な事後層化は空のセルをうまく扱えないので、交互作用を無視して周辺分布だけ事後層化したり、セルを併合したりしてお茶を濁すのがふつうであった。MRPではマルチレベルモデルを使って群推定値を自動的にregularizeする。そのやり方はいくつかあるけれど、ふつうは群レベルの誤差の独立性を仮定する。たとえば選挙予測だったら、州レベルの切片を想定し、過去の投票率のような州レベル予測子に回帰させ、かつ州レベル誤差を独立とみなすわけである。
 しかし場合によっては、回帰の予測子では捉えきれない基盤構造を含めるのがよろしいかもしれない。たとえば、予測子が順序カテゴリカル変数のときに自己相関構造をいれるとか。この場合における順序変数の予測子のことを構造化事前分布と呼ぶ。[←ここだけ読んでもなんのことだかわかりにくいが、あとで明確になる]

 標本抽出後の調整をなぜやるかというと、標本がバイアスを持っている可能性があり、それと目標母集団との差を修正したいからである。
 モデルベースの調査推定量においては、無回答を調整するために事後層化ウェイティングを使うことが多い。正確性は増すけれど、事後層化推定値の品質は層の母集団サイズについての既知情報の品質に依存するし、それぞれの事後層化セルにおいては代表標本だという仮定にも依存する。
 事後層化の近似としてレイキングを用いることもあるけど、要因の数が多くなるとウェイトの分散が大きくなって推定値が不安定になる。
 こうした非確率標本の事後調整については、Elliott, Vallian, et al. (2017 Stat.Sci.)をみよ。[←ノーチェックだったよ…ありがとう著者のみなさん]

2. MRP概観
 そもそもMRPとはこういう手法である。(1)調査データに階層回帰モデルを当てはめる。(2)個々の事後層化セルの母集団サイズをつかって重みつき調査推定値をつくる。

 母集団を \(K \) 個のカテゴリカル変数で分割できるとしよう。\( k \) 番目の変数のカテゴリ数を \( J_k \) とする。母集団は \( J = \prod_{k=1}^K J_k \)個のセルで表現できる。それぞれのセル \( j \) の母集団サイズ \( N_j \) は既知とする。なお、もともと変数が連続的だったときいくつに分割するのが良いかというのは難しい問題なんだけど、ここでは扱わない。
 個人 \( i \) の反応変数を \( y_i \in \{0,1\} \)とする。

 ステップ1, マルチレベル回帰。
 セル \( j \) の母平均を \( \theta_j \) とする。個々のカテゴリカル共変量 \( k \)について、変動する切片 \( \{ \alpha_{k,j*} \}_{j*}^{J_k} \) を考える。[原文では\(\{ \alpha_{j*}^k \}_{j*}^{J_k}\)だけど、上添字 にインデクスを書くのがすごく苦手なので書き換えた。ランダム切片はセルじゃなくて個々の層別変数の水準に与えられているわけで、つまり層別変数間の交互作用は考えないってことであってる?]
 デザイン行列の行を \( X_j \)とする。個々の行は \( \theta_j \) に対応する [つまりデザイン行列の行はセルなわけだ。言い換えると、層別変数以外の共変量もみんなセルレベルの共変量であって、個人レベル共変量とかはないわけだ]。で、$$ Pr(y_i = 1) = logit^{-1} \left(X_i \beta + \sum_k^K \alpha_{k,j[i]} \right) $$ $$ \alpha_{k,j} | \sigma^k \sim^{ind} N(0, (\sigma^k)^2) $$ $$ \sigma^k ~ \sim N_{+}(0, 1)$$ $$\beta \sim N(0, 1)$$ [説明がないけど、ランダム切片の下添字 \( j[i] \) は「その層別変数について個人\(i\)が属している水準」という意味であろう。
 MRPの定式化をきちんと読んだのははじめてなんだけど、この定式化はちょっと意外だった。セルレベルの撹乱項をいれるのかと思っていた。この論文は\(j\)がセルを指したり水準を指したりしていて困るので、層別変数 \(k\) 上での \( i \)さんの所属水準を \( l[i] \in \{1, \ldots, L_k\}\), 所属セルを\(c[i] \in \{1, \ldots, \prod_{k=1}^K L_k \}\) と書き換えると、$$logit(Pr(y_i=1)) = X_i \beta + \sum_{k=1}^K \alpha_{k,l[i]} + \epsilon_{c[i]}$$ として、\( \{ \epsilon_c \} \)になんらか狭めの事前分布を与えるようなモデルを想像していたよ…
 話はちがうが、\(\beta, \sigma^k \) にしれっと弱情報事前分布を与えている。わたし、こういうところでモヤモヤしちゃうんですよね… デフォルトは無情報じゃないといかんでしょう?と思ってしまう。これは別になにかのはっきりした信念があるからではなくて、おそらく、ふだんMplusを使っているせいだろう]

 ステップ2, 事後層化。
 母集団をセルに事後層別し、セルサイズじゃなくて母集団サイズで重みづけた和をとる。事後層別行列で定義されるなんらかの下位母集団を \( S \) として、$$ \theta_S = \frac{ \sum_{j \in S} N_j \theta_j }{\sum_{j \in S} N_j } $$

3. 提案アプローチとモチベーション
 ここからは、ガウシアンマルコフ確率場 (GRMF) の形をとる構造化事前分布を考える。[←あっちゃー… 難しい話になってきちゃった…]
 情報ハイパーパラメータを与えられるんなら与えてもよい。また、構造化事前分布は有向の条件分布でも無向のでもよい。たとえば、離散時間インデクスと一緒に自己回帰過程とランダム・ウォーク過程をいれてもよいし、空間モデルでよくつかうCAR過程とかICAR過程とかをいれてもよい。
 一様じゃない情報借用を可能にするような複雑な事前構造をいれてもよい。たとえば、いちばん上の年代について推論するときにその下の年代から借りてくるとか。この場合、順序変数に自己回帰する事前分布を与えらればよい(回帰に年代を線形予測子とか二次予測子とかとして入れるような強い仮定を置くんじゃなくて)。
 構造化事前分布はMRPのマルチレベルの側面を改善するだけで、回帰構造はてつかずである。ふつうの回帰の代わりにスパース階層回帰とかベイジアン加法回帰木とかを使ってもよい。

4. シミュレーション
 [原文では、以下の式で \( X|\mu, \sigma \sim N(\mu, \sigma^2) \) というスタイルで表記しているんだけど、わたくし、この書き方がすっごく苦手で… 先生方には悪いけど、\( X \sim N(\mu, \sigma^2) \) と略記させてもらいます]

 次の場合について考える。事後層化変数は州(51水準), 年代(12水準), 収入(4水準)。ほかに予測子はない。[説明がないんだけど、式をみると51州は5つのリージョンに分かれるらしい]
 リンク関数は$$ Pr(y_i = 1) = logit^{-1} (\beta^0 + \alpha_{j[i]}^{State} + \alpha_{j[i]}^{Age} + \alpha_{j[i]}^{Income}) $$ 州の効果は $$ \alpha_j^{State} \sim^{ind} N(\alpha_{m[j]}^{Region} + \beta^{Relig} X_{j}^{Relig} + \beta^{VS} X_{j}^{VS}, (\sigma^{State})^2)$$ ただし、\( X_{j}^{VS} \in [0,1] \)は2004年の民主党の得票シェア, \(X_{j}^{Relig} \in [0,1]\)は保守派宗教(モルモン教と福音派)の割合である。地域の効果は $$ \alpha_m^{Region} \sim^{ind} N(0, (\sigma^{Region})^2) $$ とする。[原文では \( X_{VS,j}, X_{Relig, j} \)なんだけど書き換えた]
 話を戻す。年代はあとまわしにして、収入の効果は$$ \alpha_j^{Income} \sim N(0, (\sigma^{Income})^2)$$ ここまでの3つの分散の事前分布は $$ \sigma^{Income}, \sigma^{State}, \sigma^{Region} \sim N_+(0,1) $$ 係数の事前分布は $$ \beta^{VS}, \beta^{Relig}, \beta^0 \sim N(0, 1)$$

 さて、いよいよ年代の効果について。次の3通りを考える。[ここからが本番だ。深呼吸…]

  1. ベースライン方式。年代ごとの切片が独立に正規分布に従うと考える。$$\alpha_{j}^{Age} \sim^{ind} N(0, (\sigma^{Age})^2)$$ $$ \sigma^{Age} \sim N_+(0,1)$$
  2. 自己回帰方式。年代の順序構造を一次の自己回帰と捉える。下の4本目の式では\( \rho \) を \( (-1, 1) \) に制約している。定常でないと困るからね。$$ \alpha_1^{Age} \sim N \left( 0, \frac{1}{1-\rho^2} (\sigma^{Age})^2 \right)$$ $$\alpha_j^{Age} \sim N(\rho \alpha_{j-1}^{Age}, (\sigma^{Age})^2) \ {\rm for} \ j=2,\ldots,12$$ $$ \sigma^{Age} \sim N_+(0,1)$$ $$ (\rho+1)/2 \sim Beta(0,5, 0,5) $$
  3. ランダムウォーク方式。やりたいことは \(\rho = 1 \)に固定することなんだけど、上の定式化のままだとゼロ割になっちゃうので定式化を変えている。$$ \alpha_j^{Age} \sim N(\alpha_{j-1}^{Age}, (\sigma^{Age})^2) \ {\rm for} \ j=2,\ldots,12$$ $$ \sum_{j=1}^J \alpha_j^{Age} = 0 $$ $$ \sigma^{Age} \sim N_+(0,1) $$ [おおお… これはすごく勉強になった。そうか、年代の効果はsimplexだと定義すれば、カテゴリ1の分布は書かずに済むね]

 シミュレーションのためのデータは…[力尽きたので以下メモ省略。要するに、高年齢層の回答率と、選好に対する年齢の効果の関数を、あれこれ変えては推定を試したという話だと思う。Stanでやった由]
 結果、ランダムウォーク方式がよかった。

5. 米の調査データの分析
[National Annenberg Electoin Survey 2008の実データを分析している模様。めんどくさいからパス]

6. 結論
 […中略…]
 MRPの事後層化の正確性は、事後層化行列が目標母集団の真の表現になっておるかどうかに依存する。カテゴリを細かくしたほうが事後の分散とバイアスが小さくなるけどセルサイズが小さくなるというトレードオフがある。もうひとつの限界は共変量選択で、モデラーの知識とデータとに依存する。
 この論文で扱えなかった問題:異なる構造化事前分布のうちどれを選択するか。交互作用項の含め方。共変量の構造がはっきりしないときの構造化事前分布の入れ方。
 云々。

 … 年代の効果に構造をいれるというくだりが勉強になりましたです。
 それにしても、方法論にキャッチフレーズをつけることって大事なんだな、と思う。去年の冬の英国総選挙の際(ジョンソンさんが馬鹿勝ちした奴)、YouGovがMRPによる予測を発表したときに、MRPという最新の魔術的予測手法によれば…というような書き方で読み手をあおるような記事を載せていたメディアがあった。蓋をあけてみれば、セル別母比率の直接推定が難しいから階層回帰で一気に推定し、母集団サイズで重みづけて合計しようという話であり、ちょっと気の利いた分析者なら思いつくアプローチであろう。