« 読了：Jing & De Oliveira (2015) RのgeoCountパッケージ | メイン | 覚え書き：スタイン推定量とはなんぞや »
2017年10月27日 (金)
Bivand, R.S., Gomez-Rubio, V., Rue, H. (2015) Spatial data analysis with R-INLA with some extensions. Journal of Statistical Software, 63(20).
MCMCではない謎の方法(integrated nested Laplace approximation; 積分段階的ラプラス近似)で空間統計モデルを推定しちゃう謎のRパッケージ R-INLA の解説書。できればそういう難しい話とは関わり合いを持ちたくなかったんだけど。人生とはままならないものだ。
モデル。
観察$\mathbf{y} = \{y_i\}_{i=1}^{n}$が指数型分布族に従っていて、$y_i$の平均$\mu$が線形予測子$\eta_i$となんらかのリンク関数でつながっているとする。
$\eta_i = \alpha + \sum_j^{n_f} f^{(j)}(u_{ij}) + \sum_k^{n_\beta} \beta_k z_{ki} + \epsilon_i$
$f^{(j)}$は共変量$\mathbf{u}$のランダム効果を表す$n_f$個の関数。$\beta_k$は共変量$\mathbf{z}$の関数。
潜在効果のベクトルを$\mathbf{x} = \{\{\eta_i\}, \alpha, \{\beta_k\}, \ldots\}$と書く。
[←この論文、よくわからなくて困っているんだけど、このメモをとった翌日になってじっくり読み直してみて、混乱の根底はここだと気がついた。ここに$\eta_i$と$\alpha, \beta_k$の両方がはいる理由がよく分からない。$\mathbf{x}$は地点を要素にもつ、共変量の効果抜きの効果なのではないかと思うのだけれど...$\mathbf{x} = \{\epsilon_i\}$というのなら納得するんだけど...]
$\mathbf{x}$について、平均0, 共分散行列$Q(\theta_1)$の正規分布を仮定する。さらに、「$\mathbf{x}$はガウシアン・マルコフ確率場(GRMF)であると仮定する。これは$Q(\theta_1)$が数多くのマルコフ特性を満たすことを意味する。」[←素人を殺しに来たぞ... なにをいっているのかさっぱりだが、諦めずに読み進めることに]
$y_i$の分布は潜在効果$\mathbf{x}$に依存するが、他のハイパーパラメータ$\theta_2$にも依存しうる。[←どういうこと？ $y_i$の方程式の右側に$\eta_i$以外の項があるってこと？] 以下では$\theta=(\theta_1, \theta_2)$としよう。
観察$y_i$は所与の$x_i, \theta$の下で互いに独立である。なぜなら$\mathbf{x}$はGRMFだから。[←わ、わからん...いったい何を云っているのだ...]
$\mathbf{x}$と$\theta$の事後分布はこうなる。
$\pi(\mathbf{x}, \theta | \mathbf{y})$
$\propto \pi(\theta) \pi(\mathbf{x} | \theta) \prod_{i \in I} \pi(y_i | x_i, \theta)$
$\propto \pi(\theta) |\mathbf{Q}(\theta)|^{\frac{1}{2}} \exp (-\frac{1}{2} \mathbf{x}^T \mathbf{Q}(\theta) \mathbf{x} + \sum_{i \in I} \log \pi (y_i | x_i, \theta) )$
ただし、$I$は観察データのインデクス($1, \ldots, n$)。$\mathbf{Q}(\theta)$は$\theta$の精度行列。最後に出てくる$ \log \pi (y_i | x_i, \theta)$は$y_i$の対数尤度。
[えーと、数式2行目から3行目に行くところで完全にギブアップなんだけど、信じましょう]
なお、INLAは指数型分布族だけではなく、たとえば混合分布であっても対応できる。$\mathbf{x}$について上記と異なる定式化も可能。事前分布$\pi(\theta)$についてもとても柔軟である。
さて、$\mathbf{x}$と$\theta$の周辺分布は
$\pi(x_i|\mathbf{y}) = \int \pi(x_i|\theta, \mathbf{y}) \pi(\theta|\mathbf{y}) d\theta$
$\pi(\theta_j|\mathbf{y}) = \int \pi(\theta|\mathbf{y}) d\theta_{-j}$
となるのだが、この右辺の項をだね、うまいこと近似する方法があるのだ。おまえのようなアホに説明することは到底不可能だが、ラプラス近似という方法をさらにごにょごにょしちゃうことによってなにかがどうにかなっちゃうのだ。[←こんな風には書いてないが、私としてはこんな感じ。正直、ここの数段落は全く理解できない]
というわけで、この方法を実装したソフトINLA、そしてそれをRで使うパッケージR-INLAを作ったので、使っても苦しゅうないぞ。
パッケージの具体的な使い方。最初はラティスデータの分析の話から。
まず次の点に注意。無相関なランダム効果を入れる場合、ランダム効果は
$u_i \sim N(0, \tau_u)$
と定義され、$\tau_u$じゃなくて$\log(\tau_u)$に事前分布が与えられる(デフォルトでは対数ガンマ分布)。
空間的相関をモデル化するためには隣接エリアを定義しないといけない。空間的自己相関は、平均0, 共分散が
$\Sigma = (1/\tau) Q^{-1}$
の正規分布でモデル化される。$Q$は、エリア$i$と$j$が隣接じゃない時には$Q_{ij}=0$となる。$Q$はふつうすごくスパースになる。[←あー！それでマルコフ場っていうのか！おとなりさんとしか相関しないというモデルなのね。やっと意味が分かった...そりゃあ計算も普通のクリギングより速くなるかもな]
こういうときの$Q$の決め方は、いろいろご用意してある。いくつかご紹介しよう。
ひとつめ、intrinsic CARモデル。$Q_{ii}$を隣接領域数$n_i$とし、 $Q_{ij}(i \neq j)$は隣接していれば$-1$, でなければ$0$。これは結局次のモデルになる。空間的ランダム効果を$v_i$、精度を$\tau_v$として, $i \neq j$のとき
$\displaystyle v_i | v_j, \tau_v \sim N\left( \frac{1}{n_i} \sum_{i \sim j} v_j, \frac{1}{\tau_v n_i} \right) $
$\log(\tau_v)$には事前分布として対数正規分布を与える。
その２、proper CARモデル。上のモデルをちょっと変えて分布を正則にする。$d$は正の値。
$\displaystyle v_i | v_j, \tau_v \sim N\left( \frac{1}{n_i+d} \sum_{i \sim j} v_j, \frac{1}{\tau_v (n_i+d)}\right)$
その３、もっと一般的なアプローチ。$I$を単位行列、$\rho$を空間的自己相関パラメータ、$C$を隣接行列, $\lambda_{max}$を$C$の最大固有値として
$\displaystyle Q = (I - \frac{\rho}{\lambda_{max}} C)$
で、$\log(\rho/(1-\rho))$に事前分布として正規分布を与える。
[こうやって列挙されてもねえ... それぞれのアプローチのモチベーションとか長短とかを教えてくださいよ...]
ここからは、具体的な分析例がたくさん。
大事そうなことがいろいろ書いてあるんだけど、力尽きたのでいったん読了としておく。ま、雰囲気がわかったのでよしとしよう。あれだね、バリオグラム関数はモデル化せず、隣との相関だけをモデル化するわけね。
論文：データ解析(2015-) - 読了：Bivand, Gomez-Rubio, Rue (2015) R-INLAパッケージで楽しい楽しい空間統計モデリング