リアルX-Men

スラッシュドットの記事で発見。

指先にツメを隠しているカエルが発見されたということです。
危険を感じると皮膚を破ってツメが飛び出てくるそうな。

新種発見ということではなく、もともと知られていた種にこのような爪があることがわかったそうです。

サイエンスの記事の文頭には、

X-Men fans rejoice: Wolverine has come to life, as a frog.
（X-Menのファン大喜び。ウルヴァリンがカエルになって現世にやってきた。）

なんて書いてある。
アメリカ人はX-Menを想像するんですね。

元論文には、カエルのツメの写真が出ています。

このツメはある程度の攻撃力があるようで、引っかかれて血がにじむくらい鋭いそうです。

カメルーンにいる１２種のArthroleptidaeがこの形質をもっている。
種類は多いけど、一種しか確認されていない肺のないカエルよりもこっちの方が驚きです。

土をまさぐってみた

ちいさなミミズを見つけた。

　

フトミミズ科の一種。

TOEFでは、フトミミズ科は３種類知られており、いずれも一年生の卵越冬。この個体は雪が解けてから生まれたばかりの赤ちゃんだ。このくらいだとまだ種類の同定ができない。

体重数mgといったところか。こいつが、９月をこえて性成熟するころには、１０ｇを超える巨体になる。

TOEFは、土壌層が薄い割には非常にミミズが多いことが知られている。そのため、夏のTOEFは大型ツグミランドになる。
昼間はアカハラ、クロツグミが、夜はトラツグミがミミズを食べまくっているのだろう。
特に明け方の林道を走るとミミズを大量にくわえたクロツグミが林道をバンバン飛び立つ姿が観察できる。

　

こちらは、ツリミミズ科のサクラミミズ。複数年生きる種。TOEFには、ツリミミズ科も３種いるが、このサクラミミズが圧倒的に多い。

このサクラミミズは、日本中に分布しているが、なぜだかTOEFの個体群はかなり大きいのだそうだ。自分は、ここでしかサクラミミズを見たことがないが、他の地域ではもっと小さく、写真のフトミミズがもう少し育った程度くらいにしか大きくならないのだそうだ。なぜそうなのかは、未だに解明されていないそうです。

　

定規を持っていなかったので、フィールドワーカーがよく使う「セ-Y3」を使って比較。
大きい方が、ツリミミズで小さい方がフトミミズ。
この大きさが４ヶ月後には逆転する。

林道

TOEFの樹木たちの展葉はこれからだ。

　

ここはなんでこんなに遅いのだろう？

札幌はとっくに生い茂っているし、TOEF近くの支笏湖周辺もずいぶん展葉している。

カナヘビ

林内を歩いていたらカナヘビを見つけた。

けっこうおとなしいので、写真を撮らせてもらった。

　

　

　

　

撮ったその場では、「うん。なかなかよく撮れた」と思っていたけど、光条件が悪かったのか、ちょっと寒色系が強い。
フラッシュを焚いた方が良かったかもしれないが、スーパーマクロだとフラッシュの光が全体に広がらず、変な影ができてしまう。

写真はやはり難しい。

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

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

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

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

即注文しました。

楽しみです。

はやく来ないかな～

悲しいかな新川河口

朝早めに起きて、新川河口へ出向く。

シギチはシロチドリしかいなかった。
悲しい。

シギチ以外ではノゴマがわりと近くで見ることができた。
河口へと続くダートのロープに止まって一生懸命さえずっています。
カメラをかまえると逃げて行ってしまった・・・

記録　２００８年５月２３日　新川河口
シロチドリ、オオセグロカモメ、ウミネコ、トビ、マガモ、ヒバリ、ノビタキ、ノゴマ

新川から北大へ。
朝の鳥見会に参加する。
だんだんシメの数が減ってきた。
コルリの声が聞こえた場所ではすでに声は聞けず。
すでに抜けてしまったかな？

予言

今年はシラカバ花粉が多かった。

４月２２日の日記の写真にあるように花が例年より多かった。

花が多いということは、できる実も多いことが予想される。

で、ふとネットを徘徊していたらこんな写真を見つけた。

来シーズンの冬はシラカバの枯れた実が多いはず（ハンノキの花粉飛散も多かったので、ハンノキの実も多いかも）なので、種子食冬鳥の食糧がいっぱいあるので、大量に冬鳥飛来するかも。

２年連続の冬鳥当たり年か？

花粉症そろそろ終息か？

花粉症だいぶ楽になってきた。

４月２２日の日記と同じ場所で写真を撮ってみた。

　

　

　

花が枯れ、実ができているようです。

ということで、（すくなくともここらへんでは）既に花粉の飛散は終了している。

しかし、まだ症状がでるということは、大気中のエアロゾルにまだ相当量花粉があるのでしょう。

今回の雨で、かなり楽になったので、もうひと雨降ればもう大丈夫になるかな。

嵐が運んだ？

朝早く起きて、北大構内で鳥見というかセンサスのお手伝い。

マヒワはだんだん抜けているようで、あまり見ない。
相変わらずシメが多い。
一か所に百羽くらいとかいるし。

いつまでいるかこのブログ上にでも記録しておきたいけど、来週はTOEFなんだよな。

歩いている途中、原生林でコルリのさえずりを聞く。
研究室のある建物前の林でカッコウの声を聞く。
上空をアマツバメが飛ぶ。（自分は確認しなかったがヒメアマツバメもいたらしい。）

いづれも今シーズン北大構内初認だ。

前日の低気圧が運んできたのかな？

記録 2008年5月21日　北大構内
マガモ, キジバト, カッコウ, アマツバメ, アカゲラ, ハクセキレイ, ヒヨドリ, モズ, コルリ, ルリビタキ, ウグイス, エゾムシクイ, センダイムシクイ, キビタキ, コサメビタキ, ハシブトガラ, シジュウカラ, アオジ, カワラヒワ, マヒワ, シメ, スズメ, コムクドリ, ムクドリ, カケス, ハシボソガラス, ハシブトガラス
27種

シギチが気になるなあ。なんか運ばれてきているのかな？

いつ間にやら発売

ふと、kowaのページを見てみると、なんと、新しい双眼鏡が発売していた。最近、birder誌とかを読んでいなかったので、全然知らなかった。

BD25 GR

コーワのC3コートをプリズムに施した２５mm口径の改良機なんですね。

製品ページや上の写真を見ると、ヒモ通しが両側についているようなので、カバーをつけてやる必要はなくなるようだ。

軸は、舶来やニコンのHGの同クラスのようによりコンパクトになる二軸ではなく、扱いやすい一軸を維持したのか。
実機をみたことないから何とも言えないけれど、もともとコンパクトなものだし、二軸に比べてコンパクトにならずともそれほど問題ない気がする。

定価が３万切っているので、同クラスの舶来物が５万円台後半から１０万で、ニコンのHGが実売５万を切るくらいを考えると、なかなかいい値段設定ではないかな。コーワの製品だから、見えはおそらくいいのであろうし。

うーむ。これは、ほしいかも。
ちょっとお金ができたらほんと買っちゃおうかな。

ちょっと北の方を歩いてみる

次世代ポストゲノムセンターと牛　

ノビタキがさえずっていた
　

ヒバリも飛んでいた（ピントあわず）
　

くさっぱら。牧草なのか？
　

もうちょっと散策したっかったけど、鼻が限界。
ティッシュを携帯していなかった。

怠惰な一日

前日３時まで飲んでいたけれど、朝７時の鳥見会に行った。

眠い。

いつものコースを散策。

いまだにシメが多い。集団で飛んでいるところを撮りたいけど、難しい。

いつものコースを回ったあとは、有志で実験苗畑へ。
ここもシメが多い。

ここも一周して鳥見は終わり。

さすがに眠いので、いったん家に帰ってから一眠りしてから研究室に行こうと思い、家に帰る。１、２時間寝るつもりだったんだけど

・・・・・

起きたのは１８時だった。

ガックシ。

夕飯の準備して、もう家に引き籠ることにした。

記録　北大構内　5月17日
マガモ, アカゲラ, ハクセキレイ, ヒヨドリ, ウグイス, エゾムシクイ, センダイムシクイ, オオルリ, エゾビタキ, コサメビタキ, ハシブトガラ, シジュウカラ, アオジ, カワラヒワ, マヒワ, イカル, シメ, スズメ, コムクドリ, ムクドリ, ハシボソガラス, ハシブトガラス
22種

とにかく、種が大豊作なのでアトリ類がいっぱい。

シメ大群の原因

飯島さんの日記の5月１２日の記事にハルニレに来ているシメの写真が載っている。

今年はハルニレの種が豊作なので、ここ数年ないくらい多数のシメが北大構内に居ついているということなんだろう。

今年は、シラカバ花粉も尋常じゃないくらい飛散した。花粉の飛散量は、前年の気温に依存しているようで、去年の猛暑が効いているんだろう。

ハルニレ豊作もシラカバと同じく去年の猛暑が効いているのかも。で、その豊作が鳥類群集（シメとかマヒワとか）の飛来にも波及していると。

まさにボトムアップ効果ですね。

今シーズンは、冬鳥が当たり年だったけど、夏の気温が高くて全体的に冬の鳥の餌資源が多かったからなのかもしれない。

今年の夏はどうなのか？そして冬は？

冬鳥大当たりは歓迎だが、猛暑はいやだ。

師弟三代共同作業

殿のかつての指導教官の統計下請をしている。

殿の師匠は大学退職後、奥さんと日本全国温泉巡りをするついでに、ご当地の動物をサンプリングしていた。というか、温泉巡りは名目で、サンプリングがメインだったんだろうけど。

採集記録は、国内誌に掲載され、さらにそのDNAを使った研究は自分の学部時代の指導教官により解析され、成果が国際誌に掲載された。

そんな殿の師匠は、喜寿を2年後に控え、まだまだ現役。
時間がある時にフィールドに出てサンプリングをしているらしい。
去年は、新しく開発したトラップを自慢げに見せてくれた。
また、そのとき直々にトラップの仕掛けかたを教えていただいた。

動物への探求心は未だに衰えることなく、博物館のボランティアとして作業する傍ら、自分が採集した標本をつかった研究を続けられているスーパーおじいちゃんである。

自分も将来こういう風になれるのか？ひとつの目標とする人物である。

そんな殿の師匠が、短報ながら新しく論文を執筆された。
先ほどの日本全国温泉採集行脚のときのサンプルを用いた形態の研究である。

結果は非常にシンプルで、地域ごとに歯の突起があるかないかを調べ（よくそこに着目したもんだ！）、地域群ごとに突起の有無に違いがあり、それがDNAを用いた研究で得られた傾向と合致するというものだ。

そこで、突起の有無が地域ごとに違うことを統計的に示したいけれど、そのやり方がわからないということで、自分に統計のお役が回ってきた。

殿の師匠は、あるなしの割合データを持ってきた。割合データというとかの人が大嫌いなデータなので、これを現代統計学っぽくGLMの二項モデルで解いた。

もともと結果が非常にクリアーなので、サクッとゆ～い差をだして、解析結果をまとめ作図もちょっとして殿の師匠に見せた。殿と３人で統計のことのみならず、論文の内容のことを小一時間ほどディスカッション。GLMははじめて聞いた単語のようで、キョトンとしていたけれど、まあ結果が出たからいいとすることに。どうやったか、メールで説明することにしたけど、統計の部分は自分で書いた方が早そうなので、書いて送ることにした。

GLMは、ジャーナルや当たるレフェリーによって不運に見舞われるという話をよく聞く。殿の師匠が出す予定のジャーナルで使っている人がいるのかちょっと気になったので、調べてみると東のイノシシでGLMによるAkaike weightをつかったモデル選択な研究をあっさり見つけたので、まあ大丈夫なんだと安心。

原稿は、殿がチェックしたくさんのコメントをつけていた。殿の師匠はこれからそのコメントの対応と統計の部分を原稿に組み込むことになるのだろう。
いつか、自分も殿の原稿をチェックしてコメントをびっしり書くようになるのだろうか？
それを目標の一つにして頑張ろう。

この論文、もともと5年前に結果が出ていて放置していたらしい。
ある日突然思いついて、一気に書き上げたそうな。

この論文と別にして、あと二つネタが残っていて、ひとつはサンプルは集まっていて、あとは計測すれば結果が出る。これは、多変量解析をする必要があるので、また自分の出番がありそう。

もうひとつは、島にサンプリングをしに行かなければいけない。
やらないか？と言われたけど、時間がとれるかなあ？
今やってるネタと全く関係ないからなあ・・・
興味はあるんだが。

Tree diversity analysis

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

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

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

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

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

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

大群

朝、少し研究室のある建物前で鳥を見る。

林がざわめいていた。

シメとマヒワの群れが新芽や実のようなものを食べていた。

特にシメが多くて100羽単位でいたと思う。

構内全体がこんな感じだから、全部でどれくらい北大構内にいるのだろう？

記録　5月１４日 北大構内
キジバト、シメ、マヒワ、アオジ、センダイムシクイ、ハシブトガラ、シジュウカラ、ムクドリ、スズメ、ヒヨドリ、ハシブトガラス

想定外の種

先日の鳥記録をよーく見るとおかしな種が混じっている。

ふつーじゃありえない種をみてしまった。

　

　

　

近かった。
すぐ近くの電線にそいつはいた。

ハヤブサよりも明らかに小さい。

チゴハヤブサだと思ったけれど、赤い部分がない。

コチョウゲンを疑ったけれども、違和感があった。

クチバシの色を見て図鑑を見ると一発だった。

パ タパタと飛び出すと、地面におりてまた電線に上がった。

虫を捕まえて食べていた。

　

唯一撮れた飛翔写真。

デジ眼部隊は、すばらしい飛翔写真を撮っていた。
自分のカメラの性能と腕のなさが悔やまれる。

　

一時間ほどじっくり観察していたらどこかに飛んで行った。

2年前の天売で見た以来2回目。

次に見られるのはいつだろう。

朝の新川河口

昨日の飲み会で、新川にシマアジが来ていることを知る。

朝早く起きて行ってみた。

　

すでに抜けていた・・・・

悲しいので、恐る恐る砂浜を走らせ、シギチを探す。

　

この子らだけ・・・

あとは、

　

ヒバリがいたり、ホオアカがいたり、ショウドウツバメがいました。

前日飲み会だったのに早起きして、一日眠かった。
仕事もあまり捗らず。

なんか反省・・・・

記録2008年5月13日
マガモ, キンクロハジロ, シロチドリ, オオセグロカモメ, ウミネコ, キジバト, ヒバリ, ショウドウツバメ, ハクセキレイ, ノビタキ, ホオアカ, アオジ
１２種

ちなみに、今日は北大構内にはマヒワとシメの大群があちこちにいた。
とにかく葉の茂り始めた木の中が騒がしかった。

植物園ではこれに合わせてレンジャクが騒いでいた。
少なくてもキレンジャクはいた。

ついに渡りのピークか？

夏鳥ウォッチング

朝5時集合。ねむい。

まずは、千歳川中流域。

センダイムシクイ、オオルリ、キビタキ、クロツグミと代表的な夏鳥たちを見ることができた。朝早く起きて来てくれた鳥初心者の人たちもきれいな鳥を結構近くで見ることができて、喜んでくれたのではないかな。カワセミも見られたし。

ちなみに今日のメインは、ここ。
９時前にここの鳥見を終了。
あとは、おまけみたいなもの。

続いてウトナイ

白鳥の餌付場へ
　

ケガをして飛べないヒシクイと、白鳥のエサのおこぼれにあずかるオオジュリン。

さっさと、TOEFへ移動。

アカハラが多い。
芝生に５羽くらい一緒にいる場所もあり。

ここでは、他の人とちょっぴり別行動。
自分の仕事をするために、見物人2人連れて林内へ。
フィールド内で３０分くらい作業。

勇払へ。

　

青春の一ページみたいな一枚。

勇払周辺は、なぜだかカササギがいることで有名。
町中を探索すると小さな公園にいた。
　

以前は室蘭にもカササギの個体群がいたそうだが、今はいないらしい。
なぜ、ここにいるのかは不明。勇払周辺にしかおらず、分布拡大の気配もあまりない。
あまり話は聞かないが、ここ周辺でも電柱に巣を作っていたりしているのだろうか？

そして、鵡川へ。

カモメくらいしか鳥がいない。
乾燥化してしまって、ほとんどシギチも見当たらない。
残念だ。

ちょっと、鳥とは別にハマトビムシをサンプリングする必要があったので、
砂浜で探す。
木片を引っぺがすとピョンピョンと出てきた。
適当に捕まえて、アルコールで保存。

あまりに鳥が出ないので、さっさと温泉に行って、帰ることにする。

2008年5月11日　千歳川中流域～ウトナイ湖～TOEF～勇払～鵡川
アオサギ, ヒシクイ, オオハクチョウ, オシドリ, マガモ, コガモ, ヨシガモ, ヒドリガモ, オナガガモ, ハシビロガモ, キンクロハジロ, オジロワシ, トビ, ハヤブサ, アカアシチョウゲンボウ, オオセグロカモメ, シロカモメ, キジバト, カワセミ, ヤマゲラ, アカゲラ, コゲラ, ヒバリ, ハクセキレイ, タヒバリ, ヒヨドリ, モズ, ノビタキ, アカハラ, ウグイス, エゾムシクイ, センダイムシクイ, キビタキ, オオルリ, コサメビタキ, エナガ, ハシブトガラ, ヒガラ, シジュウカラ, メジロ, アオジ, オオジュリン, カワラヒワ, イカル, シメ, ニュウナイスズメ, スズメ, ムクドリ, カケス, カササギ, ハシボソガラス, ハシブトガラス
52種

夜は、調査バイトのあいた時間で大量にギョウジャニンニクとウドを採った人の家にいって、みんなでギョウニン餃子パーティをした。
飲んだのだが、あまり寝ていなくて、ドライバーだったので、疲れていたのであろう。
少しの間、コテっと寝てしまった。
すぐ起きて、顔に落書きを回避。危なかった。

レーダーで渡り鳥調査

うちの研究室がある、建物は、生物の研究だけでなく、大気、気象、惑星などの分野の研究室もある（というかそっちがメイン）。

そのため、建物の屋上には気象を計測するドップラーレーダーとかライダーとかが設置されている。

そのレーダーに鳥の影がノイズとして映るらしい。

それに目をつけたバードリサーチの方々が、この機器を使って渡り鳥の調査を始めた。

この気象レーダーに鳥が写る現象は、戦時中から知られていた現象で、欧米で研究されてきたそうな。イギリスでは、最初、レーダーに映る鳥影の正体がわからないとき、それを「エンジェル」と呼んでいたそうだ。

アメリカでは、ドップラーレーダーで、コウモリと虫を追っている研究もあるようです。

日本では、それに遅れること約半世紀。やっと研究が始まった段階です。

これが、レーダー

　

クルクルまわって電波を発信・受信しています。

表示画面
　

赤矢印で示した点が鳥影。
左（南）から右（北）方向へ向かっている。
高度は、この写真だと、１２５０ｍから１５００ｍの間。
さすがに種類はわからないけれど、点の大きさ、飛ぶ時間から考えてツグミ類だろうとのことです。
この画を１秒間隔で保存しています。

これを、夜７時位から朝の５時位まで、１分間隔で手動で記録。

　

大変な作業です。

現在、別の大学のグループと共同で、記録した画像から鳥の数を数えてくれるシステムを開発中だそうです。
うまくいけば、もうひたすら一分間隔で数える必要はなくなりますね。

そして、ここでは地面から高度数千メートルの気温分布、大気の流れなども同時に記録されているので、そのデータと付け合わせれば、鳥がどんな大気状態のところを飛ぶかがわかります。

うまくまとまれば、論文になると思うので、結果が楽しみですね。

とりあえず提出

学振書類のとりあえず版を事務に提出。

あとは、事務からの修正を受けたり、中身もよーく見直して修正するのみだ。
あの、学振の電子申請ページの入力は特に緊張する。
一度完了させて、PDFを作ったらもう修正利かないし。
研究タイトルは、もう変えることができない。

そして、あともう二つ研究助成の応募書類も送る。
どっちか当たってくれると、非常に楽になるんだが。

これでフィールドが始まるまで、顕微鏡作業に集中できそうだ。

ヒタキ類

朝のちょっとした鳥見。

研究棟前にある林で立ち止まり鳥を探す。

離れたところで、センダイムシクイ、エゾムシクイ、オオルリが鳴いている。

キビタキが近くに来てくれた。

　

証拠写真。

そろそろ、エゾセンニュウも来て鳴きはじめるのかな。

2008年5月7日　北大構内
ハクセキレイ, タヒバリ, ヒヨドリ, ツグミ, ウグイス, エゾムシクイ, センダイムシクイ, キビタキ, オオルリ, ハシブトガラ, シジュウカラ, カワラヒワ, スズメ, ハシブトガラス
14種

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 ***

積丹

鳥仲間たちはゴールデンウィークを鳥見で謳歌するため、そろって天売島に行ってしまった。

鳥残された自分は、積丹半島にでも行って鳥を見に行くことにする。

積丹岬への遊歩道の入り口近くで、アリスイとツツドリの鳴き声が聞こえた。アリスイの声は結構近かったので、探したけれども見つからず。
去年、この時期に来た時もアリスイを見たので、ここにはよく来るのだろう。北海道でアリスイは繁殖しているらしいけど、ここでやっているのだろうか？

海へ出る道すがら、旅鳥のシロハラが結構目につく。そして上空にはイワツバメが飛んでいた。他には、ベニマシコ、ノビタキ、アカハラ、クロツグミと夏鳥たちが顔を見せてくれた。

　

�