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

精度はいかほどですか

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

SPONSOR LINK

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
}
Spread the love

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です