system
で計測させておいたところ，
作業にちょうど一時間かかった，
とわかった．
多くの時間はコンパイルではなく，
標準 package のテストに使われてしまう，
というのが最近の R ビルド
(とくにオプション指定ナシの場合)
の傾向ですな
……
sudo /sbin/fsck.ext3 -C -y /dev/hdb[13]．
bugs
オブジェクト (bugs class オブジェクトの
概要)
をぽいっとわたしてオワり，
ではすまされないからなぁ
……
1915 研究室発．
雪．
寒い．
氷点下．
1930 帰宅．
晩飯の準備．
晩飯．
reshape()
なる関数があるのだが
……
えーい，
今回も細かいところで手こずってしまった．
列名に id
があるとエラーがでる，
とか．
これは idvar
オプションで対処する
……
といったことを
reshape(data.frame)
に追記．
いつものことながら，
登録番号： 0075 発表形式： 一般講演（ポスター） 演題： ウミガメ上陸数のベイズ統計モデリング 著者： *久保拓弥,重田麻衣 (北大・地球環境), 亀崎直樹 (日本ウミガメ協議会) 要旨 (著者名欄も含めて 800字以内): 個体群の時系列データから個体群の状態変化を推定しようとする方法の検討は，多く の応用生態学で必要とされているだけでなく，生態学の基礎的な理解にとっても重要 な問題である．この発表では，日本ウミガメ協議会が 50 年間以上のウミガメ上陸観 察によって蓄積したデータを使って，この問題をうまくあつかう方法について議論し たい． 今回の解析では太平洋側の 6 箇所のウミガメ上陸観測地点で得られた上陸数の時系 列データを解析するために新しく階層ベイズモデルのひとつである一般化状態空間モ デルを開発した．これによって，ウミガメ (メス) 上陸数という観測データの背後に あるウミガメ密度や個体群増減速度の時間変化をデータにてらしあわせて推定できる． 同時に，全観測地点に共通する年変動や上陸地点ごとの random effects も仮定する ことで，年ごとの気象変化や上陸観察海岸ごとの無方向なばらつきも考慮した，より 現実的な統計モデリングとなっている． 観測データを使った Markov chain Monte Carlo 法によって，この階層ベイズモデル の諸パラメーターの事後分布を推定した．年変動する繁殖活性に関して全国共通の年 変動らしきものが存在し，何らかの要因でウミガメ繁殖活性が広い範囲にわたって同 調している可能性も示された．ただし，観測地点ごとに大きな不確定性があることも 示された．これらの変動にくらべて対数成長速度や対数密度の確率論的なばらつきは かなり小さく，観測された期間内のウミガメの個体群密度の変動そのものは急激なも のではなかった，と推定された．しかしながら，その変化の方向は観測地点ごとに異 なっていて，地点によっては上陸数が減少傾向にある可能性が示された．
vim な生活をしてる私にとって，
たしかにこれって便利だけど
……
いろいろ設定を工夫してみても Gmail
ショートカット
との「ぶつかり」を回避できない．
glm()
のパラメーター推定値 (estimates of coefficients)
の summary()
の表にでてくる t だの z
だのがちょっと気になったのだが
……
family = gaussian
の場合の例を考えると，
> (y <- rnorm(10, mean = 0.1, sd = 1)) [1] 0.0035866 -0.1161750 -0.6117348 0.8111811 -0.5777867 1.4966526 0.2464104 [8] -0.7588385 3.2370962 -0.8850454上のように， てきとーなる例題データを作って
glm()
関数で推定してみると下のようになる．
> summary(glm(y ~ 1))
Call:
glm(formula = y ~ 1)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.170 -0.888 -0.341 0.385 2.953
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.285 0.404 0.7 0.5
(Dispersion parameter for gaussian family taken to be 1.6343)
Null deviance: 14.709 on 9 degrees of freedom
Residual deviance: 14.709 on 9 degrees of freedom
AIC: 36.24
Number of Fisher Scoring iterations: 2
このように係数の推定値 (Estimate)
と標準誤差 (Std. Error)
の比である t value
(これは t 検定の検定統計量とだいたい同じもの)
とかを使って
「t 検定の P 値とだいたい同じと思われる確率」
を評価している．
(後記: 翌日のぎょーむ日誌も参照)
Std. Error)
は標準偏差 / 標本数とだいたい同じ量になっている．
ここでの計算は
は標準偏差 / (標本数 - 1)
ではない．
family = gaussian
でない場合だ．
先ほどと同じような例題で計算させてみると，
下のようになる．
> (y <- rpois(10, lambda = exp(0.1)))
[1] 2 1 4 3 1 1 0 1 2 3
> summary(glm(y ~ 1, family = poisson(link = "log")))
Call:
glm(formula = y ~ 1, family = poisson(link = "log"))
Deviance Residuals:
Min 1Q Median 3Q Max
-1.897 -0.651 -0.253 0.648 1.410
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.588 0.236 2.49 0.013
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 8.6586 on 9 degrees of freedom
Residual deviance: 8.6586 on 9 degrees of freedom
AIC: 33.14
Number of Fisher Scoring iterations: 5
family = gaussian
の場合とのちがいは:
t value
でなく z value
になっている
(ただし算出方法は同じで Estimate/Std.Erro)
Pr(> |z|)
は pnorm(0, 0.588, 0.236) * 2
で計算されている
……
ただし，
じつは family = gaussian
の場合であっても同様では
(pt() を使ってないのでは)? なる疑惑あり
(後記: 翌日のぎょーむ日誌参照)
Pr(> |z|)
の確率が 0.05 より小さくても
***
などといった「お星様」がつかない
……
これは最近の version の R になって新しく改善されたところ
(後記 2009-08-17:
と思ったら私のうっかりまちがいで，
「星」がつくかどうかは
R の option の #{code show.signif.stars}
だけに依存している
(#{code options(show.signif.stars=FALSE)}
などで設定可能)
→ 参照)
summary(glm(..., family = poisson))
の coefficients の表をみて「係数の検定やりました」
とか誤解するヒトの数は少しばかり減るだろう．
Pr(> |z|)
は何をあらわしているのか?
これは係数の推定値 (estimate)
の Wald confidence interval (Wald CI; Pawitan 本 2.7 節など参照)
で「区間の端にちょうどゼロがくるように
設定した大きさ p の Wald CI の p」
をあらわしている
(信頼限界; confidence limits)，
ということになる．
上の例でいうと
「大きさ 98.7% (= 100 - 1.3)
もしくはそれより大きい Wald CI をとると
ゼロが含まれる」
という意味だ．
family = gaussian
ではない
summary(glm(...))
の係数の z 値や
Pr(> |z|)
をみたときの説明のしかたとしては
library(lattice)
でシュートごと・花茎ごとに図示してみると，
なかなか興味ぶかい，
とわかった．
しかしながら，
今回はあえて glmmML(.., family = poission)
だけを使う (先日は花と芽，本日は種子数)．
M1 に何もかもいっぺんにツメこみ教育をするのは不可能なので．
random effects は花差由来と仮定．
あとから気づいたんだけど，
これはシュート差にすればよかった．
ま，
いーか．
nlme
ぐらいしか使えないけどね．
fit <- glmmML(...)
したあとに，
print(fit)
で表示される
coef se(coef) z Pr(>|z|)
(Intercept) 0.307 0.518 0.592 0.550
x -0.511 0.292 -1.748 0.081
z だの Pr(>|z|)
だのは，
fit
オブジェクトのどこに格納されているんですか?
names(fit)
しても見つかりません，
という質問．
どこにも格納されていません．library(glmmML) したあとに
print(print.glmmML, envir = as.environment("package:glmmML"))
するとわかりますが，この中で
coef <- x$coefficients
se <- x$coef.sd
tmp <- cbind(coef, se, coef/se, signif(1 - pchisq((coef/se)^2, 1), digits - 1))
というふうに (print(result) するたびに) Wald 統計量 z を計算しています．
glmmML()
では Wald 統計量 z
の Pr(>|z|)
は 1 - pchisq((coef/se)^2, 1)
で評価していたのか
……
print(summary.glm, envir = as.environment("package:stats"))
してみると，
昨日
のぎょーむ日誌で私が憶測してたとーり，
family = gaussian
でない場合の
summary(glm(...))
の Pr(>|z|)
は pvalue <- 2 * pnorm(-abs(tvalue))
として評価しているな．
family = gaussian
なんかの場合は
pvalue <- 2 * pt(-abs(tvalue), df.r)
と計算しているのか．
予想どーり t 分布の計算してるんだけど，
これって絶対値にマイナスをつける必要があったんだ
……
知らなかった．
(後記: これは t 分布の「下側のスソ」の
確率を計算してるだけでしょ，
とあとで教えていただいた
…… たしかにそれだけのことですな)
http://hosho.ees.hokudai.ac.jp/~kubo/c
で切れてる (.../ce/ のつもりだった)．
何かしくじったかな?
http://hosho.ees.hokudai.ac.jp/~kubo/c
でアクセスしてきたヒトは
Apache2 がホントの
自由集会ペイジ
まで自動的にごあんなーい，
と設定してやろう!
……
.htaccess
に RewriteRule ^c$ /~kubo/ce/EcoSj2009.html [R=301,L]
と追加．
よしよし，
うまくいった．
