GIS、WinBugs、R、を使って、ベイズ法によるロジスティック回帰分析から、HSIモデルを作成する方法

営巣木のデータから、GIS、WinBugs、R、を使ってベイズ法によるロジスティック回帰分析を行い、その結果からHSIモデルを作成してみます。
手順は、以下の通りです。

1. 営巣木データの準備

  1. 営巣木の確認位置をGISで落とします。
  2. 営巣木が確認されなかった位置をダミーとして、ランダムに広域に落とします。

2. 環境データの準備

  1. 営巣木の有無を説明するための環境データを、GISで作成します。
    • 例えば、標高、植生、起伏、相対位置などを、ラスター形式で広域に作成しておきます。

3. 解析データの準備

  1. 営巣木と非営巣木の、環境データを取得します。
    1. ArcGISで、SpatialAnalystツール→抽出→サンプルで、1.と2.で作成したデータを利用します。
    2. 1.のデータにrowidの行を作成して、↑で作成したデータをテーブル結合します。
    3. ↑のデータをエクスポートして、dbfファイルをEXCELで開いてCSVに変換して保存します。

データの例 nest.csv

Y(有無) v1(500m圏森林面積) v2(起伏) v3(標高) v4(相対位置) v5(傾斜)
1 1500 100 890 0.3 33
0 1120 78 540 0.1 12
0 1350 910 230 0.6 10
1 500 100 1500 0.4 22
1 720 120 100 0.2 5
... ... ... ... ... ...

4. glmによるロジスティック回帰分析

  1. ベイズ法の前に、glmでロジスティック回帰分析をして、環境データの変数選択をします。
  2. Rを起動して、以下のように入力実行します。
    • 営巣木と環境の関係が、線形ではないと考える場合は、2次式でモデルを作成する(?)。
    • 考えられる環境変数を全部入れてみて、stepAICで変数選択
  3. stepAICの結果を考慮して、説明変数を選択します。
library(MASS)
g<-glm(Y~v1+v2+v3+I(v3^2)+v4+I(v4^2)+v5+I(v5,2),data=d,family=binomial)
stepAIC(g)

5. ベイズ法によるロジスティック回帰分析

  1. WinBugsをRから利用できるように準備します。方法は、こちら↓を参照
  2. WinBugsのためのモデルを作成します。model.bug.txt 参照:生態学のためのベイズ法 p144
  3. WinBugsをRから実行するスクリプトを作成実行します。
  4. plot(post.bugs)、post.bugs$mean$a、post.bugs$mean$b で推定された値を確認します。
model
{
	#中央化のためにデータの平均を計算(v4は、stepAICの結果、不採用)
	mv1 <- mean(v1[])
	mv2 <- mean(v2[])
	mv3 <- mean(v3[])
	mv5 <- mean(v5[])
	
	for (i in 1:N) {
         #ベルヌーイ分布
	 Y[i] ~ dbern(p[i])
         #ロジスティック関数。v3、v5は、二次式の項。すべての項を中央化してある。
	 logit(p[i]) <-  a + b[1]*(v1[i]-mv1)+b[2]*(v2[i]-mv2)+b[3]*(v3[i]-mv3)+b[4]*(v5[i]-mv5)+b[6]*(v3[i]-mv3)*(v3[i]-mv3)+b[7]*(v5[i]-mv5)*(v5[i]-mv5)
	}

 	a ~ dnorm(0, 1.0E-6)

	for(i in 1:7){
	 b[i] ~ dnorm(0,1.0E-6)
	}
}
source("R2WBwrapper.R")
d<-read.csv("nest.csv")

clear.data.param()
set.data("N", nrow(d))#標本数
set.data("Y", d$Y) #巣のあるなし 0:無 1:有
set.data("v1", d$v1)
set.data("v2", d$v2)
set.data("v3", d$v3)
#set.data("v4", d$v4) v4は不採用
set.data("v5", d$v5)

set.param("a", 0)
set.param("b",rep(0,7))
set.param("p", NA)

post.bugs <- call.bugs(
	file = "model.bug.txt",
	n.iter = 2000, n.burnin = 1000, n.thin = 5
)

plot(post.bugs)
post.bugs$mean$a
post.bugs$mean$b

6. ロジスティック回帰分析の結果からSIを作成

  1. SI値は、ロジスティック回帰モデルで、他の条件を最大値として固定したものとします。
m<- post.bugs$mean #サンプリングされた値の平均値
invlogit<-function(x) 1/(1+exp(-x)) #ロジスティック関数

#標高 0〜1500mの範囲で、営巣確率が最も大きくなるv3の値を取得。v3以外は、平均値をとることとして計算。
x<-seq(0,1500,(1500-0)/100) 
y<-invlogit(m$a + m$b[3]*(x-mean(d$v3))+m$b[6]*(x-mean(d$v3))*(x-mean(d$v3)))
plot(y~x)
max_v3<-x[y==max(y)]

#↑と同様に、max_v1〜max_v5までを取得

#v3のSIモデルは、v3以外の値が、max_v1〜max_v5をとるとしてロジスティック回帰分析したもの
x<-seq(0,1500,(1500-0)/100)
y<-invlogit(m$a + m$b[1]*(max_v1-mean(d$v1))+m$b[2]*(max_v2-mean(d$v2))+m$b[3]*(x-mean(d$v3))+m$b[6]*(x-mean(d$v3))*(x-mean(d$v3))+m$b[4]*(max_v5-mean(d$v5))+m$b[7]*(max_v5-mean(d$v5))*(max_v5-mean(d$v5)))
plot(y~x)

#↑と同様にして、v1〜v5までのSI値のモデルを作成する。

7. ロジスティック回帰分析の結果からHSIモデルをGISで図示

  1. HSIモデルは、ロジスティック回帰分析の結果を、そのまま利用します。
  2. 推定したい地域のHSI値は、ArcGISのラスタ演算で、GISで作成した環境変数(v1〜v5)、ベイズ法で推定した重みの平均値 post.bugs$mean$a,b、中央化のために算出した平均値を、ロジスティック関数に代入して、算出します。
[q] = -1.638 + 0.0018585 * ( [v1] - 409.6476 ) - 0.0000173 * ( [v2] - 409.6476 ) - 10.4250000 * ( [v3] - 0.3531808 ) + 0.7933500 * ( [v3] - 0.3531808 ) * ( [v3] - 0.3531808 ) + 0.0705900 * ( [v5] - 22.90048 ) - 0.0010540 * ( [v5] - 22.90048 )
HSI = 1 / (1 + Exp(-1.0 * [q] ))


自分用メモ

  • 営巣しないデータのサンプリング方法は、これでいいのか？ランダムサンプリングしたものにも、実際には、営巣している場所があるかもしれない。
  • 環境変数に相関関係があっても、stepAICで除去される？もしくは、事前に、相関係数を調べておく必要がある？
  • この例は、別にベイズ法じゃなくてもいい。glmでいいと思う。とりあえず、ベイズ法でやってみたかった。
  • carモデルを利用すれば、ベイズ法の意味があるかもしれないが、外挿できないので、HSIモデルとして広域に展開できない。
  • glmでpoly(v3,2)と指定したら推定結果がおかしい。v3+I(v3^2)としないといけない？

参考文献

↓これを読めば、すっきりベイズ法が理解できます。最初に読むと良いと思います。