Ads by Google
R で始める 10 行コードからの個体群動態乱数シミュレーション
#（以前、オフラインで披露したものですが、バージョンアップしています）
# 注1："#" 以降、行末まではコードとして読まれない
# 全文を R の画面にコピペして走らせれば動作します、
# コピペ & Enter を繰り返すたび、毎回異なる結果が得られます。
# 出力されるもののうち、数値は絶滅率を示します。
# 図は、全シミュレーションの結果を一枚の図に書き込んだものです。
# 最後の行は、実行例に過ぎないので、数値を入れ替えて試してもらえたらと思います。
# init: 初期個体数； surviv：生存率； repro：繁殖率； egg：卵数；
# gen：世代数； times：シミュレーション回数
pop.dyn <- function(init, surviv, repro, egg, gen, times) {
Popul <- matrix(init, gen+1, times) # シミュレーション結果を格納する入れ物
Extinct <- 0 # 絶滅率を格納する入れ物
for (k in 1:times) { # 繰り返しを意味するコード（シミュレーション回数)
Ncohort <- init # 各世代の個体数（初期値として、初期個体数を入れる)
for (i in 1:gen) { # 世代数を指令
Ncohort <- rbinom(1, Ncohort, surviv) + rbinom(1, Ncohort, repro)*rpois(1, egg)
# 二項乱数、ポアソン卵数を用い、生存数 + 繁殖個体数 × 卵数 を表記
Popul[i+1,k] <- Ncohort } # 各世代の個体数を記録（初項は init のまま）
Extinct <- if(Ncohort == 0) Extinct + 1 else Extinct }
# 絶滅回数（最終世代までに絶滅したシミュレーションの回数)
matplot(Popul, type="l") # 全シミュレーションの結果を図示
Extinct/times } # 絶滅率を出力（絶滅したシミュレーション回数/総シミュレーション回数）
# 実行例、実行する度に結果は変動します
pop.dyn(50, 0.5, 0.5, 1, 10, 50)
# 初期個体数, 生存率, 繁殖率, 卵数, 世代数, シミュレーション回数
# 注 2：コード部分、数えてみたら 11 行あるじゃないかといった突っ込みは無しの方向でお願いします。
# 注1："#" 以降、行末まではコードとして読まれない
# 全文を R の画面にコピペして走らせれば動作します、
# コピペ & Enter を繰り返すたび、毎回異なる結果が得られます。
# 出力されるもののうち、数値は絶滅率を示します。
# 図は、全シミュレーションの結果を一枚の図に書き込んだものです。
# 最後の行は、実行例に過ぎないので、数値を入れ替えて試してもらえたらと思います。
# init: 初期個体数； surviv：生存率； repro：繁殖率； egg：卵数；
# gen：世代数； times：シミュレーション回数
pop.dyn <- function(init, surviv, repro, egg, gen, times) {
Popul <- matrix(init, gen+1, times) # シミュレーション結果を格納する入れ物
Extinct <- 0 # 絶滅率を格納する入れ物
for (k in 1:times) { # 繰り返しを意味するコード（シミュレーション回数)
Ncohort <- init # 各世代の個体数（初期値として、初期個体数を入れる)
for (i in 1:gen) { # 世代数を指令
Ncohort <- rbinom(1, Ncohort, surviv) + rbinom(1, Ncohort, repro)*rpois(1, egg)
# 二項乱数、ポアソン卵数を用い、生存数 + 繁殖個体数 × 卵数 を表記
Popul[i+1,k] <- Ncohort } # 各世代の個体数を記録（初項は init のまま）
Extinct <- if(Ncohort == 0) Extinct + 1 else Extinct }
# 絶滅回数（最終世代までに絶滅したシミュレーションの回数)
matplot(Popul, type="l") # 全シミュレーションの結果を図示
Extinct/times } # 絶滅率を出力（絶滅したシミュレーション回数/総シミュレーション回数）
# 実行例、実行する度に結果は変動します
pop.dyn(50, 0.5, 0.5, 1, 10, 50)
# 初期個体数, 生存率, 繁殖率, 卵数, 世代数, シミュレーション回数
# 注 2：コード部分、数えてみたら 11 行あるじゃないかといった突っ込みは無しの方向でお願いします。
R関連記事のポリシー & forthcoming
（注：不確定要素が確定し次第、一部の記述を削除していく予定）
書きたいことがたまっているのですが、まだ自分の中でGOサインが出せません。私など統計のド素人ですが、知り合いの研究者の方々から、検索していたら見つけたとか、役に立ったといったコメントをけっこういただいており、非常にありがたく思います。
こんなつたない、しかも量も少ない情報でも私用の忘備録だけでなく皆様の役に立つのであれば、何よりです。統計は避けて通れないけれど、本当に勝負すべきは科学的な発見や理論。統計については知恵を出し合い、とっとと済ませて主題に行きたいものです！
さて、このブログのR関連項目は、ある程度自分自身で解析を実行して一応の結論が出たものだけを載せる方針でいます。しかし、自分自身にプレッシャーをかける意味でも、次に何を書こうとしているかを人目に付くところに宣言しておいた方が良いかもしれません。
当面の目標は、ケースバイケースでP値を判定基準とした旧来の統計学から、現代的な一般化線形モデルに統合し情報量基準により最適モデルを絞り込むスタイルにすべてを移行するために、両者の対応関係をあぶり出すことにあります。
これまでに掲載したのは、
1. P値で冷徹に切り分けられてきた有為水準が、GLMのモデル選択ではどれくらいの信頼性があるかを調べた話（ただしまだ突き詰め切れていない、私自身、グレーゾーンが残る）。
2. 変数変換などで対応してきた正規分布から外れたデータは一般化線形モデルの枠組みでは分布型の指定によってより適切に対処可能であること（残る不確定要素は1とも関連、random factorをモデルに含めることで解決するであろう部分もあり）
最後に、現在準備中・試行中の項目を挙げておきます。
いずれもある程度のネタは仕入れつつも、私自身の検証が済んでいない事項です：
・多重比較（Tukey, Bonferroniなど）は、GLMの枠組みではどう対応するのか？
・manova：複数種の全体的なパターンを解析するのに用いていた手法。TypeIエラーを除くために用いられてきたような…SPSSなどではmultivariate GLMなるものが実装されているとか、Rではまだらしい？
・repeated-measure ANOVA
・pseudo R2：正規分布（family=normal）の場合に多用されてきた決定係数（R^2：アール二乗）、実測データの何％までをモデルで説明できたかを示すのに便利だが（ex. R^2値は0~1をとり、1が100%）、正規分布以外の分布型では使用することができなかった。最近、Deviance値を利用して疑似的にR2を計算するpseudo R2なるものが使用されている論文を見かける。非常に便利そうだが、検証が済んでいないので未掲載。
書きたいことがたまっているのですが、まだ自分の中でGOサインが出せません。私など統計のド素人ですが、知り合いの研究者の方々から、検索していたら見つけたとか、役に立ったといったコメントをけっこういただいており、非常にありがたく思います。
こんなつたない、しかも量も少ない情報でも私用の忘備録だけでなく皆様の役に立つのであれば、何よりです。統計は避けて通れないけれど、本当に勝負すべきは科学的な発見や理論。統計については知恵を出し合い、とっとと済ませて主題に行きたいものです！
さて、このブログのR関連項目は、ある程度自分自身で解析を実行して一応の結論が出たものだけを載せる方針でいます。しかし、自分自身にプレッシャーをかける意味でも、次に何を書こうとしているかを人目に付くところに宣言しておいた方が良いかもしれません。
当面の目標は、ケースバイケースでP値を判定基準とした旧来の統計学から、現代的な一般化線形モデルに統合し情報量基準により最適モデルを絞り込むスタイルにすべてを移行するために、両者の対応関係をあぶり出すことにあります。
これまでに掲載したのは、
1. P値で冷徹に切り分けられてきた有為水準が、GLMのモデル選択ではどれくらいの信頼性があるかを調べた話（ただしまだ突き詰め切れていない、私自身、グレーゾーンが残る）。
2. 変数変換などで対応してきた正規分布から外れたデータは一般化線形モデルの枠組みでは分布型の指定によってより適切に対処可能であること（残る不確定要素は1とも関連、random factorをモデルに含めることで解決するであろう部分もあり）
最後に、現在準備中・試行中の項目を挙げておきます。
いずれもある程度のネタは仕入れつつも、私自身の検証が済んでいない事項です：
・多重比較（Tukey, Bonferroniなど）は、GLMの枠組みではどう対応するのか？
・manova：複数種の全体的なパターンを解析するのに用いていた手法。TypeIエラーを除くために用いられてきたような…SPSSなどではmultivariate GLMなるものが実装されているとか、Rではまだらしい？
・repeated-measure ANOVA
・pseudo R2：正規分布（family=normal）の場合に多用されてきた決定係数（R^2：アール二乗）、実測データの何％までをモデルで説明できたかを示すのに便利だが（ex. R^2値は0~1をとり、1が100%）、正規分布以外の分布型では使用することができなかった。最近、Deviance値を利用して疑似的にR2を計算するpseudo R2なるものが使用されている論文を見かける。非常に便利そうだが、検証が済んでいないので未掲載。
参戦 国際甲殻類学会
先の大型連休、東京海洋大学で開催の国際甲殻類学会へと参加・発表してきました。
これまで研究材料系の学会に参加したことはなかったのですが…（ヨコエビは甲殻類学会でもマイナー過ぎて）…今回、ビッグネームも何人も来られるようで参加を決めました。
やはりフクロエビ系のシンポが私にとって一番の目玉でした。Alistair Poore や Martin Thiel と、少しですが話ができたのは嬉しかったです。先方も私の論文を読んでいて認知してくれていたり。別のシンポで、あの Mark Hay も来ていて、けど捕まえ損ねました…残念。何にせよ、論文だけで知っていた近い研究分野の人達の話をLiveで聞けたのは、やはり国際学会ならではですね！
しかし…研究材料系学会の性かもですが…マテメソとリザルトだけの発表が多すぎるのは聞いていて少々シンドイです。それとも、それで通用するのが材料限定のよさということなのだろうか…
これまで研究材料系の学会に参加したことはなかったのですが…（ヨコエビは甲殻類学会でもマイナー過ぎて）…今回、ビッグネームも何人も来られるようで参加を決めました。
やはりフクロエビ系のシンポが私にとって一番の目玉でした。Alistair Poore や Martin Thiel と、少しですが話ができたのは嬉しかったです。先方も私の論文を読んでいて認知してくれていたり。別のシンポで、あの Mark Hay も来ていて、けど捕まえ損ねました…残念。何にせよ、論文だけで知っていた近い研究分野の人達の話をLiveで聞けたのは、やはり国際学会ならではですね！
しかし…研究材料系学会の性かもですが…マテメソとリザルトだけの発表が多すぎるのは聞いていて少々シンドイです。それとも、それで通用するのが材料限定のよさということなのだろうか…
干潟調査デビュー
この業界に入って10年突破した私ですが、
先日、今更のように干潟調査デビューを果たしました…。
千葉大時代は富津干潟にはよく行っていたものの、
通過して定質の安定したアマモ場にばかり行っていたし、
アマモのギャップは干潟といえばそれまでだったかもしれないけれど…。
とにかく、よい経験になりました！干潟の歩き方というのも理解しました。
今回はイメージ通りのいかにもな泥っぽい干潟の調査（藤前です）。
入ってみての感想：ちょっと変わり種の干潟だったかも…
干潟というと遠くからでも沢山のカニやヤドカリ、トビハゼといった
目立つ輩がわんさかしているイメージだけど、ちっとも見当たらないのです？
なんだか嫌気的（酸素欠乏気味、黒い泥）？珪藻がわんさか漂っていたり、
出てくる生物も変わり種な印象です。それから、泥がやわらかい、足がよく沈む、
うっかりしてるとすぐに膝まで沈んでいく…
今回、私にとっての一番のヒットかもしれない…
オキシジミ、渋い美しさがあります
大きさからしてサギの足跡かな、それと周囲の泥が妙に黄色いと思いませんか？珪藻がシート状に積もっている様子です。こんなの初めてみました！
そうこうするうちに調査が終了して、じゃぁ帰ろう！と目の前の牡蛎礁へと足を踏み出し残りあと50cmというところでズブッと太股まで一気に入りました…まずはすごいっと感動したのだけれど…その感動が恐怖に変わるまで長くはかかりませんでしたすっぽりハマった左足を抜こうと右足で踏ん張るもミイラ取りになりかけて止めるを繰り返し…
目の前ではヤマトオサガニが小さくプクプクと音を立てながら人生（蟹生）を謳歌しつつ…
しかし無情にも我が左足は泥に引きずり込まれるばかり…足は攣り始めしかも止まらない…同行の干潟のベテラン達は笑いながら近くを通り過ぎて生物観察をエンジョイしてるし写真は撮られるし（泣…）
サヨウナラミナサン……願わくば最期はもうすこし透明度のいい海がよかったなぁ………
先日、今更のように干潟調査デビューを果たしました…。
千葉大時代は富津干潟にはよく行っていたものの、
通過して定質の安定したアマモ場にばかり行っていたし、
アマモのギャップは干潟といえばそれまでだったかもしれないけれど…。
とにかく、よい経験になりました！干潟の歩き方というのも理解しました。
今回はイメージ通りのいかにもな泥っぽい干潟の調査（藤前です）。
入ってみての感想：ちょっと変わり種の干潟だったかも…
干潟というと遠くからでも沢山のカニやヤドカリ、トビハゼといった
目立つ輩がわんさかしているイメージだけど、ちっとも見当たらないのです？
なんだか嫌気的（酸素欠乏気味、黒い泥）？珪藻がわんさか漂っていたり、
出てくる生物も変わり種な印象です。それから、泥がやわらかい、足がよく沈む、
うっかりしてるとすぐに膝まで沈んでいく…
今回、私にとっての一番のヒットかもしれない…
オキシジミ、渋い美しさがあります
大きさからしてサギの足跡かな、それと周囲の泥が妙に黄色いと思いませんか？珪藻がシート状に積もっている様子です。こんなの初めてみました！
そうこうするうちに調査が終了して、じゃぁ帰ろう！と目の前の牡蛎礁へと足を踏み出し残りあと50cmというところでズブッと太股まで一気に入りました…まずはすごいっと感動したのだけれど…その感動が恐怖に変わるまで長くはかかりませんでしたすっぽりハマった左足を抜こうと右足で踏ん張るもミイラ取りになりかけて止めるを繰り返し…
目の前ではヤマトオサガニが小さくプクプクと音を立てながら人生（蟹生）を謳歌しつつ…
しかし無情にも我が左足は泥に引きずり込まれるばかり…足は攣り始めしかも止まらない…同行の干潟のベテラン達は笑いながら近くを通り過ぎて生物観察をエンジョイしてるし写真は撮られるし（泣…）
サヨウナラミナサン……願わくば最期はもうすこし透明度のいい海がよかったなぁ………
onもoffも調査？
先日、日本列島を沸かせた？日食、全般的に悪天候ではありましたが…皆さまのところではいかがだったでしょうか？？
磯調査といえば、大潮、大潮といえば月と太陽、地球の位置関係が直線的に並ぶ時…はい、日食の22日は鹿児島に居りました（焦…）休暇を取って調査…
島に渡らないとはいえど、最大食率98%の場所、
生態学者として一番気になるのはやはり動物の行動じゃないでしょうか！？
日食の方は悪天候で危ぶまれたものの、かろうじて食率90%ぐらいまでは見れました><
見えたとはいえど、雲の隙間から時折1~2秒ずつチラッと見える程度で、
林間を飛び回る小鳥を観察する状況と同じ（解り難い例え…）くらいでした。
なので、ちゃんとした写真は撮れておりませぬ、
50%くらいまでは楽に撮れる天候だったので、カメラの損傷を恐れずに撮っておけばよかったとちょっと後悔（→これは３年後の金環日食で挽回したい！）
ピーク時の前後、動物達の動きに注意していたのですが、やはり陸上の動物が正直だったのか、クロサギがえらく高い樹の上にとまってみたり（珍しい気がするのですが）、カラス達が山の方へと飛んでいってみたり、マツムシが鳴き始めたり…。一箇所でのわずかな観測なので何も言えませぬが、日中にいきなり暗くなると動物が変な行動をとるというのはあながちウソではないかもと思ってみたり（先入観な気もするのですけれど）。
オマケ：逆日食！？ウニの殻の中央に薄曇りの向こうの太陽を合わせてみました。管足の部分の穴のパターンが見えて綺麗です！
次のページ
磯調査といえば、大潮、大潮といえば月と太陽、地球の位置関係が直線的に並ぶ時…はい、日食の22日は鹿児島に居りました（焦…）休暇を取って調査…
島に渡らないとはいえど、最大食率98%の場所、
生態学者として一番気になるのはやはり動物の行動じゃないでしょうか！？
日食の方は悪天候で危ぶまれたものの、かろうじて食率90%ぐらいまでは見れました><
見えたとはいえど、雲の隙間から時折1~2秒ずつチラッと見える程度で、
林間を飛び回る小鳥を観察する状況と同じ（解り難い例え…）くらいでした。
なので、ちゃんとした写真は撮れておりませぬ、
50%くらいまでは楽に撮れる天候だったので、カメラの損傷を恐れずに撮っておけばよかったとちょっと後悔（→これは３年後の金環日食で挽回したい！）
ピーク時の前後、動物達の動きに注意していたのですが、やはり陸上の動物が正直だったのか、クロサギがえらく高い樹の上にとまってみたり（珍しい気がするのですが）、カラス達が山の方へと飛んでいってみたり、マツムシが鳴き始めたり…。一箇所でのわずかな観測なので何も言えませぬが、日中にいきなり暗くなると動物が変な行動をとるというのはあながちウソではないかもと思ってみたり（先入観な気もするのですけれど）。
オマケ：逆日食！？ウニの殻の中央に薄曇りの向こうの太陽を合わせてみました。管足の部分の穴のパターンが見えて綺麗です！
