Datta, G., & Ghosh (2012) Small Area Shrinkage Estimation. Statistical Science, 27(1), 95-114.

　仕事の都合でMrPモデルを推定していて思ったんだけど、地域の母平均推定量は地域の標本サイズが小さいときにふつう全体の母平均推定値に向かってシュリンクしていく。これはMrP(つまり階層回帰モデル)の性質というより、小地域推定に用いられるいろんな方法に共通した性質だし(地域レベルの単純な方法、たとえばFay-HerriotモデルやJames-Stein推定量もそうだ)、他のいろんな推定でもそういうことが起きる(回帰の正則化とか)。
　小地域推定のいろんな手法を、縮小特性という角度からレビューしたものはないだろうか？ と思って探してみたら、そのものずばりの論文があった。きっと私の能力を超える論文だろうなとは思ったが、試しに読んでみた。仕事を離れた勉強のつもりなので、いきなりメモを取りながら読みはじめ、嫌になったら挫折する気満々、という作戦である。
　掲載誌のこの号は”Minimax Shrinkage Estimation: A Tribute to Charles Stein”という特集号だったようだ。WikipediaによればSteinの死去は2016年。

1. イントロダクション
　[小地域推定の重要性が増しているという話…]

　縮小推定量は小地域推定量よりも長い歴史を持っている。その正確な定義は難しいが、なんらかの標準的な推定量(MLEとか一様最小分散不偏推定量(UMVUE)とか最小二乗推定量とか)を修正して、なんらかの望ましい基準(MSEとか二乗リスクとかバイアスとか)を最小化しようとする推定量のことである。最良線形不偏予測子(BLUP)、経験的最良線形不偏予測子(EBLUP)、経験ベイズ(EB)、階層ベイズ(HB)などが含まれる。それらはふつう、標準的な推定量とか適切なモデルの下で合理的ななんらかの推定量の重み付き平均のの形をとり、ウェイトがなんらかの最適性基準によって決められている。
　縮小推定は小地域推定において自然である。小地域推定では、MLUやUMVUEのような直接推定はSEが大きくなることが多い(元の調査はもっと上の集計レベルでの正確性を目指しているからである)。そこで、他の小地域から「強さを借りてくる」ことが必要になる。

　初期の縮小推定量は地域ごとの直接推定をなんらかの全体推定へと縮小させるものであった。補助情報が手に入るようになると、直接推定量をなんらかの推定された回帰表面に縮小させるのが普通になった。この縮小過程はモデルを必要とする。
　ベイズ推定量は2世紀以上の歴史を持つ。それらは(たとえば事前平均への)縮小推定量とみなしうることが多い。BLUP, EBLUP推定量を提案したのはHenderson(1953)だが、これもなんらかの回帰推定量への縮小推定量である。しかしなんといっても、「縮小」といえばStein(1956)である。彼は3次元以上のMVNの平均ベクトルの縮小推定量が二乗誤差損失の点で標本平均ベクトルより優れていることを示した。Stein推定量はその後Efron & Morris(1973)などによってEBとして解釈された。小地域推定の分野ではFay & Herriot(1979)へと拡張された。

　本論文では１要因のランダム効果回帰モデルに焦点を当て、BLUP, EPLUBをHB, EBと結びつける。これらは小地域推定では地域レベルモデルと呼ばれているものだ(これに対するのが個体レベルモデルで、本論文ではほとんど触れない)。また、EB信頼区間の構築についても述べる。
　詳しくは以下をみよ: Rao(2003, 書籍)、Ghosh & Rao(1994), Pfeffermann(2002), Rao(1999 SurveyMethodol.; 2003 J.Iran.Stat.Soc.), Datta(2009 in “Handbook of Statistics”), Pfeffermann (2002) [Pfeffermann(2013)の前のレビューね]。

2. 釣り合いデータのための縮小推定量
　[読み進めてわかったんだけど、ここで釣り合いbalancedといっているのは、実験計画でいうところの標本サイズの釣り合いというより、むしろ地域レベル直接推定量の分散が地域間で同じだという意味のようだ。ってことは、地域レベル直接推定量が二項分布に従う場合は非釣り合いってことかなあ]

(モデル)
　\(y_i (i=1, \ldots, m)\)を地域レベル推定量とする。以下のモデルを考えよう。$$ y_i | \theta_i \sim_{ind} N(\theta_i, V) $$ $$ \theta_i | A \sim_{ind} N(\mathbf{x}^\top_i \beta, A)$$ ただし\(\mathbf{x}_i\)は計画ベクトル。次のランダム効果モデルに書き直してもよい。$$ y_i = \mathbf{x}_i^\top \beta + u_i + e_i $$ $$ u_i \sim_{iid} N(0, A) $$ $$ e_i \sim_{ind} N(0, V) $$ 行列で書けば
$$ \mathbf{y} = \mathbf{X} \beta + \mathbf{u} + \mathbf{e} $$ \(rank(\mathbf{X}) = p (\lt m)\)と仮定する。
　当然ながら識別問題が生じる。個体レベルモデルならいいですよ、\(V\)は個体レベルデータで推定すればいいから。でもこれは地域レベルモデルなので、\(V\)は既知と仮定するのがふつうである。実務的にはなんらかの平滑化推定(たとえば一般化分散関数アプローチ)を使う。

(Aが既知の場合)
　まずは\(A\)も既知としよう。$$ B = V(V+A)$$ $$ \mathbf{P}_X = \mathbf{X} (\mathbf{X}^\top \mathbf{X})^{-1} \mathbf{X}^\top $$ とする[一見ややこしいが、\(\hat{\beta}= (\mathbf{X}^\top \mathbf{X})^{-1} \mathbf{X}^\top \mathbf{y}\)だから、\(\mathbf{P}_X\)とは要するに\(y\)の前から掛けると \(\mathbf{P}_{X} \mathbf{y} = \mathbf{X} \hat{\beta}\)になるような行列のことである]。
　以下の定理が成り立つ。

定理1. 所与のモデルの下で、\(\theta\)の事後分布は $$ N((1-B)y + B\mathbf{P}_X \mathbf{y}, V((1-B) \mathbf{I}_m + B \mathbf{P}_x)) $$ となる。

　証明: […読み飛ばしたけどそんなに難しくなさそう…]。
　注意1: 上の定理より、事後平均は $$ \hat{\theta}^B = E(\theta | \mathbf{y}) = (1-B) y + B \mathbf{P}_{X} \mathbf{y} $$ つまり直接推定量と回帰推定量の重み付き平均になっていて、重みは標本分散と事前分散の逆数に比例する。つまり、直接推定量を回帰推定量へと縮小させている。その量は\(V/A\)で決まる。なお、これは正規性仮定の下で最良不偏予測量である。
　注意2: もし\(\beta\)が既知なら事後分散は\(V(1-B)\mathbf{I}_m\)となる。つまり、\(VB\mathbf{P}_{x}\)が\(\beta\)が未知であることによる不確実性の増大を表しているわけだ。

　こんどは、分布の仮定を置かず、2つの積率だけがわかっているとしよう。

定理2. \( \hat{\theta}^B\)は\(\theta\)のBLUPである。また以下が成り立つ:$$ E[(\hat{\theta}^B – \theta)(\hat{\theta}^B-\theta)^\top] = V \{ (1-B) \mathbf{I}_m + B \mathbf{P}_x \} $$[BLUPってのはつまり、線形不偏予測量のなかでMSE最小ってことね。後半はMSEが事後分散だと述べている]

　証明: […読み飛ばしたけど、たぶんそんなに難しくないんじゃないかな…]。
　注意3: \(\mathbf{u}, \mathbf{e}\)の正規性のもとでは, BLUPなだけじゃなくて最良不偏予測量でもある。

　定理1, 2は、BLUPとHB予測量との等価性を示している。また、パラメータ\(A\)が既知の場合の、釣り合い1要因ランダム効果モデルの不確実性指標との等価性も示している。

(Aが未知の場合, EB推定量)
　しかし、\(A\)が未知だと話が変わってくる。
　[以下、途中で混乱したのでメモしておくけど、\(m\)は地域数、\(p\)は計画行列のランクね]

　上のランダム効果モデルから $$ \mathbf{y} \sim N(\mathbf{X}\beta, (V+A) \mathbf{I}_m) $$ である。ベイズ的に言えばこれは\(\theta\)を積分消去した\(\mathbf{y}\)の周辺分布である。この周辺PDFに基づくと、\( (\hat{\beta}, S = || \mathbf{y} – \mathbf{X} \hat{\beta} ||^2) \)は\((\beta, A)\)の最小十分統計量である[はい深呼吸！ 線形回帰モデルの回帰係数と誤算分散(のうち未知の部分)を推定しようとするとき、回帰係数の推定量と残差平方和がその最小十分統計量だってことね]。
　\(S \sim (V+A) \chi^2_{m-p} \)だから、\(B = V(V+A)\)のUMVUE [一様最小分散不偏推定量] はこうなる。\(m \ge p+2\) として$$ \hat{B}^{EB} = \frac{V(m-p-2)}{S} $$ これに対応するEB推定量ないしEBLUP推定量はこうなる [\(B\)をデータで決め打ちしてるのでEBなわけね]: $$ \hat{\theta}^{EB} \equiv \hat{\theta}^{EBLUP} = \left[ 1- \frac{V(m-p-2)}{S} \right] \mathbf{y} + \frac{V(m-p-2)}{S} \mathbf{X} \hat{\beta} $$ これはJames-Stein推定量である。
　この推定量に対する批判として、\(\hat{B}^{EB} \)が1を超え、回帰推定量とは逆の方向に動いてしまうかもしれないという点がある。しかし、\(P(\hat{B}^{EB} \ge 1)\)は\(m\)と共に0に近づくので大丈夫。

(Aが未知の場合, HB推定量)
　[このへんから話が難しくなる。ほぼ逐語訳]
　今度はフルベイズアプローチ。事前分布を\(\pi(\beta, A) = 1, \pi(\beta, B) = B^{-2} \)とすると以下となる。$$ \pi(B|\mathbf{y}) \propto B^{\frac{m-p}{2}} \exp\left(-\frac{1}{2V} BS \right) B^{-2} I[0 \lt B \lt 1] = B^{\frac{m-p-4}{2}} \exp\left( -\frac{1}{2V} BS \right) I[0 \lt B \lt 1] $$ [なんでこうなるのかさっぱりわからん。えええええん]
　ここまででは話を簡単にするために\(A\)の事前分布を一様だとしていたんだけど、事後分布が正則になるのならなんでもよい。\(\pi(\beta, A) = A^{-k}\)としたら事後分布は\(k \lt 1, m \gt p-2k+2\)のときに正則になる。つまり、\(\pi(\beta, A) = 1\)ならば\(m \gt p+2\)のとき事後分布は正則になる。\(A^{-1}\)とか\(A^{-2}\)とかだと事後分布は常に非正則になる。
　一様事前分布の場合、\(\theta\)の事後平均は定理1の\(B\)を\(E(B|\mathbf{y})\)に置き換えれば得られる。事後分散は、単に\(B\)を\(E(B|\mathbf{y})\)に置き換えるだけではすまず、\(B\)の推定に伴う不確実性が追加されて次式になる:$$ V(\theta(\mathbf{y})) = V[(1-E(B|\mathbf{y}))\mathbf{I}_m + E(B|\mathbf{y}) P_X] + V(B|\mathbf{y})(\mathbf{y}-P_X \mathbf{y})(\mathbf{y}-P_X \mathbf{y})^\top$$
　この推定量の性質についてみていこう。\(m\)が大きければ以下が成り立つ。$$ E(B|\mathbf{y}) \approx \frac{V(m-p-2)}{S} $$ $$ V(B|\mathbf{y}) \approx \frac{2V(m-p-2)}{S^2} $$ $$ E(\theta|\mathbf{y}) \approx \hat{\theta}^{EB} $$ $$ V(\theta|\mathbf{y}) \approx V \left( 1-\frac{V(m-p-2)}{S} \right) \mathbf{I}_m + \frac{V^2(m-p-2)}{S}\mathbf{P}_X $$ $$ + \frac{2V^2(m-p-2)}{S^2}(\mathbf{y}-\mathbf{X}\hat{\beta})(\mathbf{y}-\mathbf{X}\hat{\beta})^\top$$
　損失関数を二乗誤差\(L(\theta, \mathbf{a}) = || \theta – \mathbf{a} ||^2\)としたときの\(\hat{\theta}^{EB}\)のベイズリスクについて考えてみよう。

定理3 \(m \ge p+2\)とする。すべての\(i\)について\(h_{ii} = \mathbf{x}^\top_i(\mathbf{X}^\top \mathbf{X})^{-1} \mathbf{x}_i\)とする。
$$ E[(\theta_i – \hat{\theta}_i^{EB}]^\ = V(1-B) + VB h_{ii} + \frac{2VB(1-h_{ii})}{m-p} $$ $$ E||\theta – \hat{\theta}^{EB}||^2 = V[m-(m-p-2)B] $$

証明: [パス。なんだか難しそうだ]
注意4: 定理1,2をみるとわかるが、\(A\)が既知ならベイズリスクは\(V[m(1-B)+pB]\)である。つまり、未知の分散成分\(A\)を推定したことによるベイズリスクの超過は\(2VB\)である。[へー]
注意5: [なんか込み入った話なのでパス]

3. 非釣り合いデータのための縮小推定量
(モデル)
　残念ながら、ふつう抽出分散は地域間で等しくない。そこで出てくるのがFay-Herriot(1979)モデルで $$ y_i | \theta_i \sim_{ind} N(\theta_i, V_i) $$ $$ \theta_i \sim_{ind} N(\mathbf{x}_i^\top \beta, A) $$ 彼らは分析にあたってEBアプローチを採用した。以下とする: $$ \mathbf{G} = Diag(V_1, \ldots, V_m) $$ $$ \mathbf{D} = \mathbf{G}+A\mathbf{I}_m$$ $$ B_i = V_i/(V_i + A) $$ $$ \mathbf{B} = \mathbf{G} \mathbf{D}^{-1} = \mathbf{D}^{-1} \mathbf{G} = Diag(B_1, \ldots, B_m)$$

(Aが既知の場合)
　まずは\(\beta, A\)が既知だとしよう。\(\theta\)のベイズ推定量は $$ \hat{\theta}^B = (\mathbf{I}_m – \mathbf{B})\mathbf{y} + \mathbf{B}\mathbf{X} \beta $$ \(\beta\)のGLS推定量は、\(rank(\mathbf{X})=p (\lt m)\)として $$ \tilde{\beta}(A) = (\mathbf{X}^{-1} \mathbf{D}^{-1} \mathbf{X})^{-1} \mathbf{X}^\top \mathbf{D}^{-1} \mathbf{y} $$ $$ = [\mathbf{X}^\top (\mathbf{I}_m-\mathbf{B})\mathbf{X}]^{-1} \mathbf{X}^\top (\mathbf{I}_m – \mathbf{B})\mathbf{y} $$ \(\theta\)のBLUP推定量は \( (\mathbf{I}_m – \mathbf{B}) \mathbf{y} + \mathbf{B} \mathbf{X} \tilde{\beta}(A) \)である。

(Aが未知の場合, EB推定量)
　以下が成り立つ。$$ E \left[ \sum_{i=1}^m \frac{ (y_i – \mathbf{x}_i^\top \tilde{\beta}(A))^2}{V_i +A} \right] = m-p $$ この期待値記号を取っ払った等式を考えよう。左辺(つまり期待値記号の内側)は\(A\)の非増大関数になっているから、\(A=0\)のとき左辺は\(m-p\)より小さい。よって期待値を取っ払った等式に解はない。この場合、推定値を0とする[??? どういうこと？]。そうでない場合は\(A\)に初期値をいれて、\(\hat{\theta}^B
\)の式と等式を繰り返し解くと、\(\hat{A}\)と\(\tilde{\beta}(\hat{A})\)が得られる。こうして以下のEB推定量(EBLUP推定量)が得られる。[途中がわからんかったがまあいいや] $$ \hat{\theta}^{EB} = (\mathbf{I}_m – \hat{\mathbf{B}}) \mathbf{y} + \hat{\mathbf{B}} \mathbf{X} \hat{\beta} $$ $$ \hat{B} = Diag \left( \frac{V_1}{V_1 + \hat{A}}, \ldots, \frac{V_m}{V_m + \hat{A}} \right) $$

(Aが未知の場合, HB推定量)
　以下のHBモデルを考える。$$ y_i | \theta_i, \beta, A \sim_{ind} N(\theta_i, V_i) $$ $$ \theta_i | \beta, A \sim_{ind} N(\mathbf{x}_i^\top \beta, A) $$ $$ \pi(\beta, A) = 1 $$ 同時事後密度は $$ \pi(\theta, \beta, A | \mathbf{y}) \propto A^{-m/2} \exp \left[ -\frac{1}{2} \{ (\mathbf{y} – \theta)^\top \mathbf{G}^{-1} (\mathbf{y} – \theta) + A^{-1} || \theta – \mathbf{X} \beta ||^2 \} \right] $$ \(\theta | \beta, a, \mathbf{y}\)はこうなる[…略…]。\(\beta | A, \mathbf{y}\)はこうなる[…略…]。\(\pi(A | \mathbf{y})\)はこうなる[…略…]。というわけで、以下が得られる。$$ E(\theta | \mathbf{y}) = [\mathbf{I}_m – E(\mathbf(B|\mathbf{y}) ] \mathbf{y} + E[\mathbf{B} H_X \mathbf{y} | \mathbf{Y}] $$ $$ V(\theta | \mathbf{y}) = E[\{\mathbf{I}_m – \mathbf{B}\} \mathbf{G} | \mathbf{y}] + E[\{\mathbf{B}(\mathbf{I}_m – H_X)\} \mathbf{G}|\mathbf{y}] + V[\mathbf{B}\{\mathbf{y} – \mathbf{X} \tilde{\beta}(A)\} | \mathbf{y}] $$ ただし$$ H_X = \mathbf{X}[\mathbf{X}^\top (\mathbf{I}_m – \mathbf{B})\mathbf{X}]^{-1} \mathbf{X}^\top(\mathbf{I}_m-\mathbf{B}) $$である。[完全に写経である。南無南無]
　\(V(\theta | \mathbf{y})\)の第1項は\(\beta, A\)が既知だった場合の事後分散、第2項は\(\beta\)が未知であることによる追加分、第3項は\(A\)が未知であることによる追加分である。

　ベイジアンの場合、事後分散が自然な不確実性指標になっている。頻度主義の場合は、BLUPの平均二乗予測誤差のなかの\(A\)をなんらかの適切な推定量に入れ替えるという手が考えられる。しかし[…中略…]、結局、EBLUPのMSEは過小評価される。[以下、EBLUPのMSEについての議論が延々と続く。正直、力尽きました。まるごとパス]

　実例。[まるごとパスするけど、実データに基づき小地域の母平均推定量をいろんな方法で出して比べ、性質の違いについて議論しているようだ]

　本節の締めくくりとしてひとこと。
　本節では\(V\)を既知としてきた。識別の問題を避けるために必要な仮定である。\(V_i\)について\(Y_i\)とは独立に推定できれば、そして\(V\)が有限のパラメータに依存しているなら、ここまでの結果をモデル・ベースの小地域平均推定とその不確実性の推定へと拡張できる(EBLUPでもHBでも)。また個体レベルモデル(6節)でもそういうことをしている。
　問題は\(V_i\)が有限個のパラメータに依存していない場合である。その場合、ここまでに述べた近似は維持されない。そうした場合についてはWang & Fuller(2003 JASA)をみよ。

4. 拡張
　[\(\mathbf{V}\)が対角行列でない場合のFHモデル、多変量FHモデル。パス]

5. 小地域推定における信頼区間
[パス]

6. 小地域推定におけるそのほかの重要な発展
　本節では個体レベルモデルについて述べる。
　[半分くらい目を通したんだけど、説明が凝縮しすぎててついていけない。これ Rao & Molina (2015) を読み直したほうがよさそう… というわけであえなく断念]

7. その他の小地域推定量
[測定誤差モデル、一般化線形モデル、釣り合い損失関数というのが挙げられている。パス]

8. 要約と今後の課題
　[…]
　今後の課題: モデル・ベース推定値の欠点として、大きな地理的地域についての推定値が直接推定値と全然違ってて、後者のほうがあてになる、というようなことが起きるという点が挙げられる。これに対して、大きな地域については推定値が直接推定値と同じになるように推定値を修正するという方法がある(いわゆる「ベンチマーキング・アプローチ」)。もっとも有名なのはレイキングである。すなわち、すべての小地域推定を定数倍して重みづけ合計が直接推定値と一致するように修正することである。しかしこれはアドホックな修正で統計的基盤がない。そうではなくて、制約付きベイズ小地域推定値を使うという方法がある。すごく複雑なモデルになる。Shen & Louis(1998 JRSS:B)をみよ。
　云々。
—————
　そうではないかと思っていたが、やはり、私の能力と根気を上回る論文であった… おとなしく山に帰るよ…