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

2008年4月

獣毛観察法

29日の日記で獣毛の同定の資料を紹介したが、今日は獣毛の観察について

リンク先にあるように獣毛の観察には毛の内部構造を見る『透過標本法』と毛の表面(いわゆるキューティクル)を見る『スンプ法』があり、これを組み合わせて同定を行う。

上のリンク先では、キットなどを用いて永久標本を作る方法が紹介されているが、ただ観察するだけならもっとキットを用いず簡単に観察が可能だ。

まず、透過標本法
これは、普段のプレパラートを観察する方法と一緒で、毛を水で封入してプレパラートを作れば一時標本として十分な観察できる。

次にスンプ法
スンプ法は検索すればわかるようにSuzuki's Universal Micro-Printing Methodの略なのだが、これは毛の表面を観察するためだけではなく、繊維表面とか、金属表面とか、皮膚表面の構造を観察する目的でもつかわれているみたいだ。

いったいもともとは何の目的でこの方法が考案されたのだろうか?

以前生物技術者連絡会の29回(だったと思う)講演を聞きに行ったとき、最後の演者の方(記憶が正しければ大野氏か?)が、考案者の鈴木さんとお知り合いで、「鈴木とスタンプを合わせてスンプとした。」とおっしゃっていた気がする。

で、話を戻してと。

スンプ法の一時標本の作り方は、

http://www.kahaku.go.jp/exhibitions/vm/resource/scicom/bank/002/data1/data1p1.html
http://homepage3.nifty.com/ymorita/omnis7.htm

上記のリンクに詳しい。

まず、スライドグラスに木工用ボンドを小さいカバーグラスを使って、薄く均等に伸ばす(血液塗沫標本を作る要領)。
次に、接着剤が乾かないうちに毛を載せる。
で、乾いたらゆっくりはがす。

そうすると、ボンドが鋳型になって毛のキューティクルを写しとっているので顕微鏡で観察できる。

で、上のリンクによると、木工用ボンドよりも水絆創膏の方が使いやすいそうです。

普通に観察するだけでも楽しいので是非お試しあれ。

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

獣毛同定マニュアル

生物技術者連絡会から、『獣毛による哺乳類の同定』の資料vol4が発表されました。

哺乳類の調査は発見の困難さから調査精度の向上が難しいとされてきた分野です。それをできるだけ簡便&安価に打開するための策を練っていた同連絡会の邑井さんが出した案が獣毛の同定です。同氏は地道に研究をつづけ『獣毛による哺乳類の同定』資料を作成し公表してきました。

自分も卒論のとき獣毛の同定をする必要があったため、邑井氏の資料が非常に役に立ちました。

獣毛はフィールドサインとして地面に落ちていたり、フンに含まれていたりするので、地域に生息する動物の生息調査や食性調査には獣毛の同定法は非常に有用です。

最近はDNAをつかった手法による獣毛の種判別だけでなく、個体判別も可能でヘアトラップを用いた応用研究もされていますが、DNAを使った方法はコストがかなりかかります。

獣毛同定法は、顕微鏡があれば同定可能なので、DNAを用いた方法よりも簡便&安価です。

この資料は、出版社を通して出版されていないので、手に入れるには直接邑井氏に連絡する必要があります。連絡先は、『獣毛による哺乳類の同定 vol.4(概要版)』の中にあります。

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

更新・vegan1.11-4

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

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

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

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

とのことです。

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

肺のないカエル

2週間以上前のニュースですが、肺のないカエルがカリマンタン島で見つかったそうです。

http://www.afpbb.com/article/environment-science-it/science-technology/2377002/2820885

http://jp.reuters.com/article/oddlyEnoughNews/idJPJAPAN-31246520080410

世界初の発見だそうな。

記事には「進化プロセスの逆行」とか書かれているけど、まあ進化は複雑なものに発展していくことではないので、それは置いておいて、そもそも、「肺のない陸上の脊椎動物は珍しいのか?」ということに引っかかってしまった。

そもそも肺のない両生類は日本にも生息している。

ハコネサンショウウオは、肺のないことで有名な動物だ。
しかも、日本でかなり広い範囲に生息している動物で、昔から人に食べられているし。

と、ここまで来てall about japanで詳しく解説している記事を見つけた。

「なぜボルネオバーバーガエルには肺がないのか」という理由です。

今回の発表をしたデビッド・ビックフォード氏が言うには「浮力を捨て去るため」であろう、ということでした。

つまり、生息環境が山地の渓流で流れが速い場所だから、浮力があると流されたり、うまく泳ぐことができない。そこで、もともと両生類なわけで、皮膚呼吸ができるのだから「浮力の原因になる肺」は不要どころか邪魔。それなら捨て去ってしまえ。ということで肺がなくなってしまったのだろう、と。

なるほど、このカエルは肺を無くして浮力を無くして水に沈んで皮膚呼吸だけしているわけか。

肺呼吸を捨てて皮膚呼吸だけに頼ろうと思うのなら、水にずっと沈んでいてなおかつ酸素交換効率を上げるために表面積を上げた方がいい。そのことも上の記事で解説されているので引用。

記事にもあるのですが、これらの肺がない両生類たちは酸素と二酸化炭素の交換つまり呼吸のほとんどを「皮膚呼吸」に頼っています。

もちろん、両生類は肺がある種類でも、多くを皮膚呼吸に頼っていますが、肺のない種類はその依存度が大きいわけです。
気体は、水に溶けるときに「低温ほど溶けやすい」というルールがありますから、今回のボルネオバーバーガエルもハコネサンショウウオも冷水の環境で生活していますし、両生類が高温に弱い原因の一つでもあるわけです。

また皮膚呼吸を能率良く行うためには、皮膚の「表面積」を大きくする必要もあります。
写真を見る限りでは、ボルネオバーバーガエルもやや皮膚がたるんだ印象を受けます。
ただし、この特徴は日本をはじめ、世界中のカエルでよく知られた現象で、関東地方などの山地に生息する「ナガレタゴガエル」は繁殖期になると、ほとんどを水中で生活するため、皮膚呼吸を能率良く行う工夫として、著しく皮膚が伸びて垂れ下がって異様な印象のカエルになります。

上のリンク先にあるナガレタゴガエルの写真ではあまりブヨブヨ具合はわからないので、前に自分が撮った写真を載せておく。

Nagaretago  

 

 

 

 

 

 

 

 

 

 

上がオスで、下がメス。
普段は陸上にいるのだが、繁殖のときは水の中でオスがひたすらメスが来るのを待つために陸に上がらなくていいように皮膚呼吸をがんばるために皮膚をブヨブヨにして表面積を上げていると考えられる。

そんなわけで、この新発見のカエルこと、Barbourula kalimantanensis は、今までの発見されている動物の事実をまとめて考えるとそれほど驚かずに「さもありなん」。と思ってしまう。

で、このカエルのニュースが報道で流れたのは、4月9日。研究チームが発表した。ニュースの記事によると「カレントバイオロジー誌」つまりCurrent Biologyなんだと思うんだが、発表した当時最新にあたるVolume 18, Issue 7, Pages 471-552 (8 April 2008) に該当記事がないんだけれども。

これはつまり、論文がパブリッシュされたからではなく、アクセプトされたから報道発表したということなのかな?

開発でいつ絶滅してもおかしくないから、いち早く発表したかったのかもしれない。

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

動物に遭遇

数日前TOEFの林道を走らせていたらシカを発見。

Sikamesu1 

 

 

 

 

 

 

 

 

 

 

道脇にはいってこちらを警戒。

Sikamesu2  

 

 

 

 

 

 

 

 

 

 

そして逃げて行った。

Sikamesu3  

 

 

 

 

 

 

 

 

 

 

シカのオスも見れた。

Sikaosu  

 

 

 

 

 

 

 

 

  

角が生え始めている。

シカ程度なら大したことないのだが、林道を走らせているとその先になんかよくわからない物体があるなあ~と思ってみていたらタヌキだった。

しかも2度見た。同じ個体かもしれない。

Tanuki  

 

 

 

 

 

 

 

 

 

 

斜面の上の方に逃げて行き、尾根向こうで見えなくなった。

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

アイヌネギ

研究室のM2と二人でサンプリングついでにギョウジャニンニクもサンプリングしようということになった。

ギョウジャニンニクが採れる場所は知らないが、湿った日当たりの良い所に生えているそうなので、水辺周辺を攻めてみた。

Hitobiro

 

 

 

 

 

 

 

 

 

 

 

小さな群落を見つけた。

スーパーでよく見るギョウジャニンニクと違いすでに葉が開いている。もう時期が遅いようだ。

一部シカに食われたような跡があったが、ほとんど人の手は入っていないようだ。

ホクホクと採集。根本から抜くと次の年に生えてこないようなので、爪を立てて引きちぎるように収穫をしていると手がニンニクくさくなってくる。

ねこそぎとってしまうのは気がひけたので、次の場所に移動する。

ここら辺にはさらにありそうなので、もうひと群落あてて収穫したい。

 

 

 

川沿いをちょっと歩くと見つかった。

Hitobiro2  

 

 

 

 

 

 

 

 

 

 

あたり一面ギョウジャニンニク

Hitobiro3  

 

 

 

 

 

 

 

 

 

 

もう、ウハウハ。

Hitobiro4  

 

 

 

 

 

 

 

 

 

 

袋いっぱいになったところで、まだまだいっぱいあったけど、今回はこんなもんかということで収穫終了。

二人で山分けしたうえで、さらに人におすそ分け。醤油漬けにするともつみたいなので保存食にしておこうかな。

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

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

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

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)線形方程式系は正確に特異です

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

花粉症

19日から花粉症が発症した。

目がかゆい上に鼻水が止まらない。

関東にいた時と比べて期間も短く、症状も軽いのでまだ楽なのだが、つらいものはつらい。目が重たく眠い感じだし、花がムズムズするので、集中力がわかない。

シラカバ花粉の悲惨飛散状況は、ここで提供されている。これを見る限りあと1か月はこの症状が続くことになりそうだ。

屋上に上がってシラカバの写真を撮ってみた。

Sirakaba  

 

 

 

 

 

 

 

 

 

Sirakaba2  

 

 

 

 

 

 

 

 

 

Sirakaba3  

 

 

 

 

  

 

 

 

 

見るからに花粉だしますっていう。邪気を放っていた。

実際はたくとモアモア~と黄色い粉を放った。

ううう・・・・・ずっと室内にいたから鼻が落ち着いていたのに、写真を撮りに外に出たばっかりにまた鼻水がぶり返した。

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

シギチを求め新川河口へ

新川河口へシギチがいないか見に行く。

新川は車で砂浜を走って、車上で鳥が見られる便利な場所である。
しかしながら、油断するとスタックしてひどい目にあう。
実際自分は、スタックして痛い目にあったことがある。

なので、もう2度と一人で新川河口へ車で行くものかと思っていたが、やっぱり時期になると誘惑に勝てなくなる。

複数人で行けばいいのだろうが、ぱっと行ってぱっと帰ろうと思うと一人の方がなにかと都合がいい。

車を流したかんじだとシギチはいなかった。
まだ、時期が早いのかもしれない。

河口にはまだカモメは残っていた。

Kamome まんなかのカモメは違和感を感じた。他の個体よりも嘴が長めだし、目が赤っぽいし。なにか珍しい種なのかなあと追って写真をたくさん撮っておく。

Kamome2 同じ個体(右)。

よくわからんが、まあ個体差の範囲内なんでしょう。

Kaiju なぞの海獣が打ち上げられていた。

アザラシではないなあ。アシカかな?と思ったけど、帰ってよくよく考えてみれば、日本には「アシカ」と名前のつく哺乳類は現存しないはずで、もっとよく調べればよかったと後悔。

顔もちゃんと写真撮っておけばよかった。

トドの子供に一票。

記録 2008年4月20日 新川河口
シロカモメ、カモメ、ウミネコ、ヒバリ、ハクセキレイ

30分くらいで撤退。
あまりまじめに鳥見ていません。

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

フクジュソウ

今年度初めてTOEFに行く。

TOEF利用者がみんなで集まって、ゴミひろいとか消火訓練とか。

雪は溶けているが、まだまだ、葉っぱは出ていない。

自分の調査プロットを見に行くが、雪でつぶれるなど、特に異常はないようだ。

Hukuju  

 

 

 

 

 

 

 

  

フクジュソウが咲いていた。

植物の専門の先輩たちが、動物専門の後輩に

「フクジュソウはおいしいよ。」
「葉っぱが甘いよ。」
「気持ち良くなるよ。」
「天国見れるよ。」

なんていっていた。

さすがにだまされません。

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

Crawley本の翻訳

共立出版から来月Crawley本

Statistics: An Introduction Using R Statistics: An Introduction Using R

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

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

統計学:Rを用いた入門書 統計学:Rを用いた入門書

著者:Michael J.Crawley
販売元:共立出版
Amazon.co.jpで詳細を確認する

やるな共立出版

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

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

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

AICの調節

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

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

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

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

D/2 = LI - Ls
D/2 + Ls = LI
Ls = -D/2 + LI (式1)

となるわけです。

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

式1に代入して

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がそうなわけで、以前からデビアンスに関する記述を読んでもわかったふりして放っておいたんだけど、やっと理解できた。

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

生物地理学

研究室のお茶部屋にこんな本がおいてあった。

地球と生命の進化学 地球と生命の進化学
販売元:セブンアンドワイ
セブンアンドワイで詳細を確認する
 

地球の変動と生物進化 地球の変動と生物進化
販売元:セブンアンドワイ
セブンアンドワイで詳細を確認する
 

動物地理をやっている人の発表を聞いていると、
「日本列島は○○代には大陸と陸続きで・・・」
みたいな話をよくする。そういう地殻変動とかの話はよくしらないので、どんな本で勉強すればいいんだろう?と前々から思っていた。この本にはそこら辺のことも書いてありそうなのでなかなかよさげだ。

それにしても北大出版は生物地理学の本が好きだな。

動物地理の自然史―分布と多様性の進化学 Book 動物地理の自然史―分布と多様性の進化学

著者:増田 隆一,阿部 永
販売元:北海道大学図書刊行会
Amazon.co.jpで詳細を確認する

 

こんな本が2005年に出版されているし(これは自分も購入した)、来年には生物地理の『○○版』も出版される予定みたいだし(まだ言っちゃいけないかもしれないので、いちおー伏せ字にしておく)。生物地理はブームですね。

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

sem研究(3)線形方程式系は正確に特異です

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

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

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

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

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

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

system is exactly singular

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

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

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

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

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

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

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

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

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

East_scrofaさんのところでまとめられているように、GLMおよびGLMMの関数群にはAICを求める対数尤度の基準が2通りあるようで、違う求め方で求めた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)
}
a1-rd1 #AICとresidual devianceの差

#betabin
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))
y2<-floor(runif(6,min=1,max=20))
tdata1<-data.frame(x1,ID1,y1,y2)
res1<-betabin(cbind(y1,y2)~1+x1,random=~ID1,data=tdata1)
a1[i1]<-AIC(res1)@istats$AIC
rd1[i1]<-deviance(res1)
}

a1-rd1 #AICとresidual devianceの差

#glm.nb
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<-glm.nb(y1~1+x1,data=tdata1)
a1[i1]<- res1$aic
rd1[i1]<-res1$deviance
}
a1-rd1 #AICとresidual devianceの差

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

宮島沼:落雁

土曜日は午後から宮島沼へ珍ガン探しと塒入り(落雁)観察。

Gan車でガンがいる田圃を探してガンを観察。特にアイリングがあるかどうかをよーく見る。
ネイチャーセンターの情報だと、カリガネ、シジュウカラガンが入っているそうだが、見つけることは出来なかった。

Binihiwa今年大当たりだったベニヒワがこんなところにもいた。
曇天だったのが残念です。

デジスコすればよかったとちょっと後悔。
逃げてもあまり遠くには逃げなかったんだよね。

Hibari  

 

 

 

 

 

 

 

 

上空でひたすらさえずっているヒバリ。
この子はわりと丸っこく地面でじっとしてた。メス?

Ootaka  

 

 

 

 

 

 

 

 

運がいいことにオオタカを間近で見ることができた。
時間的に暗いのでISO感度高め。
車の中で見られるギリギリの距離。
車中からスコープ使って観察。
とってもかっこいい。
10分くらい観察できた。
重ね重ねデジスコすれば・・・

時間が迫ってきた。
ガンが帰る時間だ。

Ganko 

ガンが宮島沼に向けて帰っていく。
車を走らせながら撮影。ガンとほぼ同じスピード(30km/hくらい)。

宮島沼で待機。

どんどんガンが帰ってくる。

 Ganko2 

 

 

 

 

 

 Ganko3  

落雁

この日はめちゃくちゃ寒かった。
風が冷たくて冷たくて・・・
明るいうちにスコープでシジュウカラガンとカリガネを探すも見つからず。
とにかく冷たい風で集中ができない。
もうちょっと暖かい格好をしてくればよかった。
6日の朝がそんなに寒くなかったからなめてた。

塒入りはバラバラと群れが入ってくるので迫力は朝の塒発ちにはかなわないが、夕方の落雁はなんとも趣がある。曇天で夕日がでていないかったのが残念だけれども、それでも十分いい。この感覚は平安時代から変わらないわけです

落雁の動画も撮りました。風が強いのがわかると思います。

動画を再生する(3.84MB)

記録  2008年4月12日宮島沼とその周辺
アオサギ、マガン、ヒシクイ、オオハクチョウ、マガモ、コガモ、ヒドリガモ、オナガガモ、カワアイサ、オオタカ、ノスリ、ハヤブサ、アカゲラ、ヒバリ、ハクセキレイ、カワラヒワ、ベニヒワ、ムクドリ、ヒヨドリ、ハシボソガラス、ハシブトガラス
21種 ガンは54000羽(ネイチャーセンター情報)

第一目標のカリガネが見つけられなかったけれど、ベニヒワも見られたし、オオタカを近くで長時間観ることができたから満足。

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

lmer研究(4)AICについて

East_Scrofaさんのブログ問題提起され粕谷先生答えで一応の解決をみたglmとglmmMLのAIC問題

はてさてではこれまたGLMMでよくつかわれるlmer()ではどうなんでしょ?

ということで粕谷先生と同じ実験を試みた。

n1<-100
a1<-numeric(n1)
a2<-numeric(n1)
rd1<-numeric(n1)
x1<-seq(6)
ID1<-seq(6)
for (i1 in 1:n1){
y1<-floor(runif(6,min=1,max=20))
y2<-floor(runif(6,min=1,max=20))
tdata1<-data.frame(x1,ID1,y1,y2)
res1<-lmer(cbind(y1,y2)~1+x1+(1|ID1),family=binomial,data=tdata1)
res2<- summary(res1)
a1[i1]<-res2@AICtab$AIC
rd1[i1]<-res2@AICtab$deviance
res.ML<-glmmML(cbind(y1,y2)~1+x1,family=binomial,cluster=ID1,data=tdata1)
a2[i1]<- res.ML$aic
}
a1-rd1 #AICとresidual devianceの差
a1-a2  #lmerとglmmMLのAICの差

こんなふうに100回繰り返し、lmerのAICとresidual devianceの差とlmerとglmmMLのAICの差を調べてみた。

結果は

>a1-rd1 #AICとresidual devianceの差
  [1] 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
[35] 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
[69] 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
> a1-a2 #lmerとglmmMLのAICの差
  [1]  1.324157e-05  5.179822e-05  2.441203e-04  1.464408e-04
  [5]  4.384840e-06  1.745354e-03  1.205713e-03  3.590084e-05
  [9]  2.175130e-05  1.936024e-04  8.448943e-04  1.256820e-03
[13]  3.406060e-05  1.891679e-03  1.235944e-04  2.394006e-05
[17]  5.712714e-05  5.474305e-05  3.625226e-06  2.454966e-06
[21]  1.863478e-04  4.860338e-06  3.574641e-05  4.384132e-03
[25]  2.333348e-05  1.103509e-05  8.370281e-04  9.335986e-04
[29]  3.744454e-04  2.467052e-05  1.299870e-04  2.425289e-04
[33]  4.815803e-05  2.758913e-05  7.678939e-05  3.804263e-04
[37]  1.090735e-05  3.360833e-03  1.225639e-06  1.036787e-03
[41]  1.183095e-03  5.349143e-04  2.272831e-08  5.549182e-09
[45]  3.773663e-05  1.222656e-03  1.175915e-04  1.981372e-04
[49]  1.217419e-03  1.497944e-05  1.573509e-03  3.522962e-04
[53]  9.569729e-08  1.575468e-05  1.065551e-04  6.734617e-06
[57]  1.196103e-06  4.256138e-05  5.318053e-06  2.478905e-04
[61]  7.835918e-05  3.917400e-04  1.291786e-04  9.974890e-09
[65]  3.838827e-05  2.277255e+02  4.037893e-08 -6.601926e-08
[69]  9.562574e-05  3.439671e-05  6.205722e-06  1.291586e-05
[73]  1.528715e-04  4.164637e-05  1.105234e-04  1.005523e-03
[77]  1.025349e-06  2.781931e-03  9.924978e-05  6.234150e-06
[81]  1.031122e-09  7.213863e-04  5.700328e-04  4.373709e-04
[85]  3.168673e-04  2.794296e-04  2.029505e-04  2.377327e-04
[89]  1.202936e-03  6.089783e-05  4.066697e-05  6.186540e-04
[93]  1.336513e-04  3.112440e-09  1.673464e-04  2.553641e-03
[97]  1.418265e-03  8.741700e-05  1.260778e-03  1.839569e-04

つまりlmerとglmmMLは同じ方法でAICを求めているらしい。AICもほぼいっしょ。たまーに推定に大きなズレが生じる場合があるらしい(66番目の結果とか)。

つまりlmer()とglm()のAICを単純に比べることはできないようです。

ちなみに、これはロジスティック回帰の場合のみかもしれないので、同様にpoissonでもやってみる。

n1<-100
a1<-numeric(n1)
a2<-numeric(n1)
rd1<-numeric(n1)
x1<-seq(6)
ID1<-seq(6)
for (i1 in 1:n1){
y1<-floor(runif(6,min=1,max=20))
tdata1<-data.frame(x1,ID1,y1,y2)
res1<-lmer(y1~1+x1+(1|ID1),family=poisson,data=tdata1)
res2<- summary(res1)
a1[i1]<-res2@AICtab$AIC
rd1[i1]<-res2@AICtab$deviance
res.ML<-glmmML(y1~1+x1,family=poisson,cluster=ID1,data=tdata1,start.sigma=2)
a2[i1]<- res.ML$aic
}
a1-rd1 #AICとredisual devianceの差
a1-a2  #lmerとglmmMLのAICの差

> a1-rd1
  [1] 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
[35] 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
[69] 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
> a1-a2
  [1] 5.945945e-06 1.481231e-07 6.689094e-06 3.610293e-07 1.303709e-05
  [6] 1.443701e-05 2.868542e-03 1.651486e-07 7.309863e-05 1.724372e-08
[11] 3.623184e-05 7.338124e-06 1.716323e-03 9.544965e-04 1.597265e-05
[16] 1.373632e-03 2.274728e-05 2.113061e-04 3.111850e-06 2.044724e-03
[21] 6.613532e-05 1.055137e-03 1.613774e-06 2.248374e-06 1.667302e-04
[26] 1.415890e-04 1.534181e-04 4.010243e-05 1.126477e-04 1.665226e-05
[31] 2.344202e-08 7.540830e-06 3.225858e-04 1.712079e-04 6.658772e-04
[36] 1.575257e-04 1.089677e-05 2.162002e-04 1.813634e-05 6.073717e-06
[41] 1.206985e-04 1.279918e-03 9.174171e-06 2.327433e-05 4.515342e-04
[46] 8.745013e-08 1.407305e-04 3.531264e-05 6.483757e-06 1.210553e-04
[51] 2.590982e-05 2.293854e-06 2.042587e-04 5.248377e-04 6.756842e-06
[56] 1.770799e-04 4.533351e-09 1.537446e-03 2.865895e-08 5.706534e-07
[61] 6.225319e-06 1.752127e-04 4.919500e-05 6.520251e-05 3.153157e-04
[66] 1.218186e-05 1.685016e-04 1.455618e-05 3.178621e-04 5.736318e-05
[71] 9.202713e-04 1.029569e-05 2.974918e-08 7.128886e-04 6.555559e-05
[76] 2.713900e-04 5.879766e-04 1.432476e-05 1.155867e-08 8.446181e-05
[81] 3.857819e-04 5.713443e-04 4.592716e-09 1.502404e-05 2.842325e-04
[86] 1.776855e-07 6.350424e-05 4.541055e-06 4.022462e-05 2.729961e-07
[91] 3.643668e-04 2.859998e-05 9.815364e-05 4.802853e-05 1.243280e-03
[96] 6.598562e-04 1.872406e-05 5.225837e-07 2.062394e-03 9.789729e-04

つまり同じ結果になりました。

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

雪ぐされ菌

昼休みが終わるころ、微生物生態の人たちが何やら外に出る雰囲気をしていたので、
「お昼御飯ですか?」
と聞いたら、
「菌探しをする」
そうなので、同行した。

雪ぐされ病の原因となる菌の研究者が研究所にやってきたので、探し方を教えてもらうのだそうだ。

場所は研究所隣の牧場。
イネ科植物の雪でつぶされていた枯草を引きちぎってみてよ~く見ると黒や赤の小さい粒がある場合がある。

これは雪ぐされ菌の菌核なのだそうだ。

これが11月くらいになると、細長く伸びた菌体になるらしい。
身近な所にはまだ見ぬ不思議なものがたくさんあるようだ。
11月にこれを見つけるのが楽しみだ。

ちょっと場所を変えて地面を地面をほっくりかえして遊ぶ。

変形菌研究者の方が早速変形菌を見つける。

自分はダンゴムシを見つけた。

かつて北海道に分布されていないとされていたダンゴムシ。
これも温暖化の影響だろうか。

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

日本の哺乳類学

東京大学出版から『日本の哺乳類学』という本が出版される。

日本哺乳類学会が、IMC9の記念出版をすることになったらしい。

目次を見ると各方面の第一人者の方たちが執筆に当たっている。
ぜひ読みたい。出版が待ち遠しい。
三巻セットで1万円超えるのがちょっと痛いですけどね。

また、これとは別に哺乳類学者のバイブルになる図鑑も現在制作が進められている。もともと図鑑の方が先に出版される予定だったが、まだまだ予定は未定らしい。
ちなみにこっちは英語で書かれている。

一種につき膨大な文献が収集され記載されているようだ。執筆者の苦労が偲ばれる・・・

日本の哺乳類学〈1〉小型哺乳類 日本の哺乳類学〈2〉中大型哺乳類・霊長類 日本の哺乳類学〈3〉水生哺乳類

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

sem研究(2)記法とか雑多なまとめ

semの命令記述の仕方。Rjpwikiの該当ページやヘルプをよ~く見ればだんだんわかってくる。

semの操作で日本語の説明がある本は現状では、

Rによるやさしい統計学 Book Rによるやさしい統計学

著者:山田 剛史,村井 潤一郎,杉澤 武俊
販売元:オーム社
Amazon.co.jpで詳細を確認する

 

この本にすこしだけ書かれている。

自分が詰まった点を書いておく。(Rjpwikiの該当ページを例にする。)

まず、specify.model()の記法がわからなかった。

これは、順番に、「"パスの方向", "変数名", "初期値"」であった。
初期値はNAにしておけば自動的に決めてくれる。収束しない場合は初期値をいじる必要がある。

sem()関数の必須の引数は、
sem(パスのモデル式, 相関行列, データ数)  とする。

話は変わって
summary(ans)をみると

Parameter Estimates
       Estimate Std Error z value Pr(>|z|)                        
gamma  -0.22909 0.151662  -1.5105 1.3092e-01 満足度 <--- 小型軽量 
delta   0.45680 0.088107   5.1847 2.1641e-07 満足度 <--- 操作性   
kappa   0.74739 0.148927   5.0185 5.2069e-07 満足度 <--- 可搬性   
alpha  -0.41148 0.091601  -4.4921 7.0521e-06 操作性 <--- 小型軽量 
beta    0.84217 0.054192  15.5405 0.0000e+00 可搬性 <--- 小型軽量 
theta1  1.00000 0.142162   7.0342 2.0037e-12 小型軽量 <--> 小型軽量
theta2  0.83068 0.118097   7.0339 2.0079e-12 操作性 <--> 操作性   
theta3  0.29074 0.041353   7.0308 2.0537e-12 可搬性 <--> 可搬性   
theta4  0.63666 0.090519   7.0334 2.0155e-12 満足度 <--> 満足度 

赤で示したように、gammaはwald検定では有意でないので、この変数も抜いてもいいかもしれない。

実際やってみると

model2 <- specify.model()
操作性 -> 満足度, delta, NA
可搬性 -> 満足度, kappa, NA
小型軽量 -> 操作性, alpha, NA
小型軽量 -> 可搬性, beta, NA
小型軽量 <-> 小型軽量, theta1, NA
操作性 <-> 操作性, theta2, NA
可搬性 <-> 可搬性, theta3, NA
満足度 <-> 満足度, theta4, NA

ans2 <- sem(model2, r, nrow(camera))
bic<-matrix(c("ans","ans2",summary(ans)$BIC,summary(ans2)$BIC),2,byrow=T)
rownames(bic)<-c("model","BIC")
bic
      [,1]                [,2]               
model "ans"               "ans2"            
BIC   "-4.33514946688171" "-6.68460861285924"

BICではgammaを抜いたans2の方がBICが低い

尤度比検定もすると

> anova(ans,ans2)
LR Test for Difference Between Models

        Model Df Model Chisq Df LR Chisq Pr(>Chisq)
Model 1        1     0.27002                      
Model 2        2     2.52573  1  2.25571     0.1331

やはりgammaは抜いてもいいようだ。

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

マガンの塒発ち

4月6日に行った宮島沼の塒発ちの動画をアップ。

元動画は65MBくらいだったけれど、圧縮して8MB程度にした。

撮影したカメラの性能のためもともと動画の質は良くありません。

周囲の人の声も入っています。

動画を再生する(ダウンロード)

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

更新・sem

semパッケージが5日に新しいバージョンに更新されたようだ。

semパッケージ作者のサイトで特に何を更新したのかアナウンスはされていないようだ。

早速アップデートして使ってみると

> library(sem)
Warning message:
パッケージ 'sem' はバージョン 2.6.2 の R の下で造られました

と出るけど(自分のRバージョンは2.6.1)、特に問題はなさそうだ。

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

うまく計算ができない

semを使ってパス解析を試みる。

モデルを作って、いざ解析。

出てくるエラーは、

以下にエラー solve.default(C) :
   Lapack routine dgesv: 線形方程式系は正確に特異です
追加情報:  Warning messages:
1: In sem.mod(model.shrew, shrewcor, 243) :
  The following observed variables are in the input covariance or raw-moment matrix but do not appear in the model:
wall

2: In sem.default(ram = ram, S = S, N = N, param.names = pars, var.names = vars,  :
  The following variables have no variance or error-variance parameter (double-headed arrow):
shrew, huto, litter, isopoda, spider, collembola, chiropoda, sipider
The model is almost surely misspecified; check also for missing covariances.

ちょっと手を加えると(適当にfixed.xなるパラメータを入れる)

以下にエラー sem.default(ram = ram, S = S, N = N, param.names = pars, var.names = vars,  :
  The model has negative degrees of freedom = -5
追加情報:  Warning messages:
1: In sem.mod(model.shrew, shrewcor, 243, fixed.x = c("shrew", "huto",  :
  The following observed variables are in the input covariance or raw-moment matrix but do not appear in the model:
wall

2: In sem.default(ram = ram, S = S, N = N, param.names = pars, var.names = vars,  :
  The following variables have no variance or error-variance parameter (double-headed arrow):
litter, sipider
The model is almost surely misspecified; check also for missing covariances.

自由度が足りなくなる。

仕方がないのでパスを減らしてもまた前者のエラーがでる。

うーーーむ。

ちょっと先生と相談する機会があったのでいろいろ聞いてみると、データそのものが悪くて収束しないんじゃないかということになった。

もはや打つ手なしか・・・

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

宮島沼から鵡川

朝3時15分集合。総勢12名での強行行程。

まずは宮島沼。
4月1日のときと違い、氷が完全に溶けている。

Miyajima

 

 

 

 

 

 

 

 

 

 

場所にスタンバイしてから数分とガン達が飛び出しはじめた。
前回の経験からして写真がほとんど撮れないことはわかっていたので、動画撮影に専念する。まあまあ満足いく動画が撮れた。そのうちアップ予定。

マガン達は一回すべてが飛んでも半分くらい沼に戻ってそのまま休息し始めた。結構明るくなっても飛び立たない。

本当は、田圃に行ってカリガネ探ししたいんだけど、リーダーはじめその他のみなさんはあまり興味はないらしく、さっさと次の場所へ移動。

鵡川方面へ。いきなり鵡川方面に行くのもなんなので、通り道にある野幌森林公園へ。野幌はまだ雪が溶けてなく、歩きにくかったのでさっさと退散。
そんな中でもツミのメスを見ることができた。ラッキー。

千歳川に到着。ヤマセミ探しをするも見られたのはカワセミだけ。ここもさっさと移動。

ウトナイ湖着。白鳥餌付け場所で軽く鳥見。

Hagehaku  

 

 

 

 

 

 

 

 
 

こんな白鳥をみつけた。首の一部がきれいに剥げてる。気持ち悪い。 

勇払マリーナで鳥を見ながら昼食。
ウミアイサが求愛行動してた。
一瞬ツメナガホオジロが飛んだ。
眼デジ部隊の一人が飛翔写真をとって同定。それにしてもすごい瞬発力だ。 

鵡川で、ちょっと昼寝してカモメウォッチング。うーむこれと言って面白い鳥がいない。
飽きはじめる人多数。シギチ見つけるもシロチドリ・・・
Gundesi  

 

 

 

 

 

 

 
 
 

 

シロチドリ撮影に挑む眼デジ部隊 

あとは漁港に行ってカモメとカモ見て終了。
温泉に入って帰る。

家に着いたのは22時過ぎ。

走行距離300km越え。もうヘロヘロ。

記録 2008年4月6日 宮島沼~野幌~千歳川~ウトナイ~勇払~鵡川
アビ、ハジロカイツブリ、ウミウ、ダイサギ、アオサギ、マガン、ヒシクイ、コブハクチョウ、オオハクチョウ、コハクチョウ、オシドリ、マガモ、カルガモ、コガモ、ヨシガモ、オカヨシガモ、ヒドリガモ、オナガガモ、ハシビロガモ、キンクロハジロ、スズガモ、クロガモ、ビロードキンクロ、シノリガモ、ホオジロガモ、ウミアイサ、カワアイサ、トビ、オジロワシ、オオワシ、ツミ、ノスリ、チュウヒ、ハヤブサ、キジ、シロチドリ、ハマシギ、ユリカモメ、セグロカモメ、オオセグロカモメ、ワシカモメ、シロカモメ、カモメ、ウミネコ、キジバト、フクロウ、カワセミ、アカゲラ、コゲラ、ヒバリ、ハクセキレイ、ヒヨドリ、モズ、ツグミ、エナガ、ハシブトガラ、ヤマガラ、シジュウカラ、ゴジュウカラ、キバシリ、ツメナガホオジロ、カワラヒワ、イカル、スズメ、ムクドリ、カササギ、ハシボソガラス、ハシブトガラス 68種

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

sem研究(1)例題実行

semパッケージの関数sem()はexample(sem)とやっても作動しない。

そこでhelp画面をコピペして実行する。もっとも、実行結果は、すでにhelp画面にでているので、生成された結果から実際にパス図を描いてみる。

semのパス図を描くには別途Graphvizをインストールする必要がある。Graphvizの日本語の説明はここにある

Dhp1 

 

 

sem.dhp.1

Dhp2   

 

 

sem.dhp.2

 

Wh1 

 

 

sem.wh.1

 

Wh2 

 

 

sem.wh.2

 

 

Thur 

 

  

 

 

 

 

 

 

 

sem.thur 

   

 

Kerch  

 

 sem.kerch

 

 

 

Mcardle 

 

 

sem.McArdle 

 

 

 

 

こうしてkerchの例をみると、潜在変数なしのモデルも組めていることから、いわゆる「パス解析」もできるようだ。

CRANではないが、BioconductorRgraphvizというパッケージ(別途graphパッケージをインストールする必要あり)がある。これを使えばR上でパス図を書くことができるのだろうか?マニュアルを読まなきゃだな。

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

いまいちよくわからない

パッケージsemでいわゆる「パス解析」はできるのか?

共分散構造分析のように実測値に潜在変数を噛ませなくてはいけないのだろうか?

調べてもよくわからないんだよな。

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

プログラムの考え方

20×20のマス目に10%、20%・・・とランダムに色を塗ることが必要になった。

エクセルで主導でやるのは面倒だし、必ず恣意的になってしまうので、コンピュータでなんとかしたいと思った。

自分が使えるプログラム言語はRだけなので、Rでなんとかしようと策を巡らす。

確かRtipsに似たようなことが書いてあったと思い、Rtipsをめくってみると、image()で行列を可視化すれば簡単にできそうなことに気づく。

あとは、どうランダムに点を発生させるかだ。

乱数発生はRの得意分野なので、乱数を発生させてやれば簡単だ。しかもRtipsにはコイントスのプログラム例が載っているからそれを応用すればいいと判断し、プログラムを書いてみた。

denplot<- function(number,probability){
##nuber=マス目の数、probability=確率
     random<- function(n, p){
       vec<- numeric(n)

      for(i in 1:n){
       x<- runif(1)
       if (x <= p) vec[i]<- 1
       else        vec[i]<- 0
                    }
       return(vec)}
      
co<- random(number,probability)
mat<- matrix(co,sqrt(number))
image(mat,col=c(0,1))
return(sum(co)/number)
            }

で、これを実行すると、

> denplot(400,0.3) #400個の点で30%の密度
> [1] 0.3125

Denplot 

 

 

 

 

 

 

 

 

 

 

 

こんな感じになる。[1]0.3125 は実際には31.25%の密度の図になったことを示す。

実際は大体30%くらいでいいので、自分の求めるものはこれでもいいのだが、やはりきっちり30%にしたい。

そこでwhile文を用いて、30%になるまで回し続けると、計算は終わらずRが暴走。ひたすら新たな図が高速で生成され、スロットのようだった。なかなかきっちり30%にはならないようだ。

うーむ、やはり最初から0と1を決まった数生成して(400マス目で30%なら1を120個)並べ替えた方がいいなと感じるがいい方法が思いつかないので家に帰ってねることにした。(この時点で夜中の1:00)

で、自転車に乗りこむとき、ふと「あっ、sample()ってのがあった。」と思いつく。まさに「三上」だと思った(この場合「馬上」か?)。

家に帰って、rep()とsample()を使ってプログラムを書いてみた。

plotden<- function(n,p){
  cc<- c(rep(1,n^2*p),rep(0,n^2*(1-p)))
  dd<- sample(cc,n^2)
  mat<- matrix(dd,n)
  image(mat,col=c(0,1),xaxt="n",yaxt="n",
   main=paste("Matrix",n,"*",n,"p=",p))}

楽勝じゃん・・・

仕上げに手直しして、全体のマス目の数ではなく、一辺のマスの数と確率を入れるようにして、さらにグラフに題名を入れるようにした。

> plotden(20,0.3)

とすると

Plotden  

 

 

 

 

 

 

 

 

 

 

 

いとも簡単に目的の図を得ることができた。

プログラムのときは思考を柔軟にしなくてはいけませんね。

それにしてもまたしてもRの便利さをひしひしと感じた。

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

助成金

助成金の書類を書く季節になった。
生態系関連の助成金に関しては村田先生のブログにまとめられているのでリンクをはっておく。

http://blogs.yahoo.co.jp/haemoproteus_gallinulae/3193113.html

その他にある助成として、
藤原ナチュラルヒストリー振興財団

物質循環とを扱うならエスペック

北海道限定だけどノーステック

もある。

話は変わって、村田先生のブログの過去ログをあさっていたらこんなのがでてきた。

院生心得12ヶ条
http://blogs.yahoo.co.jp/haemoproteus_gallinulae/2716911.html

学部研究室生心得17ヶ条
http://blogs.yahoo.co.jp/haemoproteus_gallinulae/2657123.html

肝に銘じておきます。

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

宮島沼

午前3時起床。

同行者と3時半に待ち合わせ、出発。

空が白んできたころに宮島沼に到着。

ちょうど観察位置についたころに期待のイベントが始まった。

Ganko 

 

 

 

 

 

 

空一面がマガンの群れ(写真がブレブレで申し訳ない)

頭上を大歓声をあげて飛ぶマガンの群れに久しぶりにしびれた。

Ojiro  

 

 

 

 

 

 

 

 

 

オジロワシもいた(いい写真がない)。

6羽以上いて、ほとんどが若鳥。盛んに飛び回り魚を獲ったり、相手の持っている魚を取り合ったりと元気だった。

そして7時半には現地出発し、9時ころには何食わぬ顔で研究室にやってきた。ちょっと眠いけど。

記録 2008年4月1日 宮島沼
オジロワシ、トビ、アオサギ、マガン、ヒシクイ、オナガガモ、ヒドリガモ、カワアイサ、オオハクチョウ、ハクセキレイ、ヒバリ、スズメ、ムクドリ、ヒヨドリ、ハシボソガラス、ハシブトガラス
16種
マガンの大群の堪能でおなかいっぱい。ヒバリはこの地域ではさえずり初認。

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

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