install.packages(rjags)
がエラーにはばまれてインストールできん!
jags
(1.04) がインストールされていたためだった．
sudo apt-get remove jags
してから，
JAGS の
source forge
から JAGS-2.1.0.tar.gz
をダウンロードして，
手動でインストールすれば，
rjags
の最新版はインストールできる．
boxplot()
とかの「ヒゲ」の両端は近似的な 95% 区間を示していて
(Tukey 式)
ホントの 95% 区間ではない．
bplot()
とかでいいぢゃない
……
とかいいかげんな返信したんだけど，
これの style="quantile"
指定は 95% ではなくて 90% 区間だった!
plotQuant()
関数．
mcmc オブジェクトが主な対象だが
(ふつーの data.frame でも使えるらしい)，
まあ WinBUGS な出力の作図をしたいよーで
……
はい，
そうですね．
いやー，
ぼくはこんなのは自分で
ぱっぱっと作図関数を組み合わせて目的の図を
つくっちゃう
からねえ，
と負けおしみを述べてみる．
which.max()
ってあったよな
……
こんどからあらしコメントするとき前に，
間瀬さんの R 間数本を確認することにしよう
……
make pdf
を命じて，
その 3 ペイジぶんだけダウンロードできるようにした．
はんぱな対処ではあるけれど．
\hour
と
\minute
が使えたのだけど，
いまのやつでは使えないので
……
/usr/share/texmf/ptex/platex/base/jbook.cls
などからコードをかっぱらってきたりして，
下のごとく見苦しく実装．
\hour\time \divide\hour by 60\relax
\@tempcnta\hour \multiply\@tempcnta 60\relax
\minute\time \advance\minute-\@tempcnta
\newcommand{\zeropadding}[1]{\ifnum #1<10 0#1 \else #1 \fi}
\newcommand{
\updatetime}{\the\year--\zeropadding{\the\month}--\zeropadding{\the\day}
\ \
\zeropadding{\the\hour}:\zeropadding{\the\minute}
}
pdftk とかでは，
用紙サイズ変更は無理そう．
もちろん pdftk
つかえば異なる用紙サイズのペイジを統合した PDF ファイルが
生成できるし印刷も問題ないんだけど，
甲山さんにそれ説明すんのもめんどうだしね．
convert とか mogrify
とか使えばサイズは変更できるが「A4 サイズ」とゆー指定が難しいように思える
(解像度 + DPI 指定?)
……
しょーがないので evince
で開いて「PDF ファイルとして印刷出力」わざ，
用紙サイズはもちろん A4 で「用紙サイズにあわせて縮小」
……
うーん，
これでうまくサイズ変更はできたのだけど，
とちゅう ghostscript
の呪われ介在があるので，
もともと 1.1MB だったファイルサイズが
(これ以上の PDF 圧縮不可能な) 13MB ファイルになってしまった．
pdftk member.pdf poster.pdf output lab.pdf compress
とくっつけて
……
やっぱり圧縮できないな．
もーいーや．
evince
からの印刷出力がいまいちキレイではないので
(半透過色を使うとペイジ全体が少しグレイがかる?)，
acrobat reader から印刷しようとしたのだが
(1 ペイジに 6 枚というレイアウト)
……
なぜかしら上下反転するので，
いっそのことあらかじめ上下反転してやったら，
うまく印刷出力できた．
/var/log/rkhunter.log
のチェックとか．
librarary(glmmML)
の
glmmML()
で推定すると，
> library(glmmML)
> d <- read.csv("data7.csv")
> glmmML(cbind(y, N - y) ~ x, family = binomial, data = d, cluster = id)
Call: glmmML(formula = cbind(y, N - y) ~ x, family = binomial, data = d, cluster = id)
coef se(coef) z Pr(>|z|)
(Intercept) -4.19 0.878 -4.77 1.8e-06
x 1.00 0.207 4.84 1.3e-06
Scale parameter in mixing distribution: 2.41 gaussian
Std. Error: 0.22
となる．
> library(MCMCglmm) > fit <- MCMCglmm(cbind(y, N - y) ~ x, data = d, family = "multinomial2")と推定してみる．
family は "binomial"
ではなく "multinomial2"．
random
オプション指定の必要なし．
default ではデータの行ごとに，
R ストラクチャーなる random effects が設定されてるみたいなんだよね．
> summary(fit)
Iterations = 12991
Thinning interval = 3001
Sample size = 1000
DIC: 667.64
R-structure: ~units
post.mean l-95% CI u-95% CI eff.samp
units 6.79 3.65 10.1 480
Location effects: cbind(y, N - y) ~ x
post.mean l-95% CI u-95% CI eff.samp pMCMC
(Intercept) -4.242 -6.177 -2.523 864 <0.001
x 1.020 0.633 1.474 907 <0.001
個体差のばらつきは R-structure
(共分散ゼロの事前分布)
のところに入っていて，
これは VCV 行列の要素の事後分布を示している．
SD で表示したければ，
> summary(sqrt(fit$VCV))
Iterations = 3001:12991
Thinning interval = 10
Number of chains = 1
Sample size per chain = 1000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
2.5843 0.3304 0.0104 0.0167
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
2.01 2.35 2.56 2.78 3.34
まあ，
妥当な値．
ここでは平均値や中央値が 2.58
といった値になっているのにたいして，
glmmML()
の推定結果に示されているのは mode に近い値だろうから
(最尤推定しているので)，
まあ妥当なのでしょう．
plot(sqrt(fit$VCV))
MCMCglmm()
で random =~ units
と指定すると，
これは何をどう指定してるんだか，
事前分布が多変量正規分布になって分散共分散行列
の事後分布なんかが推定される．
うーむ?
random
引数の指定がしんどい
(prior 指定の工夫も必要なのだろう)．
ここをまちがえると，
random effects の推定はむちゃくちゃになる．
なぜかしら fixed effects のパラメーターの推定は，
それほど影響うけてないように見えるが．
例題が簡単なせいか?
