はじめに
統計学の勉強，研究を行っている大学院生です．
Chan et al. (2019) の教科書（下のリンクに紹介があります）を読んでいる際に，Chapter 14 Latent Variable Models で扱われた順序プロビットモデルのベイズ推定に関して興味を持ち，具体的な計算方法を調べてみました．
順序プロビットモデルは，被説明変数が順序尺度である場合に用いられる標準的な計量経済学モデルです．代表的な応用例としては以下のようなものがあります．
- 企業の信用格付け（AAA, AA, A, ...）の分析
- アンケート（非常に満足，満足，普通，不満，非常に不満）の分析
- 厚生経済学：主観的幸福度の決定要因分析
この記事では基本的な方法とその問題点について説明し，次の記事以降で具体的な解決方法を説明します．線形回帰モデルのベイズ推定や代表的な MCMC の方法（MH, Gibbs, データ拡大など）は既知とします．
モデル
観測される順序データを $y$ として，以下のようなカットオフ $\alpha$ を用いた潜在変数表現を考える．
z_i = x_i \beta + \epsilon_i, \quad \epsilon_i \stackrel{iid}{\sim} N(0, 1), \quad i=1,\dots,n
ここで，$z_i$ は観測されない連続的な潜在変数で，$x_i$ は $1 \times k$ の共変量ベクトル，$\beta$ は $k \times 1$ の係数ベクトルである．識別のため誤差項の分散は $1$ であるとする．観測される離散変数 $y_i \in {1, 2, \dots, M}$ は，潜在変数 $z_i$ と閾値パラメータ $\alpha = (\alpha_0, \alpha_1, \dots, \alpha_M)'$ に基づいて以下のように決定される．
y_i = \begin{cases}
1 & \text{if } \alpha_0 < z_i \le \alpha_1, \\
2 & \text{if } \alpha_1 < z_i \le \alpha_2, \\
\vdots & \vdots \\
M & \text{if } \alpha_{M-1} < z_i \le \alpha_M.
\end{cases}
識別上の制約として，$\alpha_0 = -\infty$，$\alpha_M = \infty$ とおき，$\alpha_1 = 0$ とする．このモデルの尤度関数は以下のように記述される．
観測値 $y_i=k$ となる確率は，
\begin{aligned}
\Pr(y_i = k|x_i, \beta, \alpha) &= \Pr(\alpha_{k-1} < z_i \le \alpha_k) \\
&= \Pr(\alpha_{k-1} < x_i\beta + \epsilon_i \le \alpha_k) \\
&= \Phi(\alpha_k - x_i\beta) - \Phi(\alpha_{k-1} - x_i\beta)
\end{aligned}
となる．ここで $\Phi(\cdot)$ は標準正規分布の累積分布関数である．したがって，尤度関数は次式で与えられる．
L(\beta, \alpha) = \prod_{i=1}^n \sum_{k=1}^M I(y_i = k) \left[ \Phi(\alpha_k - x_i\beta) - \Phi(\alpha_{k-1} - x_i\beta) \right]
ベイズ推定
本モデルの推定において，最尤推定ではなくベイズ推定を採用することには，順序データ特有の課題に対処する上で以下のような利点がある．
第一に，計算が単純化できることである．尤度関数が非線形であることから，多変量モデルや変量効果モデルへの拡張時には数値積分が必要となり，最適化が困難になる．後ほど確認するように Albert and Chib (1993) のデータ拡大の方法では，潜在変数 $z_i$ を未知パラメータとして扱うことで，この問題を標準的な線形回帰モデルと切断正規分布からのサンプリングという単純な問題に帰着させることができる．
第ニに，係数以外の値の不確実性評価がしやすいことである．順序プロビットモデルを用いる実証分析では係数 $\beta$ そのものよりも，限界効果や特定のカテゴリの確率に関心のあることが多い．ベイズ推定を用いると，事後分布からのサンプリングを通じて，これらの非線形な関数の事後分布や信用区間を簡単に導出できる．
AC93 の標準的なデータ拡大法
Albert and Chib (1993) は潜在変数 $z_i$ を未知パラメータとして扱うデータ拡大法を提案した．これにより，切断正規分布と標準的な線形回帰の事後分布を利用したギブスサンプリングが可能となった．
識別制約として $\alpha_0 = -\infty, \alpha_1 = 0, \alpha_M = \infty$ を設定し，推定すべきカットオフパラメータを $\tilde{\alpha} = (\alpha_2, \dots, \alpha_{M-1})'$ とする．平坦な事前分布 $p(\beta, \tilde{\alpha}) \propto 1$ を仮定すると，事後分布は以下のように書ける．
\begin{aligned}
p(\beta, \tilde{\alpha}, z | y) &\propto p(y|z, \beta, \tilde{\alpha}) p(z|\beta, \tilde{\alpha}) p(\beta, \tilde{\alpha}) \\
&\propto \left\{ \prod_{i=1}^n I(\alpha_{y_i-1} < z_i \le \alpha_{y_i}) \right\} \left\{ \prod_{i=1}^n \phi(z_i - x_i\beta) \right\}
\end{aligned}
ここで，$p(y|z, \beta, \tilde{\alpha})$ は $z$ と $\alpha$ が与えられた下で $y$ が特定の値に決まるため，定義関数 $I(\cdot)$ の積となる．
サンプリングの方法
以下のステップ 1 から 3 を繰り返すことで，事後分布からのサンプリングを行う．
1. 係数 $\beta$ のサンプリング
$\tilde{\alpha}$ と $z$ を所与とすると，標準的な線形回帰モデル $z = X\beta + \epsilon$ のベイズ推定に帰着する．
\beta | z, \tilde{\alpha}, y \sim N((X'X)^{-1}X'z, (X'X)^{-1})
2. 潜在変数 $z$ のサンプリング
各 $i$ について，$z_i$ は条件付き独立であり，観測されたカテゴリ $y_i$ に対応する区間 $(\alpha_{y_i-1}, \alpha_{y_i}]$ で制限された切断正規分布に従う．
z_i | \beta, \tilde{\alpha}, y \sim TN_{(\alpha_{y_i-1}, \alpha_{y_i}]}(x_i\beta, 1), \quad i = 1, \dots, n
3. カットオフ $\tilde{\alpha}$ のサンプリング
各カットオフ $\alpha_k$ ($k=2, \dots, M-1$) は，順序制約 $\alpha_{k-1} < \alpha_k < \alpha_{k+1}$ および潜在変数 $z$ との整合性を満たす範囲の一様分布からサンプリングされる．
条件付き分布は以下のようになる．
\alpha_k | \alpha_{-k}, \beta, z, y \sim U[\underline{\alpha}_k, \bar{\alpha}_k]
ここで，下限と上限は以下で定義される．
\begin{aligned}
\underline{\alpha}_k &= \max \left( \alpha_{k-1}, \max_{\{i: y_i = k\}} z_i \right) \\
\bar{\alpha}_k &= \min \left( \alpha_{k+1}, \min_{\{i: y_i = k+1\}} z_i \right)
\end{aligned}
計算における課題
Albert and Chib (1993) による標準的なデータ拡大法は実装が容易である反面，MCMC の収束が遅くなるという欠点が知られている．
原因は，カットオフ $\alpha$ の条件付き分布にある．$\alpha_k$ のサンプリング範囲 $[\underline{\alpha}_k, \bar{\alpha}_k]$ は，そのカテゴリの潜在変数 $z$ の最大値と，次の潜在変数の最小値で決まる．$n$ が増加すると，この区間は狭くなり，$\alpha_k$ から値が動きにくくなる．これにより，連鎖の自己相関が高くなり，ESS が低下する．
この問題を具体的に確認するため，簡単な数値実験を行う．以下の設定で順序プロビットモデルから $n=500$ 個の観測データを生成した．
\begin{aligned}
z_i &= 0.5 + 0.3x_i + \epsilon_i, \quad x_i \stackrel{iid}{\sim} N(0, 1), \quad \epsilon_i \stackrel{iid}{\sim} N(0, 1) \\
y_i &= \begin{cases}
1 & (z_i \le 0) \\
2 & (0 < z_i \le \alpha_2) \\
3 & (\alpha_2 < z_i)
\end{cases}
\end{aligned}
ここで真のカットオフの値は $\alpha_2 = 1$ と設定している．推定には Albert and Chib (1993) の標準的なギブスサンプリングを用い，初期値を $\alpha_2 = 0.5$，$\beta = [0, 0]'$ とした．2500 回のサンプリングを行い，バーンインとして最初の 500 回を除去している．
この実験におけるカットポイント $\alpha_2$ の推定結果を Figure 1 に示す．
Figure 1 は，標準的なデータ拡大におけるサンプリング効率の悪さを視覚化している． まず，（a）のトレースプロットを確認すると，連鎖がパラメータ空間内を緩慢に移動している挙動が見て取れる．実際，(b) の自己相関関数を見ると，ラグが増加しても相関は緩やかにしか減衰しておらず，ラグが 40 離れた時点でも 0.8 に近い高い相関を維持している．その結果，(c) の事後分布は真値周辺をカバーしているものの形が歪んでいる．
この視覚的な評価は，以下の Table 1 の定量的な指標によって確認できる．特にカットオフパラメータ $\alpha_2$ に注目すると，事後平均値（1.034）自体は真値（1.0）を捉えているものの，ESS は 2000 回のサンプリングに対してわずか 15 である．これは標準的なデータ拡大法ではカットオフの事後分布を適切に探索できていないことを意味している．
おわりに
このように，単純な設定であっても標準的な手法ではカットオフパラメータの計算効率に課題が残ります．
以降の記事ではこのサンプリング効率を改善するためにどのような方法が考案されてきたのかのレビューを行います．
(2) ではモデルを再パラメータ化して推定する方法
(3) では MH アルゴリズムを用いる方法と prameter expantion を用いる方法
(4) では EM アルゴリズムと漸近分布を用いた MH アルゴリズムの方法
について紹介する予定です．
(5) の記事では数値実験を行い，これらの手法の計算効率について比較します．
執筆，実装に関して一部 Gemini を使用しました．
References
- Albert, J. H., & Chib, S. (1993). Bayesian analysis of binary and polychotomous response data. Journal of the American statistical Association, 88(422), 669-679.
- Chan, J., Koop, G., Poirier, D. J., & Tobias, J. L. (2019). Bayesian Econometric Methods (2nd ed.). Cambridge University Press.