この文書は、 Rao & Molina (2015) “Small Area Estimation”, 2nd Ed. を読みながらとったノートである。 いうまでもなく、私による私のための私のノートであり、すべての誤りは私に帰属いたします。
なお、このノートをとりつつ実際にRで推定を試してみた。 その記録は、実習編として別途まとめてある。
本文書はR Notebooksによって 作成しているが、実習編を分離してしまったので、もはや Rのコードは全く含まれていない。 なんだか妙な感じですね。
更新履歴:
- 2018/06/20 作成開始
- 2018/06/22 1章, 4章, 5章, 6章を読了。疲れた。
- 2018/06/23 6章の実習(生乳支出データのFay-Harrietモデル)。
- 2018/06/24 飽きたので6章の実習の続き(hbsaeパッケージでも試してみた)。
- 2018/06/25 量が多すぎて読みにくくなってきたので、 実習を実習編に分離。 7.1節, 8.4節, 8.9節, 10.1節, 10.3節を読了。
- 2018/06/26 10.4節, 10.13節を読了。読書ノートは この辺にして、実習に注力します。
- 2018/06/27 誤字などを修正。
- 2018/07/01 ちょっぴり追記。
関心ある集団なり、その下位集団(ドメイン)なりについて推定したいとき、 標本調査が広く使われている。 ここでドメインとは、地理的な地域、デモグラフィック・ソシオグラフィック 属性で定義された集団、などである。
あるドメインについての推定量がそのドメインの標本だけに基づいているとき、 これを直接推定量と呼ぶ。 直接推定量はふつう「デザイン・ベース」である。 つまり、母集団の値は固定とし、標本抽出デザインによって生じる確率分布に基づいて推論する。
なんらかの事情で標本が十分でなく、適切な直接推定値が得られないドメインがあるとき、これを 小地域と呼ぶ。
小地域推定(small area estimation, SAE)においては、他のエリアなり時点なりから値を「借りてくる」 ことが必要になることが多い。これを間接推定量と呼ぶ。
[パス]
伝統的にはデザイン・ベースの間接推定量が用いられてきた。3章で述べる。
補足変数で説明できない地域間変動を説明するために地域ごとのランダム効果を用いるモデルを 小地域モデルという。これに基づく間接推定量をモデル・ベース推定量という。 4章で述べる。
現在では、間接推定量は明示的な小地域モデルに基づいていなければ ならないと考えられている。
本書では以下に焦点を当てる:
[パス]
[パス。構成のみメモしておく]
[パス。構成のみメモしておく]
3章の伝統的な間接推定は、補足データを通じて小地域間のリンクを与えるモデルに暗黙的に 基づいていた。
ここからは、地域間分散を特定し許容する明示的な小地域モデルを扱う。 明示的なモデルを使うことには以下の利点がある。
小地域モデルは以下に大別される。
地域\(i\)(\(i=1,\ldots,m\))における、なにかの変数の平均(たとえば一人当たり所得)を\(\bar{Y}_i\)とする。 これを変換するなんらかの関数\(g(\cdot)\)があって、 \(\theta_i = g(\bar{Y}_i)\)であるとする。
直接推定量\(\hat{\bar{Y}}_i\)は手に入っているとして、 \(\hat{\theta}_i = g(\hat{\bar{Y}}_i)\)とする。
また、地域の補足データ\(\mathbf{z}_i\)(長さ\(p\)のベクトル)も手に入っているとする。
\(\bar{Y}_i\)について推測したい。 そこで以下のように想定する。
まずlinkingモデル。 \[ \theta_i = \mathbf{z}_i^T \mathbf{\beta} + b_i v_i\] ここで\(\beta\)は長さ\(p\)。\(v_i\)が地域のランダム効果、\(b_i\)は既知の正の定数。 で、\(v_i\)は平均\(E_m(v_i)= 0\)、分散\(V_m(v_i) = \sigma_v^2\)の分布に独立に従う ものとする (以後、これを\(v_i \mathop{\sim}^{iid} (0, \sigma_v^2)\)と書く)。 正規性を仮定してもよい。
次にsamplingモデル。 \[ \hat{\theta}_i = \theta_i + e_i , \ \ e_i | \theta_i \sim (0, \psi_i)\] \(e_i\)は標本抽出誤差。 慣用的には\(\psi_i\)を既知とする。\(\hat{\theta}_i\)について正規性を仮定しても よい。
[↑上は\(\psi_i^2\)の誤植ではない。この本を通じて\(\psi_i\)は分散であることに注意]
以上まとめると \[\hat{\theta}_i = \mathbf{z}_i^T \mathbf{\beta} + b_i v_i + e_i\] これは結局、線形混合モデルの特殊ケースである。 これをFay-Herriotモデルという。
結構きつい想定をしていることに注意:
[事例4.2.1, 4.2.2, 4.2.3。略]
地域\(i\)の要素\(j\)(\(=1,\ldots,N_i\))のそれぞれについて、 補足データ\(\mathbf{x}_{ij}\)(長さ\(p\))が手に入っているとしよう。 あるいは、少なくとも母集団平均\(\bar{\mathbf{X}}_i\)は手に入っているとしよう。
次のnested error線形回帰モデルを考える。 \[ y_{ij} = \mathbf{x}_{ij}^T \mathbf{\beta} + v_i + e_{ij}, \ \ v_i \mathop{\sim}^{iid} (0, \sigma_v^2) \] \(v_i\)に正規性を仮定してもよい。さらに、 \[ e_{ij} = k_{ij} \tilde{e}_{ij}, \ \ \tilde{e}_{ij} \mathop{\sim}^{iid} (0, \sigma_e^2) \] \(k_{ij}\)は既知の定数。\(\tilde{e}_{ij}\)に正規性を仮定してもよい。
上のモデルでは\(j=1,\ldots,N_i\)なんだけど、そこから抽出したサイズ\(n_i\)の標本についても 成り立つと考える。この想定は、抽出がSRSの場合、ないし抽出が\(x_{ij}\)に基づいている場合に 成立する。[…パス]
関心があるのは、地域平均\(\bar{Y}_i\)か合計\(Y_i\)である。 もし\(N_i\)が大きければ、\(\bar{Y}_i\)は \[ \mu_i = \bar{\mathbf{X}}_i^T \mathbf{\beta} + v_i\] に近づく。つまり、\(\bar{Y}_i\)の推定は\(\beta\)ならびに\(v_i\)の実現値の推定に近づく。
[事例4.3.1, 4.3.2。略]
地域特性が\(r\)個ある場合について考える。 地域\(i\)の特性\(j\)の地域平均を\(\bar{Y}_{ij}\)とする。\(\theta_{ij}=g_j(\bar{Y}_{ij})\)の ベクトルを\(\mathbf{\theta}_i\)と書き(長さ\(r\))、 調査から得られた推定量のベクトルを\(\hat{\mathbf{\theta}}_i\)と書く。
samplingモデル: \[ \hat{\mathbf{\theta}}_i = \mathbf{\theta}_i + \mathbf{e}_i, \ \ \mathbf{e}_i \sim N_r(\mathbf{0}, \mathbf{\Psi}_i) \] 共分散行列\(\mathbf{\Psi}_i\)は既知とする。
linkingモデル: \[ \mathbf{\theta}_i = \mathbf{Z}_i \mathbf{\beta} + \mathbf{v}_i, \ \ \mathbf{v}_i \mathop{\sim}^{iid} N_r(\mathbf{0}, \mathbf{\Sigma}_v) \] \(\mathbf{Z}_i\)は\(r \times rp\)の行列で, \(j\)番目の行は、長さ\(p\)の横ベクトルを\(r\)本横に並べたもので、その\(j\)番目だけは 補足変数\(z_{ij}^T\)で、ほかのはすべて\(\mathbf{0}^T\)である。\(\mathbf{\beta}\)は長さ\(np\)。
以上まとめると \[ \hat{\mathbf{\theta}}_i = \mathbf{Z}_i \mathbf{\beta} + \mathbf{v}_i + \mathbf{e}_i\] これはFHモデルの自然な拡張になっている。特性が複数あるときは、 個別にFHモデルを使うんじゃなくて、これを使ったほうが良い。
[事例4.4.1。パス]
標本抽出誤差\(e_i\)の間に相関を考える。つまり、\(m\)個の\(\theta_i\)のベクトルを\(\mathbf{\theta}\), \(m\)個の\(\hat{\theta}_i\)のベクトルを\(\hat{\mathbf{\theta}}\)と略記して、 \[ \hat{\mathbf{\theta}} = \mathbf{\theta} + \mathbf{e},\ \ \mathbf{e}|\mathbf{\theta} \sim N_m(\mathbf{0}, \mathbf{\Psi}) \] と考えるのである。共分散行列\(\mathbf{\Psi}\)は既知とする。 実際には、調査による推定量\(\hat{\mathbf{\Psi}}\)、ないし平滑化推定量で代用してしまう。
[事例4.4.2。パス]
[パス]
基本的なFHモデルは、地域効果\(v_i\)がiidだと仮定している。 しかし相関を考えたほうが現実的であることも多い。
その1, CAR(条件つき自己回帰)空間モデル。
地域\(i\)の「近接」地域の集合を\(A_i\)とする。 他のエリア\(\{v_l:l\neq i\}\)の下での\(b_i v_i\)の条件つき分布を \[ b_i v_i | \{v_l:l\neq i\} \sim N(\rho \sum_{l \in A_i} q_{il}b_l v_l, b^2_i\sigma^2_v) \] ここで\(\{q_{il}\}\)は既知の定数で\(q_{il} b^2_l = q_{li} b^2_i\)を満たす。 未知パラメータは\(\rho\)と\(\sigma^2_v\)である。
[ちょ、ちょっとまって… ええと、地域効果\(b_iv_i\)は正規分布に従う。その平均は、近接 地域\(l \in A_i\)の地域効果を、重み\(q_{il}\)をつけて合計して\(\rho\)倍したものだ。 \(v_i\)の分散は\(\sigma^2_v\)だ。ってことね]
ここから次式が得られる。\(b^2_1, \ldots, b^2_m\)を持つ対角行列を\(\mathbf{B}\)とし、 \(m\)個の\(v_i\)のベクトルを\(\mathbf{v}\)とする。\(q_{il}\)を全部突っ込んだ\(m \times m\)行列を \(\mathbf{Q}\)とする(近接じゃないところと対角は0にする)。 \[ \mathbf{B}^{1/2} \mathbf{v} \sim N_m(\mathbf{0}, \mathbf{\Gamma}(\mathbf{\delta})) \] \[ \mathbf{\Gamma}(\mathbf{\delta}) = \sigma_v^2(\mathbf{I}-\rho \mathbf{Q})^{-1} \mathbf{B} \]
[恥ずかしながら全然理解できない。一本目は、ベクトル\((b_1 v_1, \ldots, b_m v_m)^T\)がMVNに従う ということですよね。その平均がなぜ0だといえるの？ さっきまで \(\rho \sum_{l \in A_i} q_{il}b_l v_l\)だったのに。なぜだーーー]
その2, 地球統計学のモデル。
地球統計学の文献では、共分散構造について次の2通りの形式を使っている。地域\(i, l\)間の 「距離」(ユークリッドとは限らない)を\(d_{il}\)とする。
[何の話をしているのかを整理しておこう。ある小地域の地域効果\(v_i\)は、 平均ゼロのMVNに従い、分散は\(\sigma^2_v \delta_1\)であり、地域間の共分散は 地域間の距離\(d_{ij}\)で決まる、と考える。 形式1は、地域間の共分散を\(\sigma_v^2 \delta_2 \exp(-d_{ij})\)としている。 いっぽう形式2は、地域間の共分散を\(\sigma_v^2 \delta_2 \delta_3^{d_{ij}}\)と している。あれ？ これだと距離が広くなるにつれて共分散が大きくなることにならない？ よくわからん…]
その3, SAR(同時自己回帰過程)モデル。
\(b_i\)はみな1とする。 \(\mathbf{v} = (v_1, \ldots, v_m)^T\)についてこうモデル化する。 \[ \mathbf{v} = \phi \mathbf{W} \mathbf{v} + \mathbf{u}, \ \ \mathbf{u} \sim N(\mathbf{0}, \sigma^2_u \mathbf{I}) \] [なるほど、自己回帰なわけだ]
\(\mathbf{W}\)はCARモデルの\(\mathbf{Q}\)みたいなもので、地域間の近接性を定義する行列。 \(\phi\)は空間的関係の強さを表す。 \(\mathbf{I}-\phi \mathbf{W}\)が正則であるというのが唯一の条件。
\(\mathbf{W}\)の簡単な決め方は、隣接している地域に1を立てるというもの。 ないし、それを行方向に標準化してもよい(各行の和を1にする)。後者の場合、 \(\phi \in (-1,1)\)となり、相関係数と解釈できるので、空間自己相関パラメータと呼ぶ。
CARモデルやSARモデルの欠点は、\(A_i\)次第で結果が決まるという点である。 ちょっと主観的だともいえる。
[←なるほどね… この本のなかで 空間モデルの扱いが意外に小さいのは、SAEは主に公的統計の分野で 発展していて、公的統計では主観的要素が入るのをすごく嫌うから、 なのかもしれない]
[事例4.4.6。パス]
[地域がさらに下位地域に分かれていて、下位地域の平均にも関心がある場合の話。パス]
[パス。構成のみメモする]
\(y_{ij}\)が二値で、小地域の割合 \(\bar{Y}_i = P_i = \sum_j^{N_i} \frac{y_{ij}}{N_i}\)に 関心があるとしよう。
MacGibbon & Tomberlin(1989 Survey Methodology)はこう考えた。 \(y_{ij} \sim Bernoulli(p_{ij})\)と考えて、 \[ logit(p_{ij}) = \mathbf{x}_{ij}^T \mathbf{\beta} + v_i, \ \ v_i \mathop{\sim}^{iid} N(0, \sigma_v^2) \] \(x_{ij}\)はユニットレベル共変量。
このモデルを 経験ベイズ法か経験最良法か階層ベイズ法で推定し、 \(\mathbf{\beta}\)と\(v_i\)の実現値を推定する。
で、\(P_i\)の推定量を以下とする: \[ \left( \sum_{j \in s_i} y_{ij} + \sum_{j \in r_i} \hat{p}_{ij} \right) / N_i\] \(\hat{p}_{ij}\)は、\(\mathbf{\beta}\)と\(v_i\)の推定に基づく推定値。
[わ、わからない… 困った。よく理解できない点を3点メモしておく。
うーむ、これは元論文を当たるべきかもしれない…めんどくさい…]
Malec et al.(1997)は、 各地域のユニットが\(h\)個のクラスにグループ化されている場合について…[略]。
[事例4.6.1。パス]
[パス]
[パス]
[パス]
[パス]
4章ではいくつかの小地域モデルを示した。 それらのなかには一般線形混合モデルの特殊ケースとしてみなせるものもあった。 さらに、地域の母集団サイズが大きければ、小地域平均は 固定効果とランダム効果の線形結合とみなすことができた。
パラメータを推定する方法として、最良線形不偏予測(BLUP)推定量がある。 BLUP推定量とは、線形不偏推定量のクラスのなかでMSE最小の推定量である。
BLUP推定量は、ランダム効果の正規性には依存しないが、 ランダム効果の分散(そして共分散)に依存する。これを分散成分という。
分散成分の推定には、モーメント法、ML法、REML法を使う。 分散成分の推定値をつかったBLUP推定量をEBLUP推定量という。
次の一般線形混合モデルを考える。 \[ \mathbf{y} = \mathbf{X} \mathbf{\beta} + \mathbf{Z} \mathbf{v} + \mathbf{e}, \ \ \mathbf{v} \mathop{\sim}^{iid} (0, \mathbf{G}), \ \ \mathbf{e} \mathop{\sim}^{iid} (0, \mathbf{R}) \] \(\mathbf{y}\)は長さ\(n\)の観察値ベクトル。 \(\mathbf{X}\)(\(n \times p\))と\(\mathbf{Z}\)(\(n \times h\))は既知でフルランク。 \(\mathbf{y}\)の共分散行列は\(\mathbf{V}=\mathbf{R}+\mathbf{ZGZ}^T\)である。
\(\mathbf{G}\)と\(\mathbf{R}\)は長さ\(q\)のベクトル\(\mathbf{\delta}\)が決まれば決まるものとする。 \(\mathbf{\delta}\)は\(q\)次元ユークリッド空間の特定の下位空間に属しており、 その下位空間に属しているすべての\(\mathbf{\delta}\)について\(\mathbf{V}\)は正則であるとする。
我々が\(\mu = \mathbf{l}^T \mathbf{\beta} + \mathbf{m}^T \mathbf{v}\)に関心を持っている。 ここで\(\mathbf{l}\)と\(\mathbf{m}\)は定数ベクトル。
[要するに、我々はある共変量ベクトルの下で 観察値がどうなるかを知りたいということっすね]
\(\mu\)の線形推定量について考えよう。つまり、\(\hat{\mu} = \mathbf{a}^T \mathbf{y} + b\)という形になる \(\hat{\mu}\)である。上の一般線形混合モデルに照らして期待値を求めたとき、 \(E(\hat{\mu}) = E(\mu)\)であれば、\(\hat{\mu}\)はモデル不偏であるといえる。
\(\hat{\mu}\)のMSE \(MSE(\hat{\mu}) = E(\hat{\mu}-\mu)^2\)について考える。 人はこれを平均二乗予測誤差(MSPE)とか予測平均二乗誤差(PMSE)などと呼ぶ。これを、 線形で不偏な\(\hat{\mu}\)のクラスの中で最小化したい。これをBLUP推定量という。
既知の\(\delta\)に対し、\(\mu\)のBLUP推定量は次式で与えられる。 \[ \tilde{\mu}^H = t(\mathbf{\delta}, \mathbf{y}) = \mathbf{l}^T \tilde{\mathbf{\beta}} + m^T \tilde{\mathbf{v}} \] \[ \tilde{\mathbf{\beta}} = \tilde{\mathbf{\beta}}(\mathbf{\delta}) = (\mathbf{X}^T \mathbf{V}^{-1} \mathbf{X})^{-1} \mathbf{X}^T \mathbf{V}^{-1}\mathbf{y} \] \[ \tilde{\mathbf{v}} = \tilde{\mathbf{v}}(\mathbf{\delta}) = \mathbf{GZ}^T \mathbf{V}^{-1} (\mathbf{y}-\mathbf{X}\tilde{\mathbf{\beta}}) \] \(\tilde{\beta}\)は\(\beta\)の最良線形不偏推定量(BLUE)である。証明は5.6.1節で与える。
[別の証明について。パス]
[別の定式化について。BLUP推定値が「同時最尤推定値」とも呼ばれている理由の 説明。パス]
[正規性を仮定すれば\(\mu\)の最良予測推定量はBLUPになるとかなんとか。パス]
[\(\mathbf{y}\)を変換してどうのこうの… これはREMLと関連していて… 云々。パス]
[\(\mu\)を多次元に拡張すると… パス]
[パス]
BLUP推定量\(t(\mathbf{\delta}, \mathbf{y})\)は分散パラメータ\(\mathbf{\delta}\)に依存しているが、 \(\mathbf{\delta}\)は実際には未知である。
\(\mathbf{\delta}\)を推定量\(\hat{\mathbf{\delta}} = \hat{\mathbf{\delta}}(\mathbf{y})\)に置き換えると、 二段階推定量\(\hat{\mu}^H = t(\hat{\mathbf{\delta}}, \mathbf{y})\)が得られる。 これは経験BLUP(EBLUP)と呼ばれている。
ここからは、\(t(\hat{\delta}, \mathbf{y})\), \(t(\delta, \mathbf{y})\)を \(t(\hat{\delta})\), \(t(\delta)\)と略記する。
次の条件を満たすとき、\(t(\hat{\mathbf{\delta}})\)が不偏である、つまり\(E[t(\hat{\mathbf{\delta}})-\mu]=0\) となる。
[その証明について。省略]
[Kackar & Harville(1981)の紹介。残念ながら主旨が理解できない。パス]
それでは\(\mathbf{\beta}\)と\(\mathbf{\delta}\)をどうやって推定するのか。
まずはML推定量について。
正規性を仮定すると、対数尤度関数は[…パス…]
ここから\(\mathbf{\delta}\)のML推定量は[…パス…]
\(\mathbf{\beta}\)のML推定量は[…パス…]。
この方法の欠点は、\(\mathbf{\delta}\)のML推定量が、\(\mathbf{\beta}\)を 推定したことで生じている自由度の損失を 無視しているという点だ[…パス…]。 そこでREML推定量というのがあって[…パス、パス、みーんなパスだ]。
他に、正規性を仮定しないMINQU推定量というのもあって…[とかなんとかいろいろ書いてある。 全部パスだ！やってられっか！]
[パス]
[パス]
[あとで気がついたんだけど、この節はすごく重要]
一般線形混合モデル \[ \mathbf{y} = \mathbf{X} \mathbf{\beta} + \mathbf{Z} \mathbf{v} + \mathbf{e}\] に出てくるベクトルと行列が、\(m\)個の要素(典型的には小地域)に分かれていると しよう。\(\mathbf{y}\), \(\mathbf{X}\), \(\mathbf{\beta}\), \(\mathbf{Z}\), \(\mathbf{v}\), \(\mathbf{e}\)はそれそれの小地域のベクトルなり行列なりを縦に積んだものになる。 共分散行列\(\mathbf{R}\), \(\mathbf{G}\), \(\mathbf{V}\)はブロック対角行列になる。 \(i\)について \[ \mathbf{y}_i = \mathbf{X}_i \mathbf{\beta} + \mathbf{Z}_i \mathbf{v}_i + \mathbf{e}_i\] が成り立つ。
我々は\(\mu_i = \mathbf{l}_i^T \mathbf{\beta} + \mathbf{m}_i^T \mathbf{v}_i\)に関心がある。 BLUP推定量は、さきほどと同じく、 \[ \tilde{\mu}_i^H = t_i(\mathbf{\delta}, \mathbf{y}) = \mathbf{l}_i^T \tilde{\mathbf{\beta}} + \mathbf{m}_i^T \tilde{\mathbf{v}_i} \] \[ \tilde{\mathbf{\beta}} = \tilde{\mathbf{\beta}}(\mathbf{\delta}) = \left( \sum_i^m \mathbf{X}_i^T \mathbf{V}_i^{-1} \mathbf{X}_i \right)^{-1} \sum_i^m \mathbf{X}_i^T \mathbf{V}_i^{-1} \mathbf{y}_i \]
\[ \tilde{\mathbf{v}} = \tilde{\mathbf{v}}_i(\mathbf{\delta}) = \mathbf{G}_i \mathbf{Z}_i^T \mathbf{V}_i^{-1} (\mathbf{y}_i-\mathbf{X}_i \tilde{\mathbf{\beta}}) \] \(\tilde{\mu}^H_i\)のMSEは[…略]。
EBLUP推定量は、\(\mathbf{\delta}\)を\(\hat{\mathbf{\delta}}\)に置き換えて、 \[ \hat{\mu}_i^H = t_i( \hat{\mathbf{\delta}}, \mathbf{y} ) = l_i^T \tilde{\mathbf{\beta}} (\hat{\mathbf{\delta}}) + m_i^T \tilde{\mathbf{v}_i}(\hat{\mathbf{\delta}}) \] となる。
[パス]
[パス]
[パス。構成のみメモする]
[ひとしきりSASについての説明。パス]
Rでは、nlmeパッケージのlme関数で、ブロック対角共分散行列を持つ 線形混合モデルを推定できる。ML法、REML法が使える。[…略]
lme4パッケージは、lmer関数, glmer関数, nlmer関数で、 線形混合モデル、一般化線形混合モデル、非線形混合モデルを推定できる。 nmleパッケージよりメモリが効率的だといわれている。
小地域推定に特化したパッケージとしてsaeがある。[…略]
[パス。構成のみメモする]
復習しよう。
Fay-Herriotモデル(4.2節)は \[\hat{\theta}_i = \mathbf{z}_i^T \mathbf{\beta} + b_i v_i + e_i\]
ブロック対角共分散構造をもつ一般線形混合モデル(5.3節)はこうだった。
\[ \mathbf{y}_i = \mathbf{X}_i \mathbf{\beta} + \mathbf{Z}_i \mathbf{v}_i + \mathbf{e}_i\]
前者は後者の特殊ケースである。
[頭を整理するために、以下付記すると…
ということであろう。]
後者において、\(\mu_i\)のBLUP推定量は \[ \tilde{\mu}_i^H = \mathbf{l}_i^T \tilde{\mathbf{\beta}} + \mathbf{m}_i^T \tilde{\mathbf{v}}_i \] であった。ここから、前者における\(\theta_i\)のBLUP推定量は下式になることがわかる。 \[ \tilde{\theta}^H_i = \mathbf{z}^T \tilde{\mathbf{\beta}} + \gamma_i(\hat{\theta}_i - \mathbf{z}_i^T \tilde{\mathbf{\beta}}) = \gamma_i \hat{\theta}_i + (1-\gamma) \mathbf{z}^T \tilde{\mathbf{\beta}} \] \[ \gamma_i = \sigma^2_v b^2_i / (\psi_i + \sigma_v^2 b_i^2) \] [なぜこうなるのかよくわからない… あまりに自明すぎて説明してないのだろうか…悲しい…]
ここで\(\tilde{\mathbf{\beta}}\)は\(\mathbf{\beta}\)のBLUEであり、このモデルの場合は \[ \tilde{\mathbf{\beta}} = \left[ \sum_i \frac{\mathbf{z}_i \mathbf{z}^T_i}{\psi + \sigma^2_v b^2_i} \right]^{-1} \left[ \sum_i \frac{\mathbf{z}_i \hat{\theta}_i}{\psi + \sigma^2_v b^2_i} \right] \]
ここから次のことがわかる。BLUP推定量\(\hat{\theta}_i^H\)は、直接推定量 \(\tilde{\theta}_i\)と、回帰synthetic推定量\(\mathbf{z}_i^T \tilde{\beta}\)の 重みづけ平均である。重み\(\gamma_i\) (0以上1以下)は、 モデル分散\(\sigma_v^2 b_i^2\)と全分散\(\psi_i +\sigma_v^2 b_i^2\) の比で、\(\theta_i\)のモデルの不確実性を表している。モデル分散が小さければ synthetic推定量に重みがかかるし、デザイン分散\(\psi_i\)が小さければ直接推定量に 重みがかかる。
\(\hat{\theta}_i^H\)は標本抽出デザインがどうであれ妥当である[…省略]
\(\hat{\theta}_i^H\)のMSEは…[パス]
さて、分散成分\(\sigma_v^2\)は実際には未知なので、\(\hat{\sigma}_v^2\)で 置き換える[式省略]。これがEBLUP推定量である。
ここで注意。\(\theta_i\)のEBLUP推定量\(\hat{\theta}_i^H\)が手に入ったからといって、 これを小地域平均に戻すために\(g^{-1}(\hat{\theta}_i^H)\)とすると、 それはもはやEBLUPではない。だから、\(g(\cdot)\)が非線形な 場合には、EBLUPじゃなくてEBかHBがよい。
Fay & Herriot(1979)はcompromise EBLUPというのを提案していて[…略]
標本がない小地域がある場合は、直接推定量が手に入らないから、\(\mathbf{z}_l^T \mathbf{\beta}\) を使うしかなくて[…略]
方法その1、モーメント推定量\(\hat{\sigma}^2_{vm}\)。 いま、 \[ a(\sigma^2_v) = \sum_i^m \frac{ (\hat{\theta}_i - \mathbf{Z}^T_i \tilde{\mathbf{\beta}})^2 }{ (\psi_i + \sigma^2_v b^2_i) } \] と定義する。その期待値は\(E[a(\sigma^2_v)]=m-p\)である。[←なんでええええ？ わかんないいいいい]
ということは、\(a(\sigma^2_v)=m-p\)を反復して解けばよい。 Fay & Herriot(1979)いわく、その方法としては[…略]。 たいてい10反復くらいで収束する。
方法その2、単純なモーメント推定量\(\hat{\sigma}^2_{vs}\)。[略]
方法その3、ML推定量\(\hat{\sigma}^2_{vML}\)。[略]
方法その4、REML推定量\(\hat{\sigma}^2_{vRE}\)。[略]
どの方法であっても、\(v_i\)と\(e_i\)が0のまわりで対称である限り、EBLUP推定量 \(\hat{\theta}_i^H\)はモデル不偏である。
[\(b_i=1\)で\(\psi_i\)がすべて等しいとき… 略]
\(\sigma^2_v\)の4種類の推定量の、\(m \rightarrow \infty\)の漸近的分散を 比べると、[…略…]、REMLとMLが最小、次がモーメント推定量(方法1)、最大が方法2である。 方法2は特に大きい。
[パス]
[パス. 構成のみメモする]
[パス]
[パス]
[パス]
[パス]
「ビッグ・データ」は統計学のホットトピックだ。 小地域推定においてビッグデータを地域レベルの共変量として用いた例を紹介する。
その1, Marchetti et al.(2015)。イタリア・トスカーニ地方の地域の 貧困率の推定。直接推定量はサーベイで得たのだが、共変量をビッグデータから得た。
車のモビリティに注目し、GPSのデータを集めた。で、ある車\(v\)のモビリティを 次のように定義した。
\[ M_v = - \sum_{l_1}^{L} \sum_{l_2}^{L} p_v(l_1, l_2) \log(p_v(l_1, l_2)) \]
\(L\)は地点数, \((l_1, l_2)\)はそのなかの2地点。\(p_v(l_1, l_2)\)は、 (車\(v\)が\(l_1\)と\(l_2\)の間を移動した回数)/(\(v\)の移動回数)。
で、地域\(i\)にあるGPSつきの車の集合を\(A_i\), その台数を\(V_i\)として \[ \bar{M}_i = (1/V_i) \sum_{v \in A_i} M_v\] を得た。GPSつきの車というのはその地域の車のSRS標本だと考え、 \(\bar{M}_i\) をこの地域の車のモビリティの指標とみなし、その分散を推定した。 これらを使って貧困率の推定を行った。
[おもしれえええええ]
その2, Porter et al(2014)。 USの東半分におけるスペイン語世帯の割合の相対的変化を分析。 州(小地域)の直接推定量はACSで得ておいて、共変量を「よくあるスペイン語の検索語が その州で検索された回数」とし、Google Trendから得た。
[パス]
[省略]
4.3節で述べた母集団モデルが標本について成り立つものとする。
地域\(i=1,\ldots,m\), ユニット\(j=1,\ldots,n_{i}\)について \[ y_{ij} = \mathbf{x}_{ij}^T \mathbf{\beta} + v_i + e_{ij}, \ \ v_i \mathop{\sim}^{iid} (0, \sigma_v^2), \ \ e_{ij} \mathop{\sim}^{iid} (0, k_{ij}^2 \sigma_e^2) \] \(x_{ij}\)は長さ\(p\) (\(< \sum_i n_i\))。
行列で書き換えると、 \[\mathbf{y}_i = \mathbf{X}_i \mathbf{\beta} + v_i \mathbf{1} + \mathbf{e}_i\]
上の式は、ブロック対角共分散構造をもつ一般線形混合モデル(5.3節) \[ \mathbf{y}_i = \mathbf{X}_i \mathbf{\beta} + \mathbf{Z}_i \mathbf{v}_i + \mathbf{e}_i\] の特殊ケースである。
[頭を整理するために、以下付記すると…
ということであろう。]
\(V_i\)の逆行列は…[略]
我々が関心を持っているのは地域平均\(\bar{Y}_i\)である。 モデルによればそれは\(\bar{Y}_i = \bar{\mathbf{X}}_i^T \mathbf{\beta} + v_i + \bar{E}_i\) である。\(N_i\)が大きければ[←\(n_i\)でないことに注意]、\(\bar{E}_i\)は大数法則により 0に近づく。 そこでBattese, Harter, & Fuller (1988)は、 \(\bar{Y}_i\)ではなくて、\(\mu_i = \bar{\mathbf{X}}_i^T \mathbf{\beta} + v_i\) をターゲットにしようと考えた。 さらにいえば、標本割合\(f_i = n_i / N_i\)を無視すれば、\(\bar{Y}\)のEBLUP推定量と \(\mu\)のEBLUP推定量は接近する。
というわけで、以下では\(\mu_i\)の推定に焦点を当てる。\(f_i\)が無視できない場合については 7.1.3節を参照。地域内で情報的抽出が行われていた場合、ないし測定してない地域があり 地域抽出が情報的な場合については7.6.3節を参照。
さて、もとの一般線形混合モデルでは、関心の対象は \(\mu_i = \mathbf{l}_i^T \mathbf{\beta} + \mathbf{m}_i^T \mathbf{v}_i\) であった。いっぽうここでの関心の対象は \(\mu_i = \bar{\mathbf{X}}_i^T \mathbf{\beta} + v_i\) である。 \(\mathbf{l}_i\)が\(\bar{\mathbf{X}}_i\)に、\(\mathbf{m}_i\)が1に変わったわけだ。
もとの一般線形混合モデルでは、BLUP推定量は \[ \tilde{\mu}_i^H = \mathbf{l}_i^T \tilde{\mathbf{\beta}} + \mathbf{m}_i^T \tilde{\mathbf{v}}_i \] \[ \tilde{\mathbf{\beta}} = \left( \sum_i^m \mathbf{X}_i^T \mathbf{V}_i^{-1} \mathbf{X}_i \right)^{-1} \sum_i^m \mathbf{X}_i^T \mathbf{V}_i^{-1} \mathbf{y}_i \]
\[ \tilde{\mathbf{v}} = \mathbf{G}_i \mathbf{Z}_i^T \mathbf{V}_i^{-1} (\mathbf{y}_i-\mathbf{X}_i \tilde{\mathbf{\beta}}) \] であった。ここでは、まず \[ \tilde{\mu}_i^H = \bar{\mathbf{X}}_i^T \tilde{\mathbf{\beta}} + \tilde{\mathbf{v}_i} \] となる。\(\tilde{\mathbf{\beta}}\)は変わらない。\(\tilde{\mathbf{v}}\)は以下となる。 [さあ、ここからが長い…]
まず、ユニットの標本抽出誤差の分散にかける定数\(k_ij\)を \(a_{ij}=k_{ij}^{-2}\)と書き換える。 地域内での合計を\(a_{i.} = \sum_j^{n_i} a_{ij}\)と書く。
各地域について\(a_{ij}\)で重みづけた平均 \(\bar{y}_{ia} = (\sum_j^{n_i} a_{ij} y_{ij})/a_{i.}\), \(\bar{\mathbf{X}}_{ia} = (\sum_j^{n_i} a_{ij} \mathbf{x}_{ij}) / a_{i.}\)を 求めておく。[← 下添字aはインデクスでない! まぎらわしいなあ、もう]
\(\gamma_i = \frac{\sigma_v^2}{\sigma^2_v + \sigma_e^2 / a_{i.}}\)とする。 つまり、説明されていない地域間変動\(\sigma_v^2\)と、 全変動\(\sigma_v^2 + \sigma^2_e / a_{i.}\)との比である。
おまちかねの\(\tilde{v}_i\)は、 \[ \tilde{v}_i = \gamma_i (\bar{y}_{ia} - \bar{\mathbf{x}}^T_{ia}\tilde{\mathbf{\beta}}) \] [なぜこうなるのかわからないんだけど、信じます…]
BLUP推定量\(\tilde{\mu}_i^H\)はこう書き換えることができる。 \[ \tilde{\mu}_i^H = \gamma_i[\bar{y}_{ia} + (\bar{\mathbf{X}}_i - \bar{\mathbf{x}}_{ia})^T \tilde{\mathbf{\beta}}] + (1-\gamma_i) \bar{\mathbf{X}}_i^T \tilde{\mathbf{\beta}} \] 第1項は“survey regression”推定量であり、第2項はregression-synthetic推定量である。 この2つを\(\gamma_i\)で重みづけているわけである。 [第1項がよくわからない… 2章と3章を読んでいないからであろう]
すべての\(k_{ij}\)が1のとき、第1項はSRSの下で近似的にデザイン不偏である[…このくだりもよくわからない、 とりあえず省略]
すべての\(k_{ij}\)が1でSRSのとき、\(\tilde{\mu}_i^H\)はデザイン一致性がある[…省略]
\(\tilde{\mu}_i^H\)のMSEは[…省略]
\(\mathbf{\beta}\)のBLUE \(\tilde{\mathbf{\beta}}\) はOLSで求めることができて[…省略]
そんなこんなでBLUP推定量は手に入ったが、 \(\sigma^2_v\)と\(\sigma^2_e\)が未知である。これを 推定値に置き換えたのがEBLUP推定量である。
では、\(\sigma^2_v\)と\(\sigma^2_e\)をどうやって推定するか。
方法その1, “fitting-of-constant”推定量。地域の数が多い時、 計算が大変になる。[その謎の数式を紹介… 省略]
方法その2, ML推定量ないしREML推定量。 \(v_i\)と\(e_{ij}\)の正規性を仮定する必要がある。そのMSEは[…略]。
[パス]
[パス. 構成のみメモする]
[3つの例が紹介されている。パス]
[パス. 構成のみメモする]
[パス]
[パス. 構成のみメモする]
[省略]
[パス. 構成のみメモする]
[パス]
[パス]
[パス。構成のみメモする]
[この節はほぼ逐語訳]
実際の適用時には、 隣接している地域の直接推定量は相関していることがある。 FHモデルにおける地域レベル補足変数がこの空間的相関を十分に説明していない場合、 空間的に相関している地域効果\(v_i\)をモデルに含めることによって、 最終的な小地域推定の有効性が向上するかもしれない。
4.4.4節で空間的FHモデルについて述べた。 空間的FHモデルでは、 FHモデルにおける共分散行列\(G = \sigma_v^2 \mathbf{I}_m\)が \(G=\Gamma(\mathbf{\delta})\)にかわる。 この\(\Gamma(\mathbf{\delta})\)の形式は、\(\mathbf{v}\)について想定されている具体的な空間モデル によって決まる。 たとえば、
上記の空間モデルは、5.2節の一般線形混合モデルの特殊ケースである。 従って、\(\theta_i = \mathbf{z}_i^T \mathbf{\beta} + b_i v_i\) のBLUP推定量は、パラメータのベクトル\(\mathbf{\delta}\)が決まれば、 5.2節で示した一般的な公式によって手に入る。 実際には\(\mathbf{\delta}\)は未知なので、その推定量 \(\hat{\mathbf{\delta}}\)で置き換えないといけない。 \(\mathbf{\delta}\)のML推定量ないしREML推定量については5.2.4節で述べた。 Crassie & Chan(1989)は空間モデルにおける\(\mathbf{\delta}\)のML推定について論じている。
BLUP推定量のMSEは[…略]。
時空間モデルの場合は[…略]。
EBLUP推定量のMSEのブートストラップ推定について[…略].
robust EBLUP理論について[…略]。
Chandra, Salvati, & Chambers(2007)は、 SARモデルに従った空間的相関を持つランダム効果を含めたユニットレベルモデル について研究している。これは空間FHモデルに類似している。 彼らはテイラー線形化アプローチによって、EBLUP推定量とそのMSE推定量を開発している。[…略]
[事例紹介。パス]
[パス]
[パス]
[パス]
[パス]
地域\(i\)(\(=1, \ldots, m\))のなかのグループ\(j\)(\(=1, \ldots, n_i\))において、 \(K\)カテゴリの質的変数の値の割合が \(p_{ij1}, \ldots, p_{ijK}\)であるとしよう(\(\sum_k p_{ijk} = 1\))。 観察されている人数を \((y_{ij1}, \ldots, y_{ijK})\)とし、\(\sum_k y_{ijk} = m_{ij}\)とする。
\[(y_{ij1}, \ldots, y_{ijK})^T \sim multinomial(m_{ij}; p_{ij1}, \ldots, p_{ijK})\] と仮定すれば、ロジスティック混合モデル \[ \log(p_{ijk}/p_{ijK}) = \mathbf{x}_{ij}^T \mathbf{\beta}_k + v_{ik}, \ \ k = 1, \ldots, K-1 \] \[(v_{i1}, \ldots, v_{i,K-1})^T \mathop{\sim}^{iid} N_{K-1}(\mathbf{0}, \mathbf{\Sigma}_v)\] を想定できる。 なお、\(m_{ij}=1, K=2\)とすればこれは4.6節のロジスティック混合モデルである。
[ここ関心あるので、ほぼ逐語訳]
Monila, Saei, & Lombardia (2007)は、UKの406地域[説明略]の労働力の特徴を推定 するため、次の多項ロジットモデルを考えた。
\[ \log(p_{ijk}/p_{ijK}) = \mathbf{x}_{ij}^T \mathbf{\beta}_k + v_i, \ \ v_i \mathop{\sim}^{iid} N(0, \sigma_v^2) \] つまり、すべてのカテゴリ\(k\)において同じ地域効果があると考えたわけである。 カテゴリは、非雇用(\(k=1\)), 雇用(\(k=2\)), 労働力でない(\(k=3\))であった。 グループを性年代(6水準), 共変量は、登録失業者の比率の対数、その他22個のダミー変数であった[説明略]。
penalized quasi-likelihood法(PQL法)に基づく 2段階反復アルゴリズムを使ってモデルをあてはめ、 回帰係数\(\mathbf{\beta} = (\mathbf{\beta}_1^T, \mathbf{\beta}_2^T)^T\) とランダム効果\(\mathbf{v}=(v_1, \ldots, v_m)^T\)を推定した。 また、\(\sigma_v^2\)は近似的なMLないしREMLで推定した。 このような結合型のアルゴリズムはSchall(1991)が導入したもので、 小地域推定の文脈では、Saei & Chambers(2003)がGLMMのあてはめに用いている。
PQLは、\(\sigma_v^2\)の所与の固定値の下で、 回帰係数\(\mathbf{\beta}\)とランダム効果\(\mathbf{v}\)を 同時に推定する。その方法は、 データの周辺対数尤度を最大化するかわりに、 観察\((y_{ij1}, \ldots, y_{ijK})^T\)と\(\mathbf{v}\)の同時対数尤度を 最大化するというものである。 \((y_{ij1}, \ldots, y_{ijK})^T\)が多変量正規ならば、 PQLによる\(\mathbf{\beta}\)と\(\mathbf{v}\)の推定値は、 \(\mathbf{\beta}\)のML推定量と\(\mathbf{v}\)のEBLUP推定量に一致する。 ここで同時対数尤度は閉形式をとることに注意されたい。 周辺対数尤度だとそうではない。
第2段階では、\(\mathbf{\beta}\)と\(\mathbf{v}\)を 固定するために、Shall(1991)が適用した方法を用いる。すなわち、 モデルを線形化し、周辺尤度を正規近似することで、 分散成分(ここでは\(\sigma_v^2\))の近似的なML推定値ないしREML推定値を得る。
この2段階の手続きを反復することで、 必要なすべての量の推定値、すなわち\(\mathbf{\beta}\), \(\mathbf{v}\), 分散成分の推定値が 手に入る。推定された値を\(\hat{\mathbf{\beta}}_1, \hat{\mathbf{\beta}}_2, \hat{v}_i\)として、確率は次のように推定される。 \[ \hat{p}_{ijk} = \frac{ \exp(\mathbf{x}_{ij}^T \hat{\mathbf{\beta}}_k + \hat{v}_i) } { 1 + \sum_k^2 \exp(\mathbf{x}_{ij}^T \hat{\mathbf{\beta}}_k + \hat{v}_i) }, \ \ k=1,2 \] \[ \hat{p}_{ij3} = 1-\sum_k^2 \hat{p}_{ijk} \] [一本目の数式, 添字\(k\)がダブってるけど、云いたいことはわかる]
[以下略]
[パス]
[省略]
[パス。構成のみメモする]
階層ベイズアプローチでは、モデル・パラメータ\(\mathbf{\lambda}\)の 主観的な事前分布\(f(\mathbf{\lambda})\)を指定し、 データ\(\mathbf{y}\)の下で、関心ある小地域(ランダム)パラメータ \(\mathbf{\mu}\)の事後分布\(f(\mathbf{\mu}|\mathbf{y})\) を得る。
\(f(\mathbf{\mu}|\mathbf{y})\)を得るためには、 \(f(\mathbf{y}|\mathbf{\mu}, \mathbf{\lambda}_1)\)と \(f(\mathbf{\mu} | \mathbf{\lambda}_2)\)の2段階モデルと、 \(\mathbf{\lambda} = (\mathbf{\lambda}_1^T, \mathbf{\lambda}_2^T)^T\) の主観的事前分布を、ベイズの定理をつかって結合する。
HBアプローチは\(f(\lambda)\)の指定を必要とする。それは情報的なものでもいいし、 拡散的(無情報的)なものでもよい。なんなら非正則なものでもよいが、 確実に正則な事後分布が得られるものだということが重要である。 また、無情報事前分布は 頻度主義の枠組みから見てもwell-calibratedな推論ができるようなものであることが 望ましい。たとえば、いま\(\phi=h(\mathbf{\mu})\)に関心があり、 その事後平均\(\hat{\phi}^{HB} = E[h(\mathbf{\mu}|\mathbf{y})]\)を 得たとして、 その頻度主義的な意味でのバイアス\(E(\hat{\phi}^{HB} - \phi)\)と、 事後分布の分散(それはMSEの推定量になる)の頻度主義的な意味でのバイアスが、 小さくないと困る。
さて、ベイズの定理より \[ f(\mathbf{\mu}, \lambda|\mathbf{y}) = \frac{f(\mathbf{y}, \mathbf{\mu}|\lambda) f(\lambda)}{f_1(\mathbf{y})} \] である。分母の\(f_1(\mathbf{y})\)は\(\mathbf{y}\)の周辺分布で[…式省略]。 上の式から、求めている事後密度 \[ f(\mathbf{\mu} | \mathbf{y}) = \int f(\mathbf{y}, \mathbf{\mu}|\lambda) f(\lambda) d\lambda \] が得られる。
\(f(\mathbf{\mu} | \mathbf{y})\)の評価には多次元の積分が 必要になる。 単純な問題であれば解析的に解ける場合も多い。解けない場合でも MCMC法でなんとかなる。
[パス。構成のみメモする]
6.1節の基本的な地域レベルモデル \[\hat{\theta}_i = \mathbf{z}_i^T \mathbf{\beta} + b_i v_i + e_i\] にHBアプローチを適用しよう。
まずは\(\sigma^2_v\)が既知の場合について考える。 \(\mathbf{\beta}\)についてはフラットな事前分布\(f(\mathbf{\beta}) \propto 1\)を 考えよう。 すると、次のように書き換えられる。 \[ \begin{aligned} {\rm (i)} \ \ & \hat{\theta}_i | \theta_i, \mathbf{\beta}, \sigma^2_v \mathop{\sim}^{iid} N(\theta_i, \psi_i), \ \ i = 1, \ldots, m \\ {\rm (ii)} \ \ & \theta_i | \mathbf{\beta}, \sigma^2_v \mathop{\sim}^{iid} N(\mathbf{z}_i^T, b_i^2 \sigma_v^2), \ \ i = 1, \ldots, m \\ {\rm (iii)} \ \ & f(\mathbf{\beta}) \propto 1 \end{aligned} \] \(\sigma^2_v\)が未知の場合は、その事前分布\(f(\sigma^2_v)\)を設定し、(iii)を書き換えて \[ \begin{aligned} {\rm (iii)'} \ \ & f(\mathbf{\beta}, \sigma^2_v) = f(\mathbf{\beta}) f(\sigma^2_v) \propto f(\sigma^2_v) \end{aligned} \] とする。
上のHBモデルで、\(\hat{\mathbf{\theta}} = (\hat{\theta}_1, \ldots, \hat{\theta}_m)^T\) と\(\sigma^2_v\)が所与ならば、 事後分布の平均はBLUP推定量\(\tilde{\theta}_i^H\)となり、 分散は6.1節で与えたMSEとなる。
実際には\(\sigma^2_v\)は未知なので、(iii)’式で書いたように事前分布\(f(\sigma^2_u)\)を 決めないといけない。
\(\theta_i\)の推定量\(\hat{\theta}_i^{HB}\)はどうなるか。 それは\(\hat{\mathbf{\theta}}\)と\(\sigma^2_v\)に依存しているから、 \(\sigma^2_v\)の事後分布\(f(\sigma^2_v|\hat{\mathbf{\theta}})\)を通じた 期待値を求める必要がある。言い換えると\(f(\sigma^2_v|\hat{\mathbf{\theta}})\)さえ 決まってしまえば、1次元の積分で求めることができる。 \(\theta_i\)の事後分布の分散も[…式省略…], \(f(\sigma^2_v|\hat{\mathbf{\theta}})\)さえ 決まってしまえば1次元の積分で求めることができる。
では、\(\sigma^2_v\)の事後分布\(f(\sigma^2_v|\hat{\mathbf{\theta}})\)をどうやって求めるか。 これは事前分布\(\sigma^2_u\)と、制約つき尤度関数の積として 求めることができて[…式省略]。 ここからわかるのは、事前分布を\(f(\sigma_v^2) \propto 1\)としても(非正則事前分布)、 \(m > p + 2\)であれば、事後分布\(f(\sigma^2_v|\hat{\mathbf{\theta}})\)は正則だという ことである。
\(f(\sigma_v^2) \propto 1\)としたとき、\(\sigma^2_v\)の事後分布の平均は…[略]。
[以下いろいろ書いてあるけどパス]
Gibbsサンプリングの場合は、\(\sigma_v^{-2} \sim G(a,b)\)として…[以下パス]
[パス]
地域効果\(\mathbf{v} = (v_1, \ldots, v_m)^T\)に 空間相関をいれた、次のモデルを考える。\(s^2_i\)は観察値の不偏分散。\(\mathbf{R}\)は隣接行列で、 対角には「隣接している地域数」、非対角は隣接している箇所に\(-1\)がはいる。\(U\)は一様分布を表す記号。 \[ \begin{aligned} {\rm (i)} \ \ & \bar{y}_i | \theta_i, \sigma^2_i \mathop{\sim}^{iid} N(\theta_i, \sigma^2_i/n_i) \\ {\rm (ii)} \ \ & \frac{(n_i-1)s^2_i}{\sigma^2_i} | \sigma^2_i \mathop{\sim}^{iid} \chi^2_{n_i-1} \ \ i = 1, \ldots, m \\ {\rm (iii)} \ \ & (\theta_1, \ldots, \theta_m)^T|\beta, \sigma^2_v \sim N_m (\mathbf{Z\beta}, \sigma^2_v\mathbf{D}^{-1}), \ \ \mathbf{D} = \lambda\mathbf{R}+(1-\lambda)\mathbf{I} \\ {\rm (iv)} \ \ & f(\mathbf{\beta}) \propto 1, \ \ \sigma^{-2}_i \sim G(a,b,), \ \ \sigma_v^2 \sim G(a_0, b_0), \ \ \lambda \sim U(0,1) \end{aligned} \] [そうか、このモデルでは\(\phi_i = \sigma^2/n_i\)が未知なのか]
これはCARモデルになっている。解くのが簡単。[適用例…略]
[パス]
7.1節のモデルは以下であった。 \[ y_{ij} = \mathbf{x}_{ij}^T \mathbf{\beta} + v_i + e_{ij} \] \[v_i \mathop{\sim}^{iid} (0, \sigma_v^2)\] \[e_{ij} \mathop{\sim}^{iid} (0, k_{ij}^2 \sigma_e^2)\]
HBモデルに書き換える。なお、誤差分散は等しいものとする(\(k_{ij}^2=1\)) 。
まずは\(\sigma_v^2\), \(\sigma_e^2\)が既知の場合について考える。 \[ \begin{aligned} {\rm (i)} \ \ & y_{ij} | \mathbf{\beta}, v_i, \sigma^2_i \mathop{\sim}^{ind} N(\mathbf{x}_{ij}^T \mathbf{\beta} + v_i, \sigma_e), \ \ i = 1, \ldots, m, \ \ j = 1, \ldots, n_j \\ {\rm (ii)} \ \ & v_i | \sigma_v^2 \mathop{\sim}^{ind} N(0, \sigma_v^2), \ \ i = 1, \ldots, m \\ {\rm (iii)} \ \ & f(\mathbf{\beta}) \propto 1 \end{aligned} \] \(\sigma_v^2\), \(\sigma_e^2\)が未知の場合は、その事前分布を \(f(\sigma_v^2)\), \(f(\sigma_e^2)\)とし、(iii)を書き換えて \[ \begin{aligned} {\rm (iii)'} \ \ & f(\mathbf{\beta}, \sigma_v^2, \sigma_e^2) \propto f(\sigma^2_v) f(\sigma^2_e) \end{aligned} \] 以下では話を簡単にするため、 地域の母集団サイズ\(N_i\)は十分大きく、 \(\mu_i = \bar{\mathbf{X}}^T_i \mathbf{\beta} + v_i\)は 地域平均\(\bar{Y}_i\)と同じだと考えることにする。
\(\sigma^2_v\)と\(\sigma^2_e\)が既知のとき、 \(\mu_i\)の事後分布の平均は7.1節のBLUP推定量となり、 分散は7.1節で与えたMSEとなる。
\(\sigma^2_v\)と\(\sigma^2_e\)は実際には未知なので、事前分布を 設定しないといけない。
\(\mu_i\)の推定量\(\hat{\mu}_i^{HB}\)とその分散は、 \(\mathbf{y}\)と事後分布\(f(\sigma_v^2, \sigma_e^2 | \mathbf{y})\) に依存しているから、\(f(\sigma_v^2, \sigma_e^2 | \mathbf{y})\)を通じた期待値と分散 を求める必要がある。 言い換えると、\(f(\sigma_v^2, \sigma_e^2 | \mathbf{y})\)さえ決まれば 2次元の積分で求められる。
\(f(\sigma_v^2, \sigma_e^2 | \mathbf{y})\) は事前分布の積\(f(\sigma_v^2)f(\sigma_e^2)\)と 制限つき尤度関数[…式省略…]の積で決まる。 \(f(\sigma_v^2) \propto 1, f(\sigma_e^2) \propto 1\) とすれば、\(f(\sigma_v^2, \sigma_e^2 | \mathbf{y})\)は 正則になることが示せる。
さらに次のように工夫する。 \(f(\sigma_e^2) \propto 1\)ではなく、 \(\sigma_e^{-2} \sim G(a_e, b_e), \ \ a_e \geq 0, \ \ b_e > 0\) とする。こうすると、 \(\tau_v = \sigma_v^2 / \sigma_e^2\)として、 \(f(\sigma_e^2 | \tau_v, \mathbf{y})\)を通じて\(\sigma_e^2\)を 積分消去できる。 さっきは2次元の積分が必要だったけど、 今度は\(f(\tau_v | \mathbf{y})\)を通じた1次元の積分で済む。 で、 \(\tau_v\)の事前分布を \(\tau^{-1}_v \sim G(a_v, b_v), \ \ a_v \geq 0, \ \ b_v > 0\) とする。
[む、むずかしい… ]
Gibbsサンプリングの場合は、 \[\sigma_v^{-2} \sim G(a_v, b_v), \ \ a_v \geq 0, \ \ b_v > 0\] \[\sigma_e^{-2} \sim G(a_e, b_e), \ \ a_e \geq 0, \ \ b_e > 0\] として…[以下略]
[パス]
[パス]
[パス。構成のみメモする]
[パス]
[パス]
[パス。構成のみメモする]
[目的変数がポワソン分布するというモデル。パス。構成のみメモする]
[パス]
\(y_{ij}\)が二値だったらどうするか。
\[ \begin{aligned} {\rm (i)} \ \ & y_i | p_i \mathop{\sim}^{iid} binomial(n_i, p_i) \\ {\rm (ii)} \ \ & p_i|\alpha, \beta \mathop{\sim}^{iid} beta(\alpha, \beta), \ \ \alpha > 0, \ \ \beta > 0 \\ {\rm (iii) } \ \ & \alpha \sim G(a_1, b_1), \ \ \beta \sim G(a_2, b_2), \ \ a_1 > 0, \ \ b_1 > 0, \ \ a_2 > 0, \ \ b_2 > 0 \end{aligned} \] なお\(\alpha, \beta\)は互いに独立。
Gibbsサンプリングでは…[略]
[事例紹介。パス]
ベータ二項分布は共変量をいれるのがむずかしい。共変量がある場合は 以下のモデルが考えられる。
その1, 地域レベル共変量をいれる場合。
\[ \begin{aligned} {\rm (i)} \ \ & y_i | p_i \mathop{\sim}^{iid} binomial(n_i, p_i) \\ {\rm (ii)} \ \ & logit(p_i) = \mathbf{z}_i^T \mathbf{\beta} + v_i, \ \ v_i \mathop{\sim}^{iid} N(0, \sigma_v^2)\\ {\rm (iii) } \ \ & f(\mathbf{\beta}) \propto 1, \ \ \sigma_v^{-2} \sim G(a, b), \ \ a \geq 0, \ \ b > 0 \end{aligned} \] Gibbsサンプリングでは…[略]
その2, ユニットレベル共変量をいれる場合。 \[ \begin{aligned} {\rm (i)} \ \ & y_{ij} | p_{ij} \mathop{\sim}^{iid} Bernoulli(p_{ij})\\ {\rm (ii)} \ \ & logit(p_{ij}) = \mathbf{x}^T_{ij} \mathbf{\beta} + v_i, \ \ v_i \mathop{\sim}^{iid} N(0, \sigma_v^2)\\ {\rm (iii) } \ \ & f(\beta) \propto 1, \ \ \sigma_v^{-2} \sim G(a, b) \ \ a \geq 0, \ \ b>0 \end{aligned} \] (iii)はなんなら\(f(\mathbf{\beta}, \sigma^2_v) \propto 1\)としてもよい。
Gibbsサンプリングでは…[略]
[パス]
[パス]
[パス]
[パス]
[パス]
[パス. 構成のみメモする]
おしまい！