Sugasawa, S., Kubokawa, T. (2020) Small area estimation with mixed models: a review. Japanese Journal of Statistics and Data Science. 3, 693–720.
勉強のためにめくったやつ。小地域推定という地味な分野にも、なんか新しい展開とかあるんじゃないかなと思って。
小地域推定のうち混合モデルを使う手法についてのレビュー論文(小地域推定の多くの議論は混合モデルに基づくので、ほぼ小地域推定レビューといえよう)。第二著者は「現代数理統計学の基礎」の著者の先生であろう。
2章が基礎、3章が分散推定の話。4章が本題である。
1. イントロダクション
[メモ省略。要するに、混合モデルを使った小地域推定について近年の発展を含めレビューするよという話]
2. 小地域推定のための基本的な混合モデル
2.1 Fay-Herriotモデル
小地域の真の平均を\(\theta_1, \ldots, \theta_m\)とし, 直接推定量\(y_1, \ldots, y_m\)はあるんだけど分散が大きいとする。地域\(i\)の既知の地域特性(切片込み)を\(\mathbf{x}_i\)とする。ご存じFHモデルとは$$ y_i = \theta_i + \epsilon_i$$ $$ \theta_i = \mathbf{x}^\top_i \beta + v_i$$ \(\epsilon_i, v_i\)は独立に \(\epsilon_i \sim N(0, D_i), v_i \sim N(0, A)\)とし、\(D_i\)は既知とする。[あれ？ 本来のFHモデルって撹乱項の正規性を仮定してたっけ？ Rao & Molina 本では仮定してなかったと思うんだけどな]
平方誤差損失の下での\(\theta_i\)の最良予測量は $$ E[\theta_i | y_i] = \gamma_i y_i + (1-\gamma_i) \mathbf{x}^\top_i \beta $$ ただし \(\gamma_i = A/(A + D_i)\)。
モデルの下での\(y_i\)の周辺分布は\(N(\mathbf{x}^\top_i \beta, A+D_i)\)である[これって\(\mathbf{x}^\top\)は所与だけど\(y_i\)は未知の場合ってことだよね？]。\(\beta\)のGLS推定量は$$ \hat{\beta}_{GLS} = \left( \sum_{i=1}^m \frac{\mathbf{x}_i \mathbf{x}^\top_i}{A+D_i} \right)^{-1} \left(\sum_{i=1}^m \frac{\mathbf{x}_i y_i}{A+D_i} \right) $$ となる[GLSなんだから撹乱項の正規性は要らなくない？ 結局、正規性の仮定ってどこで効くんだろうか…]。
これを上の式に放り込んだ $$ \tilde{\theta}_i \equiv \gamma_i y_i + (1-\gamma_i) \mathbf{x}^\top_i \hat{\beta}_{GLS} $$ が\(\theta_i\)のBLUPとなる。安定的な推定量\(\mathbf{x}^\top_i \hat{\beta}_{GLS}\)への縮小推定量になっている。
実際には\(A\)は未知なので、\(\gamma_i, \hat{\beta}_{GLS}\)の推定には標本を使わないといけない。頻度主義的にはEBLUP, ベイジアン的には経験ベイズ推定量ということになる。\(A\)の推定量としてはML, REML, モーメント法がある。
これに対し、近年ではHBが主流となっている。\(\beta, A\)に事前分布を与えるわけだ。
2.2 nested error回帰モデル
今度はユニットレベルのデータが利用可能な場合。
地域\(i\)のユニットレベル標本を\(y_{i1}, \ldots, y_{i n_i}\), 共変量を\(\mathbf{x}_{i1}, \ldots, \mathbf{x}_{i n_i}\)とする。次のモデルを考える: $$ y_{ij} = \mathbf{x}^\top_{ij} \beta + v_i + \epsilon_{ij}$$ \(v_i, \epsilon_{ij}\)は独立に\(v_i \sim N(0 ,\tau^2), \epsilon_{ij} \sim N(0, \sigma^2)\)とする。こういうのをnested error回帰(NER)モデルという。
地域内の\(y_{ij}\)は独立でないという点にご注目。\(j \neq j’\)について \( Cov(y_{ij}, y_{ij’}) = \tau^2\)となる。
NERモデルはふつう有限母集団モデルの枠組みで用いられる。超母集団モデルは$$ Y_{ij} = \mathbf{x}^\top_{ij} \beta + v_i + \epsilon_{ij}$$で、\(Y_{i1}, \ldots, Y_{i N_i}\)のうち\(n_i\)個の実現値\(y_{i1}, \ldots, y_{i n_i}\)だけが観察されていると考えるわけだ。真の平均は $$ \frac{1}{N_i} \sum_{j=1}^{N_i} Y_{ij} = \bar{X}^\top_i \beta + v_i + \frac{1}{N_i} \sum_{j=1}^{N_i} \epsilon_{ij} $$ と書ける。ふつう\(N_i\)は十分に大きいから第3項を無視し、\(\theta_i = \bar{X}^\top_i \beta + v_i\)と定義する。つまり、補助変数の真の平均ベクトルは既知でないといけないわけね。
先ほどのNERモデルの下で、\(v_i\)の最良予測量は、\(\bar{y}_i, \bar{\mathbf{x}}_i\)を標本平均として $$ \tilde{v}_i = \frac{n_i \tau^2}{\sigma^2+n_i \tau^2} (\bar{y} – \bar{\mathbf{x}}^\top_i \beta) $$ となる。\(\beta\)は全標本を使ってGLS推定する。\(\tau^2, \sigma^2\)は正規性を仮定して推定する。\(v_i\)のEBLUP推定量\(\hat{v}_i\)が手に入ったら、\( \hat{\theta}_i = \bar{X}^\top_i \hat{\beta}_{GLS} + \hat{v}_i \)を求めればよい。
2.3 抽出モデルとリンキングモデルがアンマッチな場合
たとえば年収みたいに、正規性が仮定できない場合はどうするか。FHモデルなら、$$ y_i = \theta_i + \epsilon_i $$ $$ h(\theta_i) = \mathbf{x}^\top_i \beta + v_i$$ という手がある(\(h(\cdot)\)はたとえば\(\log(\cdot)\))。\(y_i\)の周辺分布はこうなる。\(\phi(\cdot; a, b)\)を\(N(a, b)\)の密度関数として$$ \int \phi(y_i; \theta, D_i) \phi(h(\theta_i); \mathbf{x}^\top_i \beta, A) h'(\theta_i) d \theta_i $$ HB推定とEB推定がある。
2.4 一般化線形混合モデル
二値、カウント、カテゴリカル変数とかだったらどうするか。一般化線形混合モデルをHB推定という手もあるが、ここでは二値についての頻度主義アプローチを紹介しよう。
ユニットレベルモデルを, \(P(y_{ij} = 1 | p_{ij}) = p_{ij} \)として $$ \log \left( \frac{p_{ij}}{1-p_{ij}} \right) = \mathbf{x}^\top_i \beta + u_i, \ \ u_i \sim N(0, \tau_i)$$ とする。推定対象を \(\bar{p}_i = \frac{1}{N_i} \sum_{j=1}^{N_i} y_{ij} \)としよう。\(u_i\)の(二乗損失のもとでの)最良推定量はこうなる。\(z \sim N(0, 1)\)として、$$ E[u_i | y_i] = \tau \frac{E_z [z \phi_i(z; \beta, \tau^2)]}{E_z [\phi_i(z; \beta, \tau^2)]} $$ ただし $$ \log \phi_i (z; \beta, \tau^2) = \tau z \sum_{j=1}^{n_i} y_{ij} – \sum_{j=1}^{n_i} \log \{1 + \exp(\mathbf{x}^\top_{ij} \beta + \tau z) \} $$ であり、期待値の添え字\(z\)は\(z\)に関する期待値という意味である。このように解析的にかけず、数値積分が必要になる。
[\(u_i\)の推定量はなんでこうなるのかな… わからんけど、ゆっくり考えればきっとわかるのだろう、と信じて先に進もう]
2.5 自然指数族に基づくモデル
[この節はけっこう大事な話なんだろうけど、残念ながら理解が追い付かなかった。Ghosh & Maiti (2004 Biometrika)を読むのがよいらしい]
地域レベルモデルで離散データ[?]の場合は自然指数族に基づくモデルを組むという手もある。
[おさらいすると、指数型分布族のPDFは、パラメータがひとつなら $$ f(x|\theta) = h(x) \exp[ \eta(\theta) T(x) – A(\theta)]$$ という形で書ける。正規分布(平均のみ未知)をこの形で書けば、\(h(x) = \frac{1}{\sqrt{2 \pi \sigma^2}} \exp \left( – \frac{x^2}{2\sigma^2} \right), \eta(\mu) = \frac{\mu}{\sigma}, T(x) = \frac{x}{\sigma}, A(\mu) = \frac{\mu^2}{2\sigma^2} \)とすればつじつまがあう。でも、\(h(x)\)を消して、代わりに指数関数の内側に\(\log h(x)\)を足してもいいですよね。つまり$$ f(x|\theta) = \exp[ \eta(\theta) T(x) – A(\theta) + B(x)]$$という形でもよい]
\(y_i\)の条件付き分布を $$ f(y_i | \theta_i) = \exp [ n_i \{ \theta_i y_i – \psi(\theta_i) \} + c(y_i, n_i) ] $$ とする。\(\psi(\cdot)\)はリンク関数、\(n_i\)は既知で通常は標本サイズである。[\(\eta(\theta_i) = \theta_i, T(y_i) = n y_i, A(\theta_i) = n_i \psi(\theta_i), B(y_i) = c(y_i,n_i)\)とすれば、確かに指数型分布族ではある。うーん、でもどういう理由でこの式になったんだろうか？]
\(\theta_i\)の周辺分布を、$$ \pi(\theta_i) = \exp [ v \{m_i \theta – \psi(\theta_i) \} ] C(v, m_i) $$ ただし \(m_i = \psi'(\mathbf{x}^\top_i \beta) \)とする。これは共役事前分布である。[あ、ベイジアンなのか。それにしても、この式もなぜこうなったのかわからないので辛い]
真の地域平均は\(\mu_i \equiv E[y_i | \theta_i] = \psi'(\theta_i) \)である。\(E[\theta_i] = m_i\)である。
共役事前分布のおかげで周辺尤度が閉形式で書ける。[…] \(\mu\)の条件付き期待値は、\(\hat{m} = \psi'(\mathbf{x}^\top_i \beta)\)として $$ \hat{\mu}_i = \frac{n_i y_i + \hat{v} \hat{m}_i}{n_i + \hat{v}} $$ となる。つまり、直接推定量\(y_i\)と回帰推定量\(\hat{m}\)の重み付き平均になる。
\(v = A^{-1}, n_i = D^{-1}, \psi(x) = x^2/2\)とすればこれはFHモデルになることにご注目。\(\psi(x) = \exp(x)\)とすればポアソン-ガンマモデルに、\(\psi(x) = \log(1+\exp(x))\)とすれば二項-ベータモデルになる。
3.5 備考
[略]
3. 小地域推定量の不確実性の測定
[MSEの推定の話と信頼区間の話。4ページにわたりまるごとスキップ]
4. 小地域推定の近年の発展
[いよいよ本題、なんだけど、12節もあるぞ…]
4.1 調整尤度法
FHモデルでランダム効果の分散\(A\)を推定する際、最尤推定量では0になっちゃうことがある。するとBLUPは \(\tilde{\theta}_i = \mathbf{x}^\top_i \hat{\beta}_{GLS} \)になってしまう。これを避けるため、\(A\)の尤度に調整ファクター\(h(A)\)をかけた \( L_{ad}(A) = h(A) L(A)\)を尤度にするという方法がある。\(A=0\)のときに\(h(A)\)が0に近づくようにするわけだ。\(h(A)=A\)という提案や、ほかの提案がある。
MSE推定や信頼区間構築でもこういう方法の提案があって… [パス]
4.2 有限母集団の一般パラメータの推定
NERモデルの話。推定対象をより一般的に $$ \mu_i \frac{1}{N_i} \sum_{j=1}^{N_i} T(y_{ij}) $$ としますね。たとえば貧困率とか。超母集団モデルを$$ H(Y_{ij}) = \mathbf{x}^\top_{ij} \beta + v_i + \epsilon_{ij} $$ とする。標本を\(s_i\), 残りを\(r_i\)と書いて、\(r_i\)における\(H(Y_{ij})\)の条件付き分布を $$ H(Y_{ij}) | (y_{ij}, j \in s_i) \sim N(\mathbf{x}^\top_{ij} \beta + \tilde{v}_i, \sigma^2) $$ $$ \tilde{v}_i = \frac{n_i \tau^2}{\sigma^2 + n_i \tau^2} \left( \frac{1}{n_i} \sum_{j=1}^{n_i} H(y_{ij}) – \bar{\mathbf{x}}^\top_i \beta \right) $$ とする。
\(\mu_i\)の最良予測量は $$ E[\mu_i | y_{ij}, j \in s_i] = \frac{1}{N_i} \left\{ \sum_{j \in s_i} T(y_{ij}) + \sum_{j \in r_i} E[T(Y_{ij} | y_{ij}, j \in s_i] \right\} $$ かっこ内の第2項が解析的に書けないのでMCMCで求める。なお、\(\mathbf{x}_{ij}\)の完全情報が必要になる点に注意。
[これって、最近の発展というより、Molina&Rao本にも載ってる話なのでは…?]
4.3 反応変数のパラメトリック変換
[FHモデルとかNERモデルとかで、反応変数に(対数変換をより一般化したような)パラメトリック変換をかけるという話。パス]
4.4 分散の部分についての柔軟なモデリング
FHモデルで標本抽出分散\(D_i\)をデータから推定しているとリスクの過小評価とかが起きる。そこでこういう同時階層モデルを組もうという提案がある。$$ \theta_i \sim N(\mathbf{x}^\top_i \beta, A)$$ $$ \sigma^2_i \sim\ \pi(\sigma^2) $$ $$ y_i | \theta_i, \sigma^2_i \sim N(\theta_i, \sigma^2_i) $$ $$ D_i | \sigma^2_i \sim Ga \left( \frac{n_i-1}{2}, \frac{n_i-1}{2\sigma^2} \right) $$ ところが、\(\pi(\sigma^2)\)を無情報にすると(たとえば\(Ga(a_1, b_i)\)で\(a_i, b_i\)を小さな定数にしたりすると)、結局 \(\sigma^2_i\)は\(D_i\)とほぼ等しくなってしまったりして、縮小特性が消えてしまう。そこで \(\sigma^2_i \sim Ga(\alpha, \gamma) \)として事後分布を効率的に求める方法とかが提案されている。
NERモデルでも、ランダム効果と誤差項の分散パラメータを地域によって変えるという提案があって… [ふーん。以下パス]
4.5 ランダム効果の柔軟なモデリング
ランダム効果というのはほんとにあるのかどうかわからんわけで、FHモデルをこう拡張するという提案がある。二値潜在確率変数\(s_i\)を導入して $$ v_i | (s_i = 1) \sim N(0, A) $$ $$ v_i | (s_i = 0) \sim \delta_0 $$ $$ P(s_i = 1) = \pi $$ \(\delta_0\)というのは1点の分布である。変数選択でいうところのspike-and-slab事前分布に近い。
[なるほど、いろんなことを考えるものね… 以下パス]
4.6 誤差項とランダム効果についての歪んだ分布の仮定
FHモデルで\(\epsilon_i\)についてskew-normal分布を仮定するという提案があって… [うわあ、めんどくせえええ。以下パス]
4.7 ノンパラ・セミパラなモデリング
モデルの誤指定が怖いので、NERモデルで$$ y_{ij} = f(x_{ij}) + v_i + \epsilon_{ij} $$ として、\(f(x_{ij}\)のところをP-スプラインにしようという提案がある。すなわち、ノット\(\kappa_1 \lt \ldots \lt \kappa_K\)を決めて$$ f(x) = \beta_0 + \beta_1 x + \cdots + \beta_q x^q + \sum_{l=1}^K \gamma_l (x-\kappa_l)^q_+ $$ $$ \mathbf{\gamma} \sim N(0, \lambda \mathbf{I}_K) $$ という感じ。[…後略…]
リンク関数の誤指定が怖いよねという話もあって… [パス]
ランダム効果や誤差項の分布を指定したくないからディリクレ過程mixtureモデルを使おうという提案もあって… [パス]
4.8 頑健な手法
外れ値が怖いからロバストなBLUPを使おうという提案があって [… すいません、気力も知力も足りないものでまるごとパス]
4.9 測定誤差モデル
FHモデルで、\(\mathbf{x}_i\)の真値はわからず、\(\hat{\mathbf{x}}_i\)しかわかんないけどそのMSE行列はわかるという状況で… [なるほどね。パス。NERモデルでもそういうモデルがある由]
4.10 観察された最良予測
FHモデルで、\(\mathbf{x}\top_i \beta\)のところのモデル誤指定が怖いので、真のモデルを $$ \theta_i = \mu_i + v_i, \ \ v_i \sim N(0, A)$$ として、ほんとは運悪く\(\mu_i\)が\(\mathbf{x}^\top_i \beta\)じゃないときでも\(\theta_i\)をうまく推定できるように\(\hat{\beta}\)を推定する、という提案がある。
\(B_i = A/(A+D_i)\)とする。全平均予測誤差を $$ MSPE(\beta) = \sum_{i = 1}^m \left( B_i y_i – \theta_i + \mathbf{x}^\top_i \beta(1-B_i) \right)^2 $$ とする。その期待値を最小化する\(\beta\)の推定量は $$ \hat{\beta}_{OBP} = \left( \sum_{i=1}^m (1-B_i)^2 \mathbf{x}^\top_i \mathbf{x}_i \right)^{-1} \sum_{i=1}^m (1-B_i)^2 \mathbf{x}_i y_i $$ となる。これを観察された最良予測(OBP)推定量という。最尤推定量やGLS推定量とは異なるという点に注意。[へえええええ！ そんなのあるんだ。なるほど、GLS推定量とは式が異なる。これって普通の推定量とくらべてどうなりやすいんだろうか。直感的には、0に近くなりやすく、共変量間に正の相関があるときにはどの共変量に対しても同じような値になるような気がするんだけど]
4.11 変数選択
回帰モデルの部分の変数選択が大事だが、推定対象にランダム効果が入っているのでAICやBICによる選択は必ずしも合理的でない。そこでcAICという提案がある。[..中略…] 詳しくは Muller et al.(2013 Stat.Sci.)を読むといいよ。
NERモデルでは\(\bar{\mathbf{X}}_i\)による地域レベル予測と\(\mathbf{x}_{ij}\)によるNERモデルとで共変量が異なりうるので、cAICをさらに修正するという提案がある。[… 以下パス]
4.12 そのほかのトピック
- ベンチマーキング。すなわち、EBLUP推定量の合計を直接推定値に合わせること。
- グループ化されたデータやcirlularデータの小地域推定。
- 下位地域があるときのFHモデル。
- 空間データ、時系列データへの拡張。
5. まとめ
[略]
——————
勉強になりましたですが、最近の発展に関する話題のなかに機械学習は出てこないのね… あっ、「混合モデルによる小地域推定」のレビューだから当然か！ (いま気が付いた)