Sugasawa, S., Kubokawa, T. (2020) Small area estimation with mixed models: a review. Japanese Journal of Statistics and Data Science. 3, 693–720.
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 自然指数族に基づくモデル
地域レベルモデルで離散データの場合は自然指数族に基づくモデルを組むという手もある。
[おさらいすると、指数型分布族の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)\)とすれば、確かに指数型分布族ではある。うーん、でもどういう理由でこの式になったんだろうか？]
3.5 備考
[略]
3. 小地域推定量の不確実性の測定
[MSEの推定の話と信頼区間の話。4ページにわたりまるごとスキップ]