[FOSS4G 2012 Hokkaido
T シャツ]
何やらそういうイヴェントがあって，
そのときに作られた T シャツとのこと．
FOSS4G とは free & open source software for geospatial
(この絵のもととなった
ポスター)
……
まあ，
QGIS とか GRASS とか，
そういうソフトウェアたち．
sudo do-release-upgrade -d
で簡単にできた．
時間はかなりかかったけど．
札幌初 (?) の R 集会，
7/7 (土) SappoRo.R #1
で「統計モデリング入門」だけでなく，
これも割引販売しようかなあ
(3000 円ぐらいで)，
と画策中．
てもとにはこの二冊しか残ってないんだけどね．
じつはこの問題は単純な 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
)
HTML::Template->new(die_on_bad_params=>0)というオプション指定を教えてもらった! ありがとうございます．
今回行ったノードのメンテナンスはXenServerのサービスを含む内部の問 題です。XenServerをツールにより監視し、XenServerがダウンする前に 対応したもので、障害発生までに至っておりません。そしてあいかわらず， このノードで動かしてる web server 機の reboot 直後の状況はよくない．
sudo /etc/init.d/postfix flush
で流す
……
