« 2017年9月 | メイン | 2017年11月 »
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パッケージで楽しい楽しい空間統計モデリング
2017年10月26日 (木)
Jing, L., De Oliveira, V. (2015) geoCount: An R package for the analysis of geostatistical count data. Journal of Statistical Software. 63(11).
カウントデータに特化した空間モデリングパッケージ geoCountについての紹介。実戦投入前の儀式として読んだ。
カウントデータの空間モデリングにはふつう一般化線形空間モデル(GLSM)が使われるが、実際にはいろいろ大変である。(1)なにより計算が大変。メモリが足らん。(2)MCMCがめっちゃ遅い。(3)MCMCが収束しない。
そこで新しいパッケージをお届けしよう。このパッケージは、C++で書いてあり速い。並列計算に対応しておる。そして効率的な新アルゴリズム。それではご紹介しましょう！ geoCountパッケージです！
モデル。
$\{S(\mathbf{x}): x \in A\}$を正規確率場とする。$\mathbf{S} = (S(\mathbf{x}_i), \ldots, S(\mathbf{x}_n))^T$とする。
$Y_i | S(\mathbf{x}_i) \sim p(\cdot | \mu_i), \ \ i=1, \ldots, n$
$\mu_i = t_i g^{-1} (S(\mathbf{x}_i))$
$\mathbf{S} \sim N_n (D\beta, \Sigma)$
$(\beta, \theta) \sim \pi(\beta, \theta)$
説明しよう。
一本目。$\mu_i = E(Y_i | S(\mathbf{x}_i))$である。$p(\cdot | \mu_i)$については後述。
二本目。$g(\cdot)$はリンク関数である。
三本目。$D$はフルランク計画行列で、なかに共変量がはいっておる。$\Sigma = (\sigma_{ij})$は正定な共分散行列で、要素$\sigma_{ij}=\sigma^2 \rho(u_{ij})$はユークリッド距離$u_{ij}=||\mathbf{x}_i - \mathbf{x}_j||$で決まる(つまり等方性がある)。
四本目。$\pi(\beta, \theta)$はパラメータの事前分布。
$\rho(u)$としては、球面ファミリー、Maternファミリー、power exponentialファミリーをご用意しておる。
$p(\cdot | \mu_i)$とリンク関数$g(\cdot)$については、良く用いられる次の２つをご用意。
その１、ポワソン対数正規モデル：
$Y_i | S(\mathbf{x}_i) \sim Poisson(\mu_i)$
$\mu_i = t_i \exp(S(\mathbf{x}_i))$
その２、二項ロジット正規モデル：
$Y_i | S(\mathbf{x}_i) \sim Binomial(t_i, \mu_i/t_i)$
$\mu_i = t_i \exp(S(\mathbf{x}_i)) / (1+\exp(S(\mathbf{x}_i)))$
事前分布は次のようにご用意しております...[めんどくさいので略]
プログラミング上の工夫と並列処理...[パス]
アルゴリズムの工夫...[パス]
プログラム例... [パス。実際に使うときに読めばいいや]
他のパッケージとの比較。
- 同様のモデルはgeoRglmパッケージでも組めるけど、MCMCの際のパラメータ化のやり方が違っており...[よくわからんので中略... 何だか知らんが、うちらのほうがイケているよとのこと]。ただし、計算時間はこっちのほうがかかるかも。[←よく考えてみたらオイオイって話ですね。C++かつ効率的アルゴリズムで速いんちゃうんかと]
- INLAパッケージと比べると...そもそもINLAは潜在ガウスモデルを共分散関数として特徴づけるのではなくて、ガウス・マルコフ確率場(GMRF)として特徴づける。ベイズ推論の際、MCMCは周辺分布を確率的に近似しようとするが、INLAはラプラス近似とnumerical quadrature rules[←なにそれ]を用いて決定論的に近似する。
GMRFを使う利点は計算が速いこと。欠点は、定常性とか等方性といった特徴を持つ共分散関数を持つGRMFを構築するのが難しいこと。しかしこれには部分的解決策が提案されている。Matern族の一種である、ある族に属する共分散関数を持つガウス確率場は、ガウシアン・ホワイトノイズを持つある種のstochastic partial differential equations(SPDE)の解であることが示されているのだ。これらの解はある種のGMRFによって与えられる、ランダム係数を持つある種の基底関数へと拡張できることが示されている。その結果、Matern族の共分散関数を使ってモデリングつつ、それと対応する確率場をGMRFで表現して計算をすることができる。これがINLAで採用されているアプローチである。[←ああなるほど...って、何言ってんだかさっぱり分かんないよこの野郎]
geoCountと比べてみると、計算は確かに速い。結果は似てたり違っていたり。[←詳細略]
限界。
- データがでかいと遅くなる。たまに収束しないこともある。どちらの場合もINLAがおすすめ。
- ポワソン分布のときは対数リンク、二項の場合はロジットリンクしか選べない。実はidentityリンクのほうがいいという話もある[←へー]。
- 事前分布についても改善の余地がある。
云々。
2017/10/27追記: INLAと比較しているくだりについて、メモを追加。
論文：データ解析(2015-) - 読了：Jing & De Oliveira (2015) RのgeoCountパッケージ
空間データのモデリングではバリオグラムというのを描くけれど、たいてい従属変数は量で、誤差に正規性が仮定されている。もし従属変数が二値とかカウントとかだったら、モデルはGLMとして組めるにしても、バリオグラムはいったいどうなるのか。
これ、ここんところの疑問の一つだったのだが、このたびやっと説明を見つけたのでメモしておく。すいません、純粋に自分のための覚え書きです。
$S(x)$を平均0, 分散$\sigma^2$の定常ガウス過程とする。観察$Y_i$は互いに独立で平均$\mu_i=g(\alpha+S(x_i))$, 分散$v_i = v(\mu_i)$とする。$g(\cdot)$はリンク関数$h(\cdot)$の逆関数である。
$Y$の素直なバリオグラムについて考えよう。横軸は$u=||x_i - x_j||$、縦軸は
$\tau_Y(u)=E[(1/2)(Y_i-Y_j)^2]$
これを展開するとこうなる。
$\tau_Y(u)$
$= (1/2)E_S[{g(\alpha+S(x_i)) - g(\alpha+S(x_j))}]$
$+ E_s[v(g(\alpha+S(x_i)))]$
第2項に$i$だけ出てきているのは、$S(x_i)$の周辺分布が場所を問わず一定だから。よくよくみればこの項は定数なので、以後$\bar{\tau}^2$と書く。これ、$S(\cdot)$の条件付き分布の分散を平均したもので、ガウシアンモデルで言うところのナゲットである。
さて。第一項に出てくる$g(\alpha+S(x_i))$が、$g(\alpha) + S(x_i) g'(\alpha)$とテイラー近似できるのはわかるね。[←わかんないよ！こっちは文系なんだから！]
すると結局こうなる。
$\tau_Y(u) \approx g'(\alpha)^2 \gamma_S(u) + \bar{\tau}^2$
つまり、バリオグラムの縦軸は、潜在ガウス過程$S(\cdot)$だけを取り出したバリオグラム$\gamma_S(u)$に比例し、そこにさきほどの平均ナゲット効果を乗せたものになる。
というわけで、上記のような通常のバリオグラムをナゲットのサイズ把握のために使うのはOK。しかし、GLMは$S(x)$と$Y$のあいだに非線形的関係を仮定しているわけだから、診断のために使うのはおかしい。[←「診断のため」ってどういうことだろう。バリオグラム・モデルの診断には使えないということだろうか]
バリオグラム・モデルをどうやって推定するか。
まずは一般的に、GLMの尤度関数について考えよう。ランダム効果$S$の下での$Y$の条件付き分布を$f_i(y_i; S, \theta)$とすると、パラメータ$\theta$の条件付き尤度は
$L(\theta | S) | \prod_i f_i(u_i|S, \theta)$
混合モデルの尤度は、$S$の同時分布を$g(S; \phi)$として
$L(\theta, \phi) = \int_S \prod_i f_i(u_i|S, \theta) g(s|\phi) ds$
さて、こんな風に多重積分が掛かっている尤度が計算できちゃうのは、ひとえに$S_i$が互いに独立だからである。しかし空間データでは$S_i = S(x_i)$は互いに独立でないので、こんなもん計算できっこない。実をいうとこの積分を近似する方法も提案されているんだけど、$S$の分散が大きいときには役に立たない。
というわけで... モンテカルロ最尤法、階層尤度、GEEなどの方法が提案されておる、という話になるのだが、力尽きたのでこの辺にしておこう。
以上, Diggle & Ribeiro (2010) "Model-based Geostatistics", Chapter 4-5より。勤務先の経費とはいえ2万円以上した本だったので、元を取らねば、と読んだ次第。
実をいうと、ポワソン対数線形モデルを使ったクリギングの際に、いい感じのバリオグラム・パラメータを視認で決める方法について知りたかったんだけど、それはよくわかんなかった。先の章に書いてあるのかも。
雑記：データ解析 - GLM空間モデルのバリオグラムとはなんぞや
« 2017年9月 | メイン | 2017年11月 »