トップページ | 2008年3月 »

2008年2月

NCBIのjournal検索

学術誌の略名がよくわからない。

雑誌によっては、論文のReferenceを見ると雑誌名が略名になっていて、何のジャーナルかよくわからないときがある。

そんなときにNCBIのjournal検索が使える。
たとえばEcological Reserchの略名はEcol Resなのだが、これをNCBIで検索する。

まず、『Search』のプルダウンメニューからJournalsを選ぶ。
そして、『Ecol Res』と記入し検索。

すると候補の一番目にEcological Researchが出てくる。

Title: Ecological research Links
ISSN:  0912-3814 (Print)
ISO Abbreviation:  Ecol. Res.
Publication Start Year:  1986
Publisher:  Ecological Society of Japan,
Continuation Notes:  Continues in part: Nihon Seitai Gakkai shi. 
Language:  English
Country:  Japan
NLM ID:  101143839

こんな感じで、論文のいろいろな情報がわかる。

NCBIには他にもいろいろな情報を得ることができる。

拙者が使っているのはPubmedとJournal検索とDNAの実験をちょっとやっていた時にgenbankを使っていたくらいです。

| | コメント (0) | トラックバック (0)

pubmedのRSS

論文検索は普段、web of knowlegeを使っている。

しかし、生物系ならば、NCBIpubmedもいい。
生物系のものに限られてしまうが、無料のデータベースなので、家で調べたくなった時も使えるし。

PubmedはRSSに対応している。自分の興味ある分野のキーワードを登録しておけば最新の研究が上がってきたらすぐわかる仕組みになっている。

使い方がここでわかりやすく解説されている。

RSSリーダはいろいろあるようだけど、正直拙者は、あまり使いこなしていない。
Sleipnirを使っているので、それに標準で付いているRSSリーダを使って、
「あー便利だなー。」と思う程度。

| | コメント (0) | トラックバック (0)

東京のサンショウウオ

東京には3種のサンショウウオが生息している。

トウキョウサンショウウオ、ヒダサンショウウオ、そしてハコネサンショウウオ。

サンショウウオは大別して、田んぼなどに生息する止水性と渓流に生息する流水性がいて、トウキョウが止水性、ヒダとハコネが流水性になる。

拙者がその3種の中で拙者が一番すきなのが、ヒダサンショウウオ。
3種の中では一番稀少(無名)で、そして何より美しい。

Hida1  

水の中のヒダサンショウウオ 

 

Hida2  

顔アップ

 

この時期に産卵をする。
卵は石や岩の下に産卵し、卵嚢の片方を石にくっつける。
Hidaran1  

 

 

この卵嚢は水の中にある時は、青白い色を呈してとても美しい。
卵嚢の皮の構造色だそうだ。
Hidaran2  

 

 

こんな動物が東京の山の中でひっそりと暮らしている。

移動能力は非常に低いと思われるので、開発で山を削られればすぐに姿を消してしまうだろう。

| | コメント (0) | トラックバック (0)

サンショウウオの季節

東京では、サンショウウオの繁殖が始まっているころだろう。

もう産卵しているはずだ。

Tokyoran  

 

 形はバナナ型、これは3月撮影のため発生が進んでいる。

 

親はこんなかんじ
Tokyo

 

 

 

自分の写真ライブラリーを探してもほとんど写真がなかった。
どこでどう使うかわからないから、ちゃんと真面目に写真撮らなきゃ駄目だなとかんじる。
上の写真の撮影は2003年3月。

首都大の草野先生が積極的にトウキョウサンショウウオの研究を進めており、アマチュアの方たちと連携をとり、トウキョウサンショウウオ研究会を立ち上げています。

この研究会が主催するトウキョウサンショウウオ・シンポが今年は2月23日に行われたようです
最近は、トウキョウサンショウウオの話題だけでなく、全国のサンショウオやほかの両生類のネタを持っているプロ・アマ問わず招待して講演をしてもらい、日本産両生類の情報交換の場として、毎年活発な議論がされてきました。
今年で10年目なんですね。

| | コメント (0) | トラックバック (0)

モデル選択のときの交互作用(2)

モデル選択では、施肥・伐採の交互作用のみ効いているんだから、
biomass~ cut:fet とするのが妥当だと思う(思ってしまう)。
これを実際にやってみる。

> model2<- glm(biomass~cut:fet,data=dat)
> summary(model2)
Coefficients: (1 not defined because of singularities)
                  Estimate Std. Error t value Pr(>|t|)   
(Intercept)        42.1110     0.5216   80.73   <2e-16 ***
cutFALSE:fetFALSE -12.9255     0.7377  -17.52   <2e-16 ***
cutTRUE:fetFALSE  -12.2524     0.7377  -16.61   <2e-16 ***
cutFALSE:fetTRUE  -12.4844     0.7377  -16.92   <2e-16 ***
cutTRUE:fetTRUE         NA         NA      NA       NA   

Null deviance: 22619  on 399  degrees of freedom
Residual deviance: 10775  on 396  degrees of freedom
AIC: 2462.6

これは、デビアンス、自由度、AICすべてフルモデルと一致する。
そして、係数がちょっとヘン。

おそらく、カテゴリー変数の交互作用があっても単項がない場合は、その単項の自由度も消費しているのだと思う。

実際の世界では交互作用しか効いていないのに、モデル選択では、フルモデルが選択されることになる。
考え方を変えれば、施肥および伐採をしたことにより植物の反応があったわけで間違いではないんだけど。釈然としない。これがモデル選択の限界なのかもしれない。

粕谷先生の話だと、モデル選択と検定は別物なので、うまく使いわけないといけないんでしょうね。
そもそも、「この効果しか効いてないんだ!」という決定論的な考え方がモデル選択には合っていないのかもしれない。

もともと何でこんな話をしているかというと、lmer()の正規分布ではp値が出ないんだよね。
理由は、ここ参照。p値が出ないといろいろ困る人たちが世の中にたくさんいて、そんな中の一人が拙者だったり、一緒に議論している先輩だったりするわけだ。

今のところ拙者の理解では、「GLMMのような片方ベイズに足を突っ込んだようなものだと、もともと品質管理に使われていたネイマンピアソン主義的な検定は合わない」ということなんだと思います。

だれか、詳しいかたに解説してもらいたいです。

交互作用の解説は、粕谷先生のブログ詳しくされています
交互作用が積であらわされるのはわかりましたが、カテゴリー変数の場合はどうなるんだろう?

| | コメント (0) | トラックバック (0)

モデル選択のときの交互作用(1)

モデル選択のときの交互作用について、ドクターの先輩と議論している。

こんな実験をしたとする
森林において、伐採区、施肥区、施肥伐採区、コントロール区を設けたとして、植物の反応をバイオマスで見たとする。
施肥区、伐採区ではコントロール区とバイオマスは変わらないが、施肥伐採区でのみ植物が有意によく育ったとする。
各トリートメント100サンプルずつ採集した結果、

伐採区、施肥区、コントロール区では平均30、分散5
施肥伐採区で平均42、分散5 (当社比4割増)

のデータが取れたとする。これを人工データで作る。

biomass<- c(rnorm(300,30,5),rnorm(100,42,5))
cut<- c(rep(FALSE,100),rep(TRUE,100),rep(FALSE,100),rep(TRUE,100))
fet<- c(rep(FALSE,100),rep(FALSE,100),rep(TRUE,100),rep(TRUE,100))
treat<- c(rep("Control",100),rep("Cut",100),rep("Fertilization",100),rep("Interaction",100))
dat<- data.frame(biomass,cut,fet,treat)

これをANOVA&多重比較でやると、

>anova(lm(biomass~treat,data=dat))
Analysis of Variance Table

Response: biomass
           Df  Sum Sq Mean Sq F value    Pr(>F)   
treat       3 11843.8  3947.9  145.09 < 2.2e-16 ***
Residuals 396 10775.4    27.2                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

> TukeyHSD(aov(biomass~treat,data=dat))
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = biomass ~ treat, data = dat)

$treat
                                diff             lwr          upr        p adj
Cut-Control                0.6730700 -1.230197  2.576337 0.7982632
Fertilization-Control      0.4410732 -1.462194  2.344341 0.9326686
Interaction-Control       12.9254770 11.022210 14.828744 0.0000000
Fertilization-Cut         -0.2319968 -2.135264  1.671271 0.9892166
Interaction-Cut           12.2524070 10.349140 14.155674 0.0000000
Interaction-Fertilization 12.4844039 10.581136 14.387671 0.0000000

実際は、ブロックデザインとか、正規分布じゃないデータとか扱ったりするので、GLMやGLMMで解析し、さらにモデル選択となるとこんなにきれいにはいかない。

このデータそのままつかってGLMをする。出力結果は一部省略。

>model1<- glm(biomass~cut*fet,data=dat)
> summary(model1)

Call:
glm(formula = biomass ~ cut * fet, data = dat)

Coefficients:
                Estimate Std. Error t value Pr(>|t|)   
(Intercept)      29.1855     0.5216  55.950   <2e-16 ***
cutTRUE           0.6731     0.7377   0.912    0.362   
fetTRUE           0.4411     0.7377   0.598    0.550   
cutTRUE:fetTRUE  11.8113     1.0433  11.321   <2e-16 ***

Residual deviance: 10775  on 396  degrees of freedom
AIC: 2462.6

検定では、まあ、わかりきった結果です。伐採と施肥それぞれでは効果はないけど、交互作用のみ植物バイオマスに影響があり。となります。

これをモデル選択の立場で考えた場合ちょっと問題が生じてきます。

続きは次回。書いててだんだん言いたいことが不明瞭になってきた。

| | コメント (0) | トラックバック (0)

PRIMER-E

PDの先輩とちょっとだけANOSIMに関連する群集解析の話になった。PRIMER-Eというソフトが使えるらしい。Rではできない群集の視覚化の機能がすぐれているみたい。

解析の目的によっては、このソフトも使わなければいけないのかもしれない。

本家サイト
http://www.primer-e.com/

日本での販売サイト
http://www.software4research.com/acb/showdetl_Juca.cfm?DID=10&User_ID=259880&st=7947&st2=49646310&st3=85982634&Product_ID=629&CATID=17

| | コメント (0) | トラックバック (0)

勇払鵡川

TOEFに行ったついでに勇払鵡川に行ってきた。

勇払では、ホオジロガモを筆頭とした海ガモたちが相変わらず浮いていた。今年は、コオリガモが来ていないな。それが残念だ。アカエリカイツブリも去年は結構いたのに少ない。

林の方は、マヒワがハンノキのところに2羽いた。今年はマヒワが非常に少ない。そしてちょっと珍しい鳥を見た。
Toratu   

トラツグミ
留鳥だけどなかなか姿を現してくれない鳥くん。  

久々の鳥見なのでこれだけでも十分なのだが、鵡川ではさらにいい鳥が見られた。

Tumenaga

ツメナガホオジロ

鵡川についたあとすぐに登場し、飛んで行った。

  

  

そして陸のカモメと呼ばれるハイイロチュウヒ♂も出てきた。目の前に飛び出して来たが、カメラを取り出してあたふたしていると遠くの方へ飛んで行ってしまった。

  

Washi 漁港の方では、オオセグロやシロカモメに交じってワシカモメもいた。
  

  

  

さらに写真撮れなくて残念なのが、一瞬だけシラガホオジロを見た。
去年もここで2羽出現して、情報を聞いた次の日、直行したものの、それらしい小鳥2羽のシルエットと聞きなれぬエンベリ系地鳴きだけ確認できた。あの時と似た地鳴き、そしてあの白い模様。間違えなかった。鳥仲間に情報流しておいたので、誰かが写真を撮ってきてくれることを願う。

記録
2008年2月21日 勇払~鵡川

ミミカイツブリ、アカエリカイツブリ、ウミウ、ヒドリガモ、ホシハジロ、スズガモ、クロガモ、ホオジロガモ、シノリガモ、オジロワシ、トビ、ノスリ、ハイイロチュウヒ、オオセグロカモメ、シロカモメ、ワシカモメ、ヒヨドリ、ツグミ、トラツグミ、シラガホオジロ、ツメナガホオジロ、スズメ、カササギ、ハシボソガラス、ハシブトガラス
25種

シラガホオジロは初見。もうちょっと見たかった。

鵡川橋でカササギを見た。鵡川では初めて見た。苫小牧の個体群が勢力を伸ばしてきたか?

| | コメント (0) | トラックバック (0)

vegan研究(3)adonis()

example(adonis)をやってみる。

data(dune)
data(dune.env)
adonis(dune ~ Management*A1, data=dune.env, permutations=100)

Call:
adonis(formula = dune ~ Management * A1, data = dune.env, permutations = 100)

                          Df   SumsOfSqs  MeanSqs  F.Model     R2 Pr(>F)   
Management      3.00000   1.46859    0.48953    3.26288 0.3416  <0.01 ***
A1                   1.00000   0.44089    0.44089    2.93867 0.1026   0.62   
Management:A1  3.00000   0.58918    0.19639    1.30902 0.1370   0.97   
Residuals         12.00000   1.80036    0.15003               0.4188          
Total               19.00000   4.29902                            1.0000          
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

記法、出力結果ともまるでANOVAのようなので非常にわかりやすい。
ここでは、群集間でManagement、土壌A1層、そして二つの交互作用を要因として検定し、Managementで有意に群集間の違いが説明できることがわかる。

拙者の同期のデータでこのadonis()を試したとき、うまく計算できないデータがあったのだが、引数methodを変えることでうまく動いた。methodは13種類あり、defaultは"bray"。詳しくは?vegdist参照。おそらく、類似度指数の求め方の違いによりうまく計算ができたりできなかったりするのだと思う。ここまで書いて気づいたのが、ANOSIMのときもvegdist()は使うので同じ問題があるのだろう。

自分の今の理解ではどのmethodがどのようなデータのとき向いているかはわかりません。

| | コメント (0) | トラックバック (0)

久々のTOEF

久々にTOEFに来たので、30分ほど林内を散策。

これと言って珍しい鳥はいなかったけれど、久しぶりに見たシマエナガには癒された。

冬になると、鳥たちが人に近づいてくるようになる。

池の近くにいたマガモたち

0220kamo   

  

  

腕を下から上に振り上げるとわれさきにとやってくる。0220kamo2

  

  

  

そして何もせずにカラ類はやってくる

0220hasibutogara

ハシブトガラ

  

  

  

最近パソコンばかり見ていたから、ちょうどいい目の保養だった。

記録
マガモ、コゲラ、ミゾサザイ、シメ、シジュウカラ、ゴジュウカラ、ヤマガラ、ヒガラ、ハシブトガラ、エナガ、ヒヨドリ、カケス

エゾライチョウ見たいな。

| | コメント (0) | トラックバック (0)

vegan研究(2)ANOSIM実践編

example(anosim)とやってみる。

まずは、今回使われているデータのduneとdune.envが何なのか?

str(dune)
'data.frame':   20 obs. of  30 variables:
$ Belper: num  3 0 2 0 0 0 0 2 0 0 ...
$ Empnig: num  0 0 0 0 0 0 0 0 0 0 ...
$ Junbuf: num  0 3 0 0 0 0 0 0 0 0 ...
$ Junart: num  0 0 0 3 0 0 4 0 0 3 ...
以下略

str(dune.env)
'data.frame':   20 obs. of  5 variables:
$ A1        : num  3.5 6 4.2 5.7 4.3 2.8 4.2 6.3 4 11.5 ...
$ Moisture  : Ord.factor w/ 4 levels "1"<"2"<"4"<"5": 1 4 2 4 1 1 4 1 2 4 ...
$ Management: Factor w/ 4 levels "BF","HF","NM",..: 1 4 4 4 2 4 2 2 3 3 ...
$ Use       : Ord.factor w/ 3 levels "Hayfield"<"Haypastu"<..: 2 2 2 3 2 2 3 ...
$ Manure    : Ord.factor w/ 5 levels "0"<"1"<"2"<"3"<..: 3 4 5 4 3 5 4 3 1 1 ...

duneは20地点での30種の植生データで、dune.envはその地点の環境データのようですね。A1は土壌のA1層だと思います。そして、duneはdune.envは行で対応しているようで、1行が一つの地点の種および環境データになっています。

(dune.dist <- vegdist(dune))
           2        13         4        16
13 0.6000000                                                                                          
4  0.3563218 0.5128205                           
16 0.8933333 0.6060606 0.6666667
6  0.5111111 0.7530864 0.6344086 0.8518519
以下略

と、次に総当たりで類似度を求めて、

attach(dune.env)
dune.ano <- anosim(dune.dist, Management)
summary(dune.ano)

最後にanosim()の結果を出力しています。
第一引数が類似度行列、第二引数が環境条件。この場合は、management(4因子)によって違いがあるか?を調べたい。
出力結果は、

Call:
anosim(dis = dune.dist, grouping = Management)
Dissimilarity: bray

ANOSIM statistic R: 0.2579
      Significance: 0.006

Based on  1000  permutations

Empirical upper confidence limits of R:
  90%   95% 97.5%   99%
0.122 0.165 0.192 0.235

Dissimilarity ranks between and within classes:
        0%   25%    50%     75%  100%   N
Between  4 58.50 104.00 145.500 188.0 147
BF       5 15.25  25.50  41.250  57.0   3
HF       1  7.25  46.25  68.125  89.5  10
NM       6 64.75 124.50 156.250 181.0  15
SF       3 32.75  53.50  99.250 184.0  15

"Significance"がP値で、この例だと0.006なので、有意ということなんでしょう。
で、0%から100%まであるのが信頼区間なのでしょう。

さらにplot(dune.ano)で信頼区間の箱ひげ図が確認できて、betweenからずいぶん離れていると異常なんだそうです。notch=Tになっているので、箱のへこみが重ならないなら、その群集同士は違うものととらえていいんでしょう。

| | コメント (0) | トラックバック (0)

vegan研究(1)群集間をANOVAしたい

地点1、地点2、地点3で群集を調べて、この地点間で群集が違うと統計的に言うにはどうしたらいいか?これって結構難しい。
地点1ではA種がx1個体、B種がy1個体、C種がz1個体・・・
地点2ではA種がx2個体、B種がy2個体、C種がz2個体・・・
地点2ではA種がx3個体、B種がy3個体、C種がz3個体・・・
となると、種ごとに検定かけると、A種では地点1が一番多いけど、B種では違って・・・みたいにごちゃごちゃになってしまうし、いる種といない種もいる。

全種をまとめてみる方法として、多様度指数類似度指数があるけど、どのくらいの値の差で違いがあると言っていいのかもよくわからない。

クラスター解析だと、全体の傾向が視覚的にわかるけど、指標となる値がない。

一か月くらい前、拙者の同期がこの問題にぶち当たったので、一緒にあーだこーだ方法を話あった。おそらく、この中で一番いい方法なのは、クラスター解析だと思うのだが、情報が多くなればなるほど、視覚的にわかりにくい結果となるし、なんといってもP値がでない。ちなみにマフィアボス(当研究グループの教授)はクラスター解析を推奨していた。

このP値が曲者で、P値にこだわることはいけないことはわかっちゃいるが、なかなか抜け出せない。特に先生やPD以上の先輩方がそういうのにこだわって・・・、このP値というのは、一種の麻薬のようなもののような気がする。わかりやすいもんね。5%水準で白黒つけられるから。

久保さんたちのような統計わかっている人たちに言わせれば、「ANOVAのようなネイマン・ピアソン決戦主義の・・・」と言われるんでしょうが・・・ ちょっとこの話はパス。

そんなこんなで群集もANOVAのような手法でお手軽にP値を出せたらいいな思う人が世界中にいるようで、いろいろ検討され、手法が開発されているようです。

まず、一番単純な方法として、χ二乗独立性検定
これは、割合データなので、久保さんたちのような人たちに言わせると・・・以下略。
とにかくあまり良くないと思います。

もうちょっと、データを殺さないで検定する方法として開発されたのが、ノンパラなANOSIM (Analysis of Simiralities)という方法。内部でどのような計算をやっているかはよくわかりませんが、この方法によりP値が求められ、P<0.05ならば、「比べている群集間に違いがある。」と言えるようです。

ネット(日本語限定)で、検索するとANOSIMを使った研究がいくつか出てきます。

この解析手法はRのveganパッケージに含まれています。関数はanosim()です。veganのtutrialを見ると、ANOSIMのほかにも「群集のちがいの検定」があるようです。

ひとつが、mrpp()で計算することができるMultiresponse permutation procedureという方法があるようです。

で、このanosim()とmrpp()はtutrialを読む限りロバストではない方法のようであまり推奨されておらず、その代りadonis()(Multivariate Analysis of Variance Using Distance Matrices)がロバストな方法として推奨されています。しかしながら、このtutrial(バージョン February 13, 2008)には、anonis()の解説はしておらず、「そのうち、解説するよ」としか書いていない。このadonis()が何者かを知るには原著を当たるしかないのでしょう。(Anderson, M.J. 2001. A new method for non-parametric multivariate analysis of variance. Austral Ecology, 26: 32–46.)

| | コメント (0) | トラックバック (0)

lmer研究(3)gaussian modelのp値

lmerでは正規分布でのモデルではいわゆるwald検定によるp値が表示されない。
言い換えるとsummary(model)とやった時にcoefficientsの右の方にあるp値がでない。

以前のlme4のバージョンでは表示されていたのだが、新しいバージョンでは表示されなくなった。これに関して、製作者であるDouglas Bates氏の公式見解をネットを徘徊しているときに見つけた。

https://stat.ethz.ch/pipermail/r-help/2006-May/094765.html

拙者は、英語が非常に苦手なうえに統計自体もそれほど理解している方でないので、正確な意味はわからないけれど、解読してみる。

まず理由は、

Because they feel that the denominator degrees of freedom and the
corresponding p-values can easily be calculated they conclude that
failure to do this is a sign of inattention or, worse, incompetence on
the part of the person who wrote lmer (i.e. me).

だそうな。
つまり、「出力されたアヤシイP値を多くの人が安易に信頼しているから」ということかな?耳の痛い話です。

その下の節には、さらに詳しく、なぜP値を出していないかが、統計的な説明を交えてされています。要約すると、

「階層化や並列化されたランダム効果を仮定した場合、denominaterの自由度は不正確なものになるため、P値も信用できない」

ということかな。Bates氏はSASを意識していて、行間を読むと、SASも決して正確な値を返しているわけではないと言いたいのかもしれない。Bates氏も周囲から「なんでSASでP値がでるのにlmerではでないの?」なんて聞かれているのかもしれない。

んで、Bates氏の推奨する方法は、MCMC法を使っているmcmcsamp()とcodaパッケージ。要はベイズってことですね。

話はちょっと変って、
lmerではモデルの自由度も表示されない。
モデルの自由度は(観測数)-(パラメータ数)なので、

model<- lmer(・・・)
nrow(model@frame)-attr(logLik(model),"df")

で求められる。
ちょっとめんどいけどね。
このモデルの自由度もBates氏は意図的に出していないのかもしれない。

| | コメント (0) | トラックバック (0)

aod研究(1)aodパッケージのakaike weight

aodパッケージにはdelta AICとakaike weightを求めてくれる便利な関数があります。
これはパッケージ内に含まれるGLMM関数であるnegbin()とbetabin()で記述されたモデル専用ですが、便利なので覚えておくと便利です。aodはS4なので、lmer()でも応用できるかと思ったらダメでした。残念。

aodではAIC(model)とやるとAICとAICcの両方を返してくれます。しかも、これはlmer()で記述されたモデルでもAICを返してくれます。ただし、lmer()でAICを求める場合は、AIC(logLik(model))とする必要あり。これを利用して、自動akaike weight プログラムではAICcを求めています。

さて求め方は下の例を参照のこと

> library(aod)
> model1<-  negbin(y ~ log(dose + 10) + dose, ~ 1, salmonella)
> model2<-  negbin(y ~ log(dose + 10) , ~ 1, salmonella)
> model3<-  negbin(y ~ dose, ~ 1, salmonella)
#aodに含まれるデータsalmonella使用

> AIC(model1,model2,model3)
           df AIC     AICc
model1  4 133.7792 136.8561
model2  3 137.3354 139.0497
model3  3 141.6565 143.3708

> summary(AIC(model1,model2,model3))
          df     AIC    Delta          W          Cum.W
model1  4 133.7792 0.000000 0.84144392 0.8414439
model2  3 137.3354 3.556198 0.14216958 0.9836135
model3  3 141.6565 7.877324 0.01638649 1.0000000

> summary(AIC(model1,model2,model3),which="AICc")
          df  AICc      Delta          W         Cum.W
model1  4 136.8561 0.000000 0.72863188 0.7286319
model2  3 139.0497 2.193560 0.24332270 0.9719546
model3  3 143.3708 6.514686 0.02804542 1.0000000

まあ、こんな感じで、summary()って便利ですねってことです。

そもそもnegbin()とbetabin()ってなにに使うの?って方は、こちらあたりを参照し、ランダム効果の入れ方の注意はこちらあたりを参照のこと。

| | コメント (0) | トラックバック (0)

lmer研究(2)自動akaike weight プログラム2

<2009.2.19追記>
このコードは使えなくなりました。
新しいコードはこちら
________________________________________________________

lmer()のakaike weight自動計算プログラムのランダム効果複数指定可能版。

使い方は、

model<- lmer(Z~ A + B + C + (1|D/E) + (0 + C|F),・・・)
model.average.lmer2(model,random="(1|D/E) + (0 + C|F)")

とすれば、akaike weightを求めてくれます。
random= にランダム要因部分を記入すれば計算してくれます。

lmer()はS4のため、step()が使えないので、この関数はずいぶん便利だと思います。
ちなみに、step()もlmer用に移植しようと試みていますが、いまのところまだできてません。step()およびstepAIC()のプログラムを眺めていてわかったのが、この関数はもともとlmer()の原型であるlme() (in nlme package)のオブジェクトもモデル選択できることを想定しているようだけど、実際試すとできない。今後対応するのだろうか?いまだにVRパッケージは更新しているようなので、そのうちやってくれるだろうか?いずれにせよ、lme()に対応するならlmer()に対応してくれた方が数倍便利だと思うので、maintainerのRipleyさんにがんばってほしいです。

negbin()のakaike weight自動化プログラムも同じ要領で作りたいんだけど、致命的なことにnegbin()は、update()が使えないためにできない。飯島さんのサイトglmとglm.nbのakaike weight プログラム(下の方)を参考にして作ってみるかな。

一人でイチからプログラムを作る力はまだ自分にはないです。

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

##Depends lme4, aod
#複数のランダム効果指定用
model.average.lmer2 <- function( object, criterion = "AIC",random, ... ) {
    # extract model variables
       terms <- object@terms;
       label <- attr(terms, "term.labels");
       kmax <- length( label );

    # define AIC function
  myAIC <- function(model, ...) {
    aic   <- AIC(logLik(model))
    if( criterion == "AIC" ) {
      aic[1];
    } else if (criterion == "AICc") {
      aic[2];
    }
  }

  # null model
  object.null <- update(object, paste(c(".~.-.",random),collapse="+") )
  model <- "(Intercept)";
  aic <- myAIC(object.null);
  coefficients <- c(summary(object.null)@coefs[1], rep(NA, kmax));
  coef.sd <- c(summary(object.null)@coefs[2], rep(NA, kmax));
  names(coef.sd) <- names(coefficients) <- c("(Intercept)", label);

  # updating models
  for ( n in 1:(2^kmax - 1) ) {
    m <- numeric();
    k <- n;
    for ( i in kmax:1 ) {
      m[i] <-  k %/% 2^(i - 1);
      k <- k %% 2^(i - 1);
    }
    m <- ifelse( m == 1, TRUE, FALSE);
    result <- update(object.null, paste(c(" .~.", label[m]), collapse = "+"));
    model[n+1] <- paste(label[m], collapse = "+");
    aic[n+1] <- myAIC(result);
    sd <- coef <- rep(NA, kmax + 1);
    coef[c(TRUE,m)] <- summary(result)@coefs[,1];
    sd[c(TRUE,m)] <- summary(result)@coefs[,2];
    coefficients <- rbind(coefficients, coef);
    coef.sd <- rbind(coef.sd, sd);
  }

  # calculating Akaike's weight w
  if( criterion == "AIC" ) {
    out <- data.frame(model, AIC=aic);
    out$dAIC <- aic - min(aic);
    out$w <- exp( -0.5 * out$dAIC ) / sum(exp( -0.5 * out$dAIC ));
    }
  else if (criterion == "AICc") {
    out <- data.frame(model, AICc=aic);
    out$dAICc <- aic - min(aic);
    out$w <- exp( -0.5 * out$dAICc ) / sum(exp( -0.5 * out$dAIC ));
    }

  # formatting coefficents and sds
  rownames(coef.sd) <- rownames(coefficients) <- 1:nrow(coefficients);

  # calculating z-statistics and p-value
  z <- coefficients / coef.sd;
  p <- pnorm(z);

  # averaging models
  ave.coef <- apply( coefficients * out$w, 2, sum, na.rm = TRUE );
  ave.var <- apply( ( coef.sd^2 + (t(t(coefficients) - ave.coef))^2 ) * out$w, 2, sum, na.rm = TRUE );
  ave.sd <- sqrt(ave.var);
  ave.z <- ave.coef / ave.sd;
  ave.p <- pnorm(ave.z);

  # importance of variables
  tmp <- ifelse(!is.na(coefficients[,-1]), matrix(out$w, nrow=nrow(out), ncol=ncol(out) - 1), 0);
  iv <- colSums(tmp);
  names(iv) <- label;

  # output
  if( criterion == "AIC" ) {
    o <- order(out$AIC);
    list(summary = out[o,], coefficients = coefficients[o,], sd = coef.sd[o,], p = p[o,],
       ave.coef = ave.coef, ave.sd = ave.sd, ave.p = ave.p, iv=iv);
    }
  else if (criterion == "AICc") {
    o <- order(out$AICc);
    list(summary = out[o,], coefficients = coefficients[o,], sd = coef.sd[o,], p = p[o,],
       ave.coef = ave.coef, ave.sd = ave.sd, ave.p = ave.p, iv=iv);
    }
}

| | コメント (2) | トラックバック (0)

lmer研究(1)自動akaike weightプログラム

<2009.2.19 追記>
このプログラムは使えなくなりました。
書き直したコードはこちら
_________________________________________________________________

東北大の富田さんが、lm()&glm()glmmML()用のakaike weight、model averaging用関数を作ってくれている。これがとっても便利なので、自分が良く使うlmer()やnegbin()でも使えないかと思い、富田さんのコードをlmer用に書きなおした。

現状では、交互作用項は非対応、そして、ランダム要因は一個のランダム切片のみ。
複数のランダム要因を指定したモデルの場合はこのプログラムは使えない。
複数ランダム要因が指定できるプログラムは別個に作っているので、後日アップします。

このプログラムは富田さんの同様、AICとAICcで計算してくれ、デフォルトはAICで、criterion="AICc"とすれば、AICcで計算してくれます。
ただし、AICcで計算する場合、別途aodパッケージを読み込む必要があります。
そしてoffsetにも対応しています。

使い方は、
model<- lmer(Z~ A + B + C + (1|D)・・・)
といつものようにモデルを指定してあげたら
model.average.lmer(model)
とコマンドするだけ。
出力には、富田さんの日記にあるテーブル整形関数も使うと便利だとおもいます。

あと、このプログラムは前述したとおり、顔も知らない富田さんのコードを勝手にぱくったものなので、ご本人様から文句がついたら削除します。

ダウンロードは「model_average_lmer.r」をダウンロード

##Depends lme4, aod
#一つのランダム切片指定用
model.average.lmer <- function( object, criterion = "AIC", ... ) {
    # extract model variables
       terms <- object@terms;
       label <- attr(terms, "term.labels");
       kmax <- length( label );

    # define AIC function
  myAIC <- function(model, ...) {
    aic   <- AIC(logLik(model))
    if( criterion == "AIC" ) {
      aic[1];
    } else if (criterion == "AICc") {
      aic[2];
    }
  }
 
  # null model
  object.null <- update(object, paste(c(".~.-.",
     paste("(1|",names(summary(object)@flist),")",sep="")),collapse="+") )
  model <- "(Intercept)";
  aic <- myAIC(object.null);
  coefficients <- c(summary(object.null)@coefs[1], rep(NA, kmax));
  coef.sd <- c(summary(object.null)@coefs[2], rep(NA, kmax));
  names(coef.sd) <- names(coefficients) <- c("(Intercept)", label);

  # updating models
  for ( n in 1:(2^kmax - 1) ) {
    m <- numeric();
    k <- n;
    for ( i in kmax:1 ) {
      m[i] <-  k %/% 2^(i - 1);
      k <- k %% 2^(i - 1);
    }
    m <- ifelse( m == 1, TRUE, FALSE);
    result <- update(object.null, paste(c(" .~.", label[m]), collapse = "+"));
    model[n+1] <- paste(label[m], collapse = "+");
    aic[n+1] <- myAIC(result);
    sd <- coef <- rep(NA, kmax + 1);
    coef[c(TRUE,m)] <- summary(result)@coefs[,1];
    sd[c(TRUE,m)] <- summary(result)@coefs[,2];
    coefficients <- rbind(coefficients, coef);
    coef.sd <- rbind(coef.sd, sd);
  }

  # calculating Akaike's weight w
  if( criterion == "AIC" ) {
    out <- data.frame(model, AIC=aic);
    out$dAIC <- aic - min(aic);
    out$w <- exp( -0.5 * out$dAIC ) / sum(exp( -0.5 * out$dAIC ));
    }
  else if (criterion == "AICc") {
    out <- data.frame(model, AICc=aic);
    out$dAICc <- aic - min(aic);
    out$w <- exp( -0.5 * out$dAICc ) / sum(exp( -0.5 * out$dAIC ));
    }

  # formatting coefficents and sds
  rownames(coef.sd) <- rownames(coefficients) <- 1:nrow(coefficients);

  # calculating z-statistics and p-value
  z <- coefficients / coef.sd;
  p <- pnorm(z);

  # averaging models
  ave.coef <- apply( coefficients * out$w, 2, sum, na.rm = TRUE );
  ave.var <- apply( ( coef.sd^2 + (t(t(coefficients) - ave.coef))^2 ) * out$w, 2, sum, na.rm = TRUE );
  ave.sd <- sqrt(ave.var);
  ave.z <- ave.coef / ave.sd;
  ave.p <- pnorm(ave.z);

  # importance of variables
  tmp <- ifelse(!is.na(coefficients[,-1]), matrix(out$w, nrow=nrow(out), ncol=ncol(out) - 1), 0);
  iv <- colSums(tmp);
  names(iv) <- label;

  # output
  if( criterion == "AIC" ) {
    o <- order(out$AIC);
    list(summary = out[o,], coefficients = coefficients[o,], sd = coef.sd[o,], p = p[o,],
       ave.coef = ave.coef, ave.sd = ave.sd, ave.p = ave.p, iv=iv);
    }
  else if (criterion == "AICc") {
    o <- order(out$AICc);
    list(summary = out[o,], coefficients = coefficients[o,], sd = coef.sd[o,], p = p[o,],
       ave.coef = ave.coef, ave.sd = ave.sd, ave.p = ave.p, iv=iv);
    }
}

| | コメント (0) | トラックバック (0)

ブログはじめました

このブログは生態学を研究している院生のブログです。
とある殿さまの下で研究をしているので「足軽」と呼ばれたりするので、この名前にしました。
いい名前が思い浮かんだらまた変えるかも。

日記というより、自分の研究ノートの代わりに書いていこうと思います。

三日坊主にならないようにします。

| | コメント (0) | トラックバック (0)

トップページ | 2008年3月 »