RiboSortの使い方

まずダミーデータ を用意。
１列目がフラグメントの長さ、これが１OTUにあたる。
2列目から７列目がそれぞれのサンプルに見られたフラグメントサイズとピークの高さ。
これをサンプルごとにテキストファイルに分割する。

dat<- read.csv("RiboSortsamples.csv")
library(RiboSort)
for(i in 1:6){
write.table(na.omit(dat[,c(1,i+1)]),paste("sample",i,".txt",sep=""),col.names=F,row.names=F)}

本来は、このサンプルごとのテキストファイルをひとつずつ用意して解析を始めることになると思う。
テキストファイルは、１列目がフラグメント長、2列目がピークの高さ。ラベルはいらない。

これでデータの準備ができた。
次にRiboSortを使ってソーティング。

まずは、データのリストをつくる。

file.list<- character(6)
for(i in 1:6) file.list[i]<- paste("sample",i,".txt",sep="")
file.list
[1] "sample1.txt" "sample2.txt" "sample3.txt" "sample4.txt" "sample5.txt"
[6] "sample6.txt"

これがRiboSortが他のR関数にない特有な操作。
データそのものを事前に呼び出してオブジェクト化するのではなく、先にリストを作って後でデータを呼び出すことになる。

そしてRiboSort()を実行する。
まずは引数の説明。
dataformat:データの出所。手動でつくったファイルなら"standard"。ABIのGeneMapperの出力ファイルなら"abisingle"または"abimultipel"、Beckmanなら"beckman"。デフォルトは"standard"

dye：プライマーの蛍光色素の種類？datafomartが"standard"なら関係なし。

output:abundanceの出力形式。個体数というか量のデータそのものなら"abuncances"、割合データなら"propotions"。デフォルトは、"propotions"。（MDSでは割合データを使うため）

repeats：１サンプルにつき何回繰り返しシークエンサーに流したか。繰り返したサンプルのデータの扱いは、mergerepeatsで設定

mergerepeats:繰り返しデータの扱い方。"none"は繰り返しをひとつのサンプルとして扱う。"presentinall"は、すべての繰り返しで現れたフラグメント長のものだけ採用。"presentintwo"は二つ以上の繰り返しで出たものを採用。"presentinone"は全部採用。

今回は、繰り返しなしの手動データなので

sorted<-RiboSort(file.list)

ちょっと時間がかかった後にsortedにソーティングされたデータが組み込まれる。
さらに、作業ディレクトリに４つのエクセルファイルが生成される。
Information File.xlsと、Abundances Output.xlsと、Proportions Output.xlsと、Ribotypes output.xls。内容はファイル名から大体想像ができると思う。

しかし、このエクセルファイルが厄介で、ファイルを開いてもうまく形になっていないので、区切りを変える必要がある。
Infomation File.xlsは、まず、シート左上の[A]をクリックして選択した状態で、メニューのデータ→区切り位置→「カンマやタブなどの・・・」を選択して次へ→その他を選択して":"を入力して→完了
さらにシート左上の[B]を選択して、上と同じ要領で今度は、「スペース」を選択して完了
それ以外のファイルは、[A]を選択して、上と同じ要領で「スペース」を選択して完了
こんな面倒なことしないで、そもそもCSV形式で出力すればよいものを・・・

で、ソーティングが終わったデータは、目的の解析に応じたデータ形式にしてやればいい。
RiboSortには座標付の手法であるMDSがある。
SampleをグルーピングするsamplesMDS()とOTUごとをグルーピングするspeciesMDS()。
各サンプルが一つの調査点だったとすると

samplesMDS(sorted)

で調査点の傾向が視覚的にわかる。

RiboSort

微生物群集の分析方法として、ARISAやTRFLPなどがある。

これらのデータをシークエンサーから抜き出し統計解析するとき、データをソーティングしてあげないと解析ができない。

しかし、このソーティングが結構面倒で、１サンプルに２００OTU以上あるなかでサンプルどうしを比べて、重なっているOTUがあるかどうか見てあげて、リストをつくらなくてはいけない。

これを自動的にやってくれるプログラムがRで公開されている。
それがRiboSort。

このプログラムは、ABIのGeneMapperかBeckman CoulterのCEQ 8000 Genetic Analysis Systemから吐き出されたデータをソーティングしてくれるほか、手動で入力したデータもソーティングしてくれる。

手動の場合は、一列目にフラグメントの長さ（OTU）、二列目にピークの高さ(abundance)を入れる。
ARISAの場合、abundanceを波形の面積で計算する場合もあるが、手動で二列目に入れればよいのだろう。

あとは、その気になればDGGEのデータに関しても応用可能だろう。
画像解析でバンド間の距離とバンドの濃さを計算すれば、それがフラグメントの長さとabundanceに相当するはずなので。

また、RiboSortパッケージにはmulti-dimensional scalling (MDS)による座標付けをする関数も含んでいる。これをつかわずとも、PCAやクラスター解析、ANOSIM、adonis()等のRで解析可能な解析も当然できるようになる。

8月・9月のＲ本

Rjpwikiに新しいR本情報が載っていますね。

『一般化線形モデル入門　第2版』(R 利用のガイド付き)（共立出版）8月発売

ＧＬＭの理論的な解説をしてくれている教科書ですね。

An Introduction to Generalised Linear Models (Texts in Statistical Science.) 著者：Annette J. Dobson 販売元：Crc Pr I Llc Amazon.co.jpで詳細を確認する

↑の本の和訳のようです。

アマゾンにレビューが書かれています。

この本、原著では第三版もあるようです（なぜだかアフィリエイトが貼れない）。

一般線形モデルによる生物科学のための現代統計学―あなたの実験をどのように解析するか 著者：Alan Grafen,Rosie Hails 販売元：共立出版 Amazon.co.jpで詳細を確認する

共立出版からはこっちも出てますよね。
こっちは、正規分布のモデルのみの話ですが、基礎から書かれているようなので、この本を読んでから上の本を読むと理解が一層深まるのではないでしょうか。

『ベイズ統計データ分析 ―R & WinBUGS― （統計ライブラリー）』（朝倉書店）9月発売。

朝倉書店からまたまたＲ＆ベイズ本が出るんですね。

この手の本が出てくるのは非常にうれしいのですが、
すでに自分はベイズ関連の本は５冊もっているのに、
まったくベイズが使えない（あまり一生懸命勉強していない）ので、この本はいったん保留だな。
出てから評判を聞いてから買うかどうか考えよう。

ZIPな動物

Zero Infrated な分布をする動物にはどんなのがいるだろうか？
まず、考えられるのは、数が少ない希少種。探してもなかなか見つからなければ、０ばかりが並び、Zero Infratedになるだろう。
あと、遠洋研のPDFでは、混穫ってのも挙げられている。

他にもなにかないかと考えると意外と身近にそんな動物たちがいる。
全体から見ればたくさんいて、多様だけど、非常に分布の偏りが大きい動物群集。
別名貧者の熱帯林。土壌動物群集だ。

ちょうど、土壌動物学の国際誌であるPedobiologiaにそんな論文が出版された。
The excess-zero problem in soil animal count data and choice of appropriate models for statistical inference Volume 52, Issue 1, 30 June 2008, Pages 1-17

この論文では、ザンビアのMiomboとAgroforestryで土壌動物を採ったデータを用いて、その動物たちがどんな分布をしているかを、
log変換して最小二乗法(Ordinary least squares regression)、線形混合モデル(Linear Mixed Model)、ポアソンモデル(Poisson)、擬似ポアソンモデル(Poisson with correction for Overdispersion←quasi-likelihoodをつかうやつ)、負の二項分布(Negative Binomial)、ZIP(Zero Inflated Poisson)、ZINB(Zero Inflated Negative Binomial)で比べてAICcとBICで当てはまりを調べている。

MiomboとAgroforestryでそれぞれ、ミミズ、ヤスデ、ムカデ、甲虫、アリ、シロアリの数を調べた結果、
MiomboのミミズとAgroforestryのヤスデ以外はみんな負の二項分布が一番あてはまりがよく、MiomboのミミズはZIP、AgroforestryのヤスデがZINBに当てはまっていた。

という内容です。

やはり土壌動物はポアソンにはあまり当てはまらなくて、負の二項分布の方が当てはまりがいいわけね。
これは、自分の解析の経験と一致する。

ZIP、ZINBなどは、大別した分類群ごとではなりにくいけど、種ごと（特にササラダニとかトビムシ）になると結構この分布に従うものがいるのではないかな。

あと論文の中では、それぞれの確率分布で説明変数のＰ値や信頼区間がどう変わるかとか推して測るべしな内容が述べられています。

最終的ににこの論文の言いたいことは、
「土壌動物はヘテロに分布するから変数変換＆ANOVAやノンパラじゃだめですよ」
ということです。

惜しむらくは、この論文はSASを使って解析をしている。
Appendixでご丁寧にSASコードまで載せている。
これがRコードだったらもっとおもしろかったのになあと思います。

個人的には、なにか適当にいいデータがあれば、同じようなことをＲでやってEdaphologiaに日本語で投稿したら載るんじゃないかと思ってしまった。

hurdle model

昨日紹介した遠洋水研外洋資源部鯨類管理研究室の岡村寛さんのありがたいＰＤＦにhurdle modelというのも紹介されています。

遠洋水産研究所といえば、資源選択検定のレビューとＲコードを書いてくれた人達ではないですか。
これもかなりありがたいです。

話をもとに戻して、
Hurdle modelもpsclパッケージに扱われています。

以前KuboLogにて0の切れたポアソン分布の話題が出ていました。

これもHurdle modelの一種のような気がします。

KuboLogではVGAMパッケージのvglm(,family=pospoisson)で解いているが、

psclパッケージならhurdle()で解けるのかな？
hurdle()は、0の部分とそれ以外の部分を別々に推定してくれるので。

というわけでKuboLogコードをパクって検証。

v.sample <- rpois(1000, exp(-1))
y <- v.sample[v.sample > 0]
cat("# length(y) =", length(y), ", mean(y) = ", mean(y), "\n")
# length(y) = 323 , mean(y) =  1.235294
fit1 <- vglm(y ~ 1, family = pospoisson())
fit2 <- hurdle(y ~ 1)
以下にエラー hurdle(y ~ 1) :
  invalid dependent variable, minimum count is not zero
fit3 <- hurdle(v.sample ~ 1)
summary(fit1)

Call:
vglm(formula = y ~ 1, family = pospoisson())

Pearson Residuals:
                 Min       1Q   Median       3Q    Max
log(lambda) -0.46948 -0.46948 -0.46948 -0.46948 7.5117

Coefficients:
               Value Std. Error t value
(Intercept) -0.82411    0.11102  -7.423

Number of linear predictors:  1

Name of linear predictor: log(lambda)

Dispersion Parameter for pospoisson family:   1

Log-likelihood: -136.0555 on 322 degrees of freedom

Number of Iterations: 5

summary(fit3)

Call:
hurdle(formula = v.sample ~ 1)

Count model coefficients (truncated poisson with log link):
            Estimate Std. Error z value Pr(>|z|)   
(Intercept)  -0.8241     0.1110  -7.423 1.15e-13 ***
Zero hurdle model coefficients (binomial with logit link):
            Estimate Std. Error z value Pr(>|z|)   
(Intercept) -0.74002    0.06762  -10.94   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Number of iterations in BFGS optimization: 6
Log-likelihood:  -823 on 2 Df

ごめんなさい。まったくの別物でした。

で、ここまで書いて、統計学なんでもあり掲示板の５２７２番からその一連のレスを読むともうすでにこのことについて議論されたあとでしたね。

完全に恥さらしなので、全部消そうかと思ったけど、ここまで書いてしまったのでこのままアップすることにします。

で、Hurdle modelっていったい何なんだろう？

０とそれ以外を別個に計算するのはZero Inflated modelも同じなわけだし。

雰囲気としてHurdle modelとZero Inflated modelはお仲間で計算方法が違うものという認識でいいのかな。

overdispersionが0に引っ張られるときに、0がすごくたくさんある時は、Zero Inflated Modelで、そうでなさそなときはHurdle Modelということか？

zero inflated model

ZIP&ZINBなGLMのまとめ

ネットでGLMの記述といったら、まずはKuboWebですね。
KuboWeb内を検索すると、zero inflatedなRパッケージにzicountsというのがあるのですね。
しかし、このzicountsパッケージCRANで調べてみると・・・

Package ‘zicounts’ was removed from the CRAN repository.

Formerly available versions can be obtained from the <A HREF="../../../src/contrib/Archive/zicounts>archive.

すでにCRANから削除されているようだ。
archiveには残っているようなので、そこに行ってみると、ソースがおいてあるのですが、自分はコンパイルのやり方がわからないので使用不能。

で、次にZIGP
これ、インストールしてexample()使って試してみたのですが、イマイチ使い方がわかりません。

他に日本語でZIP modelを解説しているものがないかと探すと、こんなＰＤＦ(40-43ページ目)がありました。

このＰＤＦで紹介されている、現在使用に耐えるパッケージのpscl。
このパッケージには説明PDFがあります。
このPDFの2章Models and softwareは、カウントデータを扱うGLM, GAMのRパッケージの説明が書かれているので、参考になります。

ちなみにこのPDFにはzicountsとZIGPについても言及されており、「使いづらいからオススメしません。」とのことです。

また、上記のＰＤＦにはZero infrated modelが扱えるパッケージがあと二つ紹介されている。

VGAMパッケージとgamlssパッケージだ。
二つとも、非常に盛りだくさんなパッケージで、線形モデルだけでなく非線形モデル（GAM）も扱えるみたいです。名前からしてGAMがメインなのでしょうが。
二つのパッケージとも奥が深いので名前を挙げておくだけにします。

まとめると、線形でZIPを使うならとりあえず、psclパッケージのzeroinfl()を使えばよさそうだ。

非線形まで考えるならVGAMかgamlssのどちらかのパッケージを使えばいいのかな。

関数を実際に走らせて日記書くのはもうちょっと勉強してから出直そう。

lmer研究（6）lme4 Version: 0.999375-17

lme4パッケージが更新されました。

今回は、大幅に仕様がかわったようです。

・計算アルゴリズムを改善して、より頑健になり、計算が早くなった。
・PQLが使えなくなった。
・引数のmethodがなくなった。
・glmer()という関数ができた。
・mcmcsamp()のサンプリング方法が変わった。

今回の更新で我々エンドユーザーが特に注意しなければいけないのは、methodでしょう。

尤度を求めるmethodは、正規分布の場合"REML"と"ML"、それ以外の分布では、"Laplace"と"PQL"でした。
しかしながらPQL法は、問題があり、あまり正確な推定をしないことが指摘されています。方法がそれしか使えなかった時代ならいざしらず、現在はglmmMLやlmerのようにLhaplaceやghq(glmmMLのみ)のようにPQLより正確で尤度を求める方法が実装されているので、あえてPQLを使う必要はないと思います。そのため、今回のバージョンでPQLが使えなくなったのだと思います。

で、正規分布のときREMLとMLはどうしていするかというと、
REML=TRUE
で指定するようです。デフォルトは、TRUE。FALSEにするとMLになります。

現状では、まだ引数methodは使えますが、warning messageがでます。

fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy,method="ML")
Warning message:
In lmer(Reaction ~ Days + (Days | Subject), sleepstudy, method = "ML") :
  Argument ‘methood’ is deprecated.  Use ‘REML’ instead

おそらく、次に出るバージョンでは引数methodは使えなくなるので、もう使わない方向で行きましょう。

で、もうひとつglmer()について

結論からいうと現状では、あまり気にする必要はないと思います。

lmer()はLinear Mixed Model、つまり、正規分布のGLMM、
glmer()はGeneralized Mixed Model、正規分布以外のGLMMで、両者の関数を使い分けるようです。

しかしながら、lmer(,family=poisson)とglmer(,family=poisson)は全く同じ結果になります。

test1 <- lmer(round(Reaction) ~ Days + (Days|Subject),family=poisson, sleepstudy)
test2 <- glmer(round(Reaction) ~ Days + (Days|Subject),family=poisson, sleepstudy)
summary(test1)
summary(test2)
（出力結果略）

では、なんでglmer()がでてくるかというと、引数nAGQの存在です。

これは、よくわかりませんが、Gauss-Hermite approximationにおいて設定する値のようです。整数値を入力、デフォルトは1。

Gauss-Hermite approximationは、glmmML（そういえば、最近更新したようですね。なにがどうなったかわかりませんが）のmethod="ghq"でやる方法ですね。glmmMLのヘルプには、"Both methods are now fully adaptive."だそうです。以前のバージョンではちょっと書いてあることが違ったのですが。

で、もともと以前のlmerのバージョンでは、「Lhaplaceよりも正確な、method="AGQ"が選べるようになりますよ。まだ実装してないけど」と書いてあった気がするけど、現バージョンで、引数nAGQという形で実装されたのでしょう。

ということでnAGQをいじってみる。

test3<- glmer(round(Reaction) ~ Days + (Days|Subject),family=poisson,nAGQ=2, sleepstudy)
以下にエラー mer_finalize(ans, verbose) : Code not yet written

まだ実装されていないようです。

結論

現バージョンでは、lmer()を使って引数methodについて気にすれば（多分）問題なく使える。

過去のlmer研究
（１）自動akaike weightプログラム
（２）自動akaike weightプログラム２
（３）gaussian modelのp値
（４）AICについて
（５）lmerで多重比較

aod研究（２）negbinシラミつぶし

Taglibro de Hにおいて、説明変数シラミつぶし関数というのが公開されています。

これは便利な関数なのでstepAICが使えないGLMM関数に応用しようと思う。
lmerに関しては、akaike weight自動化関数を作っているので、aodパッケージのnegbin()でやってみる。

このページのデータをちょっと改変し、negbin仕様にコードを改変し、さらにoffset項までつけてみた。

library(aod)
set.seed(13)
n <- 40
x1 <- rnbinom(n, 10, mu=5)
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<- c(rep("A",10),rep("B",10),rep("C",10),rep("D",10))
y  <- x1 * 10 + floor(rnorm(n, 0, 1)*5)
dat<- data.frame(x1,x2,x3,x4,x5,y,z)

indeps <- c("x1", "x2", "x3", "x4", "x5")
tm <- terms(indeps)
lst.negbin <- vector("list", 2^length(indeps))
vec.aic <- vector("numeric", 2^length(indeps))
for (i in 1:2^length(indeps)) {
f <- as.formula(paste("y ~", tm[[i]],"+offset(x1)"))
lst.negbin[[i]] <- negbin(f, random=~z,data=dat)
vec.aic[i] <- AIC(logLik((lst.negbin[[i]])))[1]
}
summary(lst.negbin[[which.min(vec.aic)]])

あえて結果は載せませんが、ちゃんと動きました。

AICが低いモデル順に並び替える。

res<-data.frame(AIC=vec.aic, term=t(as.data.frame(tm))[,1])
rownames(res)<- 1:nrow(res)
result<- res[order(res$AIC),]
rownames(result)<- 1:nrow(result)
result$dAIC<- result$AIC-min(result$AIC)
result

        AIC                       term        dAIC
1  295.7017                     x1 + 1   0.0000000
2  296.3513                x1 + x2 + 1   0.6495105
3  296.9963                x1 + x3 + 1   1.2945497
4  297.4092                x1 + x4 + 1   1.7074249
5  297.4898                x1 + x5 + 1   1.7880303
6  297.6380           x1 + x2 + x3 + 1   1.9362689
7  298.0365           x1 + x2 + x5 + 1   2.3347773
8  298.0886           x1 + x2 + x4 + 1   2.3868343
9  298.7274           x1 + x3 + x4 + 1   3.0256804
10 298.9093           x1 + x3 + x5 + 1   3.2075662
11 299.0251           x1 + x4 + x5 + 1   3.3234029
12 299.4915      x1 + x2 + x3 + x5 + 1   3.7897122
13 299.6038      x1 + x2 + x4 + x5 + 1   3.9020230
14 299.7911      x1 + x2 + x3 + x4 + 1   4.0893903
15 300.8139      x1 + x3 + x4 + x5 + 1   5.1121448
16 301.1139 x1 + x2 + x3 + x4 + x5 + 1   5.4121685
17 307.4005           x2 + x3 + x4 + 1  11.6987170
18 307.5981                x2 + x4 + 1  11.8963951
19 308.2918           x2 + x4 + x5 + 1  12.5900286
20 309.1244      x2 + x3 + x4 + x5 + 1  13.4226783
21 310.8429                x2 + x5 + 1  15.1411847
22 311.0418           x2 + x3 + x5 + 1  15.3400820
23 311.4495                x2 + x3 + 1  15.7477459
24 312.6703                x3 + x4 + 1  16.9685468
25 313.4375           x3 + x4 + x5 + 1  17.7357443
26 314.3119                x4 + x5 + 1  18.6101821
27 316.7691                     x2 + 1  21.0673725
28 318.1372                     x4 + 1  22.4354743
29 318.6671                x3 + x5 + 1  22.9653169
30 321.8581                     x5 + 1  26.1563971
31 325.7885                     x3 + 1  30.0867493
32 451.5478                          1 155.8460312

過去のaod研究
（1）aodパッケージのakaike weight

モデルからインパクトを予測する

このエントリーがちょうど１００エントリー目です。

研究室の先輩の統計コンサル

GLMsでモデルが得られたあとにそれをどうするか？という話になった。

モデル選択で

y~ 0.33 A + 0.2 B - 0.11 C +  2

というようなモデルが得られたら、
「ああ、AとBがプラスに効いててCがマイナスに効くんだ。」

というようなところで終わっていたのだが、保全を主目的にしている先輩は、

「実際に、ある条件（変数）が変動したときはどうなるか？」

という、予測までしたいのだそうだ。
保全の現場では、よくありそうなことで、B,Cの要因（もう撹乱域が決まっているとか）はどうしようもないけど、Aの要因はなんとか変えることができる（植物植えるとかエサ撒くとか）からAを変化させることによってどれだけ保全効果を得られるかを予測したい場合こんなことをしたいわけだ。

そうなると、予測値を求めることになるので、predict()なり、fitted()なりを使えばいい。

しかし、今回は事情が異なり、

「B,Cは固定して、Aだけが変動したときの予測値を得たい」

このような場合は、predict()では対応できない（すくなくとも自分が試した限りでは）。

先輩の話では、JMPやエクセルならできるらしい。

けれども、Aをたとえば１から５０まで動かしたとすると、マウス操作とかを５０回繰り返さなきゃいけない（もしくは、ひたすら長い行をドラッグ）わけで、そんなことをしていたらどこかでミスるか、呪われると思う。

じゃあ、ちょっと遊びでRを使ってやってみるかと思い立ち、仮想データを作ってやってみることにした。

長いので、続きに入れます。

続きを読む "モデルからインパクトを予測する"

マルコフ連鎖モンテカルロ法

マルコフ連鎖モンテカルロ法 販売元：朝倉書店 Amazon.co.jpで詳細を確認する

期待の本が発売しました。

East_Scrofaさんによると、もともとこの本はアマゾンで扱う予定はなかったようですが、発売日の２５日になって、いきなり購入可能になりました。

午前中にアマゾンで確認したときは、検索でヒットしなかったのですが、午後になって注文可能に。

East_Scrofaさんがアマゾンに問い合わせてくれたり、自分がアマゾンでこの本の検索をかけたことが効いているのだろうか？
いずれにしてもよいことです。

即注文しました。

楽しみです。

はやく来ないかな～

Tree diversity analysis

RjpWikiを眺めてみたら、パッケージ更新情報にBiodiversityRがあったので覗いてみると、これ関係の本が出版されていることに気づいた。

RjpWikiのR史、R本リストにも載っていないし、他のページでも見たことがなかったので、この存在には気付かなかった。

調べてみると2005年に出版されているにも関わらず、amazonの日本版にも本家にも検索しても出てこない。
当然、ココログのアフィリエイトも全滅だった。

不思議だなと思ってよくよく見ていると、この本は、PDFで内容のすべてがダンロードできるようだ。

対象は樹木群落だけれども他の群集にも応用可能だと思うので、BiodiversityRやveganを使うときはこの本を参照するといいかも。

特にこのブログにANOSIMと検索してやってくる人がわりかしいるようなので、そんな人はこの本の８章をよんでみると参考になるかもしれません。

lmer研究（５）&multcomp研究（２）lmerで多重比較

まず、TukeyのHSDなどの多重比較は、正規分布のモデルしか使えないはず。それについてmultcompに関しての論文のp7にgeneralized linear modelとmixed modelについて書かれてある。

読む限りだと、他の分布でも使っていいようなことが書いてある。

また、mixed modelの場合は、REMLならばよい。と書いてあるように見えるけれど、上の論文にlmerで二項分布を使った例がある（検定はしていない）ので、使えないことはないのかもしれない。

一応、Lhaplace近似して尤度の近似値を求めているので、使ってもいいんではないのかなと思う。これは、lmerでAIC水準のモデル選択ができることと同じ根拠。

実際に動かしてみる。

ポアソン分布のデータを作る。

set.seed(813)
count<- c(rpois(30,10) + floor(rnorm(30,0,3)), rpois(30,10) + floor(rnorm(30,0,3))
         ,rpois(30,14) + floor(rnorm(30,0,3)) )
treat<- factor(c(rep("A",30),rep("B",30),rep("C",30)))
site<- rep(c(rep("alpha",10),rep("beta",10),rep("gamma",10)),3)
dat<- data.frame(count,treat,site)

昨日のデータとほとんど変わりません。
AとBが平均１０で、Cが平均14
siteはブロックを想定したランダム変数（平均０分散３）。ブロックは３つ。

正規分布なデータではないので、ほんとはできないのだが、試しにmultcompBoxplot()を使ってみる。

library(multcompView)
multcompBoxplot(count~treat,data=dat)

Aが重なるデータができました。

ここからが本題

library(lme4)
mod.lmer<- lmer(count~ treat+(1|site), family=poisson)
summary(mod.lmer)

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept)  2.27899    0.05842   39.01   <2e-16 ***
treatB      -0.05267    0.08373   -0.63   0.5293   
treatC       0.18938    0.07897    2.40   0.0165 *

これを、Tukeyの方法で多重比較する。

library(multcomp)
tukey1<- glht(mod.lmer,linfct=mcp(treat = "Tukey"))
summary(tukey1)

Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts
Fit: NULL

Linear Hypotheses:
           Estimate Std. Error z value p value   
B - A == 0 -0.05267    0.08373  -0.629  0.8040   
C - A == 0  0.18938    0.07897   2.398  0.0435 *
C - B == 0  0.24205    0.08013   3.021  0.0073 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)

Warning message:
In x$model$call : $ operator not defined for this S4 class, returning NULL

ポアソンモデルでやることにより、きれいにトリートメントA, BグループとCに分けることができた。

上で紹介したmultcompの論文の１８ページにlmer()の二項分布の例があり、glht()を用いて、予測値と９５％信頼区間を推定している。モデル選択の場合、やっていることは推定だと思うので、こっちの方が求めたいもの次第によってはいいのかもしれない。
今回の例でやってみる。

mod.lmer2<- lmer(count~ treat-1+(1|site), family=poisson) #切片なしモデルをつくる
K <- diag(length(fixef(mod.lmer2)))
ci2 <- confint(glht(mod.lmer2, linfct = K))
ci2$confint<- poisson()$linkinv(ci2$confint)
ci2


         Simultaneous Confidence Intervals

Fit: NULL

Estimated Quantile = 2.3878
95% family-wise confidence level

Linear Hypotheses:
       Estimate lwr     upr   
1 == 0  9.7524   8.4817 11.2134
2 == 0  9.2499   8.0147 10.6755
3 == 0 11.7992  10.3929 13.3959

こうして、予測値と信頼区間がもとまりました。
ここで、ひとつ注意。
lmerの場合、fitted()を使って予測値を求めることができる。しかし、glm()の場合fitted()で、リンク関数を外した状態で得られるのに対し、lmer()の場合、出てくる予測値はリンク関数が外されていない。そのため、リンク関数の逆関数（この場合はexp()）をしてあげる必要がある。その行程が、上のpoisson()$linkinv()の関数だ。これで、リンク関数を外した値を得ることができる。

これを図にする。

cci<- ci2$confint
plot(1:3, cci[,1], xaxt="n",ylab="",xlab="",ylim=c(min(cci[,2]),max(cci[,3])))
axis(1,1:3,c("A","B","C"))
arrows(1:3,cci[,2],1:3,cci[,3],length=.05,angle=90,code=3)
text(1:3,13,c("a","a","b"),col=c("blue","blue","red"),cex=2)

　

棒グラフと情報量は変わらないけれど、論文に棒グラフで平均と標準誤差を載せるよりましな気がする。

最後にひとつ注意を。

Tukeyの方法など、p値を求める検定がRのバージョン2.7.0でできなくなりました（今までの結果は2.6.2で求めた）。これは、lmerでやってはいけないからではなく、推測ですが2.7.0では、S4メソッドの何かしらの仕様が変わったからだと思います。
S4のlmer()にS3のglht()を無理やりなんとかしていた感が2.6.2からあった（warningがよく出る）ので、仕方ないのかと思います。
そこのところは、将来的にR側でなんとかするのか、multcompパッケージ側で対応するのだろうと思います。

過去のlmer研究
（１）自動akaike weightプログラム
（２）自動akaike weightプログラム２
（３）gaussian modelのp値
（４）AICについて

過去のmultcomp研究
（１）とりあえず正規分布でやってみる

multcomp研究（１）とりあえず正規分布でやってみる

まずはデータを用意。

set.seed(813)
count<- c(rnorm(30,10) + floor(rnorm(30,0,3)), rnorm(30,10) + floor(rnorm(30,0,3))
         ,rnorm(30,14) + floor(rnorm(30,0,3)) )
treat<- factor(c(rep("A",30),rep("B",30),rep("C",30)))
dat<- data.frame(count,treat)
head(dat,3)
     count treat 
1 10.290013     A
2  9.909496     A
3 15.290336     A

トリートメントがA, B, Cでそれぞれ、平均AとBは平均１０、Cが１４で分散１の乱数を３０個ずつ発生。

まずはmultcompViewのmultcompBoxplot()を使ってみる。

multcompBoxplot(count~treat,data=dat)

　

multcompBoxplot()は、箱ひげ図を描いてくれた上に、TukeyのHSD多重比較の結果も載せてくれる便利な関数。ちなみにこの関数、有意差がない場合、

以下にエラー plot.window(...) :  有限な 'xlim' の値が必要です
追加情報:  Warning messages:
1: In min(x) :  min の引数に有限な値がありません：Inf を返します
2: In max(x) :  max の引数に有限な値がありません：-Inf を返します
3: In min(x) :  min の引数に有限な値がありません：Inf を返します
4: In max(x) :  max の引数に有限な値がありません：-Inf を返します
5: In plot.window(...) :
   "horizontal"はグラフィックスパラメータではありません
6: In plot.window(...) :  "add"はグラフィックスパラメータではありません
7: In plot.window(...) :
   "fontsize"はグラフィックスパラメータではありません
8: In plot.window(...) :
   "fontface"はグラフィックスパラメータではありません

エラーになって図が表示されません。（一瞬でるようですが）
p値を知りたい場合は、

TukeyHSD(aov(count~treat,data=dat))

次にmultcompを使って同じことをしてみる。
glht()関数で第一引数にモデル、引数linfctに比較水準を入れてあげる。linfct=mcp()としてかっこの中に水準を入れる。この場合は、treatを"Tukey"で比較。

library(multcomp)
mod<- glm(count~ treat,family=gaussian, data=dat)
mul<- glht(mod,linfct=mcp(treat = "Tukey"))
summary(mul)
Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts

Fit: glm(formula = count ~ treat, family = gaussian, data = dat)
Linear Hypotheses:
           Estimate Std. Error z value p value   
B - A == 0   0.8179     0.7597   1.077   0.529   
C - A == 0   4.3376     0.7597   5.710  <1e-04 ***
C - B == 0   3.5196     0.7597   4.633  <1e-04 ***

で、multcompはDunnetの方法もできる。
トリートメントAをコントロールとして、BとCを実験群として考える。
まず、Aを基準とした比較ための行列を作ってから、先ほどと同じようにglht()を使う。
c()の並び順に注意。今回の場合はA, B, Cとトリートメントが並んでいるので、その順番を間違えないように気をつける。

cont<- rbind("B-A" =c(-1,1,0),
             "C-A" =c(-1,0,1))
dun<- glht(mod,linfct=mcp(treat=cont))
summary(dun)

Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: User-defined Contrasts

Fit: glm(formula = count ~ treat, family = gaussian, data = dat)

Linear Hypotheses:
         Estimate Std. Error z value  p value   
B-A == 0   0.8179     0.7597   1.077    0.453   
C-A == 0   4.3376     0.7597   5.710 2.26e-08 ***

多重比較

GLMやGLMMで、ANOVAのように多重比較はできるのかと調べてみると、どうやらできるようだというのがわかってきた。

今まで、GLMやGLMMでは多重比較はできないものと思い、調べたい要因が有意かどうか、モデル選択や尤度比検定をしていた。
自分としては、この方法は理解しやすくて結果の解釈も自然にできるので、良かったのだけれども、長年ANOVA&多重比較で論文を通してきたPDの先輩はなかなか理解を示してくれず、うーむと思っていた。

今回、別の大学の先輩と統計の話をしていた時に、もう一度多重比較のことを調べてみようと思い、久保先生のサイトに行ってみて、多重比較とモデル選択のスライドを見直してみることにした。

そこで、multcompなるパッケージの存在に気づく。

いままで、このパッケージはボンフェローニとかシェフェとか旧来の多重比較をすることができるパッケージだと思ってスルーしていたのだが、今回調べてそれは間違いだった。

このパッケージはもともとsimtest()という関数を使って従来の多重比較ができたらしいのだけれども、最新版ではその関数は削除され、かわりにglht()を用いるようだ。

このmultcompパッケージのマニュアルをみると、suggestにlme4とあるので、lmerのオブジェクトも使えるようだ。実際、使った例が示されている。

multcompのヘルプとかマニュアルとか読むと、このパッケージは多重検定をするというより信頼区間をもとめて推定をしている。最近の統計学では、こういう方法が主流なんだろう。

そういえば、Crawley本に「Contrast」という章があった。ここでの話の行きつく先がmultcompによる多重推定なんだろうなと今になって思う。

R（または、近年の統計学）において、「ANOVA&多重比較」のお決まりコンボで使われていた、ボンフェローニやシェフェなどの多重検定はどうやら排除する方向にあるようだ。
そのかわり、モデル選択やmultcompパッケージのような多重推定をすることが、主流になるんだろう。

まとめると、GLMやGLMMをANOVA&多重比較のように使う場合（久保先生の言い方をつかうと）

ヤクザなデータを取っている人は
情報量基準（AICなど）を用いたモデル選択

カタギな人たちは
multcompパッケージを使った多重推定

ということになるのだろうか。

ちなみに、multcompパッケージの拡張として、multcompViewという多重比較の結果を箱ひげ図で表現するパッケージもある。

更新・vegan1.11-4

veganパッケージが20日に更新されたそうだ。

初めて知ったのだが、veganは何を更新したのかがちゃんとwebに上がっている。

今回は、以前このブログでも取り上げたadonis()の致命的なバグを直したようだ。

「adonis()を使っている人はすぐに新しいバージョンを入れて再解析してください」

とのことです。

sem研究(5)AICとCAIC

SEMでは、BICの求め方が公式と異なるので、AICも普通の公式と違うのではないかと思い調べてみる。

と、やはり違うようだ。

AIC= χ二乗値 + 観測変数 * (観測変数-1) - 2 * df

なので（sem.wh.1の例使用）、

summary(sem.wh.1)$chisq + sem.wh.1$n * (sem.wh.1$n - 1) - 2 * summary(sem.wh.1)$df
[1] 25.48505

となる（はず）。

あと、AICとは別にサンプル数を加味したCAICというのもあるようだ。

CAIC= χ二乗値 + (1 + log(サンプル数)) * (観測変数 * (観測変数-1) - 2 * df) / 2

sm1<- summary(sem.wh.1)
sm1$chisq + (1 + log(sem.wh.1$N)) * ((sem.wh.1$n)* (sem.wh.1$n-1) - 2*sm1$df)/2
[1] 60.50905

となる（はず）。

semによるAICとCAICを求める関数を作っておく。

AIC.sem<- function(model){
          sm<- summary(model)
          aic<- sm$chisq + model$n * (model$n - 1) - 2 * sm$df
          return(aic)}

CAIC.sem<- function(model){
          sm<-summary(model)
          caic<- sm$chisq + (1 + log(model$N)) * ((model$n)* (model$n-1) - 2*sm$df)/2
          return(caic)}

使い方

AIC.sem(sem.wh.1)
[1] 25.48505
CAIC.sem(sem.wh.1)
[1] 60.50905

一気に二つの情報基準量を求めたいときはこっち（「sem.aic.r」をダウンロード ）

sem.aic<- function(model){

          AIC.sem<- function(model){
          sm<- summary(model)
          aic<- sm$chisq + model$n * (model$n - 1) - 2 * sm$df
          return(aic)}

          CAIC.sem<- function(model){
          sm<-summary(model)
          caic<- sm$chisq + (1 + log(model$N)) * ((model$n)* (model$n-1) - 2*sm$df)/2
          return(caic)}
         
          arr<-c("AIC"=AIC.sem(model),"CAIC"=CAIC.sem(model))
          return(arr)}

ご参考まで。

モデル選択をAICでやりたい方はご自由にお使いください。

結果は保証しません。

間違っていたら誰か教えてください。

過去のsem研究
(1) 例題実行
(2)記法とか雑多なまとめ
(3)線形方程式系は正確に特異です
(4)BIC

sem研究(4)BIC

semでのモデル選択基準はBICだ。

しかしながら、本とかネットを徘徊する限り、AICやまたCAICなるものも情報量基準として使われている。

なぜ、semではBICしか使われていないのかは、sem開発者であるJohn Fox氏のサイトのsemページに行っても特に何も書いていない。

まず、BICはいったい何なのかを調べてみる。

semパッケージのヘルプにある、『model.wh.1』を例にして、いろいろいじくってみる。（青字はRの出力）

sem.wh.1 <- sem(model.wh.1, S.wh, 932)
summary(sem.wh.1)

Model Chisquare =  13.485   Df =  9 Pr(>Chisq) = 0.14186
Chisquare (null model) =  2131.4   Df =  15
Goodness-of-fit index =  0.99527
Adjusted goodness-of-fit index =  0.98896
RMSEA index =  0.023136   90% CI: (NA, 0.046997)
Bentler-Bonnett NFI =  0.99367
Tucker-Lewis NNFI =  0.99647
Bentler CFI =  0.99788
SRMR =  0.014998
BIC =  -48.051

Normalized Residuals
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max.
-1.26000 -0.13100  0.00014 -0.02870  0.11400  0.87400

Parameter Estimates
     Estimate  Std Error z value  Pr(>|z|) 
lamb   5.36880  0.433982  12.3710 0.0000e+00
gam1  -0.62994  0.056128 -11.2233 0.0000e+00
beta   0.59312  0.046820  12.6680 0.0000e+00
gam2  -0.24086  0.055202  -4.3632 1.2817e-05
the1   3.60787  0.200589  17.9864 0.0000e+00
the2   3.59494  0.165234  21.7567 0.0000e+00
the3   2.99366  0.498972   5.9996 1.9774e-09
the4 259.57583 18.321121  14.1681 0.0000e+00
the5   0.90579  0.121710   7.4422 9.9032e-14
psi1   5.67050  0.422906  13.4084 0.0000e+00
psi2   4.51481  0.334993  13.4773 0.0000e+00
phi    6.61632  0.639506  10.3460 0.0000e+00
                                  
lamb SEI <--- SES                  
gam1 Alienation67 <--- SES         
beta Alienation71 <--- Alienation67
gam2 Alienation71 <--- SES         
the1 Anomia67 <--> Anomia67       
the2 Powerless67 <--> Powerless67 
the3 Education <--> Education      
the4 SEI <--> SEI                  
the5 Anomia71 <--> Anomia67       
psi1 Alienation67 <--> Alienation67
psi2 Alienation71 <--> Alienation71
phi  SES <--> SES                  

Iterations =  78

SEMにおいて、対数尤度はχ二乗値に相当するのでこの場合は、13.485。
自由パラメータ数は、dfにあたるので、9。
サンプル数は、932（semのモデル式を作る時に記入する）。

ここで、BICの公式を見る限り、AICとの大きな違いはサンプル数を組み込むかのようだ。
公式にのっとって、計算してみる。

sm1<- summary(sem.wh.1)
2* sm1$chisq + sm1$df * log(sem.wh.1$N)
[1] 34.56589
sm1$BIC
[1] -48.05094

となり、まったくBICと異なってしまう。

ここで、John Fox氏のsemの解説PDFを眺める（John Fox氏の解説PDFにはこっちもある）と、BICの求め方が書いてあった。

それによると

BIC=χ二乗値　- df ×　log(N)

のようだ。

なるほど。semの場合、χ二乗値が小さいほどよいモデルなのでこの数式には納得。

なので、BICを手動で求めるには

sm1$chisq - sm1$df * log(sem.wh.1$N)
[1] -48.05094

これでBICと一致する。

今回はこれでやめて次回に続く。

過去のsem研究
(1) 例題実行
(2)記法とか雑多なまとめ
(3)線形方程式系は正確に特異です

Crawley本の翻訳

共立出版から来月Crawley本

Statistics: An Introduction Using R 著者：Michael J. Crawley 販売元：John Wiley & Sons Inc Amazon.co.jpで詳細を確認する

↑の翻訳本がでるようだ。

統計学:Rを用いた入門書 著者：Michael J.Crawley 販売元：共立出版 Amazon.co.jpで詳細を確認する

やるな共立出版。

これでやっとこさ、GLMをたくさん扱った日本語のR本がでますね。

参考リンク：Crawly本解説のまとめ

AICの調節

glmmML()とglm()のAIC調節問題のつづき。

East_scrofaさんのデビアンスと対数尤度の解説でわかってきた。

デビアンス(D)＝　2 * （飽和モデルの対数尤度(LI)　-　調べたいモデルの対数尤度(Ls)）

なわけなので、これを展開して変形するとすると

D/2 = LI - Ls
D/2 + Ls = LI
Ls = -D/2 + LI （式１）

となるわけです。

で、glmmMLの場合、対数尤度 は「residual devianceを計算する基準となっているモデルの対数尤度を0として出したもの」
つまり
LI = 0

式１に代入して

Ls = -D/2

となると。

ということで、検証。

まずは、データ（East_scrofaさんのコメントから引用）

alive<-c(2,5,6,9,3)
dead<-c(8,5,4,1,7)
x1<-c(1,2,3,4,3)
ID<- 1:5

では、モデルを作っていきましょう。

library(glmmML)
model.glmmML<-glmmML(cbind(alive,dead)~x1,cluster=ID,family=binomial) #glmmML
model.glmmML2<- glmmML(cbind(alive,dead)~x1,cluster=ID,family=binomial,start.sigma=0,fix.sigma=TRUE) #ランダム切片なしモデル
model.glm<- glm(cbind(alive,dead)~x1, family=binomial) #glm

まずは各AICの比較（青字はRの出力）

mat<- cbind(model.glmmML$aic,model.glmmML2$aic,model.glm$aic)
colnames(mat)<- c("glmmML","glmmML2", "glm")
rownames(mat)<- "AIC"
mat

      glmmML  glmmML2      glm
AIC 11.12108 9.135012 21.63863

やりたいことは、glmのAICとランダム切片なしモデルであるglmmML2のAICを一致させればいいわけです。

そのためには、glm()でもとめたモデルの対数尤度をglmmMLのように調節する。

ここで、Ls = -D/2
で、
AIC = -2 * Ls + 2 * パラメータ数
より、代入すると、

AIC = D + 2 * パラメータ数
になると。

さきほどのglm()のモデルのAICを調節するには

deviance(model.glm) + attr(logLik(model.glm),"df")*2
9.135012

となり、glmmML2のAICと完全に一致するわけです。

非常に簡単ですが、関数化しておきます。

adjustAIC<- function(model) deviance(model) + attr(logLik(model),"df")*2

使い方
model.glm<- glm(cbind(alive,dead)~x1, family=binomial)
adjustAIC(model.glm)
9.135012

glm.nb(), negbin(), betabin()でも使えます。
glmmMLの関数をいじるよりも、いろいろな意味で安定しているglm()関数のモデルをいじった方が調節しやすいですね。

注意点として粕谷先生のブログによるとデータの型によってglmmML()とglm()のAICが一致する場合もあるそうです。

何度読んでもわからないことがあるときはたとわかる時がある。
今回のresidual devianceがそうなわけで、以前からデビアンスに関する記述を読んでもわかったふりして放っておいたんだけど、やっと理解できた。

sem研究（３）線形方程式系は正確に特異です

自分もそうだったし、同じ研究室の研究員の方もsemを使って陥ったエラーに、

以下にエラー solve.default(C) :
   Lapack routine dgesv: 線形方程式系は正確に特異です

というのがある。
全く意味不明な日本語である。英語の直訳なんだとおもうけど。

ということで、だれかが取り上げているかもしれないからgoogle先生に聞いてみる。

最初にヒットしたページのみ一言一句変わらないので行ってみると、なんだかエラー文の英語と日本語の対応したものがずらずら並んでいる。

そこで先ほどの文を探してみると、これが

system is exactly singular

だとわかるわけです。といっても何言ってるのかよくわからないわけですが。

とりあえずエラー文にはsolve.defaultとあるので、方程式をとくsolve()に関係していることが予想される。

まあ、いずれにしてもうまくモデルが組めておらず計算できていないんだと思います（全く答えになってない）。

ちなみに自分と研究員の方はこの問題をどう解決したかというと、

自分の場合：誤差項を指定していなかったのでちゃんと指定してあげた。
つまり
A<->A, parameter, NA
ってやつ。

研究員の方の場合：外性変数の誤差項を修正して、パラメータにするのではなく、
A<->A, NA, 1
のように分散を指定してあげることにより解決し、動くようになった。

とりあえず、これが誰かの役に立てば幸いです。

GLM関数群のAIC算出法まとめ

East_scrofaさんのところでまとめられているように、GLMおよびGLMMの関数群にはAICを求める対数尤度の基準が２通りあるようで、違う求め方で求めたAICは単純に比べられない。そこで、よくつかわれるGLM, GLMMの関数群のAICの求め方を補足としてまとめる。

従来のAICの求め方をしている関数群。
glm(), glm.nb(MASSパッケージ), negbin(aodパッケージ), betabin(aodパッケージ)

「residual devianceを計算する基準となっているモデルの対数尤度を0として出したもの」
lmer(lme4パッケージ), glmmML(glmmMLパッケージ)

下にnegbin, betabin, glm.nbの検証プログラムを載せておく。

#negbin
n1<-100
a1<-numeric(n1)
rd1<-numeric(n1)
x1<-seq(6)
ID1<-factor(seq(6))
for (i1 in 1:n1){
y1<-floor(runif(6,min=1,max=20))
tdata1<-data.frame(x1,ID1,y1,y2)
res1<-negbin(y1~1+x1,random=~ID1,data=tdata1)
a1[i1]<- AIC(res1)@istats$AIC
rd1[i1]<-deviance(res1)