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]
と追加．
よしよし，
うまくいった．