« 覚え書き：スタイン推定量とはなんぞや | メイン | 読了：Ghosh & Rao (1994) 小地域推定レビュー »
2017年11月14日 (火)
Molina, I., Marhuenda, Y. (2015) sae: An R Package for Small Area Estimation. The R Journal, 7(1).
小地域推定のためのRパッケージsaeの解説。仕事の都合で読んだ。いっちゃなんだけど、面白くも何ともない話題だ... (すいません)
有限母集団$U$が$D$個の地域に相互排他的かつ網羅的に分割されている。サイズを$N_1, \ldots, N_D$とする。[ここで身構えたけど、有限母集団であることはこの後の話には顔を出さない。やれやれだぜ]
地域$d$における個人$j$の、ある関心ある変数の測定値の$Y_{dj}$とする。その地域の全員について$\mathbf{y}_d = (Y_{d1}, \ldots, Y_{dN_d})^t$とベクトルで書いちゃうこともある。目標パラメータを$\delta_d = h(\mathbf{y}_d)$とする(たとえば平均とか)。
各地域についてサイズ$n_d$の下位標本$s_d$が手に入っている。残りの部分を$r_d$とする。
saeパッケージが提供する小地域推定手法は以下の通り。
1. Fay-Herriotモデルに基づくEBLUP。
EPLUPってのはあれね、経験最良線形不偏予測量のことね。
Fay-Herriotモデルってのは79年にアメリカの所得推定のために提案されたモデル。$\delta_d$の直接推定量を$\hat{\delta}_d^{DIR}$とすると、それは不偏であって、
$\hat{\delta}_d^{DIR} = \delta_d + e_d, \ \ e_d \sim N(0, \psi_d)$
と書ける(これを抽出モデルという)。$\psi_d$は既知とする。いっぽう、地域レベルの補足変数ベクトルを$\mathbf{x}_d$として
$\delta_d = \mathbf{x}_d^t \mathbf{\beta} + u_d, \ \ u_d \sim N(0, A)$
とする(これをリンキング・モデルという)。$\beta$は全地域共通。$A$は当面既知とする。
これを併せると
$\hat{\delta}_d^{DIR} = \mathbf{x}_d^t \mathbf{\beta} + u_d + e_d$
という混合モデルになる。これをFHモデルという。[←意外にあたりまえな話でびっくり...これにわざわざ人名をつけるかね...]
さて、このモデル$\delta_d$のBLUP
$\tilde{\delta}_d^{BLUP} = \mathbf{x}_d^t \tilde{\beta}(A) + \tilde{u}_d(A)$
はすでに75年に得られていて、それはこれこれこうである[メモ省略]。問題はそれが$A$の関数だということ。じゃあ$A$はどうすんだと。しょうがない、$\hat{A}$をデータから推定しましょう、というのが、経験BLUP, 略してEBLUPである。
EBLUPは結局こういう形になる。
$\hat{\delta}_d^{EBLUP} = \hat{\gamma} \hat{\delta}_d^{DIR} + (1-\hat{\gamma}_d) \mathbf{x}_d^t \hat{\beta}$
ただし$\hat{\beta} = \tilde{\beta}(\hat{A})$。さて、謎の$\hat{\gamma}$というのが出て、よく見ると2つの項に重みをつけているんだけど、これは[...メモ省略...] $\hat{\gamma} = \hat{A}/(\hat{A} + \psi_d)$である。つまりこういうことだ。直接推定量$\hat{\delta}_d^{DIR}$が十分あてになる場合(=$\psi_d$が小さい場合)、EBLUPは直接推定量に近づく。そうでないとき、EBLUPは回帰による推定量に近づく。なるほど、うまいことできてんね。
$\hat{A}$の一致推定量はいろいろ提案されていて...[メモ略]
モデル比較にはAIC, BICが使えて...[メモ略]
というわけで、saeパッケージはEBLUPを求める関数 eblupFH() とそのMSEを解析的に求める関数 mseFH()をご提供しております。ただし、$\psi_d$は教えてやらないといけないので、事前になんらか推定しておくこと。surveyパッケージでもなんでも好きなパッケージでやるがよろしい。saeにもdirect()という関数がある。
計算例...[略]
[モデルとしてはごくふつうの階層回帰モデルなので、地域数が少なければMplusでどうにかなっちゃいそう。どのくらいスケールするのかが知りたいところだ...]
2. 空間Fay-Herriotモデルに基づくEBLUP。
せっかく空間の話してんのに、なぜただの階層モデルなんだよ、と御不満のみなさん。お待たせしました。地域効果$\mathbf{u} = (u_1, \ldots, u_D)^t$にSAR(1)をいれます。イエー。そう来なくっちゃー。
$\mathbf{u} = \rho_1 \mathbf{Wu} + \mathbf{\epsilon}$
$\mathbf{\epsilon} \sim N(\mathbf{0}_d, \sigma_1^2\mathbf{I}_D)$
$\mathbf{0}$は0の列ベクトル、$\mathbf{I}$は単位行列、添え字はサイズ。隣接行列$\mathbf{W}$は対角が0で隣接に1がはいっている行列を行ごとに標準化したものだと思いねえ。
このモデルのEBLUPも、その解析的なMSEも得られている。saeパッケージではeblupSFH(), mseSFH()をご提供している。ブートストラップMSEとして、パラメトリック版 pbmseSFH(), ノンパラ版 npbmseSFH()もご用意している。[どっち使えばいいんだろう？ 何も書いてない]
計算例... [パス]
3. 時空間Fay-Herriotモデルに基づくEBLUP。
[いま関心ないのでパス]
4. BHFモデルに基づくEBLUP。
今度は補足変数が単位レベルで手に入っている場合について考える。地域$d$の単位$j$について$(Y_{dj}, \mathbf{x}_{dj}^t)$が得られているわけだ。
$Y_{dj} = \mathbf{x}_{dj}^t \mathbf{\beta} + u_d + e_{dj}$
$u_d \sim N(0, \sigma_u^2)$ (iid)
$e_{dj} \sim N(0, \sigma_e^2)$ (iid)
とする。Bettesesさんという人たちの88年の論文の著者頭文字をとってBHFモデルと呼ぶ。
このモデルのEBLUPも得られていて...[メモ省略]。eblupBHF(), pbmseBHF()をご用意している。
計算例... [パスしちゃったけど、使うときには読んだほうがよさそう]
5. BHFモデルに基づく非線形パラメータの経験最良推定量。
[いよいよ大ボスの登場だ...]
$\delta_d = h(\mathbf{y}_d)$が非線形である場合について考える。以下、$\mathbf{y}_d$を並び替え、$(\mathbf{y}_{ds}^t, \mathbf{y}_{dr}^t)^t$としておく[標本を前に持ってくるってことね]。
もはやBLUPではなく、線形性も不偏性も捨てて、MSEを最小化したい。すなわち、手元の$y_{ds}$に条件づけられた、$y_d$のパラメータ$h(y_d)$の、$y_{dr}$の分布を通じた期待値
$\tilde{\delta}_d^{B} = E_{y_{dr}}[h(y_d)|y_{ds}]$
が欲しいわけだ。
こういう場合は、まず$\beta, \theta=(\sigma_u^2, \sigma_e^2)^t$の一致推定量$\hat{\beta}, \hat{\theta}$をMLなりREMLなりで推定する。次にモンデカルロ法で期待値を近似し、$\delta_d$を推定する。これを経験最良(EB)推定量と呼ぶ。
具体的な手順は次の通り。まず$\hat{\beta}, \hat{\theta}$を推定する。次に、$\hat{\beta}, \hat{\theta}$の下での$h(y_{dr}|y_{ds})$を求める(BHFモデルの場合で言えば正規分布である)。で、そこからベクトル$y_{dr}$をたくさん生成する。それぞれの頭に手元の$y_{ds}$をくっつけて、たくさんの架空のセンサス・ベクトル$y_d$を得る。それぞれの$y_d$から$h(y_d)$を得る。これを平均し、$\hat{\delta}_d^{EB}$とする。
ここで大変なのは、長ーいベクトル$y_{dr}$をたくさん生成するというところ。実際にはもうちょっと楽な方法がある。BHFモデルより、$j \in r_d$について
$Y_{dj} = \mathbf{x}_{dj}^t \hat{\beta} + \hat{u} + v_d + \epsilon_{dj}$
$v_d \sim N(0, \hat{\sigma}_u^2(1-\hat{\gamma}))$
$\epsilon_{dj} \sim N(0, \hat{\sigma}_e^2)$
なので、この単変量を生成すれば良い。
さらにこういう手もある。たとえば、いま各地域の貧困生起率$\delta_d$に関心があるとしよう。$O_{dj}$を年収、貧困線を$z$として
$\delta_d = \frac{1}{N_d} \sum^{N_d} I(O_{dj} \lt z)$
と定義しよう。$O_{dj}$はどうみても歪んでいるのでBFHモデルは仮定できないが、たとえば$Y_{dj} = \log(O_{dj}+c)$というようにして正規近似できるかもしれない。$O_{dj}=T^{-1}(Y_{dj})$とすれば
$\delta_d = \frac{1}{N_d} \sum^{N_d} I(T^{-1}(Y_{dj}) \lt z)$
というわけで、一般化していえば$O_{dj}$のBox-Cox変換で正規性を得ることでどうにかなるかもしれない。
saeパッケージではebBHF(), pbmseebBHF()を御用意しており、Box-Cox変換とかも指定できます。...[などなど。途中で疲れて読み飛ばした]
計算例[パス]。
最後に、小地域推定のためのそのほかのRパッケージをご紹介。
- rsaeパッケージ。ロバスト・フィッティングの関数を提供。[2014年で止まっている]
- JoSaeパッケージ。分散行列がブロック対角な単位レベルモデルのみ。[2015年で止まっている]
- hbsaeパッケージ。REMLと階層ベイズ推定を提供。[2012年で止まっている]
- mmeパッケージ。多項線形混合モデル。[2014年で止まっている]
- saeryパッケージ。時間効果のはいった地域レベルのEBLUP。[2014年で止まっている]
- sae2パッケージ。時系列地域レベルモデル。[2015年で止まっている]
[ちなみにこのsaeパッケージは2015年で止まっている。CRAN Task ViewのOfficial Statistics & Survey Methodologyには小地域推定というセクションがあって、載っているのはrsae, hbsae, josae, そしてnlmeとlme4。このパッケージ載ってないじゃん....]
論文：データ解析(2015-) - 読了：Molina & Marhuenda (2015) 小地域推定のためのRパッケージsae