負の二項分布モデルの精度検証:glmer.nb、glmmamb、JAGS

精度はいかほどですか

この前紹介したけれども、lme4パッケージの更新で、ランダム効果入りの負の二項分布モデルが扱えるようになっていた。しかし、他の関数(glmmadmbやJAGS)と比べてどの程度の精度なのか、気になるところ。そこで、簡単に推定の精度比較をしてみる。モデルは極めて単純な線形モデル。

Y ~ NB(pi, n)
pi = n/(mui + n)
log(mui) = a + bXi + εj
εj ~ Normal(0,σR2)

まず、上記を真のモデルとするシミュレーションデータを生成し、それをglmer.nbなどに推定させる、という作業を100回繰り返す。100回分の推定値の頻度分布と真の値を比べるということをした。パラメーターの真の値はそれぞれa=1.0, b=0.5, σR = 0.7とした。

結果

赤い曲線が100回分の推定値の頻度分布。縦の赤い実線は100回分の推定値の中央値、グレーの実線は真の値。左から切片(a)、回帰係数(b)、ランダム効果のばらつき(σR)。

glmer.nb

glmmadmb

JAGS

こうしてみると、推定精度の観点から一番難がありそうなのはglmmadmb。特にランダム効果のばらつきの推定値は、右側の裾が随分と重い。これはランダム効果が二つ以上になったら絶望的ではなかろうか。glmer.nbは比較的安定していたが、回帰係数の推定値のピークが平らになっているのがやや不安。中央値は真の値に近いけれども。JAGSは可もなく不可もなく、推定としては一番マシかな、という気がする。ランダム効果が増えたときが気になるけれども(たぶんJAGSの推定がかなりまマシなものになる)、面倒なのでやらない。でも、最尤推定の負の二項分布モデルで、ランダム効果の入れ過ぎはよくないということは間違いなさそう。確認してないけど。

シミュレーションデータ生成に使ったコード

glmer.nbの精度検証に使ったコード。そのほかは割愛。
apply系を使えばfor構文が要らなくなる気がするけど。。。

for(i in 1:100){
  #simulate 100 data points from 10 groups with one continuous fixed effect variable
  x<-runif(100,0,1)
  f<-gl(n = 10,k = 10)
  data<-data.frame(x=x,f=f)
  modmat<-model.matrix(~x,data)
  
  #the fixed effect coefficient
  p <- NULL
  p[1] <- 1; p[2] <- 0.5
  fixed<-c(p[1],p[2])
  
  #the random effect
  p[3] <- 0.7
  rnd<-rnorm(10,0,p[3])
  
  #the simulated response values
  size <- 3
  mu <- exp(modmat%*%fixed+rnd[f])
  prob <- size/(mu+size)
  data$y<-rnbinom(100,size=size,prob=prob)

  fit <- glmer.nb(y~x+(1|f),data)
  est[[1]][i] <- fit@beta[1]
  est[[2]][i] <- fit@beta[2]
  est[[3]][i] <- fit@theta
}

SPONSOR LINK

スーパー堤防

父の実家が陸前高田なので、節目で被災地に訪れている(特にボランティアというわけではないです)。そのときに聞いた興味深い話。

被災以降、津波の「減災」を目的としたスーパー堤防がガンガン作られており、それはもうすごい勢い。生態学の観点からそれが好ましくないことはもちろんそうなのだが、現地に生きる人の感情を考えると極めて難しい問題だなぁ、と感じていた。しかし、それは大きな勘違い(みんな堤防ほしいだろう、という先入観)で、現地からでる「堤防」に対する意見は住民の間でも大きく割れており、そしてそのパターンがどうにも奇妙だった。

奇妙だったというのは、海の近くに住んでいる住人(つまり被災リスクが高い)ほど堤防反対の意見をもつことが多かったこと。逆に、海から離れた場所に住む人(被災リスクは相対的に低いと思われる)ほど堤防建設を熱望しているという。

海際の人は年配の人が多く住んでいる可能性あるので、年齢層との偽相関は拭えないが、海際の住人の共通した意見として、海に住んでいるのだからそういったリスクは承知の上、海を見れないことのほうがよっぽどつらいし、仕事もしにくい、とのこと(もちろん、これが一般的な傾向というつもりはございません。あくまで聞いた範疇です)。

まったく定量的な調査ではないですが、現地の人に聞く限り、上述のような傾向があるそうです。だれかきちんと証明しているのかな、これ。

SPONSOR LINK

さまざまな情報量基準

生態学の界隈では、情報量基準なるものが大人気です。そのくせ、(筆者も含め)それぞれの情報量基準なるものがどんな意味をもつのか、よくわかっていない場合がほとんどです。わかっている範疇で整理。

AIC(赤池情報量基準)

情報量基準といったらAICと言わんばかりですが、なんだったら決定係数との違いすらよく分かっていないかったのが博士課程のころの僕です。AICが低いモデルでは、「同じ母集団から同じサンプル数の別のデータセットを取り出してきたときにも、それなりのパフォーマンス(当てはまり)を発揮する」ということになる。だが、これはサンプル数を無限大に増やしたとしても、必ずしも「真のモデル」には収束しないらしい。あくまで評価基準が予測精度だからだろうか。

AICc

サンプルサイズが小さいときはAICc!みたいな使われ方をしている。だが実は、正規分布を仮定した補正をしているので、正規分布以外のモデルに対して適用するのはご法度らしい。Burnhamさんの本を参照のこと。

BIC(ベイズ情報量基準)

AICと異なり、こちらは真のモデルを選ぶらしい(あくまでサンプルが無限大になったとき、ということだろうか)。あまり使ったことはないが、たしかにAICとは全く異なるモデルを選ぶ傾向がある。予測精度ではなく、真に重要な変数を選び出したい場合にはBICを使った方がよいのだろう。式的には、AICの罰則項にサンプルサイズが入っていないのに対し、BICには入っている。そのため、サンプルサイズが大きくなったときに、前者は対数尤度の項の影響が相対的に大きくなり、より複雑なモデルを選びがちになるようだ。

DIC(Deviance Information Criteria)

こちらはベイズモデルなどでよく使われる情報量基準。だが、広く使われている割に、その数学的根拠は随分と希薄らしい(論文でこういった表記をよく見かける)。DICの式の構造としては、モデルの当てはまりに関わる指標(Deviance)に対し、パラメーター数で罰則を与えている。AICでは2k(kはパラメーター数)という罰則項があるが、これには数学的根拠がある。DICの場合は、そういった裏付けがない、ということだろうか。時間のあるときに調べよう。
また、DIC利用の際の注意点として、パラメーターの正規性が担保された状態でしか役に立たないらしい。詳細はよくわからないが、ベイズのご利益がでる階層線形モデルや状態空間モデルでは、この「パラメーターの正規性」がほとんど担保されないらしい。パラメーターの正規性とはなんぞや。。。(おそらく、パラメーターの事後分布が正規分布で近似できる、ということか)。だとしたらDICを使える場面などそうそうないのでは。

wAIC(widely applicable information criteria)

最近考案されたらしいwAICなるものは、階層をもった複雑なモデルにおいても適用できるらしい。これまたAICに引き続き日本人が発見したので驚き。日本人って統計につおいのだろうか。開発者のワタナベさんのページはこちら(ワタナベさんが開発したからwAICかと思った。違うのかorz)。内容をまったく理解していないのですが、浅はかな理解では予測精度を問題にしているような気がするので、モデル選択の結果の解釈としてはAICと同様になるような気がする。階層モデルのようなものにも適用できる点が違いか(あくまでユーザー視点からすれば)。もう少し勉強する。

lme4パッケージの更新

しばらく見ないうちに、最新版のlme4パッケージでは劇的な進化を遂げていたようなので、それについてメモ。

glmer.nbの出現

これまで負の二項分布を仮定したモデルで、ランダム効果をいれるとしたら、glmmADMBパッケージか、もしくはベイズモデル化(WinBUGSなど)する必要があった(ADMBもベイズみたいなもん?中身はよくわからない)。しかし、最新版のlme4パッケージでは、glmer.nbという関数が表れたことにより、負の二項分布でもランダム効果を簡単に入れることができるようになった。仕様については基本的なlmerと変わらない模様。ただし、頻度主義的に扱うには尤度関数が複雑すぎる気がするので、いったいどこまで結果が信用できるかは不明。>3のランダム項をいれるのはきっとご法度。

mcmcsamp()の廃止

昔のlme4には、最尤法ではなく、MCMCでベイズっぽく推定するための関数が用意されていた。しかし、最新版のlme4パッケージでは除外された模様。Rのhelpには以下のように書いてある。
Due to difficulty in constructing a version of mcmcsamp that was reliable even in cases where the estimated random effect variances were near zero (e.g. https://stat.ethz.ch/pipermail/r-sig-mixed-models/2009q4/003115.html), mcmcsamp has been withdrawn (or more precisely, not updated to work with lme4 versions >=1.0.0)

どうやらランダム効果のばらつきがゼロ付近のときに、あまりいい推定をしてくれないということで、やめにしたということらしい。個人的には、ランダム効果が多いとき(>3)の最尤推定は不安なので、mcmcsamp()の取り下げは手痛い、、、JAGSとかはコード書くのが面倒だからなぁ、、、

lmerにおけるREMLのデフォルト化

REMLとはなんぞや、と調べてみると、REstricted Maximum Likelihoodの略らしい。最新版のlmerでは、MLではなくREMLがデフォルトになっている。REMLとMLの違いは、(筆者の浅はかな理解によると)、前者は固定効果による自由度減少を考慮したうえでランダム効果の分散を推定し、後者は考慮しないらしい。つまり、後者の場合、ランダム効果のばらつきの推定は偏りが生じているということになる(不偏推定量でない?)。ただ、REMLではモデル選択ができないため(なぜかはよくわからない)、MuMInパッケージなどでモデル選択する場合はREML=Fとする必要がある。正規分布以外ではMLしか使えないので心配の必要なし。

glmerとlmerの(よくわからない)分離

正規分布の場合はlmer、それ以外はglmerということでよさそう。

QGIS: 属性テーブルの一括編集

QGISで属性テーブルの値を一括編集したいときの処理

1.変更対象の地物を選択する
2.属性テーブルからフィールド演算をひらく
3.既存カラムを更新を選択
4.変更したいカラムを選択
5.演算スペースに一括入力したい文字や数値を入力。
※文字の場合、シングルコーテーション(’ ‘)で囲む
※選択フィーチャーだけに適用のチェックボックスにチェック

ATMの悪夢

住み慣れた国を離れると、これまでの常識が通じないことも多々あり、いろいろと苦労しております。2017年の4月からアメリカに移動し、所属がUniversity of Minnesotaとなったのですが、さっそく新生活の立ち上げでひやり(というか実際損した)とした経験です。

日本ではお金を振り込んだり、預け入れをする際に、ATMの機械がその場で紙幣を数えてくれるのが普通だと思います。しかし、アメリカではそうではありませんでした。Depositのボタンを押すと、なにやら封筒がいるか、と聞いてくるのです。こちらの常識で考えると、当然封筒はいらないので、いらないと押し、そのまま紙幣を挿入口へ。ただウイーンといいながら吸い込まれていく紙幣。まったく数える気配がありません。よくわからないので、そのままその場を後にしたのですが、そもそも紙幣を数えていなかったこと、そして挿入金額が自己申告(途中で入力させられる)だったことに違和感を覚えるべきでした。不安になり、しばらく経ってからググってみると、、、げっ、、、封筒に金額と名前を記入して入れるだと、、、?オンラインアカウントの口座残高をみると、入金されていなことになっている(ちなみに500ドル)。なんとアメリカでは、ATMの預け入れは夜間金庫みたいなものになっており、行員があとでATMからお金を運び出し、手作業で数えているとのこと。なんじゃそりゃ。最近では日本のようなATMが大多数のようですが、運悪く私の使ったATMは旧式のもので、従来の方式が採用されていた模様。

そんな馬鹿な、と思い、支店にいって事情を話すと、行員の方が丁寧に対応してくれました。しかし時すでに遅し、紙幣をそのまま挿入したためか、500ドル入金したはずが480ドルになってしまいました。とほほ。全部なくなるよりはマシか。。。