「ぎょーむ日誌」目次に戻る | KuboWeb top に戻る | twilog | atom

ぎょーむ日誌 2012-02-24

苦情・お叱りは，たいへんお手数かけて恐縮ですが，久保 (kubo@ees.hokudai.ac.jp) までお知らせください．

本日 (kubolog20120224) | 次の日 | 1 日前 | 7 日前 | 31 日前 | 365 日前 | top

2012 年 02 月 24 日 (金)

0620 起床．ラジオ体操．弁当の準備．朝飯．コーヒー． 0750 自宅発．昨日午後からずっと大雪つづいている． 0820 研究室着．

しばらくしたら，雪はやんだ．
今日のりはびり …… データ解析こんさるに由来する，ちょっとばかりめんどうな R プログラミング．
昼飯．
うーん，とりあえずこんなところかな．やったことは，
- R の library(MCMCglmm) を使った GLMM のパラメーターの MCMC サンプリングを複数回 (3 回とか) を独立にやる
- それらの結果をたばねて， library(R2WinBUGS) が定義する bugs オブジェクトに変換する
かくのごとく bugs オブジェクトに変換すれば， Gelman-Rubin の R-hat の計算とか，あるいはいろいろな図示がラクになる …… といったもくろみ．
いやはや，院生とかに「lmer() なんかは使うのヤメて， MCMCglmm() を!」とススめているんだけど (皆さん WinBUGS とか JAGS なんかを使用すると，なぜかしら健康をいささか損なうおそれがある …… といった迷信にとりつかれているようなので) …… いつものごとく「あーいえば，こーいう」といいますか，「これってどうやってサンプリング数とか決めたらいいんですか?!」といかにも院生じみた質問・詰問ばかりが返ってくる日々に，はなはだしく疲労困憊してしまったため．

とりあえず，現在の動作はこんなかんぢで．

library(MCMCglmm)
burnin <- 1000
thin <- 50
nitt <- 1000 * thin + burnin
n.chain <- 3
list.MCMCglmm <- lapply( # ここは並列化したほうが良い
	1:n.chain,
	function(chain) MCMCglmm(
		cbind(peak_no, flower_no - peak_no) ~ log(stem_Do) * sex,
		random =~ stem_id + ind_id,
		data = d, family = "multinomial2",
		pr = TRUE, # !!! これはめんどうのもと (後述) !!!
		nitt = nitt, thin = thin, burnin = burnin
	)
)

source("merge.MCMCglmm.R") # 関数定義ファイルのよみこみ
post.bugs <- merge.MCMCglmm( # 変換
	list.MCMCglmm,
	n.iter = nitt, n.burnin = burnin, n.thin = thin
)

以上のように R 内で実行して，いつものごとく plot(post.bugs) などとやると，以下のように図示されたりする．

あるいは， print(post.bugs, digits.summary = 3) すると以下のように出力される．

Inference for Bugs model at "cbind(peak_no, flower_no - peak_no) ~ log(stem_Do) * sex", fit using MCMCglmm(),
 3 chains, each with 51000 iterations (first 1000 discarded), n.thin = 50
 n.sims = 3000 iterations saved
                    mean    sd   2.5%    25%    50%    75%  97.5%  Rhat n.eff
stem_id            0.005 0.019  0.000  0.000  0.000  0.000  0.063 1.089    28
ind_id             0.012 0.042  0.000  0.000  0.000  0.000  0.139 1.008   280
units              0.502 0.102  0.332  0.429  0.491  0.565  0.727 1.001  3000
(Intercept)       -0.357 0.288 -0.930 -0.549 -0.350 -0.164  0.203 1.001  3000
log(stem_Do)       0.680 0.257  0.161  0.508  0.678  0.853  1.188 1.001  2600
sexm              -0.551 0.386 -1.310 -0.810 -0.549 -0.300  0.201 1.001  3000
log(stem_Do):sexm -1.012 0.438 -1.862 -1.303 -1.008 -0.710 -0.167 1.001  3000
r.stem_id[01]      0.007 0.077 -0.105  0.000  0.000  0.001  0.193 1.018  3000
r.stem_id[02]     -0.003 0.070 -0.149 -0.001  0.000  0.000  0.118 1.006  3000
r.stem_id[03]     -0.010 0.079 -0.217 -0.001  0.000  0.000  0.086 1.010  2800
r.stem_id[04]     -0.004 0.064 -0.143 -0.001  0.000  0.000  0.106 1.029  1500
r.stem_id[05]      0.010 0.074 -0.091  0.000  0.000  0.001  0.213 1.041   470
r.stem_id[06]      0.003 0.067 -0.122  0.000  0.000  0.001  0.174 1.038  1200
r.stem_id[07]      0.009 0.074 -0.079  0.000  0.000  0.001  0.172 1.018  3000
r.stem_id[08]     -0.004 0.066 -0.168 -0.001  0.000  0.000  0.114 1.042  1400
r.stem_id[09]      0.000 0.067 -0.124  0.000  0.000  0.000  0.128 1.010  2600
r.stem_id[10]      0.003 0.070 -0.109  0.000  0.000  0.000  0.145 1.007  3000
r.stem_id[11]      0.004 0.079 -0.118  0.000  0.000  0.001  0.154 1.041   670
r.stem_id[12]     -0.004 0.071 -0.162 -0.001  0.000  0.000  0.114 1.011  3000
r.stem_id[13]      0.002 0.069 -0.126  0.000  0.000  0.001  0.142 1.033  1700
(略)

また「どわすれ」するときのために書いておくと …… この MCMCglmm → bugs オブジェクト変換で面倒なのは， pr = TRUE と指定したときに MCMC サンプリングされる，いわゆる random effects なパラメーターの名前を書き換えねばならぬところ．
もとの MCMCglmm オブジェクト (注: pr = TRUE と指定して生成したもの) が x であるとすると， x$Sol という行列には "stem_id.i01", "stem_id.i02", "stem_id.i03", ... といった列名 (colnames) で MCMC サンプルが格納されているので，これを sub("(.+_id)(\\.)([A-Za-z_]*)([0-9]+)", "r.\\1[\\4]", ...) といったかんじの 呪術的 というしかない R 正規表現で置換した．
アタマに r. という文字列をつけたのは，これをやっとかなないと plot(post.bugs) がエラーになるため …… さらに，これらの添字が 1 から始まる通し番号でないときには， plot(post.bugs) の左パネルには正しく図示されないみたいだ．上の図でも ind_id は示されていない (右パネル一番下には図示されているが)．
まあ，自力解決できるよう院生のアタマを reprogramming するよりは，自分で R ぷろぐらみんぐやってたほうがラク …… ということで．
で，説明メイル書き …… これで一件落着 …… すればいいけど．
ふう
ちょっと休憩．
岩波書店から「統計モデル本の校正ゲラを送りました」メイルが． 3/12 までに返せとのこと．作業期間は 2 週間か ……
そして今から 3 週間後の 3/19 は生態学会大会のデータ解析自由集会なんだが，準備は中座したままだな．
R Graphics Manual で分割表 (contingency table) の図示法を調べてみる．標準 graphics pacakge で検索してみると …… mosaicplot() とか plot.table(), spineplot() とかかな．
library(lattice) で検索すると barchart.table() (総称関数 barchart) …… まだインストールしてないけど， library(vcd) にはいろいろな関数が用意されてるみたいだ． vcd とは visualizing categorical data とのこと．
インストールしてないと思ってたら，すでに入ってた． R-help HTML 版 page も生成されてるし． demo(package = "vcd") もいろいろ準備されてるな．これはタイタニック号遭難事件時の，客室等級・男女別・年齢ごとの生死．

ちょっと関係ないけど (いや，大いに関係あるか?) …… RGM でいろいろ検索してたら，量的食物網 (quantitative food web, QFW) 作図の関数が見つかった． library(bipartite) の plotweb() とか．これはすごい．以前，苫小牧にいた平尾君も同様の関数を作ってはいたのだが (ぎょーむ日誌 2005-03-17)．

1730 研究室発．今日は自分の晩飯だけでいいので， AEON 桑園によってあれこれと買いもの …… 時間かかってしまった． 1840 帰宅．晩飯の準備．かんたん晩飯．
明日は北海道若手生態学研究会なので弁当の準備とか．
[今日の運動]
- 朝のラジオ体操．
[今日の食卓]
- 朝 (0700): 米 0.5 合．ジャガイモ・タマネギ・ニンジン・シイタケ・鶏肉のクリームシチュー．
- 昼 (1230): 研究室．米 0.7 合．鶏肉・シイタケの炒めもの．レンコンさつまあげ．インゲン．セロリ．リンゴ．
- 晩 (2030): スパゲッティー．ジャガイモ・タマネギ・ニンジン・シイタケ・鶏肉のクリームシチュー．ヨーグルト．

本日 (kubolog20120224) | 次の日 | 1 日前 | 7 日前 | 31 日前 | 365 日前 | top

KuboLog | KuboWeb