累積ロジットとGLM二項分布の比較・再＆続

（うっかり、同様の検証記事を消してしまったので、ついでにアップデートします）
前回に引き続き、段階的なカテゴリーデータのモデリングに用いられる累積ロジット（cumulative logit）、
例えば、悪い、ふつう、よい、のようなデータを関連しそうな要因で解析するような場合に用いる。

しかし、悪い＝0、ふつう＝1、よい＝2、のように数値化してしまえば、二項分布型のGLMではダメなのだろうか？たぶん、間隔のいびつなデータでは累積ロジットの方が適しているのだろうが。

テストデータを用いて、両者の推定と推定値の求め方を比較する。

# まず解析用データの設定
logistic <- function(xx) 1 / (1 + exp(-xx))
N <- 3 # 最大3、つまり 0, 1, 2, 3 の値を取りうる
X <- rep(c(1:10), each=10)
Y1 <- rbinom(100, N, prob=logistic(-5 + 0.8*X)) # logisticの中身を基にして、二項乱数を発生
Y2 <- factor(Y1, ordered=T) # 累積ロジット用に、ランク化した応答変数を用意
D <- data.frame(X, Y1, Y2)

# 解析の実行

require(ordinal)
require(VGAM)
M1 <- glm(cbind(Y1, N-Y1) ~ X, family=binomial, data=D) # 二項分布のGLM
M2 <- clm(Y2 ~ X, data=D) # 累積ロジット
M3 <- vglm(Y2 ~ X, family=cumulative(parallel=T), data=D)
 # 比較用にvglm版、parallel=Tで切片のみ複数になる、Fにすると回帰係数もランクごとに推定（切片のみランク毎の場合、比例オッズモデルともいうようだ）

### GLM

summary(M1)

# Coefficients:

#             Estimate Std. Error z value Pr(>|z|)    

# (Intercept)  -5.9324     0.6644  -8.929   <2e-16 ***

# X             0.9241     0.1011   9.141   <2e-16 ***

# ---

# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#

# (Dispersion parameter for binomial family taken to be 1)

#

#     Null deviance: 287.723  on 99  degrees of freedom

# Residual deviance:  84.861  on 98  degrees of freedom

# AIC: 139.14

### 累積ロジット（clm）

summary(M2)

#  link  threshold nobs logLik AIC    niter max.grad cond.H 

#  logit flexible  100  -67.06 142.12 5(0)  2.20e-08 2.6e+03

#

# Coefficients:

#   Estimate Std. Error z value Pr(>|z|)    

# X   1.2295     0.1659    7.41 1.26e-13 ***

# ---

# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#

# Threshold coefficients:

#     Estimate Std. Error z value

# 0|1    5.955      0.913   6.522

# 1|2    8.082      1.130   7.153

# 2|3    9.716      1.283   7.575

### 累積ロジット（vglm）# 回帰係数の正負が逆になる

summary(M3)

# Coefficients:

#              Estimate Std. Error z value

# (Intercept):1   5.9548    0.92423  6.4430

# (Intercept):2   8.0821    1.14814  7.0393

# (Intercept):3   9.7158    1.31086  7.4117

# X                  -1.2295    0.16832 -7.3046 

#

# Number of linear predictors:  3 

#

# Names of linear predictors: logit(P[Y<=1]), logit(P[Y<=2]), logit(P[Y<=3])

#

# Dispersion Parameter for cumulative family:   1

#

# Residual deviance: 134.1208 on 296 degrees of freedom # 過小分散気味、clmには無い情報！

#

# Log-likelihood: -67.06038 on 296 degrees of freedom

### 推定値を得るには少し手間がかかる

pre1 <- round(3*fitted(M1)) # GLM：試行回数*確率、を整数値に丸める

pre2 <- as.numeric(predict(M2,type="class")$fit) - 1 

# 累積ロジット（clm）：predict関数で推定ランクを求め、

# ラベルに合うように 1 を引く（ランクは1,2,3,4、元の値は0,1,2,3なので）

pre3 <- apply(fitted(M3), 1, which.max) - 1

# 累積ロジット（vglm）：何番目のランクの確率が最大かを求め、

# ラベルに合うように 1 を引く

cbind(X, Y1, pre1, pre2, pre3) 

# 推定は完全ではないが、GLMと累積ロジットの推定結果は近いものになった。

     X Y1 pre1 pre2 pre3

1    1  0    0    0    0

2    1  0    0    0    0

3    1  0    0    0    0

4    1  0    0    0    0

5    1  0    0    0    0

（中略）

51   6  1    1    1    1

52   6  1    1    1    1

53   6  2    1    1    1

54   6  2    1    1    1

55   6  2    1    1    1

56   6  1    1    1    1

57   6  2    1    1    1

58   6  1    1    1    1

59   6  1    1    1    1

60   6  0    1    1    1

61   7  3    2    2    2

62   7  1    2    2    2

63   7  1    2    2    2

64   7  3    2    2    2

65   7  3    2    2    2

66   7  1    2    2    2

67   7  3    2    2    2

68   7  2    2    2    2

69   7  1    2    2    2

70   7  3    2    2    2

71   8  3    2    3    3

72   8  2    2    3    3

73   8  3    2    3    3

74   8  3    2    3    3

75   8  2    2    3    3

76   8  1    2    3    3

77   8  1    2    3    3

78   8  3    2    3    3

79   8  3    2    3    3

80   8  3    2    3    3

81   9  2    3    3    3

82   9  1    3    3    3

83   9  3    3    3    3

（後略）

### では、推定確率の曲線を描くにはどうしたらよいか？

# これがけっこう面倒くさい。まずはvglmを用いてチェックする。

# vglmをfitted()すると、ランク毎の確率が出てくる

head(fitted(M3)) # head()で頭出しする

                       0                     1                      2                  3

1   0.991209875 0.007734525 0.0008493615 0.0002062383

2   0.991209875 0.007734525 0.0008493615 0.0002062383

3   0.991209875 0.007734525 0.0008493615 0.0002062383

4   0.991209875 0.007734525 0.0008493615 0.0002062383

5   0.991209875 0.007734525 0.0008493615 0.0002062383

6   0.991209875 0.007734525 0.0008493615 0.0002062383

# clmでは、
M2@fitted.values で同様にランク毎の確率が得られる。

# まず、この出来合いの推定値を図示してみる

plot(X, fitted(M3)[,1], xlim=c(0,10), ylim=c(0,1), col="lightblue", ylab="probability") # 0

points(X, fitted(M3)[,2], col="turquoise") # 1

points(X, fitted(M3)[,3], col="royalblue") # 2

points(X, fitted(M3)[,4], col="darkblue") # 3

# モデルの構造を考えると推定確率曲線は次の計算でいいはず

#（coef(M3)[4]は回帰係数だが、正負が逆なので - を付ける）

# 累積ロジットの名前通り、累積で表されている

curve(1 - logistic(-coef(M3)[4]*x - coef(M3)[1]), 
  add=T, lwd=2, col="lightblue") # ランク0: 1 - ランク1の累積確率

curve(logistic(-coef(M3)[4]*x - coef(M3)[1]) - logistic(-coef(M3)[4]*x - coef(M3)[2]), 

add=T, lwd=2, col="turquoise") # ランク1: ランク1の累積確率 - ランク2の累積確率

curve(logistic(-coef(M3)[4]*x - coef(M3)[2]) - logistic(-coef(M3)[4]*x - coef(M3)[3]),

add=T, lwd=2, col="royalblue") # ランク2: ランク2の累積確率 - ランク3の確率

curve(logistic(-coef(M3)[4]*x - coef(M3)[3]), 
  add=T, lwd=2, col="darkblue") # ランク3: ランク3の確率（もはや累積でない）

curve(logistic(coef(M1)[1] + coef(M1)[2]*x), add=T, col="tomato", lwd=2) # 比較用にGLMも

# 色が薄い〜濃いにかけて、それぞれ、0, 1, 2, 3 になる確率、赤はGLMの確率

# ちゃんと関数による推定プロットと、自前で作った推定曲線が一致した。計算の仕方はこれでよさそうだ。

# ちなみにclmの場合はこう計算する（回帰係数の前の - が不要）

curve(1 - logistic(coef(M2)[4]*x - coef(M2)[1]), add=T, lwd=2, col="lightblue") # 0

curve(logistic(coef(M2)[4]*x - coef(M2)[1]) - logistic(coef(M2)[4]*x - coef(M2)[2]), 

add=T, lwd=2, col="turquoise") # 1

curve(logistic(coef(M2)[4]*x - coef(M2)[2]) - logistic(coef(M2)[4]*x - coef(M2)[3]), 

add=T, lwd=2, col="royalblue") # 2

curve(logistic(coef(M2)[4]*x - coef(M2)[3]), add=T, lwd=2, col="darkblue") # 3

### 次に、推定値の曲線を図示してみる（上のは推定確率でした）

# 比較対照用にGLMの曲線と比較する

# 元データ

plot(Y1 ~ X, data=D, xlim=c(0,10), ylim=c(0,3))

 # GLM

curve(3*logistic(coef(M1)[1] + coef(M1)[2]*x), add=T, col="red", lwd=2)

# 累積ロジット（全部足し合わせるような累積構造になる）

curve(
   1*(1 - logistic(coef(M2)[4]*x - coef(M2)[1])) # ランク0

+ 2*(logistic(coef(M2)[4]*x - coef(M2)[1]) - logistic(coef(M2)[4]*x - coef(M2)[2])) # ランク1

+ 3*(logistic(coef(M2)[4]*x - coef(M2)[2]) - logistic(coef(M2)[4]*x - coef(M2)[3])) # ランク2

+ 4*(logistic(coef(M2)[4]*x - coef(M2)[3])) # ランク3

 - 1, # 0, 1, 2, 3なので、1,2,3,4から1を引く

add=T, col="blue", lwd=2)

# 赤：GLM、青：累積ロジット。この単純なケースではほとんど推定値は変わらない。

#### ついでにベイズ版でも計算

# 下準備、データを累積に変換。1以上、2以上、3以上にまとめる
rank <- matrix(0, nrow=100, ncol=3) # 3はランクの数-1、100はいわゆるN数
for (r in 1:3) rank[Y1 >= r, r] <- 1 
# GLM用の用意していた数値データY1を利用、ランク1~3以上の場合に各列に1を入れる

model <- function() {
  for (j in 1:3) { # ランクの数-1
  for (i in 1:100) { # いわゆるN数
  rank[i,j] ~ dbern(p[i,j]) # それぞれベルヌーイ分布
 logit(p[i,j]) <- alpha[j] + beta*X[i] # 切片だけランク毎になってる
 }
  alpha[j] ~ dnorm(0, 1E-6) # 切片をランク毎に推定
  }
  beta ~ dnorm(0, 1E-6) # 回帰係数は共通
}

data <- list(X=X, rank=rank)
parm <- list(beta=0, alpha=rnorm(3))
source("WBUGS.R") # 自作ラッパー関数
out <- wbugs(data, parm, model)

# 3 chains, each with 11000 iterations (first 1000 discarded), n.thin = 10
# n.sims = 3000 iterations saved
#                 mean    sd    2.5%     25%     50%     75%   97.5%  Rhat n.eff
# beta         1.328 0.166   1.022   1.217   1.318   1.434   1.680 1.001  2200
# alpha[1]  -6.377 0.889  -8.217  -6.960  -6.327  -5.770  -4.703 1.002  1100
# alpha[2]  -8.645 1.142 -11.040  -9.370  -8.563  -7.851  -6.528 1.002  1800
# alpha[3] -10.518 1.336 -13.331 -11.373 -10.420  -9.585  -8.139 1.001  3000
# deviance 154.787 2.934 151.200 152.600 154.050 156.200 162.000 1.002  1500
# DIC info (using the rule, pD = Dbar-Dhat)
# pD = 4.0 and DIC = 158.8

# BUGS版では、切片が正負が逆になった

# きちんとコードを書いた結果がこれなので、符号が逆なのはvglmやclmなのだが…じつにややこしい。

累積ロジットの汎用Rパッケージ {ordinal}

累積ロジット（cumulative logit model）を使う際に、今ひとつ使い勝手のいいRパッケージがないのが気になっていた。
（cf. 累積ロジット：よい、ふつう、わるい、のような段階的な現象についての推定に用いる。等間隔に近ければ二項分布のGLMで構わないだろうけれど、そうでない場合にはこれが適しているようだ）

例えば、vglmではランダム効果が入れられない。

mixcatでは、ランダム効果が入れられるが、対数尤度までは出力できても、AICなど情報量規準は出してくれない（手計算すれば良い話ではあるが…）。

最近、ordinalというパッケージを見つけた。１つのパッケージでGLM版、GLMM版の両方が含まれているし、使用法もglm()やglmer()と同じようにしてくれていて使いやすい。

ちなみに、GLM版にstepAICは使えたが、dredgeはダメだった。
GLMM版では、モデル選択関数はdrop1が使用できた。

require(ordinal)
example(ordinal) # パッケージで用意されている例を出力

# まずはGLM版
ordinl> ## A simple cumulative link model:
ordinl> fm1 <- clm(rating ~ contact + temp, data=wine)

ordinl> summary(fm1)
formula: rating ~ contact + temp
data:    wine

 link  threshold nobs logLik AIC    niter max.grad cond.H
 logit flexible  72   -86.49 184.98 6(0)  4.01e-12 2.7e+01

Coefficients:
           Estimate Std. Error z value Pr(>|z|)    # 回帰係数部分（この例ではカテゴリカルだが）
contactyes   1.5278     0.4766   3.205  0.00135 **
tempwarm     2.5031     0.5287   4.735 2.19e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Threshold coefficients:           # 各ランクの切片
    Estimate Std. Error z value
1|2  -1.3444     0.5171  -2.600
2|3   1.2508     0.4379   2.857
3|4   3.4669     0.5978   5.800
4|5   5.0064     0.7309   6.850


# 比較対照用にvglmでの結果
require(VGAM)
vm1 <- vglm(rating ~ contact + temp, family=cumulative(parallel=T), data=wine)

Coefficients:
Estimate Std. Error z value
(Intercept):1 -1.3444 0.50850 -2.6438
(Intercept):2 1.2508 0.43908 2.8487
(Intercept):3 3.4669 0.59711 5.8061
(Intercept):4 5.0064 0.72906 6.8669
contactyes -1.5278 0.47362 -3.2258 # 回帰係数部分の正負がclmと逆なことに注意
tempwarm -2.5031 0.53199 -4.7052

# ordinalに戻って、こちらはGLMM版
ordinl> ## A simple cumulative link mixed model:
ordinl> fmm1 <- clmm(rating ~ contact + temp + (1|judge), data=wine) # glmer()同様の構造

ordinl> summary(fmm1)
Cumulative Link Mixed Model fitted with the Laplace approximation

formula: rating ~ contact + temp + (1 | judge)
data:    wine

 link  threshold nobs logLik AIC    niter    max.grad cond.H # ちゃんとAICも算出する
 logit flexible  72   -81.57 177.13 331(996) 1.05e-05 2.8e+01

Random effects:                                     # ランダム効果
 Groups Name        Variance Std.Dev.
 judge  (Intercept) 1.279    1.131
Number of groups:  judge 9

Coefficients:
           Estimate Std. Error z value Pr(>|z|)  
contactyes   1.8349     0.5125   3.580 0.000344 ***
tempwarm     3.0630     0.5954   5.145 2.68e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Threshold coefficients:
    Estimate Std. Error z value
1|2  -1.6237     0.6824  -2.379
2|3   1.5134     0.6038   2.507
3|4   4.2285     0.8090   5.227
4|5   6.0888     0.9725   6.261

（次の記事に続く）

生態学会2014広島大会で発表します

広島の生態学会に来ています。発表のネタの準備中に大きなミスを見つけて取り返したりで時間に追い詰められグッタリしてますorz
発表は明後日17日のポスター、国内温帯域の藻場の衰退とサンゴの分布拡大の話です。温暖化影響や生物多様性データベースのGIS化などに興味のある方には興味の持てる内容かと思います。来聴歓迎します。
http://www.esj.ne.jp/meeting/abst/61/PB3-083.html

（旧サイトより移行）過分散データ：GLM負の二項分布、GLMMによる推定をAICでモデル選択することは可能か？

# 2014.02.20追記：下記のlmer()は現在はglmer()に相当します

かつて、数値実験をしてみて「できないようだ」と書いたのですが、GLMM側（lmerやglmmML）のAIC計算方法をGLM側（glmやglm.nb）へと合わせるように修正すれば可能であることを確かめました。じつは、両者でAICの計算方法が異なっていたことが分かりました。道理でいつもGLMMのAICが小さ過ぎるわけです…。

★以前アップしていた数値実験を大幅に修正・加筆しておきました（cf.「過分散データ：GLM負の二項分布、GLMMによる解析の比較」）

なお、glmやglm.nb と、lmerやglmmMLのAICの計算方法の違いを検証するところまでは、検索すると複数のブログで見つかります。

ところがさらに検証してみたところ、逆にglmやglm.nbのAICの方をlmerやglmmMLのAIC計算方法へと修正すると、不当にglm、とくにglm.nbのAICが小さくなり、誤ったモデル選択となる確率が大幅に増大することも分かりました（ほぼ逆転してしまう！）。もっとも、ちゃんと数値実験をやっているとはいえ、なぜそうなるのか理論的裏付けがよく分からないので、なにか思いついたら検証したいと思います。

ちなみに、まったくのポアソン分布データのパラメータ推定でも、glm、glm.nb、lmerはいずれも真の値をほぼ推定できていました。過分散問題がけっこう厄介なことを思うと、もはや glm(, family=poisson) を紹介する意義はほとんど無いのかもしれません…。



数値実験では、乱数発生させた、ポアソン分布データと、負の二項分布データ、ポアソン分布に正規乱数ノイズを加えたデータの三者で、GLM（ポアソン分布：glm）とGLM（負の二項分布：glm.nb）、GLMM（ポアソン分布：lmer）が真の値を推定できているかをクロスチェックし、同時にAICが正しい推定をちゃんと反映できているかをもチェックしました。lmerの方は本来のGLMMの使い方とは異なり、各データを異なるid（1データ1id）としているので過分散の対処に用途は限定しています（こういう使い方は、help(lmer)や、Warton & Hui (2011) Ecology 92:3–10、和文では粕谷先生のGLM本「一般化線形モデル」で紹介されています）。

混乱しやすいのでメモ的に整理します。
まず、AIC()を用いた時の、glmやglm.nbと、lmerやglmmMLの計算方法の違い：
glmやglm.nbでは、
AIC = deviance + 2*推定するパラメータ数
lmerやglmmMLでは、
AIC = residual deviance + 2*推定するパラメータ数

また、deviance()は、いずれの場合もresidual devianceを計算する。
一方、logLik()では、glmやglm.nbではdeviance/(-2)、lmerやglmmMLではresidual deviance/(-2)


lmerやglmmMLのAICをglmの方へ合わせるには、下記のようにする（上述の理由により、逆はやめておきましょう）。

# 例えば、こんな２つのモデルがあるとき（id=1~N数）
pois <- glm(Y ~ X, family=poisson)
pois.m <- lmer(Y ~ X + (1|id), family=poisson)

glmのAICはふつうに、AIC(pois)、で求められる。
GLMMのAICは、上述の違いに基づくと、こうすればglmのと同じように求められる：
AIC(pois.m) - deviance(pois) -2*logLik(pois) 
# (residual deviance + 2*パラメータ数) - residual deviance + deviance = deviance + 2*パラメータ数
#（ランダム効果を抜く以外は同等であるglmモデルを用意してやる必要がある）

# ちなみに、AICがそのように変換できることの検証として、GLMMでidを全部 1 にしてやると（1グループのみ）、グループ間の分散は 0 と推定されます。切片と回帰係数の推定も、単なるGLMのそれと同一の推定結果となります。ただし、分散値は推定するパラメータ数としてカウントされるので、変換後のGLMMのAICはGLMのAICよりもキッカリ2だけ増加した値となります。これは、AICの計算式にある、2*パラメータ数（ここでは、2*1）、の部分に相当します。

しかし、なぜGLMMのパッケージがresidual devianceに基づいたAIC計算にしているのか、その意図が何なのかわかりません。こういう計算方法の修正をしないでくれと言っているようにも見えます。さしあたり、過分散データへの対処という限定的な使用法について数値実験をしてみた限りでは問題は無いようです。

こういう煩わしさを思うと、もはや個人的にはBUGSに全面移行してしまおうかなどと考えたりします…Stanの発展も待ち遠しいです。


# （当時いただいたコメント）
2012/11/15(木) 23:39:11 | URL | 高橋
興味深く読ませていただきました。
ちょっと気になったのですが、
＞# (residual deviance + 2*パラメータ数) - residual deviance + deviance = deviance + 2*パラメータ数
この式の２つのresidual devianceは別物だと思います。
residual deviance = (当該モデルのdeviance) - (飽和モデルのdeviance)
なので、
＞- deviance(pois) -2*logLik(pois)
は、
- ((ポアソンモデルのdeviance) - (ポアソンモデルの飽和モデルのdeviance)) + (ポアソンモデルのdeviance)
となります。
したがって、
＞AIC(pois.m) - deviance(pois) -2*logLik(pois)
は、
(混合モデルのdeviance) - (混合モデルの飽和モデルのdeviance) + 2*パラメータ数 + (ポアソンモデルの飽和モデルのdeviance)
となり、「混合モデルの飽和モデル」が「ポアソンモデルの飽和モデル」と等しい場合（多分そうなのだと思いますが）、混合モデルのdevianceに基づいてAICを求めているということになると思います。
glm.nb()の関数定義を見ると、推定されたThetaのもとで飽和モデルのdevianceを計算しているようです。この値は、ポアソンモデルにおける飽和モデルのdevianceよりずっと大きくなるので、residual devianceを用いてAICを計算すると、負の二項分布の常勝となるのでしょう。
ちなみに飽和モデルの対数尤度は、応答変数をyとすれば、
sum(dpois(y, y, log = TRUE))
で求められます。負の二項分布の場合は、glm.nb()が返すThetaパラメータを使って、
sum(dnbinom(y, mu = y, size = Theta, log = TRUE))

> 高橋さま 

コメントありがとうございます、おかげで理解が深まりました。この検証に興味を持っていただき嬉しいです。 

ポアソンモデルと、ポアソン混合モデルの飽和モデルが等しい時にのみ成り立ちうるというご指摘、たしかにその通りでした。 

これはもう数値実験を見る限りそうなるとしか答えようがないのですが、1グループのみの混合モデルの変換後AICから2を引いた値とポアソンモデルのAICが一致する、というのが一応の証拠になっていると考えています（こちらの検証コード掲載は省いてしまっていますが）。 

glm.nbの飽和モデルの求め方がそうなっているとは気づいていませんでした。言われてみて改めて関数定義をチェックするとたしかにそれらしきコードを見つけることができました。 

混合モデルならば、グループを1つだけにするなどにより、ポアソンモデルとのすり合わせを試みることができましたが、負の二項モデルで同じようなすり合わせをするとなると…glm関数でThetaを1に固定してみるか…また追々試してみたいと思います。

（旧サイトより移行）連続量を単位あたりに直してGLMする場合の対処：gaussian("log")、Gamma("log")、offsetなどなど）

# 2014.02.20 旧サイトを閉じるため、移植しました。
# なお、下記のようなケースではゴンペルツ曲線などを用いた成長モデルを使った方がいいかなと最近思ってます。

最近、実験系の統計を引き受けていて気になったので、連続量を単位あたりでGLM推定する方法について検証してみました。

単位あたりの量をGLMで解析する時、個体数や頻度などの整数の場合の対処はoffset項の利用で解決するのが常套手段ですね（密度や○×率のような割り算した値ではなく、元の値はいじらずに単位量を係数1として説明変数に加える。当ブログでもかつて紹介：http://nhkuma269.blog77.fc2.com/blog-entry-9.html）。
例えば、同じく個体密度0.1でも、1/10と100/1000とでは意味合いが違うが、割り算すると両者は同一密度として扱われてしまう。前者では1個体の増減が大きな誤差を生むが、後者ではほとんど影響なし。

では、元々が連続量のものを単位あたりの量にする場合も同様にoffset項による対処がよいのだろうか？例えば、サンゴの枝あたりのクロロフィル量とか、成長量（実験後重量 / 実験前重量）のような量は、必ず単位量（この場合は枝）のバラツキによる誤差が出てしまう（なるべく条件は揃えるように努力はしているだろうが）。これを割り算してしまうと誤差が統計結果に悪影響するだろう。

CrawleyのR統計本によると、単位量は共変量として説明変数に加えるか、log(目的の量 / 単位量)というように、割り算した上で対数変換せよ。と書かれていた。前者は要因の1つとしてカウントするならばそれでいいとして、後者はすでに古典的な対処法でしょう…（本自体がすでに一時代前のものなので仕方ないですが）。

ついでに、論文などで時々見掛ける、family=gaussian(link="log")、正規分布でリンク関数が対数という一見すると対数正規分布っぽいけれど違うらしい（対数正規分布の分散は平均と共に増大するようだ→平均や分散の式をチェックされたし）、こいつの性質もチェックしてみます。もう一つおまけに、古典統計で対数変換する際によく使うlog(x + 1)変換もチェックしてみました。ただし、これはlog(0)が計算できない問題を回避する目的に限定します。


比較するモデルは下記の6バージョン。これを応答変数が対数正規分布とガンマ分布の場合のそれぞれについて一次回帰式のパラメータ推定（切片と回帰係数）をします：
（色は下記の図とおおよそ対応）
m11：正規分布モデル、割合にした上でlog(x + 1)の対数変換という古典統計の常套手段（グラフ上では黒ライン）
m12：正規分布モデル、割合にした上でlog(x + 1e-10)変換、1の代わりにごく小さい小数(1.0 x 10^(-10))でlog(0)を回避
m13：正規分布モデル、割合にしているが、変数変換はせず、連結関数を対数に指定（以下、log(0)回避には1e-10を使用）
m14：正規分布モデル、単位量をoffset(log(off))として係数1の説明変数に加え、連結関数を対数に指定
m23：ガンマ分布モデル、割合にしているが、変数変換はせず、連結関数を対数に指定（ガンマ分布は正の値しか取れないので、1と2の変数変換のモデルは作れない）
m24：ガンマ分布モデル、単位量をoffset(log(off))として係数1の説明変数に加え、連結関数を対数に指定


結論から言うと、ケースバイケースな複雑な結果になりました…。実際の利用では、m13とm23、またはm14とm24をAICで比較するのがよいでしょう。
・やはりデータの対数変換は止めた方がよい、連結関数：logを用いるべきです。対数変換したモデルでは推定が大きく乱れる場合があり信頼できない（図1_1, 1_2のm11黒、m12赤）。
・データが対数正規分布の時、割合の正規分布モデル（m13緑）は、offsetの正規分布モデル（m14青）よりも推定の分散が小さかった（図1_1、2_1）。ただし、データがガンマ分布の時は、m14の方がm13よりも推定の分散が小さかった（図1_2、2_2）。両者の違いはさほど大きくなかった。
・ガンマ分布モデルでは、割合（m23水色）と、offset（m24ピンク）とは推定値がまったく同じ！（そのため、m23の水色は完全にマスクされている）
・つまり、この数値実験の限りでは、予想に反して連続量を割合や比率にした量をGLMで解析する際、割り算をした値を用いても大した問題はないことになる。offsetの使用による論文中での説明のややこしさや、記述可能なモデルの可塑性などを考えると…割り算していい気がしてきました。しかし、なぜ割り算にしても大丈夫だったのか、背景にある数学的なロジックは朧気なままです…たぶん、元の値と割算値とで確率分布が変わらないからだと思います（よくある整数値の割り算の場合は、小数点を取るようになるので本来のポアソンや二項分布が使用不可になるが）。例外があるとすれば、単位量の方にも確率誤差（その値が属している確率分布からのズレ）がある場合でしょう。その場合はoffsetすべきということかと。
・ちなみに、常套手段なlog(x + 1)変換（m11黒）の推定は危険です。とくに平均値の小さい推定では危険（図1_1、1_2）。log(0)の回避には1の代わりにごく小さな値（1e-10くらい）を足すのがよいでしょう。ただし、古典統計で整数を対数変換する場合は、平均と分散の関係を調整する意図があるので + 1 のままで）。

なお、ここでやっているのは入り口と出口がちゃんと一致するかを確認するだけの数値実験に過ぎません。参考にする際には自己責任でお願いします。



# 以下、こんな関数でGLM推定を1000回繰り返し、パラメータ推定の精度をチェックしました。
# 平均 exp(alpha + beta*X)という一次回帰を考え、こちらで指定したalphaとbetaを用いて、対数正規分布とガンマ分布のデータを発生させ、各モデルによってalphaとbetaを逆推定し、初めに指定した値を再現できるかどうかチェックします。
# alpha, beta, distributionは1が対数正規分布、2がガンマ分布、n.itrは繰り返し数を表します。

cont.skew <- function(alpha, beta, distribution, n.itr) {
estim <- numeric(0) # 後でデータをくっつけるためのイントロンみたいなもの
set.seed(1)
X <- rep(c(1:10)*0.1, each=10)
off <- runif(length(X), min=1, max=10) # 下記のoffsetで用います。数字に特に意味はなし
mean <- exp(alpha + beta*X)*off # 平均：log(Y/off) ~ alpha + beta*X、なので。
sd <- exp(0.5)
for (n in 1:n.itr) {
Y <- rbind(rlnorm(100, meanlog=log(mean), sdlog=log(sd)), # 対数正規分布
rgamma(100, shape = mean^2/sd^2, scale = sd^2/mean)) # ガンマ分布（cf. rgamma()のhelp）
Y <- Y[distribution, ]
D <- data.frame(X, Y, off)
m11 <- glm(log(Y/off + 1) ~ X, family=gaussian(link="identity"), D) # 正規分布モデル（割合、対数変換1）
m12 <- glm(log(Y/off + 1e-10) ~ X, family=gaussian(link="identity"), D) # 正規分布モデル（割合、対数変換2）
m13 <- glm((Y/off + 1e-10) ~ X, family=gaussian(link="log"), D) # 正規分布モデル（割合、連結関数：対数）
m14 <- glm((Y + 1e-10) ~ X + offset(log(off)), family=gaussian(link="log"), D) # 正規分布モデル（オフセット、連結関数：対数）
m23 <- glm((Y/off + 1e-10) ~ X, family=Gamma(link="log"), D) # ガンマ分布モデル（割合、連結関数：対数）
m24 <- glm((Y + 1e-10) ~ X + offset(log(off)), family=Gamma(link="log"), D) # ガンマ分布モデル（オフセット、連結関数：対数）
estim2 <- c(rbind(coef(m11), coef(m12), coef(m13), coef(m14), coef(m23), coef(m24)))
estim <- rbind(estim, estim2)
} # forループ、ここまで
par(mfcol=c(1,2))
plot(density(estim[,1]), lwd=4, xlim=c(alpha-0.5, alpha+2), ylim=c(0,3), main="alpha (intercept)")
for (i in 2:6) { lines(density(estim[,i]), lwd=4, col=i)
abline(v=alpha) } # alpha値の図示
plot(density(estim[,7]), lwd=4, xlim=c(beta-1, beta+1), ylim=c(0,3), main="beta (coefficient)")
for (j in 2:6) { lines(density(estim[,j+6]), lwd=4, col=j)
abline(v=beta) } # beta値の図示
estim.m <- matrix(apply(estim, 2, mean), ncol=2) # alpha, betaの推定値を計算
colnames(estim.m) <- c("alpha", "beta")
estim.m} # 推定値を表示（cont.skew関数、ここまで）

# 図1_1： 対数正規分布のパラメータ推定、平均値が小さい場合
cont.skew(alpha=-1, beta=1, n.itr=1000, distribution=1)
# m11：黒、m12：赤、m13：緑、m14：青、m23：水色、m24：ピンク（水色はピンクと完全に重なって表示されない）
# m12赤、m23&m24ピンク、m13緑、の順に推定がよかった。m11黒は大きくずれた。
# 推定値（1~6の順にm11～m24）、m23とm24（一番下の2つ）は完全に推定値が一致！
# alpha beta
# [1,] 0.3131626 0.3959152
# [2,] -1.0009638 0.9999089
# [3,] -0.8811454 1.0011294
# [4,] -0.8879904 1.0073653
# [5,] -0.8793282 0.9996438
# [6,] -0.8793282 0.9996438

# 図1_2：ガンマ分布のパラメータ推定、平均値が小さい場合
cont.skew(alpha=-1, beta=1, n.itr=1000, distribution=2)
# m11：黒、m12：赤、m13：緑、m14：青、m23：水色、m24：ピンク（水色はピンクと完全に重なって表示されない）
# m14青、m13緑、m23&m24ピンク、の順に推定がよかった。m11黒とm12赤は大きくずれた。
# なぜかガンマ分布の推定なのにガンマ分布モデル（m23&m24）が最強にならない…

# 図2_1：対数正規分布のパラメータ推定、平均値が大きい場合
cont.skew(alpha=2, beta=1, n.itr=1000, distribution=1)
# m11：黒、m12：赤、m13：緑、m14：青、m23：水色、m24：ピンク（水色はピンクと完全に重なって表示されない）
# 平均値が大きくなると、モデルによる推定の違いは小さくなる（中心極限定理による正規分布への収束）
# m12赤、m23&m24ピンク、m11黒、の順に推定がよかった。

# 図2_2：ガンマ分布のパラメータ推定、平均値が大きい場合
cont.skew(alpha=2, beta=1, n.itr=1000, distribution=2)
# m11：黒、m12：赤、m13：緑、m14：青、m23：水色、m24：ピンク（水色はピンクと完全に重なって表示されない）
# 平均値が大きくなると、モデルによる推定の違いは小さくなる（中心極限定理による正規分布への収束）
# m11黒以外は真の値にほぼ収束

（旧サイトより移行）過分散データ：GLM負の二項分布、GLMMによる解析の比較

## 2014.02.20 追記：旧サイトを閉じるため、このページを移植しました。
なお、現在は lmer()は正規分布専用になっており、その他の分布ではglmer()関数を使用します。

##（2012.10.17 追記：GLM関数群とGLMM関数群との間でのモデル選択について、数値実験結果を追加・修正しました）

# 生態学会で話していて、過分散データの取り扱いが話題に上った。負の二項分布を使えという意見と、GLMMを使えという意見と両方があるが、実際どっちの方がよいというのはあるのだろうか？（cf. 2014.6.24追記： GLM負の二項分布とは実態はランダム効果がガンマ分布となったGLMM）

# 過分散 overdispersion：ポアソン分布や二項分布のモデルに当てはめる場合、これらの分布型が平均と分散が表裏一体なのだが（ポアソンでは平均=分散)、実際のデータは大抵はそうでないので分散が過剰となることが多い、という問題。
# GLMのsummary出力結果で、下の方に、こんなのがあって…
# Residual deviance: 5094.7 on 71 degrees of freedom
# …ここで、5094.7/71の値が、1.5を超えているならば、何らかの対処をした方がよいだろう（cf. モデル選択ペンギン本、Zuur et al (2009) Mixed effects models and extensions in Ecology with R. Springer。これのp. 225-6には"1.5"が目安と書かれている）

# 過分散しているデータとは下記のように裾野が片方に広く延びている（上：ふつうのポアソン分布、中：負の二項分布、下：正規ノイズを加えたポアソン分布）。負の二項分布はゼロのあたりに集中、正規ノイズは正負にばらつく。

par(mfcol=c(3,1)) # グラフ領域を縦に3分割
hist(rpois(10000, lambda=exp(1)), freq=F, xlim=c(0,100), breaks=seq(0,20,by=2), col=1) # ポアソン分布（平均 exp(1)）
hist(rnegbin(10000, mu=exp(1), theta=1/4), freq=F, xlim=c(0,100), breaks=seq(0,10000,by=2), col=2) # 負の二項分布（平均 exp(1)、θ 1/4。θはガンマ分布のshapeパラメータに相当、shape=平均^2 / 分散）
hist(rpois(10000, lambda=exp(1 + rnorm(10000, mean=0, sd=4^0.5))), freq=F, xlim=c(0,100), breaks=seq(0,10000,by=2), col=3) # ポアソン分布（平均 exp(1)）＋正規ノイズ（平均0、分散4）

# では、ポアソン分布のデータと、負の二項分布をしたデータ、ポアソン分布に正規乱数ノイズを加えたデータを人工的に作成し、それぞれGLM (poisson)、GLM（負の二項分布）、GLMM（データ数分の個体差があると見なす）で推定した結果を比べてみよう（実行するためのプログラムは末尾に）。


OD(alpha=2, beta=0.5, siml=100, distribution=1) # データがポアソン分布の時
#（マシンパワーによってはけっこう時間が掛かる、アラートも沢山出るが無視）
# ランダム効果がN数と同じ数だけあるよ！？というアラートが大量に出てくるが、推定はちゃんとできているので無視
# Number of levels of a grouping factor for the random effects is *equal* to n, the number of observations

# 黒：GLM (poisson)、赤：GLM（負の二項分布）、緑：GLMM。100回分の推定結果
#（３本のラインがあるはずですが完全に重なっていて、視覚的には識別不能）

# models: glm.pois glm.nb glmm.pois
# aic.d 710.5256 712.1975 712.2007 # devianceに基づくAIC
# prob.d 0.9500 0.0400 0.0100 # 各モデルが選択される確率
# aic.rd 102.8072 101.6366 104.4823 # residual devianceに基づく AIC
# prob.rd 0.6000 0.4000 0.0000 # 各モデルが選択される確率

# GLM (poisson)、GLM（負の二項分布）、GLMMともalpha、betaの推定値はほぼ真の値となった（ということは、単なるGLM（poisson）を使う理由はないのではないか！？）。AIC（deviance）によるモデル選択は確率0.95で真のモデルであるGLM（ポアソン分布）を選択したが、AIC（residual deviance）ではAICの推定値が僅差で負の二項分布モデルが最小だし、誤って負の二項分布モデルを選択する確率が0.4もあった。とはいえ、いずれのモデルもほぼ真の値を推定できているので問題なしか！？



OD(alpha=2, beta=0.5, siml=100, distribution=2) # データが負の二項分布の時
# 黒：GLM (poisson)、赤：GLM（負の二項分布）、緑：GLMM。100回分の推定結果

# models: glm.pois glm.nb glmm.pois
# aic.d 36960.67 902.9074 922.3840 # devianceに基づく AIC
# prob.d 0.00 1.0000 0.0000 # 各モデルが選択される確率
# aic.rd 36572.47 121.0299 534.1787 # residual devianceに基づく AIC
# prob.rd 0.00 1.0000 0.0000 # 各モデルが選択される確率

# GLM (poisson)は両方向に大きく広がった。GLM（負の二項分布）はalphaもbetaも真の値に近かった。GLMMはbetaの値は真の値に近かったが分布の裾が広かった、alphaは盛大にズレた。
モデル選択はすべてGLM（負の二項分布）となった。



# ポアソン分布に正規乱数ノイズを加えたデータを同様に推定してみる
OD(alpha=2, beta=0.5, siml=100, distribution=2)

# 黒：GLM (poisson)、赤：GLM（負の二項分布）、緑：GLMM。100回分の推定結果

# models: glm.pois glm.nb glmm.pois
# aic.d 845220.9 1309.4785 1278.7899 # devianceに基づく AIC
# prob.d 0.0 0.0100 0.9900 # 各モデルが選択される確率
# aic.rd 844616.2 136.0312 674.1263 # residual devianceに基づく AIC
# prob.rd 0.0 1.0000 0.0000 # 各モデルが選択される確率

# GLM (poisson)は真の値から大きく外れた。GLM（負の二項分布）はalphaは盛大にズレた、betaは真の値に近かったが裾野が広かった。GLMMはalphaもbetaも真の値に近かった。AIC（deviance）によるモデル選択は0.99の確率で真の母集団のモデルGLMMを選択したが、AIC（residual deviance）によるモデル選択はすべてのケースで誤ってGLM（負の二項分布）を選択してしまった…。



# 結論として、過分散データと言っても一辺倒ではなく、負の二項分布、正規ノイズ、どちらに近いかによって結果が変わってくるという、ごく当然の結果となった。glm.nb、GLMMのどちらがよいかは、結局データ依存ということになりました。

# 選択基準としては…glm.nb、GLMMでAICを基にモデル選択は、一手間加えれば可能であるようだ。また、glmとglm.nbの方をresidual devianceに基づくAICに修正する計算も試してみたところ（2*logLik(pois) + deviance(pois) + AIC(pois)）、誤ってglm.nbモデルが選択される確率が大きく増大することがわかった。なので、GLMとGLMMのモデルをAICで比較する場合には、GLMMのAICをdevianceベースのAICに修正して用いる必要がある（-2*logLik(pois) - deviance(pois) + AIC(pois.mm)）ということですね。（注：poisはglm(,family=poisson)のモデル、pois.mmはlmer(,family=poisson)のモデル、両者の違いはランダム効果の有無のみ → 下記のプログラム中にモデル式あり）


# しかし、lmerやglmmMLのデフォルトがresidual devianceに基づくAICを計算していること自体が、そういう比較を推奨しないというメッセージにも思える。

# もっとも、ポアソン分布の場合は今回のような選択肢があるけれど、二項分布の場合はGLMMを使うくらいしか対処の方法がない（怪しげなQAICの使用を回避するならば）。

# 2014.06.24 追記：実際のデータのモデル化を考えた時、固定効果のみで推定値を求めるときGLMM (poisson)とglm.nbは異なる推定値になってしまうのは気になるが、もはやモデルの違いと割り切るしかない気がしている。双方ともGLMMなので、誤差はランダム効果でオーバーフィットしてしまうため、AICなどで比較して差異があっても実質的なモデルの良し悪しを反映しているのだろうか？現実のデータの場合は真の答えは分からないし、モデルを使用している以上は仮定を捨てられないので、悩んでも仕方がないと振り返って思います。

# 参照：
# lmerはresidual devianceに基づいたAICを返しているので、直接比較をするにはdevianceに基づいた計算にし直す必要がある
# deviance: -2*logLik()
# residual deviance: deviance()
# AIC(glm())の場合は、deviance + 2*推定するパラメータ数
# AIC(lmer())では、residual deviance + 2*推定するパラメータ数





# 上記の数値実験を実行するためのプログラムはこちら：
# 平均 exp(alpha + beta*x)とし、真の値はalpha = 2、beta = 0.5、あと負の二項分布の歪みのパラメータθ= 1/4、正規乱数ノイズの分散4も設定
# （ただし、θや分散の値については今回はチェックはしない）

library(MASS) # glm.nbの呼び出し
library(lme4) # lmerの呼び出し

OD <- function(alpha, beta, siml, distribution) {
# 切片、傾き、シミュレーション回数、データ分布（1：ポアソン分布、２：負の二項分布、３：ポアソン分布＋正規ノイズ）
estim <- aic.d <- aic.rd <- prob.d <- prob.rd <- numeric(0) # 後でデータをくっつけるためのイントロンみたいなもの
set.seed(36) # 乱数の種を指定（他の値にすると推定エラーが出るかもですglm.nbやlmerの問題）
for (s in 1:siml) {
Noise <- rnorm(100, 0, 4^0.5) # ノイズの分散を4とした
x <- rep(c(0:9), each=10)
y <- rbind(rpois(100, lambda=exp(alpha + beta*x)), # ポアソン分布
rnegbin(100, mu=exp(alpha + beta*x), theta=1/4), # 負の二項分布
rpois(100, lambda=exp(alpha + beta*x + Noise))) # 普通のポアソン分布にノイズを追加
y <- y[distribution,] # データ分布の選択
ID <- 1:length(y) # データの数だけ個体差（ID）を設定
D <- data.frame(x, y, ID)
pois <- glm(y ~ x, family=poisson, data=D) # GLM (poisson)
negb <- glm.nb(y ~ x, data=D) # GLM (負の二項分布)
pois.mm <- lmer(y ~ x + (1|ID), family=poisson, data=D) # GLMM
estim2 <- c(cbind(coef(pois), coef(negb), fixef(pois.mm)))
estim <- rbind(estim, estim2) # simlの数だけ推定値を重ねていく
aic.p <- AIC(pois) + 2*logLik(pois) + deviance(pois) # residual devianceに基づくGLMのAIC
aic.n <- AIC(negb) + 2*logLik(negb) + deviance(negb) # residual devianceに基づくGLM.nbのAIC
aic.m <- -2*logLik(pois)[1] - deviance(pois) + AIC(pois.mm) # devianceに基づくGLMMのAIC
Aic <- c(AIC(pois), AIC(negb), aic.m) # deviance に基づくAIC
Aic2 <- c(aic.p, aic.n, AIC(pois.mm)) # residual deviance に基づく AIC
aic.d <- rbind(aic.d, Aic)
aic.rd <- rbind(aic.rd, Aic2)
prob.d <- c(prob.d, which.min(Aic))
prob.rd <- c(prob.rd, which.min(Aic2))
}
aic.d <- apply(aic.d, 2, mean)
aic.rd <- apply(aic.rd, 2, mean)
prob.d <- factor(prob.d, levels=factor(c(1:3))) # モデルを示す番号をカテゴリー化
prob.d <- table(prob.d)/siml # AICモデル選択の結果を集計
prob.rd <- factor(prob.rd, levels=factor(c(1:3))) # モデルを示す番号をカテゴリー化
prob.rd <- table(prob.rd)/siml # AICモデル選択の結果を集計

select <- rbind(aic.d, prob.d, aic.rd, prob.rd)
colnames(select) <- c("glm.pois", "glm.nb", "glmm.pois")
par(mfcol=c(1,2))
plot(density(estim[,1]), lwd=4, xlim=c(-2,6), ylim=c(0,1), main="alpha (intercept)")
for (i in 2:3) { lines(density(estim[,(2*i-1)]), lwd=4, col=i)
abline(v=alpha) } # alpha値の図示
plot(density(estim[,2]), lwd=4, xlim=c(0,1), ylim=c(0,6), main="beta (coefficient)")
for (j in 2:3) { lines(density(estim[,(2*j)]), lwd=4, col=j)
abline(v=beta) } # beta値の図示
return(select) }

負の二項分布のパラメータ推定：R と WinBUGSの比較

過分散した整数データのモデリングによく使う負の二項分布、Rではglm.nb関数でできるが、WinBUGSの例をほとんど見ない。ちょうど使う用事があったので使い方をテストしてみた。

# サンプルデータは、平均 exp(2)、dispersion parameter 1.5 の負の二項分布から抽出した N = 1000 のデータ（後で対数リンク関数を使用するのでexp()で表記）
require(MASS)
Y <- rnegbin(1000, mu=exp(2), theta=1.5) 
# 注：乱数なので実行する度に結果は少しずつ変わる

# まずはRのglm.nb関数による推定
require(MASS)
summary(glm.nb(Y ~ 1))

Coefficients:

                  Estimate Std. Error z value Pr(>|z|)    

(Intercept)  1.97519    0.02771   71.28   <2e-16 ***

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Negative Binomial(1.5895) family taken to be 1)

# 平均 1.97519（expの中身に相当）、dispersion parameter 1.5895 と順当な推定

# 次にWinBUGS、負の二項分布だからdnegbin、と試してみるが…（thetaがdispersion parameterに相当しそう）

data <- list(N=1000, Y=Y)

parm <- list(theta=rlnorm(1), alpha=rnorm(1))

model <- function() {

theta ~ dgamma(1E-1, 1E-2)

alpha ~ dnorm(0, 1E-6)

for (i in 1:N) {

Y[i] ~ dnegbin(p[i], theta) # 負の二項分布では平均をそのまま使えない

p[i] <- theta/(theta+mu[i]) # pと平均muの関係（pは0~1の値を取る）

log(mu[i]) <- alpha # muは正の値

}

}

mcmc <- wbugs(data, parm, model, n.iter=10000, n.burnin=5000, debug=T) # 自作ラッパー関数wbugsを使用

# order of negative binomial *** must be an integer というエラー

# theta は負の二項分布の狭義の定義通り、整数でないといけないようだ

# 戸惑いつつも多項分布の事前分布で thetaが正の整数しか取らないようにし再計算

data <- list(N=1000, Y=Y, q=rep(1,20)/20)

parm <- list(theta=sample(c(1:5), 1), alpha=rnorm(1))

model <- function() {

theta ~ dcat(q[]) # 事前分布1~20の整数（多項分布使用）

alpha ~ dnorm(0, 1E-6)

for (i in 1:N) {

Y[i] ~ dnegbin(p[i], theta)

p[i] <- theta/(theta+mu[i]) 

log(mu[i]) <- alpha 

}

}

mcmc <- wbugs(data, parm, model, n.iter=10000, n.burnin=5000)

               mean        sd          2.5%      25%        50%        75%        97.5%   Rhat   n.eff

theta        2.000       0.000    2.000    2.000       2.000       2.000       2.000          1        1

alpha       1.976       0.026    1.927    1.958       1.976       1.993       2.026          1  1500

deviance 6032.409 1.446 6031.000 6032.000 6032.000 6033.000 6036.000    1  1500

# 平均 alphaの推定は1.976と正しく推定できているのだが、thetaは当然のように整数値となった…

# Rの?rnbinomなどを読み返してみると、"dispersion parameter (the shape parameter of the gamma mixing distribution). Must be strictly positive, need not be integer" とある。

# gamma mixing distributionのshapeパラメータに相当するということか。

# なので、ポアソン分布にガンマ分布を混ぜた階層ベイズとして負の二項分布の推定を試みる

data <- list(N=1000, Y=Y)

parm <- list(shape=rlnorm(1,log(1),log(2)), alpha=rnorm(1), lambda=rlnorm(1000,log(1),log(2)),F) 

model <- function() {

shape ~ dgamma(1E-1, 1E-2)

alpha ~ dnorm(0, 1E-6)

for (i in 1:N) {

Y[i] ~ dpois(lambda[i]) # 負の二項分布を階層ベイズで表現

lambda[i] ~ dgamma(shape, rate[i])

rate[i] <- shape / mu[i] # rate と shape、平均 mu の関係を表す

log(mu[i]) <- alpha 

}

}

mcmc <- wbugs(data, parm, model, n.iter=10000, n.burnin=5000)

                    mean     sd     2.5%      25%      50%      75%    97.5%  Rhat   n.eff

shape           1.587  0.087    1.422    1.522    1.587    1.647    1.755 1.002  1200

alpha           1.975  0.028    1.917    1.957    1.975    1.993    2.029 1.003  1200

# lambdaは省略

# 平均 1.975、shape 1.587 と 元の値やglm.nbの推定とほぼ一致、これが正解のようだ

# それにしてもややこしい…負の二項分布には異なる定義が多すぎる！