PDF出力の日本語文字

PDF出力で題名とかに日本語を入れているとうまく表示されないでいたのであきらめていたのだが、調べてみると対処法があるようで、一番簡単な方法は、

pdf(family="Japan1")

とすればいいだけのようだ。
実際、試したら問題なく日本語表示された。
やったね。
詳しくは、久保さんのサイトで
http://hosho.ees.hokudai.ac.jp/~kubo/ce/RgJpFont.html

RODBC

今更ながら、エクセルをRに取り込むパッケージであるlibrary(RODBC)は使えると実感。

自分のデータを解析するには、R用にデータを入力してcsvにしてしまうのだが、
最近、人様の入れたデータを扱う用件がいくつかあり、そのときに複数のワークシートをそれぞれcsv化しなくてすむので、ファイル数をコンパクトにまとめるたまま作業することが可能だ。

ただ、人様のデータはR用にはなっていないので、大抵どこかが呪われているので、ある程度エクセルでシャナクを唱えながらの作業になってしまう。

エクセルは、呪いがあって厄介な部分もあるが便利でもあるので不滅だ。

GUIでsource()

source()で、読みたいファイルを指定するのが、作業ディレクトリが違うと絶対パスを書かないといけないのが結構面倒くさい。

しかし、該当ファイルをR console上にドラッグ＆ドロップをすると自動的にsource()が実行され、ファイルが読み込まれる。

あんだけCUIなRでこんなことができるとは思わなかった。
最近知った意外な事実。

factor()の順番

カテゴリー変数の場合、カテゴリーの順番はアルファベット順になってしまうのだが、引数levelsで順番が指定できるようだ。これを使えば箱ひげ図や棒グラフが思った順番に並べられるのか。知らんかった。

(fac<- factor(c("syo","tiku","bai")))
[1] syo  tiku bai
Levels: bai syo tiku

(fac<- factor(c("syo","tiku","bai"),levels=c("syo","tiku","bai")))
[1] syo  tiku bai
Levels: syo tiku bai

lmer研究7 プログラム改訂

一年前、公開したlmer()自動Akaike weight計算プログラムですが、lmerの仕様が変わったため使えなくなっていた。

原因を探ったら、意外と簡単に解決できたので、新しいプログラムを載せておく。

「model_average_lmer_new.r」をダウンロード
ランダム要因が一つのもの用

「model_average_lmer2_new.r」をダウンロード
ランダム要因が二つ以上のもの用

これを作ってもまた仕様がかわるかもなあ。
説明変数シラミつぶし関数をつかったものに変えた方がいいのかも。

MCMCglmm

今年はじめ、MCMCglmmなるパッケージがリリースされました。

今後、kubowebでも詳しく紹介されてくると思いますが、その前に軽く触ってみる。

id <- factor(rep(1:20, rep(5, 20)))
y <- rpois(100, lambda = rep(runif(20), rep(5, 20)))
x <- rnorm(100)
dat <- data.frame(y = y, x = x, id = id)

mod<- MCMCglmm(y~x,random=~id, family="poisson",data=dat
,nitt=15000, thin=10, burnin=5000)

plot(mod$Sol)

plot(mod$VCV)

summary(mod$Sol)
Iterations = 1:1000
Thinning interval = 1
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
(Intercept) -1.01394 0.2257 0.007137        0.02015
x           -0.08527 0.1976 0.006248        0.02585

2. Quantiles for each variable:

               2.5%     25%     50%      75%   97.5%
(Intercept) -1.4389 -1.1746 -1.0222 -0.85541 -0.5704
x           -0.3692 -0.2212 -0.1500  0.04049  0.3296

こんな感じで、ほぼlme()と同じような記法で、MCMC計算ができる。
ただ、上の図であるようにいいデータじゃないとなかなか収束しない。
iteration数を増やすと、計算が落ちてpriorを指定しろと言われるので、事前分布も指定しなくてはいけないのかも。

分布形は、正規分布、ポアソン分布、カテゴリカル分布、多項分布、指数分布、zero infrated ポアソンcengaussian、cenpoisson、cenexponentialが可能。
cen~は、 censored and the prefixな分布らしい。
二項分布は多項分布の一部として扱われ、family="multinomial2"と指定してあげる。
また、バイナリデータ（0, 1）は、カテゴリカル分布で指定する。
分布形の指定は、family=　のあとは"poisson"のように""でくくらなくてはいけない。

AICやBICは出力はされない。anova()も使えないので、いわゆるp値をだすことはできないので、そっちを期待していた人は残念でした。
WinBUGSででてくるDICやRhatもでてこないので、DICはいいとしてもRhatは自力で計算コードを書くべきなのか？

チェーン数は指定できないので、３回計算して重ねてプロットするしかない？

カエルベイズ本を読んでなんとなくベイズ推定できるような気がしてきたけど、自分でモデル式を書くことができず、RのGLMMならできるという人はお手軽にMCMC計算ができるパッケージだと思います。

もっとも、今回のようになかなか収束しない場合は結局WinBUGSの道に引き込まれるのかもしれない。

自然対数・・・ではない

R（確かエクセルでも）で、数字の0が多かったり小数点以下の0が多かったりすると、2e+05とか2e-05とか表示される。
これはそれぞれ200000と0.00002なのだが、このeってどういう意味なんだろう？

殿はこれを自然対数と思い、2e+05を2*(2.72)^5と計算するのでは？と言われてしまった。

確かにeって自然対数だよなあ。

WinBUGSとRのmatrix記法

カエルベイズ本のコードを実行していてWinBUGSとRで共通のようで共通でないコードを一つ見つけた。

structure()なのだが、

WinBUGSだと、

Y1=structure(.Data = c(
1, 1, 1, 1, 0, 0, 0,
1, 1, 0, 1, 1, 1, 0,
1, 1, 0, 0, 0, 0, 0,
1, 0, 1, 0, 0, 0, 0,
1, 0, 0, 0, 0, 0, 0),
.Dim=c(5,7))

こう入力すると、５行７列の見たままの行列になる。つまり、Rでいうmatrix(,byrow=TRUE)と一緒。

これは、Rでもエラーなくコピペできるのだが、これが落とし穴。上のコードで入力すると

    [,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,]    1    0    1    1    0    0    0
[2,]    1    0    1    0    1    0    0
[3,]    1    1    1    0    0    0    0
[4,]    1    1    0    0    1    1    0
[5,]    0    0    1    0    0    0    0

となってしまうので、カエルベイズ本のコードにstructure()で始まる行列のコマンドをRに移すときは、

Y2=matrix(c(
1, 1, 1, 1, 0, 0, 0,
1, 1, 0, 1, 1, 1, 0,
1, 1, 0, 0, 0, 0, 0,
1, 0, 1, 0, 0, 0, 0,
1, 0, 0, 0, 0, 0, 0),
nrow=5,byrow=TRUE)

と書き換えてあげなくてはいけない。

Bayesian Methods for Ecology 著者：Michael A. McCarthy 販売元：Cambridge University Press Amazon.co.jpで詳細を確認する

glm() 自動Akaike weight プログラム

生態屋さんで統計やっている人が一度は見たことあるであろう有名サイト『飯島の雑記帳』

ここで、glm()のakaike weightを求めるプログラムが配布されている。

しかし、自分はこのプログラムを試してみたことがあるのだが、エラーを吐き出してうまく動かなかった。特に使わなかったので、長らく放置していたのだが、久しぶりにRネタの記事を書くためにちょっと書きなおしてみる。

といっても使うのはTaglibro de Hにある説明変数シラミつぶし関数である。

「glm_aw.r」をダウンロード

↑をクリックしてファイルダウンロードしてください。

使い方
glm.aw(y, indeps, dis, offset=NULL,data=NULL)
y=目的変数
indeps=説明変数　c("x1","x2")のように入れていく
dis=分布形
offset=offset項（なくてもよい）

例

set.seed(13)
n <- 40
x1 <- runif(n, 0, 1)
x2 <- x1 + rnorm(n, 0, 0.2)
x3 <- x1 + rnorm(n, 0, 0.2)
x4 <- x1 + rnorm(n, 0, 0.2)
x5 <- x1 + rnorm(n, 0, 0.2)
z  <- x1 + rnorm(n, 0, 1)　#offset用
y  <- x1 * 10 + rnorm(n, 0, 1)
dat<- data.frame(y, x1, x2, x3, x4, x5, z)

result<- glm.aw(y, c("x1","x2","x3","x4","x5"),gaussian,z, dat)
head(result)

   as.character.tm.  vec.aic     dAIC akaike.weight
2            x1 + 1 144.7587 0.000000    0.23329566
4       x1 + x2 + 1 145.9729 1.214172    0.12713133
10      x1 + x4 + 1 146.2318 1.473145    0.11169080
6       x1 + x3 + 1 146.7564 1.997742    0.08592162
18      x1 + x5 + 1 146.7587 1.999999    0.08582473
12 x1 + x2 + x4 + 1 147.4386 2.679901    0.06109047

2008.1.28　追記
引数dataを追加。

パッケージchron

時系列というわけではないが、日時の期間を計算するときのRのパッケージを探したら、やはり願いがかなうR。そんなパッケージがありました。

chronパッケージは、月/日/年と時間を入力してその期間の計算が簡単にできるパッケージだ。

example(chron)のコードを流用。

dts <- dates(c("02/27/92", "02/27/92", "01/14/92",
                "02/28/92", "02/01/92"))
tms <- times(c("23:03:20", "22:29:56", "01:03:30",
                "18:21:03", "16:56:26"))
x <- chron(dates = dts, times = tms)
x
[1] (02/27/92 23:03:20) (02/27/92 22:29:56) (01/14/92 01:03:30)
[4] (02/28/92 18:21:03) (02/01/92 16:56:26)
x[1]-x[2]
[1] 00:33:24
diff(x)
Time in days:
[1]  -0.02319444 -44.89335648  45.72052083 -27.05876157

hours(x)
[1] 23 22  1 18 16