Silviase/jgov_v1.4_full_ann huggingface dataset:
(4) 個体群モデルによるモニタリング手法及び生息数推定法の確立に関する研究

横浜国立大学 環境情報研究院	松田 裕之
独立行政法人 森林総合研究所	堀野 眞一
<研究協力者>
独立行政法人 国立環境研究所	深澤 圭太
横浜国立大学 環境情報学府	太田 海香
財団法人 自然環境研究センター	諸澤 崇裕

平成21～23年度累計予算額：1,511千円（うち、平成23年度予算額：514千円）

予算額は、間接経費を含む。

[要旨] 標識再捕獲法における捕獲の空間情報は対象個体群の移動距離に関する情報を持つ。このため、捕獲の空間情報を活用することで調査地・調査期間に特有な有効トラップ面積を明らかにし、客観的にクマ類の生息密度を推定することができると考えられている。しかしながら、これまでに標識再捕獲法における空間情報から生息密度を推定する手法は複数考案されており、手法によって推定密度がどの程度変化するかについては明らかでない。本研究においては、ダミーデータと本プロジェクトで得られた、岩手県北上山地におけるツキノワグマのDNA標識再捕獲データを用いた解析により、平均最大移動距離法、空間明示最尤法、ベイズ空間明示標識再捕獲モデルの3手法を用いて生息密度推定手法の精度比較を行った。ダミーデータを用いて推定精度の比較を行ったところ、ベイズ空間明示標識再捕獲モデルは設定した密度と近い数値を示し、推定精度が高いと考えられた。また、空間明示最尤法も同様に比較的高い推定精度を示した。一方で、平均最大移動距離法を用いた推定では、設定密度よりも低い値を示し、推定精度は2つの空間明示モデルと比べると劣る可能性があることが示唆された。2010年と2011年に北上山地で行ったヘア・トラップ調査から得られたデータを用いて生息密度推定を行った結果、ベイズ空間明示標識再捕獲モデルと空間明示最尤推定法の生息密度推定値は、平均最大移動距離法に比べると両年ともダミーデータの結果と同様に手法間の推定結果の差は小さかった。以上の結果、ベイズ空間明示標識再捕獲モデルと空間明示最尤法は、クマ類の実用的な生息密度推定法であることが示された。

[キーワード] ベイズ空間明示標識再捕モデル、空間明示最尤法、標識再捕獲法、行動圈、生息密度推定

### 1. はじめに

マイクロサテライトなどDNAによる個体識別の技術の向上とコストダウンにより、ヘア・トラップによる生息密度推定はクマ類の保護管理の現場においても一般的なモニタリングツールとなりつつある。ヘア・トラップによる生息密度推定は、現地調査、DNA分析、統計解析という3つの段階から成る。日本においてもこれまでに北海道や長野県を中心に20の自治体で取り組まれて
きている $ ^{1)} $ 。それにもかかわらず、いずれの段階においても方法論にはさまざまなバリエーションがあり、最適な選択肢に関するコンセンサスは得られていない。

クマ類の生息数推定においては、これまでLincoln-Petersen法などの閉鎖個体群モデルが一般的に用いられてきた。しかし、この手法で推定された生息数はヘア・トラップからどの程度の範囲に生息する個体を対象としたものであるか自明ではなく、生息密度の算出のためには別途有効トラップ面積を計算する必要がある。従来の日本におけるヘア・トラップ調査においては、有効トラップ面積としてトラップが含まれるメッシュの面積やトラップから発生させたバッファの面積が用いられてきた。しかし、これらの方法はメッシュやバッファのサイズの決め方に科学的な根拠がなく、評価者の主観による影響を強く受ける点で恋意的になるという問題点がある。

標識再捕獲法における有効トラップ面積は、対象生物の行動圏サイズとの関連が古くから議論されてきた $ ^{2} $ 。クマ類の行動圏サイズは堅果類の豊児など、時空間的に変化する環境条件によって影響をうけ、それによって有効トラップ面積も大きく変動すると考えられる。標識再捕獲法における捕獲の空間情報は対象個体群の移動距離に関する情報を持つため、それを活用することで調査地・調査期間に特有な有効トラップ面積を明らかにし、客観的に生息密度を推定することができると考えられる。しかしながら、これまでに標識再捕獲法における空間情報から生息密度を推定する手法は複数考案されており、手法によって推定密度がどの程度変化するかについては明らかでない。このため、実際にどのようなモデルを用いるべきか判断するには、モデルの推定精度を比較する必要がある。

生息密度推定の精度においては、現地調査における試料採取数や再捕獲数の確保も課題として挙げられる $ ^{3)} $ 。例えば、現地調査における再捕獲数が多くなるほど推定精度が高くなることが報告されている $ ^{4)} $ 。しかしながら、現地調査において再捕獲数をコントロールすることは不可能である。したがって、十分な再捕獲数を確保するためには、トラップ数やセッション数を増やすといった努力が必要になる。このようなトラップ設置条件については、様々なトラップ条件のもとで生息密度推定を行い推定結果の変動を調べることによって、密度推定の結果を現地調査にフィードバックできると考えられる。また、ヘア・トラップの設置条件はクマ類の生息密度や環境の時間的、空間的な違いにも依存すると考えられる。様々な地域のデータを用いて生息密度推定を行うことで、より一般性の高いトラップの設置条件が抽出できることが期待できる。

## 2 ．研究開発目的

本研究では、(i) 空間明示型標識再捕獲モデルによる生息密度推定の精度検証および、(ii) より効率的なトラップ設置方法（トラップ配置、トラップ数、セッション数）の検討を目的として以下の検証を行った。まず、空間非明示モデルと空間明示モデルを用いて、均一トラップ配置と不均一トラップ配置のもとでの密度推定をダミーデータで行い、推定精度ならびにトラップ配置について検証した。次に、密度既知のダミーデータにより複数の空間明示モデルによる推定値と真の値のずれを検証した。同様の傾向が野外データに対する推定においても見られるかどうかを確かめることで、それぞれの手法の推定バイアスを明らかにすることを試みた。さらに、ダミーデータと岩手県北上山地における大規模調査データを用いた解析により、生息密度推定の精度比較を行った。生息密度推定手法の対象とした手法は下記の3つである。

## (i) 平均最大移動距離に基づく方法 $ ^{5)} $
(ii) 空間明示最尤法 (Program “DENSITY”) $ ^{6)} $

(iii) ベイズ空間明示標識再捕獲モデル $ ^{7,8)} $

2010年度の大規模調査のデータを用いて、トラップ数、セッション数と密度推定結果のばらつきの関係を検証した。加えて、過去にヘア・トラップ調査が行われた岩手県以外のデータも用いて他の地域への推定方法の適用可能性を検証するため、生息密度推定の成功事例あるいは失敗事例から、再捕獲数やトラップ数、セッション数が密度推定の成否とどのような関係にあるかを分析した。

## 3 ．研究開発方法

## （１）生息密度推定手法の精度検証

1）空間非明示モデルと空間明示モデルの比較およびヘア・トラップ配置の検討

本サブテーマにおいては、北上山地モデル調査地における大規模ヘア・トラップ試験における2010年のトラップ設置地域（図(4)-1a）、および同じ範囲に均等に配置した同数のトラップ（図(4)-1b）という2通りの空間配置のもとでダミーデータを生成した。トラップ数は262基、試験地の面積は674.74km²であり、セッション数は6とした。この地域におけるツキノワグマの生息密度は知られていないため、過去のヘア・トラップ調査 $ ^{9)} $ に基づく北上山地の推定生息密度0.195個体/km²を真の値と仮定した。ダミーデータから推定された生息密度がこの値に近いほど望ましい結果であるということになる。

図(4)-1 (a)北上山地モデル調査地における大規模へア・トラップ試験のトラップ設置地域、および(b)均等に配置した仮想的なトラップ位置図。黒い点がトラップの位置、灰色は調査地の範囲(674.74km²)を示す。

大規模ヘア・トラップ試験のトラップ設置地域と均一トラップ配置の下で生成された2種類のダミーデータを用いて、Gardner et al.(2009) $ ^{7)} $ の空間明示モデル（詳細は後述）と従来の空間非明示捕獲率均一モデルによる生息密度の推定精度を比較した。空間非明示モデルはOtis et al.(1978) $ ^{10)} $ の $ M_{0} $ モデルを使用した。このモデルは初回の捕獲と再捕獲が同じプロセスによって行われた標識
再捕獲データの解析に適した方法であり、ヘア・トラップによる標識再捕獲データはそれに合致していると考えられる。

空間明示モデルを用いると任意の長方形の範囲における生息数が算出されるため、生息密度は自然に導出される。一方、 $ M_{0} $ をはじめとする空間非明示モデルでは、推定された個体数がどの範囲に生息する個体群を対象としたものなのかは自明ではない。これまでの生息密度の母数には、トラップ設置位置の凸包の面積やトラップ設置位置から経験的に決めた距離のバッファに含まれる面積など、さまざまなバリエーションがある。

ここでは、便宜的に調査地面積674.74 km²を個体密度計算の際の分母とした。

解析に用いるダミーデータの生成においては、Gardner et al. (2009)の空間明示型標識再捕獲モデルを用いた。このモデルはクマの各個体が円形の行動圏を持つことを仮定しており、生成されるダミーデータは毛が採取されるトラップが行動圏の中心の周りで空間的に凝集するパターンをもつ。採取率と行動圏中心・トラップ間距離の関係を決めるパラメータσ、トラップ接触強度が+∞のときの採取率$p_{0}$、各個体の行動圏の中心座標$\mathbf{s}_{i}=(s_{1i}, s_{2i})$を与えれば、乱数によりダミーデータを生成することができる。

σ は Gardner et al.(2009)がアメリカ・ニューヨーク州のアメリカクロクマのヘア・トラップデータから推定した2.77とした（図(4)-2）。距離による相対的な採取率の変化はクマの利用頻度に相当する。テレメトリー調査による既存の日本におけるツキノワグマの最外殻行動圏の報告は、円形行動圏を仮定した半径に換算すると1.27～6.26 km^{11,12}（季節移動が明らかな例は除いた）の範囲であり、図(4)-2のパターンはその範囲に含まれる。

また、$p_{0}=0.165$ とした（図(4)-3）。これは、(サンプルが採取されたトラップ・セッション・個体数 / 全トラップ・セッション数)がおよそ0.1になるように調整した値である。実際にはこの値は数

行動圏の中心とトラップ間の距離(km)

図(4)-2 相対的な採取率と行動圏の中心からトラップまでの距離の関係。Gardner et al. (2009)によるアメリカクロクマのデータからの推定値(σ=2.77)に基づく。

図(4)-3 ダミーデータ生成の際に用いた、いずれかのトラップで検出される確率とトラップ接触の強度の関係(p0=0.165)。

倍大きい場合が多い（例えば、Miura & Oka (2003) $ ^{13} $ ）の北上山地における調査では0.289）。採取率が低い条件下においても正しく個体密度が把握できるならばより採取率が高い状況においてもよい推定が得られると考えられるため、ここではあえて悲観的なシナリオを想定した。各個体の行動圏の中心座標は一様乱数により決定した。各個体の行動圏の中心座標がとりうる範囲は、ヘア・トラップ位置座標の上限と下限を囲む長方形の上下左右10 kmの内部とした。
## 2 ）ダミーデータを用いた生息密度推定精度の検証

一般に、生物の行動は2次元ランダムウォークとして記述されるが、クマの場合はその軌跡が行動圏に従うパターンを示すことが多い。それをダミーデータ生成の際に再現するため、以下に示すようなアルゴリズムを考案した。

### a. 確率論的行動シミュレーションモデルとトラップにおける捕捉

単純なランダムウォークにおいては、時点tにおける個体の位置$\mathbf{u}_{\mathrm{t}}$は、t-1の位置$\mathbf{u}_{\mathrm{t-1}}$を中心とする2次元正規分布に従うと考えることができる。この分布は、単位時間当たりの個体の移動を規定する移動距離分布である。

$$ \mathbf{p}(\mathbf{u}_{\mathrm{t}};\mathbf{u}_{\mathrm{t}-1},\sigma^{2})=\mathrm{N}_{2}(\mathbf{u}_{\mathrm{t}-1},\sigma^{2}\mathrm{I}) $$

$ \sigma^{2} $ は移動距離を決めるパラメータであり、大きな値を取るほど単位時間当たりの移動距離は長くなる。

ここで、このようなランダムウォークモデルに行動圏の制約を組み込むことを考える。一般に、行動圏は十分長い時間行動したときの個体位置の系列が従う分布としてみなされる。ここで、長い時間行動したときの個体位置が従う行動圏分布を任意の確率分布 $ \mathbf{p}_{\mathrm{h}}(\mathbf{u}) $ として定義すると、時点tにおける個体の位置 $ \mathbf{u}_{t} $ の分布 $ \mathbf{p}(\mathbf{u}_{t};\mathbf{u}_{t-1},\sigma^{2}) $ は次のように記述することができる。

$$ \mathbf{p}(\mathbf{u}_{\mathrm{t}};\mathbf{u}_{\mathrm{t}-1},\sigma_{\mathrm{m}}^{2})\propto\mathbf{p}_{\mathrm{h}}(\mathbf{u})\times\mathrm{N}_{2}(\mathbf{u}_{\mathrm{t}-1},\sigma^{2}\mathrm{I}) $$

t=0における個体の位置$\mathbf{u}_{0}$を与え、そこからこの分布に従う乱数列を逐次発生させることで、行動圈分布に制約される行動シミュレーションを行うことができる。このようなプロセスの模式図を図(4)-4に示した。図のように、時点tにおける個体の位置が行動圈分布の外側にある場合、移動距離分布と行動圈分布の積として表わされる時点t+1の個体位置分布は行動圈分布の中心に偏ったものになる。その結果、個体の軌跡は行動圈分布の確率密度が高い範

図(4)-4 行動圏分布と移動距離分布による行動シミュレーションの模式図

囲を動くことになる。なお、マルコフ連鎖の収束定理より、初期状態における個体の位置によらずこのプロセスの定常分布はホームレンジ分布 $ p_{h}(u) $ に一致する。

さらに、個体がトラップに接近したときのトラップの利用と検出についても確率モデルにより表現した。クマ個体とトラップの距離が近いほどトラップにおける検出率は高まるため、位置 $ x_{j}=(x_{1j},y_{1j}) $ にあるトラップjに時点 $ t+1 $ で検出される確率 $ p_{c} $ を個体の位置とトラップの距離の関数として以下のように表わした。

$$ \mathbf{p}_{\mathrm{c}}(||\mathbf{u}_{\mathrm{t}}-\mathbf{x}_{\mathrm{j}}||)=\exp(-(||\mathbf{u}_{\mathrm{t}}-\mathbf{x}_{\mathrm{j}}||/\sigma_{\mathrm{c}})^{2}) $$
この関数は正規分布と形状が同じであるが、距離0における検出率は1となる。シミュレーションの各時点において最も近くにあるトラップに対して $ p_{c} $ を計算し、ベルヌイ乱数により検出の有無を決定した。

### b. ダミーデータの生成

ダミーデータの生成においては、全てのトラップを1週間に1回、計6回の点検を行い、その期間内に検出された個体番号とそのトラップを記録する調査を想定したシミュレーションを行った。その結果、検出個体、トラップ、セッションそれぞれにおける検出の有無を示す3次元の行列が検出履歴として得られる。トラップは、2km格子の交点に $ 10\times10=100 $ 個を設置したと仮定し、トラップ位置の凸包は20 km四方の正方形となる。

ダミーデータ生成の際のパラメータの設定は下記のとおりである。

・真の生息密度：0.2個体/km²

・行動シミュレーションの時間ステップ：10ステップ/日

・1ステップ当たりの移動距離分布のパラメータσ：0.18

※この値は、1日の移動距離の平均を1.28 kmとすることに相当する。

・行動圈分布ph(s)：95%確率円の面積が4.85 km²となる2変量正規分布

・検出確率の距離依存性パラメータσ:0.02

※この値は、検出確率の変曲点が20 mにあることを意味する。

以上の条件を以下の二つのシナリオについて40セットのダミーデータを生成し、検証した。

・シナリオ1：すべての個体の行動圏分布は円形

・シナリオ2：個体ごとに、行動圏分布の長径/短径が1:2でランダムに変化（楕円形）

個体の行動圈の中心は、トラップ設置範囲を中心とした30×30 kmの範囲にランダムに発生させた。シミュレーションの初期状態として、t=0における個体の位置を各個体の行動圈分布からランダムに発生させ、上記のアルゴリズムに従いデータ生成を行った。1個体についてデータ生成を行った一例を図(4)-5に示した。シミュレートされた行動の軌跡は行動圈分布の中にとどまることがわかる。また、一般的な標識再捕獲調査は短期間のうちに行われるため、その間で個体は行動圈分布全体をくまなく利用できるとは限らない。シミュレーションが生成する行動パターンはそのような性質も再現している。
図(4)-5 データ生成シミュレーションの1例。青色の同心円は内側が行動圏分布の標準偏差半径円、外側が95%確率円、緑色の折れ線はシミュレートされた行動軌跡、■は非検出トラップ、■は検出トラップを示す。

### c. 密度推定手法の比較

本研究においては、(i) 平均最大移動距離法、(ii) 空間明示最尤法、(iii) ベイズ空間明示標識再捕獲モデルの推定精度の比較を行った。ダミーデータを用いた比較においては、それぞれのモデルにおいては40回の繰り返しを行い、推定値と真の値のずれを95%CI coverage（95%CIが真の値を含む割合）と平均95%CI幅により評価した。下記にそれぞれの推定手法の概要を示した。

## (i) 平均最大移動距離法(MMDM) $ ^{5)} $

この手法は、複数回検出された個体ごとに検出トラップ間の距離の最大値を計算し、その平均値の半分をホームレンジ円の半径rとみなすものである。トラップから生成した半径rのバッファに含まれる範囲の面積を有効トラップ面積とし、任意の閉鎖個体群モデル(例えば、

Otis(1978)の捕獲率均一モデル $ M_{0} $ や不均一モデル $ M_{h} $ による生息数推定値とあわせて生息密度を計算する手法である。この方法は簡便であるが、有効トラップ面積の計算方法はアドホックであり統計学的な根拠に乏しいとされる。なお、今回の評価における生息数の推定には、捕獲率均一モデル $ M_{0} $ を用いた。

## (ii) 空間明示最尤法 (Program “DENSITY”) 6

空間明示最尤法は、トラップ設置範囲を含めた十分に広い空間S内において動物の行動圏中心の空間分布を記述するサブモデルと、特定のトラップにおける捕獲確率と個体と行動圏中心からトラップまでの距離からなるモデルの二つのサブモデルから生息密度を推定する手法である。行動圏中心の位置は実際にはわからないため、行動圏中心の位置をポアソン過程にしたがってランダムに発生させた状態において、1個体ずつ行動圏中心と各トラップにおける捕獲確率の尤度計算を行う。
各個体は円形の行動圏中心を持ち、行動圏中心からの距離が遠くなればなるほど低下すると考える。このため、トラップにおける捕獲率は多くの場合ガウス分布（図(4)-6）を仮定し、捕獲率は、以下のように記述される。

$$ \mathrm{p}=\mathrm{g}0\times\exp(-\frac{\mathrm{d}^{2}}{2\sigma^{2}}) $$

pはトラップにおける捕獲確率、g0は行動圏中心における捕獲率、σが距離に対する捕獲確率の減衰、dが行動圏中心からあるトラップまでの距離を表す。つまり、g0とσの2つのパラメータによってトラップにおける捕獲確率が推定されることにより、設定した空間S内の行動圏中心の数が推定される。さらに、生息密度については、推定された行動圏中心の数を空間Sの面積で割ることで算出することができる。

行動圏中心からの距離

図(4)-6 捕獲確率と行動圏中心の関係

## (iii) ベイズ空間明示標識再捕獲モデル $ ^{7)} $ $ ^{8)} $

生息密度、各個体の行動圈の位置、トラップでの捕獲の関係を階層ベイズモデルとして記述し、検出位置情報から各個体の行動圈の中心、捕獲率と距離の関係、生息密度推定を同時に推定する手法である。行動圈の中心とトラップの距離が大きいほど捕獲される確率は低下すると考えられるため、距離rと捕獲確率pの関係を下記のように定式化している。

$$ \mathbf{p}=\mathbf{p}_{0}\times\exp(-(\mathbf{r}/\sigma)^{2}) $$

このモデルでは、捕獲確率は距離に対して釣鐘型(ガウス型)に減衰すると仮定されている。$\mathbf{p}_{0}$ は距離0の時の捕獲確率、$\sigma$ は釣鐘型曲線の変曲点までの距離で表わされる捕獲確率の減衰を決めるパラメータであり、これらはデータから推定される。Tセッションの標識再捕獲調査は、この確率を期待値とする繰り返しTの2項過程とみなすことができる。

$$ \mathrm{Y}_{\mathrm{ij}}\sim\mathrm{Binomial}(\mathrm{p},\mathrm{T}) $$

$ Y_{ij} $ は個体iのトラップjにおける捕獲回数である。

各個体iの行動圈の中心siは未知であり、これもデータから推定する必要がある。そこで、トラップを含む任意の長方形の空間Sを設定し、siに対してSの範囲に関する一様事前分布を与える。調査において1度も検出されなかった個体に対してもsiを定義するが、それによってsiの長さは生息数によって変化することになる。推定には、パラメータ値を確率的探索的に変化させながら、データに対して当てはまりのよいパラメータ値のセットを得るマルコフ連鎖モンテカルロ法(MCMC法)を用いる。それにより、少なくとも1回検出された個体のsiはトラップの周囲
を動き、1度も検出されなかった個体のsiは消失・生成を繰り返しながらSの中をほぼ自由に動き回ることになる。個体群中に一度も採取されなかった個体がどの程度存在するかはData Augmentation（データ付加）という手法により、逆算することができる。MCMC法の過程で得られるsiの数は、空間S内の生息数Nとみなすことができる。したがって、生息密度はNをSの面積で割った値となる。

また、ベイズ空間明示標識再捕獲モデルを統計解析環境RのGUI上で操作が可能なパッケージSPACECAPについても評価を行った。上記のベイズ空間明示標識再捕獲モデルと異なる点として、例えば行動圏中心を1000 m間隔といったように格子状に配置する点があげられる。本研究においては、行動圏中心は500 m間隔として密度推定を行った。

## 3 ）野外データおよび他地域のデータを用いた推定精度の検証

ダミーデータと同様に平均最大移動距離法、空間明示最尤法とベイズ空間明示標識再捕獲モデルを用いて、2010年と2011年の野外データ用いた密度推定値を比較することで推定手法のバイアスの検討を行った。さらに、2010年と2011年ではトラップ数や配置など違いがあるため、2010年と2011年で共通な3次メッシュにおける80基のトラップから採取されたデータを用いて、生息密度推定を行った。なお、2010年と2011年のデータは同じ基準で識別されたDNAデータを解析に用いた。

## (2) トラップ配置の検証

## 1 ）北上山地のデータを用いた解析

単位面積当たりのヘア・トラップ数、セッション数と生息密度推定値の関係を検証するため、ヘア・トラップ数、セッション数を変化させた場合の生息密度推定値の変動を調べた。ヘア・トラップ数については、約2.5、3、4、5、6、8 km²あたり1基のヘア・トラップ密度とした場合の生息密度推定を行った。245基のヘア・トラップによる標識再捕獲のデータを用い、設置したトラップの中から、上記のトラップ密度になるようにランダムにトラップを抽出して解析に用いた。ただし、トラップの位置に空間的な偏りがないデータを使用した。セッション数については、ヘア・トラップ数と同様に6セッションからランダムに3、4、5セッションを抽出したデータセットを作成し、生息密度推定を行った。このような解析をそれぞれのトラップ密度、セッション数において20回行い、平均値、95%上限値、下限値の平均値をそれぞれ算出した。このように、実際の調査より努力量が低い場合の推定精度については、実証的に評価することができる。

## 2 ）他の地域の事例を用いた解析

岩手県の事例以外でも推定方法の適用が可能であるかを検証するため、ヘア・トラップ先行調査事例（富山、奈良、滋賀（2007年、2008年）、神奈川、京都、福井、山梨）のデータを用いて、生息密度の推定と分析成功率に影響するパラメータの評価を行った。生息密度推定には空間明示最尤法 $ ^{6)} $ とHartstack et al.(1971) $ ^{14)} $ の空間明示モデルを用いた。なお、再捕獲数が0または、同じトラップのみでしか再捕獲がなかったデータについては解析を行わなかった。次に、上記の7地域8サンプルに北上山地における2010年と2011年のデータを加えて生息密度推定の成功率（分析成功率）を推定した。データは生息密度推定が成功したか否か、トラップ数、再捕獲数、セッション
数、トラップ間の距離を用いた。生息密度推定の成否を目的変数、それぞれのパラメータを説明変数としたロジスティック回帰を行い、生息密度推定の成否とトラップ数、再捕獲数、セッション数、トラップ間の距離との関係を評価した。

### 4. 結果及び考察

## (1) 生息密度推定手法の精度検証

1）空間明示モデルと空間非明示モデルおよびトラップ配置ごとの推定値の比較

均一配置と北上山地モデル調査地における実際の配置である2種類のトラップの配置、および空間非明示モデルと空間明示モデルの組合せによる過去のヘア・トラップ調査データに基づく生息密度推定結果を図(4)-7に示した。トラップの配置によらず、空間非明示モデルは実際よりも高い生息密度が推定された一方、空間明示モデルの推定値は真の値に近い推定が得られた。また、空間非明示モデルを用いた場合、均一トラップ配置ではトラップ設置予定地よりも大きな生息密度が推定された。従来の手法における生息密度の算出においては、母数となる調査地面積の決め方によって誤った生息密度を導いてしまう可能性があることが明らかとなり、空間非明示のモデルは推定精度が低い可能性が示された。一方で、空間明示モデルは真の値に近い推定値を示し、推定精度が高い可能性が示された。

図(4)-7 トラップ設置予定地、均ートラップ配置における、空間明示モデルと空間非明示モデルの生息密度推定値。黒点は事後分布中央値、エラーバーは95%信用区間を示す。赤破線は真の個体密度(0.195/km²)を示す。

## 2 ）ダミーデータを用いた推定精度の検証

行動圈を円形と仮定したダミーデータによる推定結果を図(4)-8、それぞれのモデルの95%CI幅、相対バイアスと95%CI coverageを表(4)-1に示した。平均最大距離法は生息密度を過小評価する傾向があったが、少数のケースで真の値よりも高い値を示していた。平均的に95%CIの幅は小さく、coverageは24%であった。空間明示最尤法の推定値は真の値に近かったが、平均95%CIの幅が広く、真の値よりも大きい方へ偏りがみられた。ベイズ空間明示標識再捕獲モデルは4つのモデルの中で95%CI coverageが高くかつ、平均95%CI幅も狭く、精度の高い推定手法であることが明らかとなった。Coverageでは空間明示最尤法とベイズ空間明示標識再捕獲モデルではほぼ同等の評価であったが、空間明示最尤法は95%CI幅が4つのモデルの中で最も大きかった。SPACECAPは、25%点が真の値を上回っており、過大評価の傾向がみられた。さらに、95%CI幅は比較的狭かったものの、coverageはやや低い結果となった。また、真の値よりもやや高い推定結果となる傾向があった。行動圈を楕円形としたダミーデータによる推定結果では、空間明示最尤法、ベイズ空間明示標識再捕獲法ともに若干高い値を示したものの、円形行動圈のシナリオと同様の傾向が見られた。
表(4)-1 円形を仮定したダミーデータ解析における各密度推定手法の95%CI幅と95%CI coverage

手法	平均95%CI幅(95%区間)	95%CI coverage
平均最大距離法 (n=40)	0.06 (0.02-0.15)	0.225
空間明示最尤法 (n=40)	0.26 (0.15-0.42)	1.00
ベイズ空間標識再捕獲モデル (n=40)	0.13 (0.09-0.15)	0.975
SPACECAP (n=40)	0.19 (0.06-0.32)	0.375

図(4)-8 ダミーデータによる生息密度推定モデルの推定結果。箱ひげ図は、内側から中央値、四分位点、四分位点×1.5の範囲の最大値を示し、○はそれより外側の値を示している。赤点線は真の推定密度0.2頭/km²を示す。MMDMは平均最大移動距離法を示す。

## 3 ）野外調査データを用いた精度検証

北上山地における野外調査を4つのモデルで解析した結果を表(4)-2に示した。ベイズ空間明示標識再捕獲モデルを用いて生息密度推定を行なったところ、2010年は0.400（0.340-0.489）個体/km²となった。2011年のデータを用いた推定結果は、0.188（0.147-0.235）個体/km²であった。空間明示最尤法による生息密度推定では、2010年と2011年の推定値は、それぞれ0.403（0.329-0.493）、0.185（0.139-0.245）個体/km²であった。SPACECAPを用いた生息密度推定値は2010年と2011年の
それぞれで0.374（0.298-0.437）、0.173（0.139-0.213）個体/km²であった。一方で、平均最大移動距離法を用いた推定結果は、2010年が0.309（0.261-0.381）、2011年が0.191（0.165-0.242）となった。2010年と2011年では様々な条件の違いがあったものの、ベイズ空間明示標識再捕獲モデルと空間明示最尤法、それぞれの年で近い推定結果を示した。SPACECAPについてはやや推定値が低くなる傾向があったが、ベイズ空間明示標識再捕獲モデルと空間明示最尤法と近い値を示した。一方で、平均最大距離法では、2011年は他の手法との差がほとんどなかったが、2010年においてはダミーデータほど極端ではないものの、他の手法よりも低い推定値を示した。野外データを用いた結果においてもダミーデータによる推定結果と類似した結果が得られた。以上の結果から、野外において標識再捕獲データを用いた生息密度推定においてもベイズ空間明示標識再捕獲モデルと空間明示最尤法、SPACECAPは生息密度推定に有用な手法である可能性が示唆された。

表(4)-2 北上山地のデータにおける各手法による推定生息密度

手法	2010年推定生息密度(95%CI)（頭/km²）	2011年推定生息密度(95%CI)（頭/km²）
平均最大距離法	0.309 (0.261-0.381)	0.191 (0.165-0.242)
空間明示最尤法	0.403 (0.329-0.493)	0.185 (0.139-0.245)
ベイズ空間明示標識再捕獲モデル	0.400 (0.340-0.489)	0.188 (0.147-0.235)
SPACECAP	0.374 (0.298-0.437)	0.173 (0.139-0.213)

2010年と2011年の密度推定結果の違いを検証するため、2011年と同じ3次メッシュ内に位置する80基のトラップを用いて、ベイズ空間明示標識再捕獲モデルにより生息密度推定を行った。その結果、2010年の密度推定値が0.323(0.232-0.489)個体/km²、2011年の密度推定値が0.188(0.147-0.235)個体/km²であり、同じ3次メッシュに設置されていた80基のトラップデータを用いても2010年の方が2011年よりも推定値が大きくなった。推定されたパラメータ値については2010年と2011年では2011年に捕獲確率がほぼ2倍となった（表(4)-3）。さらに、2011年はσが2010年に比べ1.4倍ほど大きい値となっている。つまり、1個体あたりの行動圏が2011年の方が広く推定され、生息密度が低く推定されたことになる。2010年の245基と80基のデータでは、80基のデータの方がややσが大きかったが、その他のパラメータはほぼ同様の値を示した。

表(4)-3 ベイズ空間明示標識再捕獲モデルによるパラメータ推定値。中央値および括弧内に95%CIを示す

2010年	2010年80トラップ	2011年
密度	0.400 (0.340-0.489)	0.323 (0.232-0.489)	0.188 (0.147-0.235)
捕獲確率（p）	0.044 (0.033-0.059)	0.044(0.025-0.066)	0.090 (0.064-0.124)
行動圈中心から変曲点までの距離（σ）	1.76 (1.59-1.98)	2.15(1.82-2.63)	2.39 (2.09-2.81)
トラップでの検出個体の割合（Ψ）	0.498 (0.419-0.610)	0.496 (0.352-0.756)	0.181 (0.136-0.236)

ダミーデータと野外データによる推定精度比較の結果においても、ベイズ空間明示標識再捕獲モデルと空間明示最尤法は平均最大移動距離法よりも推定精度が高いと考えられた。平均最大移動距離法による有効トラップ面積の計算は再捕獲の回数に影響を受けやすく、再捕獲が少ないと有効トラップ面積は小さくなり、密度を過大評価することがあると知られている $ ^{5} $ 。また、生息数の推定に用いた $ M_{0} $ モデルは、個体ごとの行動圏の位置などに起因する検出率の違いがある状況
では、生息数を過小評価する傾向がある。ダミーデータ解析の結果から明らかになった平均最大距離法の推定生息密度の過小評価およびばらつきの大きさは、これらのバイアス要因が複合的に作用した結果であると考えられる。空間明示最尤法は、距離に依存する捕獲率を明示的に扱っているモデルであるが、そのような空間構造は空間明示最尤法においては平均移動距離として集約されてしまい、各個体の捕獲位置に関する情報は考慮されない。そのことがダミーデータ解析における95%CI幅の広さにつながったと考えられる。その点で、捕獲位置に関する情報を最大限に活用しているベイズ空間明示標識再捕獲モデルはより優れていると考えられる。

空間明示最尤法のメリットとしては、計算時間の短さとAICによるモデル選択が可能である点が挙げられる。北上山地のデータを用いた場合、ベイズ空間明示標識再捕獲モデルは推定値を得るまでに数日を要した（Core2 DUO E8600搭載のPCを使用）。一方、空間明示最尤法では10分程度で結果を得ることができた。多数の場所、または広域を対象とした標識再捕獲調査を行い、推定に際してさまざまな試行錯誤を行う必要があるときには、ベイズ空間明示標識再捕獲モデルの適用が計算時間の面で非現実的になる場合があるかもしれない。その場合には、空間明示最尤法は比較的よい代替になりえると考えられる。さらに、AICによるモデル選択を行うことで、例えば雌雄差を考慮に入れたモデルと、雌雄差を入れないモデルで推定を行い、推定結果の比較をモデル選択で行うことも可能である。

2010年と2011年の推定値が大きく異なった原因についてははっきりしないものの、すべての手法において2011年の方が低い傾向を示したことから推定手法の問題というよりもデータに起因している可能性が高いと推測される。2010年と2011年ではトラップを設置した空間スケールの違いなど調査デザインの違いばかりでなく、体毛サンプルの試料採取率や再捕獲率（ヘアトラップ班報告）、さらには、捕獲確率やσのパラメータ推定値にも違いが見られた。調査動物の行動圏に対する調査面積 $ ^{15} $ や再捕獲個体の不均質性（性別、齢、経験など） $ ^{16} $ によって生息密度推定値がバイアスを受けることが報告されている。2010年と2011年の間で行動圏の違いがある可能性が示唆されたので、上記の要因が複合的に影響したのかもしれない。推定値の違いについては以上のようなことが推察できるものの、今後さらなる精査が必要である。

## (2) トラップ配置の検証

## 1 ）トラップ数、セッション数の検証

単位面積あたりのヘア・トラップ数とセッション数を変化させ密度推定を行なった結果を図(4)-9に示した。トラップ数においては、 $ 4km^{2} $ に1基まではトラップ密度の変化にも関わらず、推定値の平均値は0.40個体/ $ km^{2} $ 程度の値を示し、変動は小さかった。ただし、95%CIはトラップ数が少ないほど広くなる傾向があった。 $ 5km^{2} $ あたり1基のトラップ数では中央値が0.320個体/ $ km^{2} $ とやや推定値が小さくなった。 $ 6km^{2} $ に1基のトラップ密度では、中央値が0.301個体/ $ km^{2} $ 、 $ 8km^{2} $ に1基では、0.355個体/ $ km^{2} $ となり $ 5km^{2} $ と同様に $ 4km^{2} $ に1基までに比べるとやや低い推定結果となった。セッション数を変えて生息密度推定を行なった結果では、3、4、5セッション共に6セッションでの密度推定値である0.40個体/ $ km^{2} $ 程度と大きくは異ならなかった。ただし、推定値の95%CIはセッション数が少なくなるほど広くなる傾向があった。単位面積あたりのトラップ数は地域ごとのクマの生息密度によって異なると考えられるが、トラップ数が生息密度推定の精度に影響しやすいことが示唆された。
図(4)-9 トラップ数とセッション数と密度推定結果の関係。エラーバーは95%CIを示す

2）他地域のデータを用いた生息密度推定結果および分析成功率とトラップの関係

他地域のデータを用いて生息密度推定を行った結果を表(4)-4に示した。空間明示最尤法では再捕獲数が5個体以下と少ない地域（奈良、滋賀2007年）では生息密度を推定することが出来なかった。Hartstack et al. $ ^{[14]} $ の方法では全ての地域で密度を推定することができたが、再捕獲数が少ないほど95%CIが広くなる傾向がみられた。特に奈良県の事例では95%CIの下限と上限の値が100倍近く異なり、推定精度が低いと考えられた。以上のことは、再捕獲数の重要性を示唆しており、個体数が少ない地域や再捕獲数の少ない地域では密度推定を行う際に注意が必要である。推定できた地域においては、空間明示最尤法とHartstack et al.の方法で推定結果に大きな違いは見られなかった。

表(4)-4 他地域のデータを用いた個体数推定結果。95%CIを括弧内に示した。

地域	セッション数	再捕数	トラップ数	空間明示最尤法	Hartstack et al.の方法
富山	5	9	53	0.059(0.029-0.122)	0.087(0.049-0.164)
奈良	8	2	47	-	1.04(0.11-10.4)
滋賀 2007年	10	4	44	-	0.001(0.0003-0.007)
滋賀 2008年	11	7	44	0.1(0.045-0.26)	0.087(0.034-0.234)
丹沢	8	9	66	0.016(0.008-0.031)	0.022(0.015-0.032)

表中の──は生息密度推定ができなかったことを表す

生息密度推定の成否とパラメータの関係を評価したところ、サンプル数、再捕獲数、トラップ数が生息密度推定の成否と関連することが明らかとなった。ロジスティック回帰に用いたサンプル数が少なかったため、推定値の標準誤差が非常に大きいことを考慮に入れる必要があるが、ヘア・トラップでの採取サンプル数の少なさ、再捕獲数やトラップ数の少なさが生息密度推定の失敗に関係することが明らかとなった。また、調査を行う際に考慮ができるトラップ数およびセッション数の関係と分析成功の関係を図(4)-10に示した。セッション数は少なくても密度推定が出
来ている事例はあるが、トラップ数が少ないと密度推定は失敗する傾向が見られた。さらに、再捕獲数とセッション数、トラップ数の関係を図(4)-11に示した。トラップ数の違いはあるものの、セッション数は再捕獲数と負の相関を示しているのに対し、トラップ数は再捕獲数と正の相関を示した。特に、セッション数が7から11のデータではほとんどトラップ数も変わらないため、セッション数を増やすことが必ずしも再捕獲数の増加にはつながらないかもしれない。トラップ数を可能な限り多くすることで、再捕獲数が確保され、標識再捕獲法を用いた生息密度推定の成功につながると考えられる。

図(4)-10 分析成功とトラップ数およびセッション数の関係。○×が各地域における事例をそれぞれ示す

図(4)-11 再捕獲数とセッション数、トラップ数の関係
## (3) まとめ

生息密度推定手法間の精度検証を行った結果、標識再捕獲法による生息密度推定には、ベイズ空間明示標識再捕獲モデルないしは空間明示最尤法を用いることが望ましいと考えられた。今後、標識再捕獲データを用いた生息密度推定において標準的な手法になることが期待される。ベイズ空間明示標識再捕獲モデルはフリーソフトWinBUGS (http://www.mrc-bsu.cam.ac.uk/bugs/)により推定が可能であり、フリーの統計解析環境Rからデータおよびモデルコードの読み込み、推定結果の出力までの一連の操作を行うことができる。また、検証を行ったRのパッケージ"SPACECAP"¹は、GUI上での操作で推定を実行することも可能である。空間明示最尤法については、Program "DENSITY"²で実装されており、これについてもGUI上での操作が可能である。また、同様の機能をもつRのパッケージ"secr"も公開されている。このように、すでにソフトやパッケージが開発されていることも空間明示型標識再捕獲モデルの敷居を低くしている。

トラップ設置については、本研究の結果から、不均質なトラップ配置においてもベイズ空間明示標識再捕獲モデルを用いることにより頑健な推定結果が得られる可能性が高いこと、トラップ数が推定精度や推定の成功に重要であることが明らかとなった。これまで、生息密度推定の精度を高めるには、再捕獲数が多いこと $ ^{4} $ やジェノタイピングエラーを減らすこと $ ^{17} $ などが指摘されてきた。本研究の成果から推定精度や推定の成功にはトラップ数を増やすことが有効であることが示唆された。ただし、実際のヘア・トラップのデータ解析においては、常に空間明示モデルの推定が良い結果を与えるとは限らない。実際の調査においては地形的な要因や季節変化などのさまざまな要因が資料採取率に影響を与えると考えられ、それらの影響については現実を単純化したダミーデータによる検証は不可能である。実際に調査を行う際には、あらかじめ予備調査を行った上で、行動圏サイズ、再捕獲率を考慮に入れ、トラップ数やセッション数を決定する必要があると考えられる $ ^{18} $ 。その上で、複数年にわたって調査を行うことで、このような不確実性にも対応できるかもしれない。

## 5 ．本研究により得られた成果

## （1） 科学的意義

捕獲の空間情報を用いることにより、個体の移動距離の情報を活用し、より客観的な密度推定を行う手法を確立できた。精度が高い生息密度推定法の確立により、クマ類の個体群動態の把握や堅果類の豊児との関係の解明にも貢献できると考えられる。

## （２）環境政策への貢献

本研究によって、生息密度推定の精度や成功にはトラップ数や再捕獲数が重要であることが明らかとなった。自治体等が密度推定を行う際には、予算や調査面積など状況に応じてトラップ設置方法を決定しなければならない。本研究の結果は精度の高い密度推定を目指してヘア・トラップ設置における予算配分を行うための資料として役立つと期待される。さらに、精度が高い密度推定手法は保護管理において有用な手法であり、クマ類の保護管理政策への貢献が期待される。
研究成果を自治体等が利用する際の参照とするため、下記のマニュアルを作成し、都道府県鳥獣保護部局へ配布した。

クマ類の個体数を調べる：ヘア・トラップ法とカメラトラップ法の手引き（統合版）（印刷配布）URL：http://www.bear-project.org/pdf/Tebiki/tebiki-tougou.pdf

## 6 ．国際共同研究等の状況

特に記載すべき事項はない

## 7 ．研究成果の発表状況

## （1） 誌上発表

<論文（査読あり）>

1) Ohta, U., Jusup, M., Mano, T., Tsuruga, H. and Matsuda, M. 2012. Adaptive Management of the Brown Bear Population in Hokkaido, Japan. Ecological Modelling, DOI:10.1016/j.ecolmodel.2012.05.011 (in press, online first)

<査読付論文に準ずる成果発表>

特に記載すべき事項はない

くその他誌上発表（査読なし）>

特に記載すべき事項はない

## （2） 口頭発表（学会等）

1) 深澤圭太・太田海香・松田裕之：第16回野生生物保護学会・日本哺乳類学会2010年度合同大会、自由集会F-10「クマ類の個体数を推定する：DNAマーカ個体識別の有効性と課題」(2010)「空間明示型Capture-Recaptureモデル：個体密度はどのように計算すべきか」

2）深澤圭太：2011年度日本哺乳類学会自由集会（2011年9月、宮崎）、自由集会「クマ類の個体数推定法の開発：環境省究総合推進費プロジェクトから」。（2011）空間明示型標識再捕獲モデルの密度推定精度の比較検討.

3) Ohta, U., Mano, T., Tsuruga, H., Jusup, M. and Matsuda, H.: Session "Wildlife Management Challenges - Bears around the world", Pathways to Success 2010 Conference, (September, 2010, Colorado, USA). 'Adaptive management plan of brown bear ( $ Ursus\ arctos $ ) in Oshima peninsula, Hokkaido, Japan', (2010).

4) 太田海香・間野勉・釣賀一二三・Marko Jusup・松田裕之：第26回個体群生態学会大会(2010年9月、横浜). 北海道渡島半島地域におけるヒグマ個体群の順応的管理手法の開発，(2010)

5) 山内貴義・鞍懸重和・深澤圭太・米田政明：2011年度日本哺乳類学会（2011年9月、宮崎）.

ヘア・トラップ法を用いた岩手県北奥羽地域個体群に生息するツキノワグマの生息数推定（2011）.

6) 山内貴義・鞍懸重和・深澤圭太・米田政明：第17回日本野生動物医学会発表要旨（2011年9月、東京）．ヘア・トラップ法を用いた岩手県北上高地北部地域個体群に生息するツキノワグ
マの生息数推定. （2011）.

## （3） 出願特許

特に記載すべき事項はない

## （４）シンポジウム、セミナーの開催（主催のもの）

太田海香, 今野建志郎, シンポジウム「現場への適用を目指した個体群管理モデルの展望と課題」,『第26回個体群生態学会大会』,横浜, (2010年)

## （５）マスコミ等への公表・報道等

特に記載すべき事項はない

## (6) その他

特に記載すべき事項はない

### 8 \. 引用文献

1）森光由樹. 2008. 各都道府県のヘア・トラップ調査の実施状況と長野県における実施例．哺乳類科学48:133-138.

2) Dice, L. R. 1938. Some census methods for mammals. Journal of Wildlife Management 2:119-130.

3) Lukacs, P. M. and Burnham, K. P. (2005) Estimating population size from DNA-based closed capture-recapture data incorporating genotyping error. The Journal of Wildlife Management 69:396-403.

4) Efford, M. G., Borchers, D. L. and Byrom, A. E. 2009. Density estimation by spatially explicit capture-recapture::likelihood-based methods. Pp259-269 in D. L. Thomson, E. G. Cooth, and M. J. Conroy, editors. Modelling demographic processes in marked populations. Springer, New York, USA.

5) Wilson, K. R. and Anderson, D. R. 1985. Evaluation of two density estimators of small mammal population size. Journal of Mammalogy 66:13-21.

6) Borchers, D. L. and Efford, M. G. 2008. Spatially explicit maximum likelihood methods for capture-recapture studies. Biometrics, 64:377-385.

7) Gardner, B., Royle, J. A. and Wegan, M. T. 2009. Hierarchical models for estimating density from DNA mark-recapture studies. Ecology 90:1106-1115.

8) Royle, J. A., Karanth, K. U., Gopalaswamy, A. N. and Kumar, S. N. 2009. Bayesian inference in camera trapping studies for a class of spatial capture–recapture models. Ecology 90: 3233–3244.

9）岩手県．2007．第2次ツキノワグマ保護管理計画（巻末参考資料；ツキノワグマの生息数の推定について）．岩手県.

10) Otis, D. L., Burnham, K. P., White, G. C. and Anderson, D. R. 1978. Statistical inference from capture data on closed animal populations. Wildlife Monographs 62:1-135.

11）米田一彦. 1990. 秋田県太平山地域におけるツキノワグマの生態・テレメトリー調査. In 人間活動との共存を目指した野生鳥獣の保護管理に関する研究 (eds. 財団法人日本野生生物研究
センター). 環境庁自然保護局, 東京.

12) Izumiyama, S. and Shiraishi, T. 2004. Seasonal changes in elevation and habitat use of the Asiatic black bear (Ursus thibetanus) in the Northern Japan Alps. Mammal Study 29:1-8.

13) Miura, S. and Oka, T. 2003. Evaluation of apple bait hair-traps for genetic tagging of Asiatic black bears in the Kitakami Highland, northern Honshu, Japan. Mammal Study 28:149-152.

14) Hartstack, A. W., Hollingsworth, J. P., Ridgway, R. L. and Hunt, H. H. 1971. Determination of trap spacings required to control an insect population. Journal of Economic Entomology 64, 1090-1100.

15) Marques, T. A., Thomas, L. and Royle, J. A. 2011. A hierarchical model for spatial capture-recapture data: comment. Ecology, 92:526-528.

16) Ebert, C., Knauer, F., Storch, I. and Hohmann, U. 2010. Individual heterogeneity as a pitfall in population estimates based on non-invasive genetic sampling: a review and recommendations. Wildlife Biology, 16:225-240.

17) Knapp, S. M., Craig, B. A. and Waits, L. P. 2009. Incorporating Genotyping Error Into Non-Invasive DNA-Based Mark-Recapture Population Estimates. Journal of Wildlife Management 73:598-604

18) Foster, R. J. and Harmsen, B. J. 2012. A critique of density estimation from camera-trap data. Journal of Wildlife Management, 76:224-236.