2011年09月18日

標識再捕法のお勉強と気候変動によるニッチシフト

標識再捕法についてちょっとわからんことが出てきたので、簡単なシミュレーションでまず理論を確認。

#################
###標識再捕のお勉強###
#################
> set.seed(1)
> Estimated_N <- Estimated_Pab <- NULL #計算に使うオブジェクトの用意
>
> for(i in 1 : 5000){
+ N <- 100 #真の個体数
+ Pa <- 0.8 #観察者aの真の発見率
+ Pb <- 0.5 #観察者bの真の発見率
+
+ ID <- seq(1, N) #個体番号
+ Na <- sample(ID, N*Pa) #aの発見した個体番号
+ Nb <- sample(ID, N*Pb) #bの発見した個体番号
+
+ Nab <- Na[Na%in%Nb] #aとbで共通の発見した個体番号
+ Estimated_N[i] <- (length(Na)*length(Nb))/length(Nab) #推定個体数
+ Estimated_Pab[i] <- length(Nab)/Estimated_N[i] #aとbが同時に発見する確率(推定値)
+ }
>
> par(mfrow=c(2, 1))
> hist(Estimated_N, breaks=13) #推定個体数の分布
> hist(Estimated_Pab, breaks=13) #推定同時発見確率の分布
> mean(Estimated_N) #ほとんど100
[1] 100.2966
> mean(Estimated_Pab) #ほとんど0.4(Pab=Pa*Pb=0.4だから当然)
[1] 0.4006456
Estiamted_N_Pab.jpeg
うむ、とりあえず確認出来た。ここからもう少し掘り下げて考えねば。

それと今日の論文。ある外来生物は、侵入前にいたニッチと、侵入後のニッチを変えるのか、というお話。

Evidence of climatic niche shift during biological invasion. Broennimann O, Treier UA, Müller-Schärer H, Thuiller W, Peterson AT, Guisan A.Ecology Letters. 2007 10(8):701-9.

侵入先と侵入前でニッチを変えないという仮説(niche conservatism)という考えと、変えるという仮説(niche shift)の両方がある。従来、niche conservatismに従って、外来生物が侵入先ではどのような分布になるのか、という予測がなされてきた。一方で、それらの研究は、niche shiftが明示的な対立仮説として強調されていなかった。両方とも、どのように種形成をしてきたのか、という議論に重要な示唆を与えてくれる。例えば、niche conservatismは異なる場所での種形成をうまく説明し、niche shiftは同所的にいる場合の種形成の議論に重要である。また、niche conservatismは過去の気候変動に応える形で、種の移住を引き起こしたかもしれないが、そのような気候変動は、同時にニッチの区分化や進化を引き起こした可能性がある(例えば空いたニッチを埋める形で、ニッチの占める場所を広げるなどの、つまりniche shift)。

このニッチの区分化のカギとなるのが、基本ニッチと実現ニッチ。基本ニッチは、占めることのできる可能性のあるニッチで、実現ニッチは実際に生物が占めているニッチ。前者が他種との関係を考えていないのに対し、後者は考えている(と僕は理解している)。niche shiftが起こる原因としては、例えば、侵入することが競争者や捕食者からの解放を意味する場合は、実現ニッチが変わるために、それがniche shiftを引き起こすかもしれない。一方で、侵入先という新しい環境に適応すれば、それで基本ニッチが変化するかもしれない。あるいは両方とも起こるのかもしれない。

そこで、ある侵入種を使って、それがniche conservatismとniche shiftのどちらを支持しそうか、研究している。ここで著者は、過去の研究でniche conservatismが起こっているかどうかテストしたほうが良いと提案されている地域でのデータを使って、解析を行っている(このへんがモチベーションか)。結果としては、niche conservatismではなく、niche shiftが起こっていた。面白いのがFig 2で、侵入前のデータに従って予測をした結果を使って侵入後の分布を予測したときと、侵入後のデータを使って侵入前にどういう分布をしていたのかという予測をした結果がまるで違う。niche conservatismが起こっていれば、どっちからどっちを予測しても結果は同じになるはずだが、全くそうなっていない。このように、著者らはniche conservatismに従って予測をすると、完全には予測できないと結論付けている。

まぁ当たり前の結論と言えば結論だが、すべての種がniche conservatismにならないとも限らない。両方起こっていて、他種や気候との関係で、どちらが尤もらしそうに見えるのかは変化するだろうなぁ、と考えた。でも結論まではしっかりした手法で導かれていて、見ていて安心できた論文でした。
posted by しばきん at 12:31| 神奈川 ☀| Comment(0) | TrackBack(0) | R | このブログの読者になる | 更新情報をチェックする

2011年05月24日

evenかoddか

さっきボスと打ち合わせして、シミュレーションがやり直しになったので久々にRに触った。高速化のため、rnorm()で別々の分散を持つ乱数を一気に発生出来ないかなぁ、と思ったら出来たので
それの確認のためにもメモ。偶数と奇数を分ける方法を例に。
> a <- rnorm(100000, 0, c(1, 100))
> b <- 1:100000
> sd(a[which(b%%2==1)])
[1] 0.9966524
> sd(a[which(b%%2==0)])
[1] 99.66625
#ベクトルの長さの確認。丁度半分になっている。
> length(a[which(b%%2==1)])
[1] 50000
> length(a[which(b%%2==0)])
[1] 50000
ラベル:R
posted by しばきん at 14:05| 神奈川 ☔| Comment(0) | TrackBack(0) | R | このブログの読者になる | 更新情報をチェックする

2011年02月02日

Sumの使い方で怖かったこと

さっき計算してて、一瞬冷や汗かいたのでメモ。
ベクトルの長さが同じならカンマで区切っても+でも一緒の結果。

> sum(1:3, 1:3)
[1] 12
> sum(1:3+1:3)
[1] 12

で、さっき長さの違うものを足したらこうなった。

> sum(1:3+2)
[1] 12
> sum(1:3, 2)
[1] 8

そりゃそうだ。上は3 + 4 + 5をやっているのに対し、下は1 + 2 + 3 + 2をやってるわけだから。気をつけよう。
posted by しばきん at 20:11| 神奈川 ☁| Comment(0) | TrackBack(0) | R | このブログの読者になる | 更新情報をチェックする

2010年11月10日

作図時に任意の数字や文字に使用する軸を変更する方法

後輩からよくわからないという連絡を受けたので、自分でも再確認のためにメモ。
確かに検索しても見つけにくいような気もする。
plot(seq(1, 10), xaxt = "n", yaxt="n")
axis(1, at=1:10, LETTERS[1:10])
axis(2, at=1:10, seq(100, 1000, length=10))

出来た図
change_axis.pdf

ラベル:R
posted by しばきん at 22:16| 神奈川 ☀| Comment(0) | TrackBack(0) | R | このブログの読者になる | 更新情報をチェックする

2010年10月12日

Rで並列計算


#################
###Rで並列計算###
#################

###並列計算用に関数作成###
start_simulation <- function() {

###データの生成###
a <- rnorm(400)
b <- rnorm(400)

###回帰###
c <- lm(a ~ b)

###とりあえず、係数でも記録します###

return(c(a = coef(c)[1],
b = coef(c)[2])
)
}

###繰り返し数###
k <- 10

library(snow)
###乱数の都合でこいつも使います。ないと面白いことになります。###
library(rlecuyer)


#Core 2 Duoなので2つです。クアッドコアとかなら4ですかね。開きすぎると閉じさせられます。
#Firewallは無視です。ブロック解除です。
cl <- makeCluster(2, type = "SOCK")

###乱数種はここでセット###
clusterSetupRNG (cl,seed=123)

###並列計算するfunctionをここで明示###
clusterExport(cl, "start_simulation")

###スタート。ついでに時間も計っちゃいましょう###
###係数はpに、時間はtimesに入るようになっています###
times <- system.time(p <- parSapply(cl, 1:k, function(x){start_simulation()}))

###applyするときは転置してからであります###
apply(t(p), 2, mean)
apply(t(p), 2, sd)

#################
###Rで並列計算###
#################


しかし、あんまりにも複雑で長い時間がかかるものを並列計算させると、終わってるはずの計算がいつまで経っても終わらないこともある(経験談)ので、少し注意が必要みたいです。

posted by しばきん at 17:20| 神奈川 ☁| Comment(0) | TrackBack(0) | R | このブログの読者になる | 更新情報をチェックする
×

この広告は1年以上新しい記事の投稿がないブログに表示されております。