午後はもうメイルなんぞ書いてたまるかー
……
と思っていたのだけど，
自分のデータを「GLM にかける」方法がわかりませーん
(ついでに「統計モデリング入門」は買わずに借りてまーす，とか)
……
というこんさるメイルを
(大学院生からではなく中堅研究者から)
いただいて「いまどき『GLM にかける』とか言ったりするやつがいるのかよー」
とやや逆上ぎみとなり
(これは某大統領を連想させるいいぐさで，
しかも内容が「アリ」データであったのもよくなかったのかも)
……
一時間以上をこの返信かきについやしてしまった．
うう，
これは敗北なのかも
……
じつはこの問題は単純な GLM ではダメなのですが，ともあれ R の glm()
から始めたいと思います．ふつうのデータなので，なぜ GLM あてはめが
できないのか，その点がまったく不明なのですが，以下のように操作すれ
ばよいでしょう．
まず，よみこみから作表までです．データを目にみえるかたちにしないと
データ解析は始まりません．
d <- read.table("Ant_Data.txt", header = TRUE)
table1 <- table(factor(d$NEST))
table2 <- tapply(d$SSR, d$NEST, sum)
tab1 <- rbind(table1, table2)
rownames(tab1) <- c("total", "SSR")
print(tab1)
ここで作った tab1 を利用することもできますが，以下ではもっと簡単に
やります．R の glm() 関数を使って，GLM の一部であるロジスティック
回帰 (このあたりの区別はきちんとできてきるでしょうか?) を実施し，線形
予測子の回帰係数を推定します．
fit1 <- glm(SSR ~ NEST, data = d, family = binomial)
print(fit1)
これで一見するともっともらしい推定結果が表示されます．査読者氏が
指摘するのはこういうことでしょうが，これでは地区間の差はいえません
ここでは地区ごとのパラメーターが出てきますが，これらを比較して「この
地区は他とは違う」ときちんと言うための方法がないからです．
またここでは，「地区ごとに SSR 発生しやすさは，完全に独立にきまる」
といったことを仮定していますが，この実験のわくぐみではそれは，強す
ぎる仮定でしょう．同種の**を同地域とかで見ているわけですから．
こういう強すぎる仮定のために，そもそも glm() で得られたパラメーター
の推定値はちょっと偏った値になっています．
まず，モデルを改良して地区をランダム効果のようにあつかう必要があり
ます．このあたりの考えかたは「統計モデリング入門」の第 7 章とかを
参照してください．
さて，第 7 章で紹介している glmmML() 関数を使った GLMM 版のロジ
スティック回帰をすれば問題が解決するのかというと，これでもまだ
不足です．glmmML() では地区ごとの効果を正確に推定できても，それら
を比較するのが難しいためです (めんどくさいプログラミングをやるつ
もりなら，glmmML() を使ってもだいじょうぶです)．
けっきょく，この問題は第 10 章とかで紹介している階層ベイズモデル
を使うのがらくで，階層ベイズ版のロジスティック回帰をやった結果は
mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff
a 0.161 0.644 -1.135 -0.218 0.141 0.524 1.635 1.006 360
r[1] 2.368 0.852 0.858 1.796 2.305 2.854 4.257 1.002 1100
r[2] -1.044 0.696 -2.587 -1.452 -1.013 -0.617 0.296 1.005 380
r[3] 0.032 0.695 -1.411 -0.377 0.048 0.468 1.376 1.004 460
r[4] -0.329 0.732 -1.876 -0.748 -0.308 0.115 1.103 1.009 220
r[5] 0.144 0.684 -1.431 -0.277 0.156 0.555 1.490 1.004 540
r[6] -1.028 0.687 -2.474 -1.426 -1.008 -0.586 0.251 1.006 350
r[7] -0.332 0.671 -1.756 -0.744 -0.311 0.074 1.034 1.007 280
sigma 1.550 0.696 0.736 1.090 1.384 1.798 3.451 1.001 1500
このようになります．r[1] から r[7] までが Cj, .., T に対応した地区の
効果の大きさになっています．「パラメーターの分布の 95% 区間に
ゼロが含まれていない地区の効果については，『平均的な地区』とは
異なっているとみなす」という規準を適用したとすると，地区の 1 番
つまり Cj だけが他の地区より SSR の発生率が高いと言えます．ただ
単に高いのではなく，格段に高いと言えます．
以上のような推定計算をするためには，下のような R のコードを
実行する必要があり，さらに R から WinBUGS (第 9 章を参照) を
呼び出して使用するので，WinBUGS がインストールされている
必要があります．
source("http://hosho.ees.hokudai.ac.jp/~kubo/ce/r/R2WBwrapper.R")
clear.data.param()
set.data("SSR", d$SSR)
set.data("Nest", as.numeric(d$NEST))
set.data("N.sample", nrow(d))
set.data("N.nest", length(unique(Nest)))
modelf <- function()
{
for (i in 1:N.sample) {
SSR[i] ~ dbern(p[i])
logit(p[i]) <- a + r[Nest[i]]
}
a ~ dnorm(0.0, 1.0E-4)
for (j in 1:N.nest) {
r[j] ~ dnorm(0.0, tau)
}
tau <- 1 / (sigma * sigma)
sigma ~ dunif(0, 100)
}
write.model(modelf, "model.bug.txt")
set.param("a", 0)
set.param("r", rep(0, N.nest))
set.param("sigma", 1)
post.bugs <- call.bugs(
file = "model.bug.txt",
n.chains = 3,
debug = FALSE,
n.iter = 5500, n.burnin = 500, n.thin = 10
)
