密度指標トレンドに基づく外来生物防除の評価
2011/10/21 深澤圭太
0. はじめに
外来生物の防除に伴う密度指標[1]の時間変化から、防除努力量と除去率の関係や自然増加率を推定し、さらにその推定値を用いて根絶の将来予測を行うRスクリプトを作成した。
ダウンロード
プログラム eradestim1.r
サンプルデータ example1.csv
1. モデルの構造と仮定
1.1. 1つの個体群を対象とする場合
防除の影響を受けている外来生物の個体群動態は、自然増加と除去のバランスによって決まると考えられる。たとえば、防除努力量[2]が大きければ外来生物の密度は減少するが、不十分である場合には除去による減少よりも自然増加が上回り、外来種の増加をくい止めることができないだろう(図１)。
図１ 防除努力量と個体密度の関係
このようなプロセスは、次のような簡単な個体群動態モデルで表すことができる。
相対密度t = 相対密度t-1×残存率t-1×自然増加率 (1)
なお、「相対密度t」のような添え字tは「t年目における値」を意味する。残存率t-1×自然増加率 < 1のときに個体密度は前年度から減少する(図1)。残存率は防除努力量によって変化するため、ここでは以下のモデルを仮定する。
残存率t = exp(- 防除係数×防除努力量t) (2)
この式は、防除努力量が0のときは残存率が1で、防除努力量が増えるほど残存率が0に近づくという関係に相当する。防除係数(> 0)が大きいほど、1防除努力量が個体群に与えるインパクトが大きくなる。なお、防除係数が十分小さいときに限り、この値は近似的に「1防除努力量を投入したときの除去率(= 1 – 残存率)」としてみなすことができる[3]。式(1)と(2)を１つにまとめると、次のような式が得られる。
相対密度t = 相対密度t-1×exp(- 防除係数×防除努力量t-1)×自然増加率
防除係数、自然増加率、そして初期状態における相対密度が決まれば、防除努力量に対する毎年の相対密度が予測できる。
密度指標は相対密度に比例する値である。ここでは単純化のためそれらの比例係数は1であるとする。密度指標 = 発見数 / 調査努力量である事から、下記の関係が成り立つ。
密度指標の平均t = 平均発見数t / 調査努力量t = 相対密度t
⇔ 平均発見数t = 相対密度t×調査努力量t
発見数はカウントデータであるため、負の2項分布にしたがうことを仮定する。
上記をまとめると、下記の３つの式からなる構造方程式が得られる。
相対密度t = 相対密度t-1×exp(- 防除係数×防除努力量t-1)×自然増加率 (3a)
平均発見数t = 相対密度t×調査努力量t (3b)
発見数t ~ Negative-Binomial(平均発見数t, 分散パラメータ) (3c)
自然増加率、防除係数、分散パラメータ、初期相対密度のパラメータ推定は最尤法により行う事ができる。
なお、モデルには下記の仮定が含まれる事に留意すべきである。
・ 自然増加率は時間によらず一定、密度効果は無視できるほど小さい
・ 防除効率は密度によらず一定
もし推定したい対象に上記を仮定できない場合、自作コードによる最尤推定か、一般化状態空間モデルなどのより洗練された手法を使うことをおすすめしたい。
1.2. 複数の個体群を対象とする場合
異なる防除圧を受ける複数のサブ個体群の時間変化が得られている場合、式(3a)～(3c)にサブ個体群のインデックスiを付加することで同じようにパラメータ推定が可能である。ただし、サブ個体群間の移出入は無視できると仮定する。
相対密度it = 相対密度it-1×exp(- 防除係数×防除努力量it-1)×自然増加率 (4a)
平均発見数it = 相対密度it×調査努力量it (4b)
発見数it~Negative-Binomial(平均発見数it, 分散パラメータ) (4c)
初期状態の相対密度はサブ個体群の個数存在する事になり、防除係数、自然増加率、分散パラメータとともにすべて最尤推定する事になる。
2. 使い方
2.1. データの準備
“example1.csv”を例に解説する。このデータは、環境省「奄美大島におけるジャワマングース防除事業」における2001-2010年度までの防除記録である[4](複数のわなをプールしたデータ)。
|
occasion |
site |
nobs |
surveyeffort |
removeeffort |
|
1 |
1 |
3375 |
165367 |
165367 |
|
2 |
1 |
2191 |
147353 |
147353 |
|
3 |
1 |
2565 |
221403 |
221403 |
|
4 |
1 |
2524 |
318359 |
318359 |
|
5 |
1 |
2591 |
630822 |
630822 |
|
6 |
1 |
2713 |
1051026 |
1051026 |
|
7 |
1 |
783 |
1380751 |
1380751 |
|
8 |
1 |
945 |
1899238 |
1899238 |
|
9 |
1 |
598 |
2174339 |
2174339 |
|
10 |
1 |
311 |
2101116 |
2101116 |
解析に際しては、下記のような5列のデータセットをcsvファイルで作成する必要がある。
各列の意味は下記のとおりである。
occasion･･･時間ステップ(年、月など。最初のステップは1とする)
site･･･サブ個体群のインデックス
複数のサブ個体群がある場合、1から順に番号をつける。
nobs･･･モニタリングにおける検出個体数
surveyeffort･･･モニタリング努力量
removeeffort･･･除去努力量
上記の例ではわなによる防除を行っており、CPUEを密度指標としているため、モニタリング努力量 = 除去努力量となっている。
2.2. Rのインストール
CRAN(http://cran.r-project.org/)のミラーサイトからRをダウンロードし、インストールする。Rの基本的な使い方については教科書やWebサイトに情報が集積されているためそちらを参照されたい。
2.3. ソースファイルの読み込み
Rを起動したら、下記のコードにより推定に必要な関数群を含むRスクリプト”eradestim1.r”を読み込む。
source(“パス名/eradestim1.r”)
なお、同様の処理は、メニューバー[ファイル]→[スクリプトを開く…]によっても可能である。ソースファイルをテキストエディタなどで開くと、それぞれの関数の引数や返り値の意味を参照する事ができる。
2.4. パラメータ推定
下記のRコードにより、推定計算、要約統計量の計算、除去努力量と残存率の図化が可能である。以下のコードはRにコピーペーストすることで実行する事ができる。
#データの読み込み example<-read.csv("パス名/example1.csv") #パラメータ推定 res1<-estim.pd(data=example,noc=10,nsite=1) #nocは時間ステップ数、nsiteはサブ個体群数 #返り値の詳細については”eradestim1.r”のコメント行を参照。 #(.rファイルはテキストエディタで開ける) #要約統計量の出力 summary.pd(res1) #捕獲努力量と捕獲率の関係 removalrate.pd(res1)
※ 現実的な自然増加率が推定できない場合
推定の結果、自然増加率が現実的でない値をとる事がある。これは、捕獲努力量のばらつきが小さく、自然増加の効果と捕獲の効果をうまく分離できない場合によく起こる。マングースの例においても推定された自然増加率は1未満であり、このことは自然に個体密度が減少してしまうことに相当し、当然ながらこれは非現実的な結果である。そのような場合は、自然増加率を事前の知識から決定し、それ以外のパラメータをデータから推定するのが次善の策となるだろう。ただし、自然増加率の値に対して防除係数は敏感であることがあるため、現実的な範囲で自然増加率の値を何通りか試して解の反応を試す必要がある。
#パラメータ推定 res2<-estim.pd(data=example,noc=10,nsite=1, fixr=1.43) #fixrは事前に自然増加率を決める際に指定する引数。ここでは1.43とした
自然増加率をデータから推定するには、場所または時間ごとの捕獲努力量のばらつきを大きくとる必要がある。最初の一定期間は防除せずにモニタリングのみを行うことや防除圧をかけない対照区を設定することは、パラメータ推定上は有効である。
2.5. 将来予測
防除の将来予測のための関数として、①個体数に関するモンテカルロシミュレーション、②初期の相対密度を1としたときの相対密度の時間変化の計算、の2つを用意した。前者は希少種の絶滅リスク評価にも用いられる存続可能性分析(PVA)であり、パラメータセットおよび年ごとの翌年の平均個体数は
個体数t = 個体数t-1×exp(- 防除係数×防除努力量t-1)×自然増加率
により計算され、実際の個体数はポアソン分布に従ってばらつくと仮定している。後者は平均的な相対密度の時間変化を計算するのみであり、年ごとのばらつきは発生しない。コードの例を下記に示す。
##10年間の予測 #10年後までの防除努力量 effort1<-rep(2174339,10) #10年間の個体数シミュレーション #(100反復、初期個体数を300-1500の一様乱数で決定) npred1<-pred.nind(res2$fixr,exp(res2$result$par[1]),effort1,10, nind.init=round(runif(100,300,1500))) #10年間の相対密度計算 rdpred1<-pred.rdens(res2$fixr,exp(res2$result$par[1]),effort1,10)
個体数の情報が得られる場合はモンテカルロシミュレーションを用いることにより期待根絶年数が計算できる。そうでない場合は危険率を見込んで相対密度の閾値を設定し、関数pred.rdensにより予測される相対密度がそれを下回った時点を便宜的に根絶達成とみなすのが現実的と思われる。