兵庫ワイルドライフレポート 1: 44-55. 2012

原著論文

イノシシの個体群動態の推定（兵庫県本州部 2011 年）
坂田宏志 1,2 *・岸本康誉 1,2・関香菜子 11 兵庫県森林動物研究センター・2 兵庫県立大学自然・環境科学研究所

要 点
・ 2011 年までに入手されたデータから、兵庫県本州部のイノシシの自然増加率や個体数を、階層ベイズモデルを構築し、マルコフ連鎖モンテカルロ法によって推定した。
・ 推定は、銃猟時の目撃効率、狩猟捕獲数、有害捕獲数のデータを基に、それぞれデータの誤差変動を組み込んだモデルを構築した。
・ 自然増加率は、66.5%(90%信頼限界で 15.2~110.4%)と推定された。
・ イノシシは増加率が高いため、推定モデルの中で設定する自然増加の前後で大きく推定値が異なるため、それぞれの個体数を推定した。
・ 推定個体数は、単純な増加傾向にあり、2010 年の自然増加前の段階で、20609.1 頭（90%信頼限界では 11385.6～70761.9 頭程度）
、その自然増加後の最大個体数で35236.5 頭（90%信頼限界では 20878.0 ～84918.3 頭程度）と推定された。
key words：個体数管理 自然増加率 ベイズ推定 マルコフ連鎖モンテカルロ法 個体数推定

1. はじめにこの論文では、兵庫県本州部におけるイノシシ（Sus scrofa）の保全と管理に資するため、自然増加率や生息個体数の推定を行う。
推定には、兵庫県で体系的に収集している 2002 年から 2010 年までのデータを用いる。具
体的には、兵庫県森林動物研究センターが収集している狩猟登録者の報告に基づく銃猟時の平均目撃数（目撃効率）、狩猟による捕獲数、有害捕獲許可による捕獲数である。
自然増加率や個体数の推定は、上記のデータと時系列的な関係を記述する階層ベイズモデルを構築し、パラメータの推定にはマルコフ連鎖モンテカルロ法によるランダムサンプリングを用いる。これらの作業の中では、自然増加率、個体数の他に、捕獲率や目撃効率と個体数の関係を表す係数に加え、観測データに含まれる誤差変動の大きさなどを構築したモデルの中で推定する。

* 連絡先：〒669-3842 兵庫県丹波市青垣町沢野 940
e-mail: sakata@wmi-hyogo.jp

44

兵庫県森林動物研究センター

兵庫ワイルドライフレポート 1 ,2012

2. 方法
対象とする期間と地域
対象期間は、2002 年から 2010 年までとした。対象地域は兵庫県本州部とする。
用いたデータ
推定には以下のデータを使用した。

有害捕獲数（年度）i_yugai[i]：i 年度の有害許可による捕獲数。個体数を反映する指標として用いる。
目撃効率 spue[i]：i 年度の狩猟期間中に、狩猟者登録者から得られた銃猟時の目撃効率。個体数を反映する指標として用いる。
狩猟捕獲数 r_ca[i]：i 年度に狩猟による捕獲数。個体数を反映する指標として用いる。
有害捕獲数（年） y_ca[i]：i 年度の 1 月から i+1 年度の 12 月までの有害許可による捕獲数。
森林面積 f_area ：兵庫県本州部の森林面積。生息密度の期待値を計算する際に用いる。
以上の方法で収集したデータセットを表 1 に示す。

表１ 入力データセット
year

i_yugai

spue

r_ca

y_ca

f_area

2002

1302

0.25

4061

1888

5452.07

2003

1704

0.21

4522

2337

5452.07

2004

2394

0.26

6456

2015

5452.07

2005

2006

0.26

4877

2278

5452.07

2006

2186

0.18

4460

2170

5452.07

2007

2264

0.19

3583

3753

5452.07

2008

3828

0.21

6543

3159

5452.07

2009

2973

0.17

4204

6139

5452.07

2010

6383

0.20

10580

5452.07

推定するパラメータ
以下の考え方に基づいて、 lire、lr_spue、pr、py、lnNins、v_spue、v_ryo、v_yugai の 8の変数について推定し、目的である自然増減率や個体数を推定する。推定変数の初期値と事前分布は表２のとおりで、各推定変数の定義と事前分布設定の際の考え方は、以下のとおりである。
1.

自然増減率の対数値

lire[i]：出生と自然死亡の結果としての雌雄合わせた全個体数に対する増減の比率とする。lire については、環境省の特定哺乳類生息動向調査の個体数推定（環境省生物多様性センター 2011）に採用された事前分布を用いる。また、exp(lire)
を自然増加率 ire とする。
2.

生息密度と目撃効率の比率を示す係数の自然対数値

lr_spue： 事前分布は正規分布を仮定し、事前の情報は十分にないため、その分散は大きめに設定した。また、exp(lr_spue)を rs とした。
3.

狩猟時の捕獲率pr[i]：狩猟の捕獲数の生息個体数に対する比率を表す。この係数は、0 か45

兵庫ワイルドライフレポート 1 ,2012

ら 1 の間で変動すると考え、pr=1/(1+exp(-prp))とし、prp を推定する。また、pr は事前の情報は十分にないため、その分散は大きめに設定した。
4.

有害時の捕獲率 py[i]：有害の捕獲数の生息個体数に対する比率を表す。この係数は、0 から 1－ pr[i]の間で変動すると考え、py[i]=(1-pr[i])/(1+exp(-pyp))とし、pyp を推定する。
また、py は事前の情報は十分にないため、その分散は大きめに設定した。

5.

2009 年の生息個体数の自然対数値

lnNins：2009 年度の個体数に関する事前の情報はないため、分散を大きめに設定した。
6.

目撃効率、狩猟捕獲数、有害捕獲数の期待値からの誤差分散 v_spue v_ryo v_yugai：
それぞれ、観測モデルで示す確率分布の誤差分散として観測データから推定する。これ
らの誤差分散の事前分布は、それぞれ、形状母数、尺度母数ともに 0.01 の逆ガンマ分布を用いた。

7.

各推定変数の初期値は、事前分布の期待値とした。捕獲率の尤度関数の変動部分 v_spue、v_ryo、v_yugai については、それぞれ初期値を 0.01 とした。

表２推定した変数とその初期値および事前分布

正規分布は（期待値，分散）を、逆ガンマ分布は（形状、尺度）をそれぞれ示す。
ブロック

推定変数

初期値

事前分布

1lire

0.0865 正規分布 ((log(1.4)-0.5*0.5),var=0.5)

1prp-0.8473 正規分布(log(0.3/(1-0.3)),var=3)1

pyp-1.3863 正規分布(log(0.2/(1-0.2)),var=3)

1

lr_spue-2.3026 正規分布((log(0.1)),var=5)

1

lnNins10.0000 正規分布(10,var=10)2v_spue0.0100 逆ガンマ分布(0.01,scale=0.01)

2

v_ryo

0.0100 逆ガンマ分布(0.01,scale=0.01)

2v_yugai

0.0100 逆ガンマ分布(0.01,scale=0.01)

個体群動態の過程モデル
個体群動態の過程モデルは、全生息個体数は 2009 年を起点とし、翌年の 2010 年までの変化をN[2010] = ire ×N[2009]－caa [2009]2000 年までの変化を
N[i] = (N[i+1]+caa[i])/ ireのように変化するものと仮定する
ここで、N[i]は、i 年の生息個体数を示す。また、caa[i]は、i 年の捕獲数であり、i 年の狩猟捕獲数 r_ca[i]と有害捕獲数 y_ca[i]の合計値である。
2009 年の個体数は N2009=round(exp(lnNins))とした。なお、生息個体数は、年末時点で46兵庫ワイルドライフレポート 1 ,2012の個体数を想定している。また、自然増加を踏まえて推定できる最大個体数 Nmax[i]を、Nmax[i]= ire ×N[i]として計算した。
観測モデル
推定する個体数と観測されるデータとの関係を示す観測モデルは以下のとおりとする。
1.

目撃効率に関する観測モデル
log(SPUE[i]) = log(rs×N[i]/f_area) -0.5×v_spue + e_spue[i]2.

狩猟捕獲数に関する観測モデル
log(r_ca [i]) = log(pr[i]×N[i]) -0.5×v_ryo + e_ryo[i]3.

有害捕獲数に関する観測モデル
log(i_yugai [i]) = log(py[i]×N[i]) -0.5×v_yugai + e_yugai [i]e_spue[i]、e_fun[i]、e_ryo[i]、e_yugai [i]は、誤差変動を示し、それぞれ期待値 0、分散がv_spue 、v_fun 、v_ryo 、v_yugai の正規分布に従うものとする。
マルコフ連鎖モンテカルロ法
これまで述べたデータとモデルおよび事前分布の設定にもとづいて、マルコフ連鎖モンテ
カルロ法（Gilks et al. 1996）による推定を行った。この推定は SAS/STAT9.3 の MCMC
Procedure を用いた（SAS Institute Inc. 2011）。
サンプリング
推定変数を表２のとおり２つのブロックに分けて、メトロポリス法と conjugate サンプリングによる独立サンプラーを用いて事後分布をサンプリングした。サンプリング回数については、最初の 500 万回はサンプリングせず、次の 2000 万回のうち 2,000 回に 1 回サンプリングし、計 1 万回のサンプリングを行った。
提案分布は、正規分布とし、実際のサンプリング回数に合わせて 5 万回のサンプリングによる事後分布にもとづいて、Roberts et al.（1997）の示した最適な採択率 23.4％を目標に±7.5％の範囲の採択率になるように、スケールと共分散行列のチューニングを行った。
収束判定
収束判定は、有効サンプルサイズ（Kass et al. 1998）と Geweke 検定（Geweke 1992）の2 つの基準で確認した。有効サンプルサイズによる判定では、これが 1,000 以上であることを基準とした。Geweke 法では、サンプリングされたデータのうち、最初の 1,000 回と最後の 5,000 回の期待値の差を検定し、棄却水準が 0.05 にならないことを基準とした。

3. 結果
収束状況
いずれの推定変数についてもサンプリングの際の自己相関はほとんどなく、有効サンプル数は 5,000 を超え、良好なサンプリングができたと判断された。Geweke 検定の結果も基準に達しなかった変数はなかった。
推定値
推定した変数の事後分布は表３の通りであった。また、事前分布と事後分布の形状を図１に示した。
47

兵庫ワイルドライフレポート 1 ,2012

表３の結果に基づいて計算した自然増加率(ire)と、目撃効率の係数(rs)、狩猟捕獲率(pr)、有害捕獲率(py)は表４のとおりであった。また、得られたデータの観測値と期待値との関係を図２、図３に示した。
有害捕獲率以外で事後分布の幅は事前分布の幅より狭まった。有害捕獲率を表す変数は、絞られる幅が少なく、事前分布の設定が推定に影響を与えていた。
自然増加率は、66.5%(90%信頼限界で 15.2~110.4%)となり、推定幅は広かった。また、捕獲率も狩猟捕獲率が 39.3%(90%信頼限界で 9.2~64.5%)、有害捕獲が 18.9%(90%信頼限界で4.6~31.7%)となり、推定幅は広かった。
さらに、これらの結果に基づいて計算した個体数と最大個体数、増加個体数を表５、それらの動向を図４、図５、図６に示す。個体数は、2002 年以降単調な増加傾向にあり、2010年末の段階で、中央値で 20609.1 頭（90%信頼限界では 11385.6～70761.9 頭程度）と推定され、増加個体数についても、個体数の増加に伴い増加していると推定された。また、最大個体数は、中央値で 35236.5 頭（90%信頼限界では 20878.0 ～84918.3 頭程度）と推定された。
表３ 事後分布の統計量
変数

平均

標準偏差

5%点

中央値

95 %点

lire

0.48140.18820.14170.50960.7439

prp-0.5821

0.9031-2.2809-0.4308

0.5981

pyp-0.7110

1.4109-2.9160-0.76691.7388lr_spue-2.5850

0.6732-3.9172-2.4185-1.7949

lnNins

9.96820.57859.31349.819611.1317v_spue0.17270.15780.03030.12850.4584v_ryo

0.13110.11450.04700.10460.2907v_yugai0.13600.13690.02900.09300.3803表４推定された自然増加率(ire)と、目撃効率の係数(rs)、狩猟捕獲率(pr[i])、有害捕獲率(py[i])
変数

平均

標準偏差

5%点

中央値

95 %点

ire

1.64650.2982

1.1522

1.66462.1041rs
pr0.09020.38400.04690.17510.01990.09270.0891
0.39390.16620.6452py

0.18560.08490.04620.18960.3170lnN201010.04960.56469.3401

9.9335

11.167148

兵庫ワイルドライフレポート 1 ,2012

自然増加率(log)

目撃係数(log)

有害係数(logit)

狩猟係数(logit)

個体数 t-1(log)図１

パラメータの事前分布と事後分布との関係

左上図

自然増加率

右上図

生息密度と目撃効率の比率を示す係数の自然対数値

左中図

狩猟による捕獲率(ロジット変換値)

右中図

有害による捕獲率(ロジット変換値)

左下図

1 年前(2009)年の生息数個体数の自然対数値

実線は事後分布を破線は事前分布をそれぞれ示す。

49

兵庫ワイルドライフレポート 1 ,2012

目撃効率

狩猟捕獲数

図２

観測値と期待値との関係

上図

目撃効率の観測値と期待値

下図

狩猟捕獲の観測値と期待値

中央値と 50%信頼限界、90％信頼限界を示す。

50

兵庫ワイルドライフレポート 1 ,2012

有害捕獲数

図３

観測値と期待値との関係
有害捕獲の観測値と期待値

中央値と 50%信頼限界、90％信頼限界を示す。

51

兵庫ワイルドライフレポート 1 ,2012

図４

兵庫県のイノシシの推定生息個体数の動向
中央値と 50%信頼限界、90％信頼限界を示す。

図５

兵庫県のイノシシの推定最大生息個体数の動向
中央値と 50%信頼限界、90％信頼限界を示す。

52

兵庫ワイルドライフレポート 1 ,2012

図６

兵庫県のイノシシの推定増加個体数の動向
中央値と 50%信頼限界、90％信頼限界を示す。

53

兵庫ワイルドライフレポート 1 ,2012

表５ 推定された生息個体数 N[i]、最大生息個体数 Nmax[i]、増加個体数 inc[i]
中央値

95％点

N2002

変数

平均18092.3

標準偏差25550.8

5％点
6016.210497.453296.5

N2003
N2004

19273.2
20184.7

25816.526021.26710.1
7252.511509.1
12299.2

55372.957120.1

N2005

19984.2

26201.1

6788.211948.5

57574.2

N2006

20906.626563.77123.4

12673.5

59827.7

N200722826.5

26986.9

8329.414514.3

63267.7

N2008
N2009

25157.5
26586.027341.027590.5

10099.2
11086.5

16897.4
18390.566760.5
68299.0

N201028409.2

28007.2

11385.6

20609.1

70761.9

Nmax2002
Nmax2003

25222.2
27043.7

25816.526021.2

12659.1
14111.5

17458.1
19158.2

61321.963979.1

Nmax2004
Nmax2005

28455.2
28061.6

26201.1
26563.7

15259.2
14278.420419.5
19828.566045.2
66982.7

Nmax2006
Nmax2007

29456.532493.5

26986.927341.014959.4
17435.2

21144.3
24233.4

69897.7
74096.5

Nmax2008
Nmax2009

36288.038752.2

27590.5
28007.2

20788.521728.6

28092.530952.1

78001.081104.9

Nmax2010

42027.828951.820878.035236.5

84918.3

inc2002inc20037130.07770.5

888.3
937.6

6628.47378.16927.07567.48366.89009.0inc2004
inc2005

8270.58077.51041.21260.17658.7
7356.28074.87783.09606.1
9939.3inc2006
inc20078549.99667.01486.1
1772.37482.07732.68224.69430.710799.7
12289.7

inc2008
inc2009

11130.5
12166.2

2298.03191.17817.17396.711069.9
12092.7

14413.017098.3

inc201013618.64828.46587.813328.521765.8

4. 考察
イノシシは、平均して毎年 4～5 頭出産するため、個体数の年内変動も年次変動も大きい。
本論文では、モデルの中で想定される年内の最大個体数を推定したが、これは、出産率から自然死亡率を差し引いた自然増加率から捕獲数を差し引く前の数値である。実際の出産時期直後で自然死亡の起こる前には、さらに多くの頭数が想定される。
このように、イノシシは出産率や死亡率が高いという特徴があるため、多数生まれた幼獣を捕獲して間引いても、最終的に生き残る幼獣の頭数は捕獲しなかった場合と変わらないといった、補償作用が働いている可能性も高い（Sandrocock et.al. 2010）。また、自然増加率の年次変動も激しいことが想定される。さらに、ニホンジカの推定で用いた糞隗密度のような、捕獲以外の情報がない。このため、現時点では、今回のモデルと推定結果をもとに、捕獲頭数を設定して予測シミュレーションは困難だと考える。
また、今回の推定値の事前分布と事後分布を比較すると、狩猟や有害捕獲の捕獲率の事前分布の設定が推定の制限になっていることが示唆される。今後の経過をモニタリングして、必要に応じてモデルの修正や事前分布の設定を検討する必要がある。また、イノシシではニ54

兵庫ワイルドライフレポート 1 ,2012

ホンジカの糞隗密度に対応するような密度指標がないため、推定を捕獲情報だけに頼らざるを得ないのが現状である。適切な密度指標の開発が望まれる。
このような現状の制限はあるが、イノシシの捕獲と個体群動態の傾向について、推定の結果から一定の判断をすることは可能である。データを収集し始めてから９年間の平均的な自然増加率や推定個体数と、現在の捕獲規模や捕獲効率を考え合わせると、毎年の捕獲数にばらつきはあるが、2002 年から 2010 年までの捕獲数は平均して自然増加頭数以下で、生息個体数は増加傾向にあると判断できる。モニタリングを行いながら順応的管理を行う限りでは、短期間でイノシシに絶滅の危惧が生じる状態ではないと考えられる。

謝辞
本研究の一部は、環境省の環境研究総合推進費（D-1003）により実施された。

引用文献
Geweke J 1992 Evaluating the Accuracy of Sampling-Based Approaches to theCalculation of Posterior Moments. In Bayesian Statistics 4 (Bernardo JM,BergerJO, Dawid AP, Smith AFM, eds), pp.169-193, Oxford Univ Press, Oxford.
環境省自然環境局生物多様性センター 2011 平成 22 年度自然環境保全基礎調査特定哺乳類
生息状況調査及び調査体制構築検討業務報告書. 411pp.
Kass RE, Carlin BP, Gelman A, Neal R 1998 Markov Chain Monte Carlo in Practice: A
Roundtable Discussion. The American Statistician 52:93–100.
Roberts GO, Gelman A, Gilks WR 1997 Weak convergence and optimal scaling of random
walk Metropolis algorithms. Annals of Applied Probability 7:110-120.
Sandercock BK, Nilsen EB, Broseth H, Pedersen HC 2010 Is hunting mortality additive
or compensatory to natural mortality? Effects of experimental harvest on the survival
and cause-specific mortality of willow ptarmigan. Journal of Animal Ecology
80:244-258.
SAS Institute Inc. 2011 SAS/STAT® 9.3 User’s Guide. SAS Institute Inc., Cary, NC.

55