« 読了：Molina & Marhuenda (2015) 小地域推定のためのRパッケージsae | メイン | 読了：Frank, et al. (2017) ベイジアン自白剤 in クラウド・ソーシング »
2017年11月16日 (木)
Ghosh, M., Rao, J.N.K. (1994) Small Area Estimation: An Appraisal. Statistical Science, 9(1), 55-93.
タイトルの通り、小地域推定についてのレビュー。仕事の都合で泣きながら徹夜で読んだ。辛い。というかさ、小地域の推定なんて諦めようよ... 標本が小さいんだから...
1. イントロダクション
小地域推定の歴史とか、活用場面の紹介とか。全部すっ飛ばして...
2. 人口学的方法
人口統計学で用いられてきた、一番古典的な方法。行政からもらってきた、出生とか死亡とか新築とか入学とかのデータ(これをsymptomtic変数という)を、センサスのデータと併せて使う。
そのひとつがVital Rates(VR)法。直近のセンサス年$t=0$について、小地域の粗出生率を$r_{10}$, 粗死亡率を$r_{20}$、その上の親地域の粗出生率を$R_{10}$, 粗死亡率を$R_{20}$とする。年$t$の小地域の粗出生率$r_{1t}$と粗死亡率$r_{2t}$をそれぞれ次とする：
$r_{1t} = r_{10}(R_{1t}/R_{10})$
$r_{2t} = r_{20}(R_{2t}/R_{20})$
[えーと、つまり今年の出生率なり死亡率なりのセンサス年に対する比が小地域と親地域で同じだとするわけね]
すると、年$t$の小地域について粗出生数$b_t$、粗死亡数$d_t$がわかっているとき、人口を次のように推定できる：
$P_t = (1/2) (b_t / r_{1t} + d_t / r_{2t})$
[ああそうか、粗出生率とか粗死亡率ってのは出生数とか死亡数を人口で割ったものだからね]
この手法の弱点は、いうまでもなく、$r_{1t}/r_{10}$とかを$R_{1t}/R_{10}$とかで近似しちゃうところ。
これを拡張した手法として、composite法とか、米センサス局のCM-II法とか、Administrative Records法とかHousing Unit法とかがあって...[なんとなく雰囲気がわかったので省略]
こういう手法はすべて重回帰の特殊ケースとみることができる。そのほかに、symptomatic変数を独立変数にとった重回帰という手もある。
そのひとつがratio-correlation法。前々回センサス年を$\alpha=0$, 前回センサス年を$\alpha=1$、今年を$\alpha=t(>1)$としよう。小地域$i$の年$\alpha$における人口を$P_{i \alpha}$, $j$番目のsymptomatic変数を$S_{ij\alpha}$とする。
まず全部割合にする[うわ...気色悪い...]。
$p_{i\alpha} = P_{i \alpha} / \sum_i P_{i \alpha}$
$s_{ij\alpha} = S_{ij \alpha} / \sum_i S_{ij \alpha}$
センサス年との比、センサス年同士の比をとる。
$R'_{i} = p_{i1} / p_{i0}$
$R_{i} = p_{it} / p_{i1}$
$r'_{ij} = s_{ij1} / s_{ij0}$
$r_{ij} = s_{ijt} / s_{ij1}$
まず、センサス年の間の人口の伸びを目的変数にして回帰する。[めんどくさいので略記するけど$\beta$にはみんなハットがついている]
$R'_{i} = \beta'_0 + \beta'_1 r'_{i1} + \ldots + \beta'_j r'_{ij}$
この係数でセンサス年との比を当てに行く。
$\tilde{R}_{i} = \beta'_0 + \beta'_1 r_{i1} + \ldots + \beta'_j r_{ij}$
で、今年の全地域人口(これは別途わかっているとしよう)に$p_{i1}$を掛け、さらに$\tilde{R}_{i}$もかけて、小地域人口の推定値とする。
この手法の弱点は、センサス年の間で作った$\beta'_j$が、今年と直近センサス年の間でも変わらないという仮定である。
そこで出てきたのが標本回帰法。
$m$個の小地域のうち$k$個については、センサス年に対する比$R_i$の標本推定値$\hat{R}_1, \cdots, \hat{R}_k$があるとしよう。[←なるほどね、ありそうな話だ。小地域のなかにもでかい奴あるしね]
これを説明する回帰モデルを組む。[例によって$\beta$の上のハットは略記]
$\hat{R}_i = \beta_0 + \beta_1 r_{i1} + \ldots + \beta_p r_{ip}$
このモデルを全ての小地域にあてはめる。
3. 合成推定量(synthetic estimators)と関連する推定量
要するに、大地域で得た推定量を使っちゃえという方法。この路線も歴史が古い。
人口が大きな領域$g$に分割されている[えーと、性別みたいなもんですかね]。で、各領域の合計$Y_{.g}$については信頼できる推定量$\hat{Y}'_{.g}$がある。で、それとは別に人口が小地域$i$に分かれてて、セル$(i,g)$の合計$Y_{ij}$を小地域を通じて足しあげると$Y_{.g}$になる。で、なんらかの補足情報$X_{ig}$も手に入るとしよう。
いま、$Y_{i} = \sum_g Y_{ig}$を推定したい。そこで
$\hat{Y}^S_i = \sum_g (X_{ig}/X_{.g}) \hat{Y}'_{.g}$
とすればいいんじゃないかしら。というのが最初期の合成推定量である。
[以下、推定量の分散とかバイアスとかの話が続くので端折って...]
Purcell & Kish (1980)とかがSPREE (構造保存推定)というのを提案していて、これは多元クロス表における反復比例フィッティングを使っているのだけれど、このSPREEは合成推定量を一般化したものだと捉えることができる。[←うおおお。SPREEってKishの論文に出てきて意味不明だったんだけど、そういう意味だったのか...]
[ここから関心ある箇所なので細かくメモ]
合成推定量には潜在的バイアスがあり、いっぽう直接推定量は不安定なので、重みづけ平均を取ろうという路線もある。
つまりこうだ。小地域$i$の直接推定量を$\hat{Y}_{1i}$, なんらかの間接推定量を$\hat{Y}_{2i}$として、適切な重み $w_i$ (0以上1以下)を用意して
$\hat{Y}^C = w_i \hat{Y}_{1i} (1-w_i) \hat{Y}_{2i}$
こういうアプローチは多いんだけど、ここでは直接推定量が不偏推定量$\hat{Y}_i$、間接推定量が合成推定量$\hat{Y}^S_i$であるときの最適な重み$w_i$(opt)の決め方について考える。
$w_i$(opt)を$\hat{Y}^C$のMSEを最小化される重みと捉えると、ふたつの推定量の共分散が0だとして、
$w_i$ (opt) $= MSE(\hat{Y}^C) / [MSE(\hat{Y}^C)+V(\hat{Y}_i)]$
となる。これは推定できるんだけどすごく不安定である。
Purcell & Kish (1979)は、MSEの平均を最小化する共通の$w$を使うというのを提案している。すると次の形になる：
$w_i$ (opt) $= 1 - \sum v(\hat{Y}_i) / \sum (\hat{Y}_i^S - \hat{Y}_i)^2$
もし$\hat{Y}_i$の分散がだいたい等しかったら、$v(\hat{Y}_i)$を$\bar{v} = \sum v(\hat{Y}_i)$で置き換えて、
$w_i$ (opt) $= 1 - m\bar{v} / \sum (\hat{Y}_i^S - \hat{Y}_i)^2$
というわけで、James-Steinタイプのウェイトになる。[←ちょ、ちょっと待って？！ これがなぜJames-Steinタイプなの？ というかタイプってなに？ どこまでがタイプなの？]
もっとも、個々の分散$V(\hat{Y}_i)$が大きく変動するとき、共通の$w$を決めるのが難しくなる。また、プールされる小地域のなかに、他と似ていない小地域があるとき、Jame-Stein推定量は非効率になる。
もっと単純に$w_i$を決めようという路線もある。２つ紹介しよう。以下、地域の人口$N_i$の直接的な不偏推定を$\hat{N}_i$とする。
ひとつめ、Drew, Singh, & Choudhry (1982)。もし$\hat{N}_i \geq \delta N_i$だったら$w_i=1$、そうでなかったら$w_i=\hat{N}_i / \delta N_i$とする。$\delta$は適当に決める。
ふたつめ、Sarndal & Hidiroglou (1989)。もし$\hat{N}_i \geq N_i$だったら$w_i=1$、そうでなかったら$w_i=(\hat{N}_i / N_i)^{h-1}$とする。$h$は適当に決めていいんだけど、著者らは$h=2$を推奨。
なお、$\delta=1, h=2$とするとこの2つは等しくなる。
例として、人口$N$から$n$を単純無作為抽出する場面を考えよう。$\hat{N}_i = N(n_i/n)$である。
もし、$n_i$がその期待値$E(n_i)=n(N_i/N)$と同じくらい大きければ、($\delta=1$とすれば)どちらのやりかたでも$w_i=1$となる。つまり、他の小地域から値を借りてくることができなくなる。$E(n_i)$がいくら小さくてもそうなるわけである。
もし、$\hat{N}_i \lt N_i$だったら、($h=2$として)どちらのやりかたでも$w_i$は$n_i$の現象と共に減少する。合成推定量の重みがどんどん増すことになる。
もうひとつ欠点を挙げると、こういう方法は小地域以内の変動に対する小地域間の変動の相対的サイズを無視しているわけで、小地域間の等質性がどうであれ、重みは同じになる。
[うぐぐ... 途中から混乱しちゃって写経状態になってしまった。あとでゆっくり考えよう]
Holt, Smith & Tomberlin(1979)はこういうのを考えた。セル$(i,g)$の単位$l$について、
$y_{igl} = \mu_g + e_{igl}$
とする。大地域$g$の固定効果$\mu_g$のBLUP(最良線形不偏推定量)は$\hat{\mu}_g = \bar{y}_{.g}$である。ここから$Y_i$のBLUPは下式となる：
$\hat{Y}_i^B = \sum_g \hat{Y}_{ig}^C$
ただし$\hat{Y}_{ig}^C$は、直接推定量$\hat{Y}_{ig} = N_{ig}\bar{y}_{ig}$に重み$w_{ig} = n_{ig}/N_{ig}$、合成推定量$\hat{Y}_{ig}^S=N_{ig}\bar{y}_{ig}$に重み$1-w_{ig}$を与えて足し上げた推定量である。この方法もまた小地域間変動を無視しているわけだけど、モデルに小地域のランダム効果をいれることで改善できる。これが次節以降の手法のもとになっている。[なるほど...]
4. 小地域モデル
ここからは小地域のランダム効果をいれたモデル。大きく２つにわかれる。
ひとつめ。小地域の補足データ$\mathbf{x}_i = (x_{i1}, \ldots, x_{ip})^t$をいれる。関心あるパラメータを$\theta_i$として
$\theta_i = \mathbf{x}_i^t \beta + v_i z_i$
$z_i$は既知の正定数。$v_i$は$E(v_i)=0, V(v_1)=\sigma^2_v$、なんなら正規性を仮定してもよい。
ここで、直接推定量$\hat{\theta}$について
$\hat{\theta} = \theta_i + e_i$
かつ$E(e_i | \theta_i) = 0$(つまりデザイン不偏)、$V(e_i | \theta_i) = \psi_i$という風に仮定し、$\psi_i$は既知とみることが多い。実際には$\theta_i$が非線形な関数になっていたりして、不偏という仮定には無理があったりするんだけど。まあとにかく、
$\hat{\theta}_i = \mathbf{x}_i^t \beta + v_i z_i + e_i$
ということになる。線形混合モデルの特殊ケースである。
[あとで出てくるけど、こういうのをFay-Herriotタイプのモデルという由]
ふたつめ。測定単位の補足データ$\mathbf{x}_{ij} = (x_{ij1}, \ldots, x_{ijp})^t$をいれる。
$y_{ij} = \mathbf{x}_{ij}^t \beta + v_i + e_{ij}$
ここで$e_{ij} = \tilde{e}_{ij} k_{ij}$, $k_{ij}$は既知の定数で、$\tilde{e}_{ij}$は平均0, 分散$\sigma^2$、なんなら$v_i$と$\tilde{e}_{ij}$に正規性を仮定してもよい。こういうのをnested error 回帰モデルという。関心があるのは合計$Y_i$、ないし平均$\bar{Y}_i = Y_i / N_i$である。
このモデルは選択バイアスを考慮してないことに注意。したがって標本が単純無作為抽出でないとそのままでは使えない。
このモデルはいろんなところで使われていて[...中略...]、最後に触れるけど、拡張の提案も沢山ある。
5.1 EBLUP(分散成分)アプローチ
BLUPというのは、固定パラメータのBLUEみたいなもので、正規性を仮定することなく、線形不偏推定量のなかでMSE最小の奴を探すアプローチ。50年代に遡る。
さっきの$\theta_i$のほうのモデルでいうと、そのBLUPは...[中略] そのMSEは... [中略]。ところがこれらは分散成分$\sigma^2_v$を既知ととっている。実際にはふつう未知なわけで、これをまずMLとかREMLとかで推定する、というのがEBLUP。[... 中略... ]
では、nested error 回帰モデルのほうはどうなるか ... [...中略...]
5.2 EBアプローチ
まずデータの周辺分布からモデルのパラメータを推定し、それを使って関心あるパラメータの事後分布を得る、というアプローチ。
[...時間が無くなってきたのでメモ省略...]
5.3 HBアプローチ
モデルのパラメータの事前分布を決め、関心あるパラメータの事後分布を得る、というアプローチ。[... 心底面倒くさいので読み飛ばした...]
Datta & Ghosh (1991) はHB, EB, EBLUPを比較している。[...中略。要するに、点推定についていえば大差ない]
6. 事例
[略]
7.1 モデル診断
[略。適合度の話ではなくて、q-qプロット描きましょう的な話]
7.2 制約付き推定
[小地域の推定値の合計がなにかの既知の値にぴったり合わないと困る、というような場面の話。あるある、そういうの。時間がないので飛ばしたが、結構切実な話だ。いずれ読み直そう]
7.3 拡張
Fay-Herriotタイプのモデル：
- 標本抽出誤差に相関を入れようという話。
- 繰り返し測定をやっている場合の話。
nested error回帰モデル：
- 多変量化するという話。[あああ... そうか... これは結構切実だ... Fuller & Harter (1987)]
- 小地域のなかで2段階抽出しているときのモデル。
- $\mathbf{x}_{ij} = 1$のときのモデル。[？ それってどう特別なのだろうか。Kleffe &: Rao(1992)]
- 反応が二値の場合のモデル。ロジスティック回帰モデルを使う。
- ポワソン回帰を使うモデル。死亡力とか疾病割合とかに使う。
- 二変量のモデル。2つのガンの同時死亡力について使う。
8. 結論
最後にひとこと。小地域の間接推定量を使うときには十分に気をつけるように。それはもともと、直接推定量の代替に過ぎないんだからね。
やれやれ、やっと終わった... 疲れた...
論文：データ解析(2015-) - 読了：Ghosh & Rao (1994) 小地域推定レビュー