« マルコフ連鎖モンテカルロ法 | トップページ | カナヘビ »

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

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

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

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をたとえば1から50まで動かしたとすると、マウス操作とかを50回繰り返さなきゃいけない(もしくは、ひたすら長い行をドラッグ)わけで、そんなことをしていたらどこかでミスるか、呪われると思う。

じゃあ、ちょっと遊びでRを使ってやってみるかと思い立ち、仮想データを作ってやってみることにした。

長いので、続きに入れます。

ストーリー

とある稀少なウサギさんを保全する条件を検討するため、いろいろなデータ(n=80)をとってモデル選択したら、”生息地面積”、”餌の量”、”キツネの個体数”が残ったとする。

set.seed(10)
number<- c(rpois(20,9),rpois(20,20),rpois(20,15),rpois(20,27))
area<- c(rep(10,20),rep(25,20),rep(18,20),rep(32,20))
food<- c(rnorm(20,14), rnorm(20,33),rnorm(20,25),rnorm(20,45))
fox<-  c(rnorm(20,6),rnorm(20,4),rnorm(20,4),rnorm(20,3))

numberがウサギさんの個体数、areaが生息地面積、foodが餌の量、foxがキツネの個体数とします。(発生データは事実的な根拠のない適当乱数)

個体数データなので、ポアソン回帰で求める

mod<- glm(number~area+food+fox,family=poisson)

ここでstep(mod)とやるとareaしか変数が残らないのですが、ここでは無視します。(いいデータ例を用意できなくてごめんなさい)

で、得らた係数は、

(coe<- coef(mod))
(Intercept)        area        food         fox
1.93672138  0.07372001 -0.01937107 -0.03605800

となる。

つまり予測値は(係数を四捨五入)、

個体数=exp( 0.07*area -0.02*food - 0.04*fox + 1.94)

となる。ポアソン回帰のデフォルトのリンク関数はlogなので逆関数のexponentialをかけてやる。

これを手動でもとめて、predict()と比較してみる。

pre<- exp(coe[1] + coe[2]*area + coe[3]*food + coe[4]*fox)

unique(round(pre-predict(mod,type="response")))
[1] 0

で、面積が変わったらウサギさんの個体数がどうなるかを予測する。

まずは、面積が10のときを計算してみると、

pre2<- exp(coe[1] + coe[2]*10 + coe[3]*food + coe[4]*fox)
mean(pre2)
[1] 7.192213
quantile(pre2,c(0.025,0.0975))  #95%信頼区間
    2.5%    9.75%
5.217238 5.487587

こんな要領で、1:30だったらどうなるかやってみる。

length(10:30)
means<- numeric(21)
qtl<- numeric(21)
qtu<- numeric(21)
for(i in 1:21){
pred<- exp(coe[1] + coe[2]*(i+9) + coe[3]*food + coe[4]*fox)
means[i]<- mean(pred)
qtl[i]<- quantile(pred,0.025)
qtu[i]<- quantile(pred,0.975)
}

#予測値と95%信頼区間
plot(10:30, means,ylim=c(min(qtl),max(qtu)),ylab="予想個体数",xlab="面積")lines(10:30, means)
lines(10:30,qtl,lty=2)
lines(10:30,qtu,lty=2)

Reg 

 

きれいに予測値と信頼区間の線が意外と簡単に引けました。

ちなみに、データを全部プロットするとこうなる。

mat<- matrix(rep(0,80*21),nrow=21)
means<- numeric(21)
qtl<- numeric(21)
qtu<- numeric(21)
for(i in 1:21){
pred<- exp(coe[1] + coe[2]*(i+9) + coe[3]*food + coe[4]*fox)
mat[i,]<- pred
means[i]<- mean(pred)
qtl[i]<- quantile(pred,0.025)
qtu[i]<- quantile(pred,0.975)
}
plot(10:30, means,ylim=c(min(mat),max(mat)),ann=F,col="red",type="n")
for(i in 1:21){points(rep(i+9,80),mat[i,])}
lines(10:30, means,col="red")
lines(10:30,qtl,lty=2)
lines(10:30,qtu,lty=2)
points(10:30,means,col="red",cex=2,pch=15)

Reg2 

 

 

 

 

 

 

 

 

 

 

こんなこんな感じで、80×21=1680点生成されるわけで、エクセルでやるとどこかで確実にミスるか呪われそうですね。

 

 

本当は、

y= 0.2A + 0.1B + 2

になるような乱数を発生させてもうちょっとスマートなデータでやりたかったけれど、そのやり方がわかりませんでした。

|

« マルコフ連鎖モンテカルロ法 | トップページ | カナヘビ »

統計/R」カテゴリの記事

コメント

この記事へのコメントは終了しました。

トラックバック


この記事へのトラックバック一覧です: モデルからインパクトを予測する:

« マルコフ連鎖モンテカルロ法 | トップページ | カナヘビ »