ベイズデータ解析７

ベイズ階層モデル

Bayesian
Statistics
Author

司馬博文

Published

12/12/2024

Modified

12/13/2024

概要

データが自然な階層構造を持つ場合，これを取り入れた自然な事前分布を，一つ上の階層に回帰モデルを付け加えることで構成できる． このようなモデルをベイズ階層モデルという． 本稿ではベイズ階層モデルの縮小効果を概観する． 事前知識を構造に関する知識としてモデルに取り入れることでデータによりフィットする尤度構造を獲得することは，データ解析の一つの目標として，（線型）回帰モデルの自然な拡張と理解できる．

Keywords

変動係数モデル, Mixed-effects model, NER model, 分散成分モデル, Fay-Herriot model, Two-way (crossed) classification model, 項目応答モデル, Rasch model

関連記事

1 はじめに

1.1 階層モデル

データが自然な 階層構造 を持つとする．例えば標本がクラスター抽出された場合，パネルデータや経時的繰り返し観測による標本，共変量統制とマッチングなど，自然な階層構造を持つデータは多い．

このようなデータに対して，個々のサブグループに対して全く同じ回帰モデルを繰り返し適用し，結果のプロットを小窓に分割して並べることで，グループごとの効果の違いを比較することができる．

この方法はナイーブながらも有力で，(10.9 節 Gelman et al., 2020, p. 148) では “secret weapon” と呼んでいる．

この別々の回帰モデルはベイズ的に統合して推論することができる．このグループごとに変動する係数を許したモデルが（ベイズ）階層モデル である．

1.2 いつ使うか？

線型回帰の枠組みで交差項を入れることで，グループごとに緩やかに異なる回帰係数を許すことができる．

さらには一般化線型モデルを用いることで，連続・離散データの簡単な非線型関係まで扱うことができる．

それでも階層モデルに進む動機として次のような点が挙げられる：

階層モデルの美点1
  1. 説明変数が多すぎる場合は馬蹄事前分布などの構造を持った事前分布が使用可能である．しかしデータが階層構造を持つ場合は，その構造に関する事前知識を用いて，自然な事前分布を組織的に導入することができる．これが階層モデリングに他ならないとも見れる．その結果大量の説明変数がある場合でも，その関連度をモデルに知らせて適切に推定する方法にもなる2．
  2. 単一のグループではデータ数が小さい場合に，グループレベルの回帰係数 \(\alpha_{j[i]}\) を安定して推定する手法とも理解できる．他のグループの回帰モデルと組み合わせた大きな階層モデルをベイズすることで，確度が低い（が関連する）モデルも捨てずに推定に利用するため (partial pooling)，より好ましい統計的推定を実行できる．これは 縮小推定 のキーワードの下でも追及されており，階層モデルは自然な縮小推定が実行される典型的な設定である．3
  3. グループ毎の回帰係数の違いに興味がある場合，グループレベルの変動の分散 \(\sigma_j^2\) を推定する最も自然なモデルである．実際，一階層の回帰分析は \(\sigma_j=0\) の場合，グループ毎に別々の回帰分析を実行する場合は \(\sigma_j=\infty\) という（場合によっては非現実的な）仮定を置いた縮退した階層モデルと見れる．

全体の係数のみに興味がある場合，グループレベルの変動やデータの階層構造は局外構造である．このような状況では，誤差が相関を持つ場合にも頑健な点推定手法として，一般化推定方程式法なども用いることができる．4

2 線型変動係数モデル

2.1 変量効果モデル

階層モデル \[ y_i=\alpha_{j[i]}+\beta x_i+\epsilon_i,\qquad\epsilon_i\sim\mathrm{N}(0,\sigma^2_y), \tag{1}\] \[ \alpha_j=\mu_\alpha+\eta_j,\qquad\eta_j\sim\mathrm{N}(0,\sigma^2_\alpha), \tag{2}\] は２つの別の階層にある回帰モデルが，１つに統合された形をしている．5

(1), (2) は 変量効果モデル，または NER (Nest Error Regression) model (Battese et al., 1988) と呼ばれる．6

2.2 中庸としての階層モデル

全く等価だが \(D_i\) を \(i\) 番目の成分のみが \(1\) のベクトル，\(\alpha=(\alpha_1,\cdots,\alpha_J)^\top\) として \[ y_i=D_i\alpha+\beta x_i+\epsilon_i \] と表すことができる．

こうみると \(J\) 個のクラスそれぞれに関する指示変数 \(D_1,\cdots,D_J\) で定数項を置き換えた単一の回帰モデルと見ることもできる．仮に複数の変動係数が存在しても共線型性が問題にならないのは，\(\alpha\) の事前分布に構造が入っているためである．

階層モデル (1), (2) において \(\sigma^2_\alpha\to0\) の極限を考えると，グループ毎に変動がないという単一の線型回帰モデルに帰着する．

他方で実は \(\sigma^2_\alpha\to\infty\) の極限を考えると，古典的な点推定の文脈では，\(J\) グループのそれぞれに別々の線型回帰を実行するというモデルに帰着する．

これは「グループ毎に全く別々で互いに関係がない」というモデルであり，真実は２つの中庸にあると思われる．

グループ毎の変動と，その間の緩い関係性の双方を許し，１つのモデルに取り込んだものが 階層モデル である．グループレベルの変動 \(\sigma^2_\alpha\) を推定するための自然なモデルでもある．

2.3 BLUE

線型最良不偏推定量 (BLUE) は (Henderson, 1950) によって提案された，\(\beta,\alpha_{j[i]}\) の推定量である．

\(\beta\) の推定量は計量経済学の文脈で変量効果推定量と呼ばれる GLS 推定量に他ならない．7

\(\alpha_{j[i]}\) の BLUE は，モデル (1), (2) においては \[ \widehat{\alpha}_j=\frac{\frac{n_j}{\alpha_y^2}}{\frac{n_j}{\alpha_y^2}+\frac{1}{\sigma^2_\alpha}}(\overline{y}_j-\beta\overline{x}_j)+\frac{\frac{1}{\sigma^2_\alpha}}{\frac{n_j}{\alpha_y^2}+\frac{1}{\sigma^2_\alpha}}\mu_\alpha \] と表せる．8

2.4 分散成分モデル

2.4.1 一般的な定義

(Henderson, 1950) に始まる分散成分モデルは，一般的には \[ y_{ij}=\beta x_i+\alpha_{j_1[i]}+\cdots+\alpha_{j_k[i]}+\epsilon_{ij} \] \[ \alpha_{j_l}\overset{\text{i.i.d.}}{\sim}\mathrm{N}(0,\tau_l^2), \] と表されるモデルと理解される (Sugasawa and Kubokawa, 2023, p. 73)．\(k=1\) の場合は (Fay and Herriot, 1979) モデルともいう．

変量効果モデルにおいて誤差には完全な階層関係があったが，分散成分モデルはその non-nested な場合への拡張とみなせる．

2.4.2 分散成分モデルの例

この場合データ \(y_{ij}\) の変動を成分ごとに分解しようという志向が，分散分析の発展と見れる．例えば \[ y_{ijk}=\mu+\mu_{1i}+\mu_{2j}+\epsilon_{ijk} \] という形の場合は特に 二元分類モデル，\(\mu_{3ij}\) も加えた場合は二元交差モデルと呼ばれる．

さらにグループ階層のモデルに関して，分散成分モデルではグループレベル変数 \(\alpha_{j_l}\) は互いに独立であるとしているが，空間モデルでは互いに相関を持つと仮定することもある．

さらにこの階層は潜在変数の階層とみるとコピュラモデルや理想点モデルに近くなる．項目反応モデルは，個人毎の変動と項目毎の変動の，non-nested な２つの変動を持ったロジスティック階層モデル 3 になる．

2.4.3 分散成分の点推定法

変量効果モデルや分散成分モデルにも適用可能な，一般の階層モデルに適用可能な分散成分の推定法として，一般化推定方程式法がある．

その特別な荷重を取った場合に REML (Restricted Maximum Likelihood) 推定量がある (2.3 節 Sugasawa and Kubokawa, 2023, p. 13)．

こうして得た分散成分，特に分散比 \(\rho:=\sigma^2_v/\sigma^2_e\) の推定量を得たあと，これを BLUE \(\widehat{\beta}\) に代入して得る２段階推定量を 経験 BLUE という (久保川達也, 2006, p. 143)．

その縮小効果については次稿も参照：

2.5 ベイズ推定

ベイズ階層モデルの基礎は (Lindley and Smith, 1972) が敷いたと言われる．

3 階層ロジスティックモデル

3.1 項目応答モデル

3.1.1 １母数ロジットモデル

\[ g(x):=\operatorname{logit}(x)=\log\frac{x}{1-x} \] \[ g(\mu_{i})=\alpha_{j[i]}-\beta_{k[i]},\qquad\mu_{i}=\operatorname{P}[Y_{ik}=1], \] というモデルを １母数応答モデル または (Rasch, 1960) モデルという．

このモデルは \((\alpha_j,\beta_k)\) 上の平行移動に関して識別不可能であるため，点推定の文脈では追加の制約条件が導入される．しかしベイズの文脈では不適な事前分布を避けることで自然に回避される．

階層ベイズモデルでは，自然に関心が能力パラメータ \(\alpha_j\) と難易度パラメータ \(\beta_k\) の分散に向けられる： \[ \alpha_j=\mu_\alpha+\gamma_\alpha X^\alpha_j+\epsilon_j,\qquad\epsilon_j\overset{\text{i.i.d.}}{\sim}\mathrm{N}(0,\sigma^2_\alpha), \] \[ \beta_k=\mu_\beta+X^\beta_k+\epsilon_k,\qquad\epsilon_k\overset{\text{i.i.d.}}{\sim}\mathrm{N}(0,\sigma^2_\beta). \]

3.1.2 ２母数ロジットモデル

さらに項目毎の 識別力母数 (discrimination parameter) \(\gamma_k\) を導入したモデル \[ g(\mu_i)=\gamma_{k[i]}\biggr(\alpha_{j[i]}-\beta_{k[i]}\biggl), \] を ２母数ロジットモデル という．

3.1.3 多次元の潜在空間

読解力と数学力の別々の能力を要求するテストなどのデータに対して，多次元の潜在空間を持つモデルを考えたい．

両方が要求される AND の論理の場合は \[ \mu_i=g^{-1}\biggr(\gamma^{(1)}_{k[i]}(\alpha_{j[i]}-\beta^{(1)}_{k[i]})\biggl)g^{-1}\biggr(\gamma^{(2)}_{k[i]}(\alpha_{j[i]}-\beta^{(2)}_{k[i]})\biggl) \] というモデルが考えられるかもしれない．OR の論理の各因数を \(1-\mu_i, 1-g^{-1}(-)\) で置き換えれば良い．

一方で潜在空間からの関数を設定しても良いだろう．最も直感的には能力値の和で \[ \mu_i=g^{-1}\biggr(\gamma^{(1)}_{k[i]}(\alpha_{j[i]}-\beta^{(1)}_{k[i]})+\gamma^{(2)}_{k[i]}(\alpha_{j[i]}-\beta^{(2)}_{k[i]})\biggl) \] と表すものである．

3.1.4 応答曲線

リンク関数 \(g\) の選択は，潜在変数 \(\alpha_j,\beta_k,\gamma_k\) の変化が応答確率の変化にどう関係するかを規定する．

\(g\) のプロットは ICC (Item Characteristic Curve) または trace line と呼ばれる (Fox, 2010, p. 6)．

これは \(g(0.5)\) を中心に対称になっているが，この対称性が不適切な場合も多い．

(Bafumi et al., 2005) などは \[ \mu_i=\pi_1+(1-\pi_0-\pi_1)g^{-1}\biggr(\gamma_k(\alpha_j-\beta_k)\biggl) \] として，最低限の \(y_i=0,1\) への反応確率 \(\pi_0,\pi_1\) をあらかじめ設定し，\([\pi_0,\pi_1]\subset[0,1]\) の範囲だけに応答曲線をフィッティングすることを考えた．

\(\pi_0,\pi_1\) も推定すべきパラメータとする．

3.2 選択モデル

3.2.1 はじめに

ロジスティックモデルは１次元の潜在空間を持つ 選択モデルと見ることができる．

これをさらに精緻化し，価値関数や効用関数などを通じてアクターの意思決定と行動をより詳細にモデリングすることも考え得る．

このようなモデルはミクロ計量経済学で歴史的に考えられたほか，政治科学における 理想点モデルはこのような選択モデルを通じて生まれた．

3.2.2 ロジットモデルになる場合

個人 \(i\in[N]\) ごとのパラメータ \(a_i,b_i,c_i\) が存在して， \[ \operatorname{P}[Y_i=1]=\operatorname{P}[a_i>b_i+c_iX_i]=\operatorname{P}\left[\frac{a_i-b_i}{c_i}>X_i\right] \] によって意思決定が決まるとした場合，\(d_i:=\frac{a_i-b_i}{c_i}\) と潜在変数 \(x_i\) との大小によって応答が変わると要約できる．

すると \(d_i\) の分布関数 \(F_i\) がロジスティックであるか正規分布であるか \(t\)-分布であるかに依って，ロジスティック回帰・プロビット回帰・ロビット回帰というモデルが得られる．

3.2.3 Tobit モデル

さらに \(y_i\) の観測値が \(x_i\) の値に依存するというモデルが (Tobin, 1958) によって提案されている： \[ Y=Y^*\lor0. \] \[ Y^*_i=\beta X_i+\epsilon_i, \]

このモデルは（１型）トービットモデルの他に，打ち切り回帰 (censored regression) とも呼ばれる．

ある閾値を超えた場合にプライバシーの問題で公開を制限する top-coding データなどに用いられる．

3.2.4 Heckman モデル

(Heckman, 1979) は標本に入るかどうかのメカニズムをモデリングする，選択バイアスのためのモデル \[ Y=\begin{cases}Y^*&S=1\\\text{missing}&S=0\end{cases} \] \[ Y^*=\beta X+\epsilon, \] \[ S=1_{\mathbb{R}_+}(S^*),\qquad S^*=\gamma Z+\eta. \] を提案した．これは２型トービットモデルとも呼ばれる．

そしてこの共分散構造 \[ \begin{pmatrix}\epsilon\\\eta\end{pmatrix}\sim\mathrm{N}\biggr(0,\begin{pmatrix}\sigma^2_\epsilon&\sigma_{\epsilon\eta}\\\sigma_{\epsilon\eta}&1\end{pmatrix}\biggl) \] を局外母数として \(\beta\) を一致推定する．

4 文献紹介

(12 章 Gelman and Hill, 2006) が筆者の知る限り最も丁寧な階層モデルへの導入である．

変動係数モデルに関して (Sugasawa and Kubokawa, 2023) が極めて見通しが良い．小地域推定を主なテーマとしている．

項目応答モデルと理想点モデルは (14.3 Gelman and Hill, 2006) において一般化線型階層モデルの文脈で扱われている．

選択モデルは (27 章 Hansen, 2022) を参照．

References

Bafumi, J., Gelman, A., Park, D. K., and Kaplan, N. (2005). Practical Issues in Implementing and Understanding Bayesian Ideal Point Estimation. Political Analysis, 13(2), 171–187.
Battese, G. E., Harter, R. M., and Fuller, W. A. (1988). An error-components model for prediction of county crop areas using survey and satellite data. Journal of the American Statistical Association, 83(401), 28–36.
Carlin, J. B., Wolfe, R., Brown, C. H., and Gelman, A. (2001). A case study on the choice, interpretation and checking of multilevel models for longitudinal binary outcomes. Biostatistics, 2(4), 397–416.
Fay, R. E., and Herriot, R. A. (1979). Estimates of income for small places: An application of james-stein procedures to census data. Journal of the American Statistical Association, 74(366), 269–277.
Fox, J.-P. (2010). Bayesian Item Response Modeling. Springer New York.
Gelman, A. (2005). Analysis of variance—why it is more important than ever. The Annals of Statistics, 33(1), 1–53.
Gelman, A., and Hill, J. (2006). Data analysis using regression and multilevel/hierarchical models. Cambridge University Press.
Gelman, A., Hill, J., and Vehtari, A. (2020). Regression and other stories. Cambridge University Press.
Hansen, B. E. (2022). Econometrics. Princeton University Press.
Heckman, J. J. (1979). Sample selection bias as a specification error. Econometrica, 47(1), 153–161.
Henderson, C. R. (1950). Estimation of genetic parameters. The Annals of Mathematical Statistics, 21(2), 302–315.
Lindley, D. V., and Smith, A. F. M. (1972). Bayes estimates for the linear model. Journal of the Royal Statistical Society. Series B (Methodological), 34(1), 1–41.
Rasch, G. W. (1960). Studies in mathematical psychology: I. Probabilistic models for some intelligence and attainment tests. Nielsen & Lydiche.
Sugasawa, S., and Kubokawa, T. (2023). Mixed-effects models and small area estimation. Springer Singapore.
Tobin, J. (1958). Estimation of relationships for limited dependent variables. Econometrica, 26(1), 24–36.
久保川達也. (2006). 線形混合モデルと小地域の推定. 応用統計学, 35(3), 139–161.

Footnotes

  1. (11.5 節 Gelman and Hill, 2006, p. 246) も参照．↩︎

  2. 一般に階層モデルは，全ての説明変数を組み込んだモデルから，特定の係数の部分集合を括り出し，その部分集合に super-population model を仮定したものと見れる．(13.6 節 Gelman and Hill, 2006, p. 296) の議論も参照．↩︎

  3. 例えば小地域推定においては極めて自然なモデルである (久保川達也, 2006), (Sugasawa and Kubokawa, 2023) も参照．↩︎

  4. この点については (11.6 節 Gelman and Hill, 2006, p. 248)，(Carlin et al., 2001) などを参照．↩︎

  5. この第二階層の構造を \(\alpha_j\sim\mathrm{N}(\mu_\alpha,\sigma^2_\alpha)\) という事前分布とみて，この未知パラメータの推定から始める点推定法を経験ベイズという (Sugasawa and Kubokawa, 2023, pp. 11–12)．↩︎

  6. (Sugasawa and Kubokawa, 2023) では，\(\beta\) が \(j[i]\in[J]\) 毎に異なるとき，これを変量係数モデルと呼んでいる．ここでは (Gelman, 2005) などに倣って，全てまとめて変動係数 (varying coefficient) モデルと呼ぶことにする．↩︎

  7. (Sugasawa and Kubokawa, 2023, p. 10) に GLS 推定量が BLUE であることの証明がある．変量効果推定量については (17.6 Hansen, 2022, pp. 601–603) を参照した．↩︎

  8. 正規性の仮定の下では一般に経験ベイズ推定量と一致することが (Sugasawa and Kubokawa, 2023, pp. 11–12) で示されている．↩︎