R ネタです。一般化線形モデル（Generalized linear mixed model : GLMM）をベイズ推定する MCMCglmm というパッケージを使ってポアソン回帰をするためのメモです。もしかすると MCMCglmm の入門的、あるいはチュートリアル的なものになっている、かもしれません。
MCMCglmm について
- 基本的にギブスサンプラー。一部（狭い意味の）メトロポリス・ヘイスティング、スライスサンプラー。
- パラメータの事前分布はデフォルトで平均ゼロ、分散 108 の正規分布。分散は逆ウィッシャート。デフォルトパラメータは \(V=1, \nu=0.002\) 。この逆ウィッシャート分布のモードは \(\nu V/(\nu+2)\) で平均は \(\nu V/(\nu-2)\) （つまりデフォルトでは平均値は存在せず、モードはおよそ 0.001 でゼロの周りでピーキーな分布）。
- lmer とかと大体同じ使い方だけど、random effect の設定はかなりちがう。指定の仕方は豊富だけども、そのぶん理解するのが大変。
- 使える family がかなり豊富。
- 罰則スプライン回帰（additive model）もできるっぽい
- Fixed effect は一番最初の引数に、random effect は random=～ で指定、それ以外の誤差構造は rcov=～ で指定。rcov をきちんと書けばマルチレスポンスな回帰もできるみたい。空間統計もできるかな？
- イテレーション回数は nitt、バーンインは burnin、サンプリング間隔は thin で指定します。つまり、サンプリングされる系列の数は (nitt-burnin)/thin となります。
ということで、使いこなせればいろいろできそうなのだけれども、やってみると結構つまずいたりする。それが MCMCglmm。一番簡単と思われるポアソン回帰でもかなりハマりました（４時間！！）。この記事はその点について書いてみます。
ポアソン回帰
ポアソン回帰はこんな感じのモデル式です。random effect はとりあえず考えません。
\begin{align*}&y\sim Poisson(\exp(\eta))\\ &\eta=\beta’x\end{align*}
テストデータ
超簡単なテストデータ（実はこれがハマる原因なのですが…）を準備。
\begin{align*}Y_i\sim Poisson(\exp(2x+1))\end{align*}
とりあえず、glm （ふつうの一般化線形モデル、mixed じゃないなやつ）をつかえば何の問題もなく推定できるたぐいのデータです。
library(MASS) x <- rnorm(200) y <- sapply(x,function(z)rpois(1,exp(2*z+1))) d <- data.frame(y=y,x=x) g <- glm(y~x,data=d,family="poisson") summary(g) ## Coefficients: ## Estimate Std. Error z value Pr(>|z|) ## (Intercept) 1.09239 0.04224 25.86 <2e-16 *** ## x 1.95143 0.02524 77.31 <2e-16 ***
これはベイズ推定ではなく、点推定ですが、うまく推定できていますね。
ベイズ推定してみる
MCMCglmm を使ってみましょう。推定結果（サンプリング結果）は MCMCglmm オブジェクトの Sol に格納されます。
library(MCMCglmm) x <- rnorm(200) y <- sapply(x,function(z)rpois(1,exp(2*z+1))) d <- data.frame(y=y,x=x) g <- MCMCglmm(y~x,data=d,family="poisson") posterior.mode(g$Sol) ## (Intercept) x ## 0.9828072 2.0057581
最後は事後分布のモードを表示させています。この結果を見る限り、いいんじゃないの？と思うかもしれません。
しかし、サンプル時系列を表示させてみると…
plot(g$Sol)
あぁ、これはひどい！となるわけです。
実際、自己相関をとってみると、次のようにものすごく高い自己相関をもっていることがわかります。
autocorr(g$Sol) ## , , (Intercept) ## ## (Intercept) x ## Lag 0 1.0000000 -0.8345695 ## Lag 10 0.9647718 -0.8168551 ## Lag 50 0.8683077 -0.7808086 ## Lag 100 0.7777437 -0.7071574 ## Lag 500 0.3318530 -0.3564110 ## ## , , x ## ## (Intercept) x ## Lag 0 -0.8345695 1.0000000 ## Lag 10 -0.8232463 0.8830354 ## Lag 50 -0.7382554 0.6859523 ## Lag 100 -0.6551926 0.6397873 ## Lag 500 -0.2632603 0.2820885
いまは問題が単純だからなんとなく推定できてしまっているけど、もっと高次元になったらヤバイ感じがする[1]ので、早急に原因を解明しておかなければなりません。
いろいろ試してみる
とりあえず事前分布などを変更してみましょう。
prior <- list(R=list(V=1,nu=200))
事前分布はこんな感じのリストを指定します。R の他に G, B が指定できます。R は誤差項の分散の事前分布。G は random effect の分散の事前分布。B は正規分布を仮定するパラメータ（係数）の平均、分散を事前分布を指定するリストです。とりあえず、いまは random effect を指定していないので、G は関係しません。また、ポアソン分布は平均＝分散なので、誤差項の分散の事前分布 R も関係なさそうです。しかしこれを上のように変更してみると「なぜか」結果が劇的に変化します。事前分布は先程のリストを prior 引数に渡します。
g2 <- MCMCglmm(y~x,data=d,family="poisson",prior=prior) autocorr(g2$Sol) ## , , (Intercept) ## ## (Intercept) x ## Lag 0 1.00000000 -0.51934518 ## Lag 10 0.26322261 -0.24273577 ## Lag 50 -0.02917192 0.04835053 ## Lag 100 -0.07810344 0.10765773 ## Lag 500 0.03006962 0.01106049 ## ## , , x ## ## (Intercept) x ## Lag 0 -0.51934518 1.000000000 ## Lag 10 -0.27520432 0.288995514 ## Lag 50 0.01301393 -0.032719865 ## Lag 100 0.02160729 -0.045266254 ## Lag 500 -0.01611708 0.001679549 plot(g2$Sol)
まだ、自己相関は残っていますが、それでも劇的に改善しています。
しかし、事後分布のモードを求めてみると今度はバイアスっぽいものが現れます。明らかに推定精度が落ちていますね。
posterior.mode(g2$Sol) ## (Intercept) x ## 0.7320478 2.1568932
一度ここまでをまとめてみると
- glm ではふつうに解くことができるポアソン回帰モデルを MCMCglmm で解こうとすると、おかしなことがおこる。
- ポアソン回帰では存在しないはずの誤差項の分散を指定するとなぜか、MCMCglmm の挙動が変わり、自己相関が減少。かわりにバイアスが現れる
という感じでしょうか。ちなみに prior に fixed=整数[2] という変数を加えると、分散を V で指定した定数とみなすことができますが、これをやってみても傾向は変わりません。
結局原因はなんなの？
鍵はポアソン回帰では関係なさそうな事前分布 R を変更したら挙動変わった、というところにありそうです。実は
plot(g$VCV)
とするとこの「なぞの誤差項」の事後分布をプロットすることができます。次のグラフは事後分布がデフォルトの場合のものです。
実際、この謎の事後分布はマニュアルを読むとその正体がわかったりします。
実は MCMCglmm で仮定されているモデルは
\begin{align*}y\sim\exp(l)\qquad l:=\beta’x+e\end{align*}というモデルだったりします。\(l\) は潜在変数、e は平均ゼロ、分散 units の正規分布だそうです。これが上の units の正体です。
ポアソン分布が平均＝分散となるのにたいして、このモデルでは平均＜分散となります。つまり、over dispersion（過分散）を暗黙的に仮定していることになります。最初の例では units をゼロと仮定したことになり、over dispersion が全く存在しない、という理想的すぎる状況を作ってしまったため、うまくいかなかった（不安定になった）わけです。
逆にいえば、現実的なデータがほとんど over dispersion であることを思えば、現実的なデータに関してはうまくいく可能性が高くなります。
over dispersion でないデータに対して推定を行う方法は、まだよくわかりません。この units を仮定するモデル化は、サンプリングスキームの中でおそらく本質的な役割を担っているようなのでオプションなどで無効化することができなさそうです。サンプリングスキームはだいたい次のような感じ。
MCMCglmm のサンプリング・スキーム
- メトロポリス法で潜在変数 \(l\) をサンプリング。これはいい感じの提案分布が存在するらしい。
- ギブスサンプラーで fixed effects と random effects（の実現値）をサンプリング。これは単なる正規分布からのサンプリング。
- ギブスサンプラーで units および random effects の分散をサンプリング。逆ウィッシャート分布。
ということで、最初の over dispersion なしのデータセットに対してはどのように推定したらいいのかわからないので、あきらめることにして新しいデータセットでリトライしてみることにしましょう。
新しいデータセット
こんなデータを準備しましょう。
\begin{align*}y\sim Poisson(\exp(2x+1+\sqrt2e))\end{align*} e は標準正規分布に従う確率変数です。
x <- rnorm(200) e <- rnorm(200) y <- rpois(200,exp(2*x+1+sqrt(2)*e)) d <- data.frame(y=y,x=x) m <- MCMCglmm(y~x,data=d,family="poisson")
こんどはうまくいっているはずです。プロットしてみましょう。
よさそうですね。units が正しい分散の値＝２の周辺でサンプリングされていることがわかります。最初の例で見られたような自己相関も見られません。念のため、事後分布と自己相関の数字を出しておきましょう。
posterior.mode(m$Sol) ## (Intercept) x ## 1.011231 1.960540 ## -> OK! posterior.mode(m$VCV) ## units ## 1.922535 ## -> OK! autocorr(m$Sol) ## , , (Intercept) ## ## (Intercept) x ## Lag 0 1.000000000 -0.38218637 ## Lag 10 0.174746578 -0.23397818 ## Lag 50 0.030926348 -0.01045872 ## Lag 100 0.008258976 -0.03562686 ## Lag 500 -0.034086677 0.04925067 ## ## , , x ## ## (Intercept) x ## Lag 0 -0.382186365 1.000000000 ## Lag 10 -0.216486193 0.252820118 ## Lag 50 -0.056063463 -0.039203250 ## Lag 100 -0.000795888 -0.005915734 ## Lag 500 0.011238684 -0.021532814 autocorr(m$VCV) ## , , units ## ## units ## Lag 0 1.000000000 ## Lag 10 0.225311785 ## Lag 50 -0.039435254 ## Lag 100 -0.008188887 ## Lag 500 0.048201071
これでも自己相関が高すぎる！と思う場合は thin 引数を大きな値に変更しましょう。
まとめ
というわけで本質的な解決はしていないような気がしますが、一応、個人的には納得したのでこのあたりで終わりにします。
- MCMCglmm を使ってポアソン回帰をしてみた
- MCMCglmm はポアソン分布などでも線形予測子の後に正規誤差を仮定する。これは現実的なデータに対して over dispersion を回避するという意味では有用かもしれないが、テスト問題や一部の non-over dispersion なデータでは収束性が悪くなるかもしれない（もしかすると収束性を改善するワザがあるかもしれない）。
- すこし書き方が冗長だったかもしれない。最後まで読んでくださった方、ありがとうござます。
- 実は MCMCglmm で罰則付きスプライン回帰なども簡単に実装できるので、そのあたりのこともいつか書いてみたいとおもいます。
