読了：Gao et al.(2019,2021) 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とかも使ってるっぽい.. . くそう…

　2025/04/05追記: 事情により、公刊された版を再読した。
Gao, Y., Kennedy, L., Simpson, D., Gelman, A. (2021) Improving multilevel regression and poststratification with structured priors. Bayesian Analysis, 16(3), 719-744.
　再読に伴い、このメモも一部追記・修正した。

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

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

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}\)だけど、上添字 にインデクスを書くのがすごく苦手なので書き換えた。ランダム切片はセルじゃなくて個々の層別変数の水準に与えられているわけで、つまり層別変数間の交互作用は考えないってことであってる？]
　これらは\( \theta_j \)を大域的にフィットする回帰モデル\(X_j \beta\)へと部分的にプーリングする効果を持つ。スパースなセルはこの正則化からほとんどを得ることになる。[つまり、票数の小さなセルの\( \theta_j \)の推定は他のセルからの情報に依存するってことね]
　デザイン行列の行を \( 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を使っているせいだろう]
　[→ 上のコメントは2020年に書いたものである。たった5年前なのに、ずいぶん正直に初々しいことを書いているなあ、と感心した。最近ではこういう新鮮な驚きがなくなってしまった。疲れているのかも…]

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

3. 提案アプローチとモチベーション
　ここからは、ガウシアンマルコフ確率場 (GMRF) の形をとる構造化事前分布を考える。[あっちゃー… 難しい話になってきちゃった…]

　ある共変量について、次の2つの場合を考える。

  • Case 1. もしあるカテゴリカル共変量についてなんらかの構造をモデル化したくないのなら、その変動切片は独立正規だとモデル化する。すると上記のモデルになる。
  • Case 2. ある共変量についてなんらかの構造をモデル化したい場合、そして、この構造を用いて空間的に平滑化するのが関心あるアウトカムに関してsensibleであるように思えるならば、適切なGRMFを事前分布とする。

　後者の場合、GMRFに情報的なハイパーパラメータを与えられるんなら与えてもよい。Simpson et al.(2017 Stat.Sci.)をみよ。
　構造化事前分布は有向の条件分布でも無向のでもよい。たとえば、離散時間インデクスと一緒に自己回帰過程とランダム・ウォーク過程をいれてもよいし、空間モデルでよくつかうCAR過程とかICAR過程とかをいれてもよい。
　一様じゃない情報借用を可能にするような複雑な事前構造をいれてもよい。たとえば、いちばん上の年代について推論するときにその下の年代から借りてくるとか。この場合、その順序変数に自己回帰的な事前分布を与えらればよい(回帰に年代を線形予測子とか二次予測子とかとして入れるような強い仮定を置くんじゃなくて)。

　構造化事前分布はMRPのマルチレベルの側面を改善するだけで、回帰構造はてつかずである。ふつうの回帰の代わりにスパース階層回帰とかベイジアン加法回帰木とかを使ってもよい。しかし、回帰のステップにおいてなんらかの正則化を入れて、調整因子と交互作用(それはすごくたくさんになりうる)を説明する能力を失わないようにしておくことである。回帰に構造化事前分布を入れるのは、機械学習的な正則化法に比べて解釈しやすい。

　GMRFはガウシアン過程(GPs)と密接な関係を持っている。MRPにおける構造化事前分布とは、構造化された共変量に対するGP事前分布の離散的な近似である。GPsがユニバーサルな関数近似であるのと同じように、MRPにおける構造化事前分布、データ生成過程を説明すると同時に解釈可能な正則化を提供する、ひとつの柔軟なモデリング方略だと考えることができる。

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

4.1 有向の構造化事前分布の例: 群レベル誤差を部分プーリングするモデル
　事後層化変数は州(51水準), 年代(12水準), 収入(4水準)。ほかに予測子はない。[説明がないんだけど、式をみると51州は5つのリージョンに分かれるらしい。なお、原文では年代はわざわざAge Cat.と表記されているが、見ずらいのでAgeと書き換える]
　$$ Pr(y_i = 1) = logit^{-1} (\beta^0 + \alpha_{j[i]}^{State} + \alpha_{j[i]}^{Age} + \alpha_{j[i]}^{Income}) $$ $$ \beta^0 \sim N(0, 1)$$

  • 州の効果は [原文では　\( X_{VS,j}, X_{Relig, j} \)なんだけど書き換えた] $$ \alpha_j^{State} \sim^{ind} N(\alpha_{m[j]}^{Region} + \beta^{Relig} X_{j}^{Relig} + \beta^{VS} X_{j}^{VS}, (\sigma^{State})^2)$$ $$ \beta^{Relig} \sim N(0, 1)$$ $$ \beta^{VS} \sim N(0, 1)$$ $$ \sigma^{State} \sim N_+(0,1)$$ $$ \alpha_m^{Region} \sim^{ind} N(0, (\sigma^{Region})^2) $$ $$ \sigma^{Region} \sim N_+(0,1)$$ ただし、\( X_{j}^{VS} \in [0,1] \)は2004年の民主党の得票シェア, \(X_{j}^{Relig} \in [0,1]\)は保守派宗教(モルモン教と福音派)の割合。
  • 年代はあとまわしにして…
  • 収入の効果は$$ \alpha_j^{Income} \sim N(0, (\sigma^{Income})^2)$$ $$ \sigma^{Income} \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の分布は書かずに済むね]

上記3種類の事前分布はどうちがうか。年代のランダム効果における、隣接カテゴリが共有する情報の量が異なる。\(\alpha^{Age}_j\)と\(\alpha^{Age}_{j-1}\)について考えよう。ベースライン方式ではなにも共有されず、自己回帰方式の場合は部分的な情報が共有され、ランダムウォーク方式の場合はすべてが共有される。情報の共有というのは、一方の事後分布が他方にシュリンクする程度と言い換えてもいい。自己回帰方式とランダムウォーク方式の場合は、\(\alpha^{Age}_j\)の事後分布は\(\alpha^{Age}_{j-1}\)の事後分布へとシュリンクする。その強さが自己回帰係数\(\rho\)で決まる。
　1次自己回帰仮定は、ある年代のある人の意見は、同じデモグラで年代だけちょっと若い人の意見に似ている、という事前の想定を組み込んでいるわけである。

　シミュレーション。[パス。高年齢層の回答率と、目的変数に対する年齢の効果の関数を、あれこれ変えては推定を試したという話だと思う。Stanでやった由。ランダムウォーク方式がよかった、とかなんとか]

4.2 無向の構造化事前分布の例: 空間的MRP　[この節は2019年のプレプリントにはなかった]
　マサチューセッツ州の52地域(PUMA)について推定する場面を考える。データはカウントで、\(i = 1, \ldots, n\)について可能な生起数が\(T_i\)、データが\(t_i\)だとする[\(i\)さんがある課題を\(T_i\)回やって\(t_i\)回成功した、というような話であろう]。$$ t_i \sim \mathrm{Binomial}(\mu_, T_i) $$ $$ \mu_i = \mathrm{logit}^{-1} \left( \beta^0 + \alpha^{PUMA}_{j[i]} + \alpha^{Education}_{j[i]} + \alpha^{Ethnicity}_{j[i]} \right) $$ $$ \beta^0 \sim N(0,1)$$

  • PUMAの効果は後回し。
  • 教育の効果は $$ \alpha^{Education}_j \sim N(0, (\sigma^{Eduction})^2 ) $$ $$ \sigma^{Eduction} \sim \mathrm{PC-Prior}(1,0.1) $$ PC(複雑性罰則)事前分布についてはSimpson et al.(2017 Stat.Sci.)をみよ。
  • エスニシティの効果は $$ \alpha^{Ethnicity}_j \sim N(0, (\sigma^{Ethnicity})^2 ) $$ $$ \sigma^{Ethnicity} \sim \mathrm{PC-Prior}(1,0.1) $$

さて、問題はPUMAの効果である。

  • Besag-York-Mllie (BYM2) specificationモデル。PUMA \(j\)の1次近隣の集合を\(A_j\), その濃度[個数のことね]を\(d_j\)として、$$ \alpha^{PUMA}_j = \frac{1}{\sqrt{\tau}} \left( \sqrt{1-\rho} \theta^*_j + \sqrt{\frac{\rho}{s}} \phi^*_j \right)$$ $$ \theta^*_j \sim N(0, I_{52}) $$ $$ \phi^*_j \sim N \left( \frac{\sum_{k \sim A_j} \phi^*_k}{d_j}, \frac{1}{d_j} \right), \ \ \sum_{j=1}^{52} \phi^*_j = 0$$ $$ \tau \sim \mathrm{PC-Prior}(1, 0.1) $$ なにをやっているのかというと、まず\(\theta^*_j\)は等方な多変量正規事前分布である。\(\phi^*_j\)はICAR事前分布になっている [そうなの？ まあここは信じよう]。\(s\)はスケーリング・ファクターで、すべてのPUMAにおいて、\(\sqrt{\frac{\rho}{\tau s}} \phi^*_j\)の分散がおよそ1になるようにしている。
  • IID specificationモデル。\(\alpha^{PUMA}\)は等方な多変量正規分布で精度は\(\tau \sim \mathrm{PC-Prior}(1, 0.1)\)。

　シミュレーション。[パス。R-INLAでやった由]

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

6. 結論
　[…中略…]
　MRPの事後層化の正確性は、事後層化行列が目標母集団の真の表現になっておるかどうかに依存する。カテゴリを細かくしたほうが事後の分散とバイアスが小さくなるけどセルサイズが小さくなるというトレードオフがある。もうひとつの限界は共変量選択で、モデラーの知識とデータとに依存する。
　この論文で扱えなかった問題：異なる構造化事前分布のうちどれを選択するか。交互作用項の含め方。共変量の構造がはっきりしないときの構造化事前分布の入れ方。
　云々。
—————-
　… 年代の効果に構造をいれるというくだりが勉強になりましたです。
　それにしても、方法論にキャッチフレーズをつけることって大事なんだな、と思う。去年の冬の英国総選挙の際(ジョンソンさんが馬鹿勝ちした奴)、YouGovがMRPによる予測を発表したときに、MRPという最新の魔術的予測手法によれば…というような書き方で読み手をあおるような記事を載せていたメディアがあった。蓋をあけてみれば、セル別母比率の直接推定が難しいから階層回帰で一気に推定し、母集団サイズで重みづけて合計しようという話であり、ちょっと気の利いた分析者なら思いつくアプローチであろう。

　[上記も2020年に書いたコメント。それから5年経ち、MRPはすっかり標準的手法になってしまい、もはやこういう素朴な感想は書きにくい。でも、これって個々の分野の実践家がそれぞれの現場で思いつくような手法なんじゃないかな、という思いはいまだに残っている]