μ t

2014年3月16日
第61回日本生態学会大会
企画集会T13 生態学における状態空間モデルの利用
状態空間モデルによる推定の考え方
農業環境技術研究所，生物多様性研究領域
山村光司
具体的なイメージを与えるため，
ここでは害虫を念頭に置いて解説（後の講演ではシカとクマ）
ニカメイガ
ツマグロヨコバイ
茨城県病害虫防除所HPより
ウンカヨコバイ図鑑HPより
1/13
ライトトラップによる年間誘殺数（5～9月合計値，水戸）
4
6
常 3.5
用
3
対
2.5
数
2
個
体 1.5
数
1
5.5
0.5
0
1950
5
4.5
4
3.5
3
ニカメイガ
1960
1970
2.5
1980
1990
2000
ツマグロヨコバイ
2
1950
1960
1970
1980
1990
2000
状態空間モデルの使用目的
（個体数変動に関する問題の場合）
• 個体数の推定
• 個体数の変動要因の推定
• 個体数を決定する環境要因の推定
今のデータの場合はこれが主目的：
将来に地球温暖化で害虫が多発
→気温の影響を推定
三重県病害虫防除所HPより
2/13
個体数変動をモデル化
• 第t年の個体数（対数値）をDt
• 第t年の個体数に伴う変動成分をet
et−1
1
Dt−1
et
1
Dt
et+1
1
Dt+1
（DtはDt−1の影響を受けて決まり，そこに変動成分が加わる）
• しかし，私らの観測しているDtは真の個体数そのものではない。
→ 実際の図式はもっと複雑。
3/13
• 第t年の真の個体数（対数値）をµt
• 第t年の真の個体数に伴う変動成分をet
• 第t年の観測個体数（対数値）をDt
• 第t年の観測個体数に伴う観測誤差をεt
et−1
et
1
μt−1
1
μt
1
Dt−1
1
εt−1
1
Dt
実際のモデル
はもっと複雑！
状態空間モデル
et+1
1
μt+1
1
Dt+1
1
εt
状態方程式
（プロセスモデル）
（システムモデル）
（遷移方程式）
1
観測方程式
（データモデル）
（測定方程式）
εt+1
（□は観測可能な変量，〇は観測できない変量）
4/13
複数の並行する時系列がある場合（複雑化その１）
• 今の場合は２種（ニカメイガとツマグロヨコバイ）が時系列が並行
• 複数の場所の並行的な時系列を扱う場合も同様
矢印を省略
et−1
μ't−1
状態方程式
観測方程式
1
μt−1
1
Dt−1
• A種個体数 µt
• B種個体数 µ't
1
εt−1
et+1
et
μ't
1
μt
1
Dt
1
εt
μ't+1
1
μt+1
1
Dt+1
1
εt+1
5/13
複数のタイプの観測値がある場合（複雑化その２）
• 今の場合は１種類のトラップデータだけ。
• 複数のタイプのデータがある場合はそれらをすべて活用できる。
状態方程式
観測方程式
et−1
et
1
1
μt−1
1
• トラップ捕獲数 Dt
εt−1
• その観測誤差 εt
• 捕虫網すくい取り捕獲数 Gt
• その観測誤差 ψt
1
μt
1
Dt−1
et+1
μt+1
1
1
Gt−1
1
ψt−1
Dt
Gt
1
1
εt
ψt
Dt+1
1
εt+1
Gt+1
1
ψt+1
6/13
環境変数（外生変数）がある場合（複雑化その３）
• 今の場合は温暖化を問題にしているので気温を組み込む
et−1
Tt−1
状態方程式
観測方程式
• 第t年の気温Tt
Tt は「外生変数」
1
μt−1
1
Dt−1
1
εt−1
et+1
et
Tt
1
μt
1
Dt
1
εt
Tt+1
1
μt+1
1
Dt+1
1
εt+1
（「観測時間」等の外生変数の場合は，状態方程式内の矢印は論理的に不要。）
7/13
どのように推定するか
• 実際にはソフトウエアで計算。中では何が行われているのか？
• [最尤推定法] データが生じる確率（尤度）
et
を計算して，その尤度が最大になるパラ
1
メーターを決める。
状態
方程式
• 左図の場合はデータはDtだけ。
µt
仮に Dt =5 のとき，さまざまな可能性。
1
μt =3, εt =2 や μt =4, εt =1 など…。
観測
「データが生じる確率」を計算するにはそ
Dt
方程式
れらの可能性をすべて積分する必要。
1
• もしet と εt の両方が正規分布に従う線形
モデルの場合には，積分結果は数式で求
εt
まり，それは正規分布となる。
→カルマンフィルターと呼ばれる計算へ。
• 積分が簡単に求まらない場合にはMCMCベイズ推定などへ。
8/13
具体的な計算例：まず時代の影響を除去する (害虫の場合は重要)
4
6
3.5
5.5
常
3
用
対 2.5
2
数
個 1.5
体
1
数
0.5
（防除技術等が変化するので）
5
4.5
4
3.5
3
2.5
ニカメイガ
0
1950
1960
1970
1980
1990
2000
ツマグロヨコバイ
2
1950
1960
1970
1980
1990
2000
→ 平滑化線からの偏差Dtを計算 (標準偏差5年のLOWESS平滑化)
（農業技術の普及には±約5年かかるので）
平
滑
値
と
の
差
Dt
1.5
1.5
1
1
0.5
0.5
0
0
-0.5
-0.5
-1
-1.5
-1
ニカメイガ
1950
1960
1970
1980
1990
2000
-1.5
（ただし，本来は状態空間内で行うべき）
ツマグロヨコバイ
1950
1960
1970
1980
1990
2000
9/13
モデルを図ではなく式で表現すると…
（ただし時代成分は除去済み）
状態方程式（2年前の個体数の影響まで考慮，気温は冬季と夏季に分ける）
µt = a1 µt−1 + a2 µt−2 + a'1 µ't−1 + a'2 µ't−2 + bµW TWt + bµS TSt+ eµt
真の対数個
体数
2年前の個体数
の影響
夏季気温の影響
2年前の他種
個体数の影響
誤差変動
1年前の個体数の影響
1年前の他種個体数の影響
観測方程式
前冬季気温の影響
Dt = µt + c µ't + bDW TWt + bDS TSt + eDt
観測された
対数個体数
他種個体数
の影響
真の対数個体数
夏季気温の
影響
前冬季気温の影響
誤差変動
• ノイズを拾わないように，重要なパラメーターだけを選択する。
（モデル選択基準を使用：AIC, RD など）
10/13
AIC基準で選択されたモデル [Yamamura et al. (2006)]
1.5
平
滑
値
と
の
差
Dt
減衰振動型
1.5
1
1
0.5
0.5
0
0
-0.5
-0.5
-1
-1.5
ニカメイガ
1950
真の個体数
1960
ランダム変動型
ツマグロヨコバイ
-1
1970
1980
前年個体数
1990
2000
観測された個体数偏差
1950
２年前の個体数
µt = 0.50 µt−1 − 0.28 µt−2 + eμt
Dt = µt + 0.087 TWt
誤差変動
真の個体数
-1.5
前冬季気温
観測方程式の誤差変動は
採択されていない。
1960
1970
1980
誤差変動
1990
2000
前冬季気温
µt = eμt
Dt = µt + 0.70 µ't + 0.17 TWt
真の個体数
真のニカメイガ個体数
間接効果：未知の環境要因があって，ニカメ
イガの真の個体数の増加率とツマグロヨコバ
イの年内増加率の両方に寄与している。
11/13
将来の害虫発生量の予測
R を発生量の倍率とする（50年後の発生量/現在の発生量）
R<2
2≤R<3
3≤R<4
4≤R
R < 1.4
1.4 ≤ R < 1.6
1.6 ≤ R < 1.8
1.8 ≤ R
ニカメイガ
50年後に本州では1.6～1.8倍の増加
ツマグロヨコバイ
50年後に本州では3～4倍の増加
（気象庁・気象研究所のモデルによる気温予測値を用いた場合）
12/13
（まとめ）状態空間モデルの効能
 複数タイプの観測データをフル活用して，それらに伴う観測誤差
を除去し，個体数等を正しく推定する。
 「モデル選択」を行うと，個体数変動の仕組みを特定できる。
ニカメイガの場合 µt = 0.50 µt−1 − 0.28 µt−2 + et どちらの式が優勢か
は対象生物に依存
Dt = µt + 0.087 TWt
• 冬季気温の影響は観測方程式だけ。翌年への持ち越しは無い。
 これは1950年代の密度調節論争の意味を式で示している。
Nicholson (1954)の主張：生物の個体数は密度依存的に制御。
Andrewartha and Birch (1954)の主張：密度依存的な制御は意
味をもたず，その年の環境条件が好適かどうかで個体数が決定。
 状態方程式（ µt の式）： Nicholson (1954)が見ていた視点
観測方程式（ Dt の式）： Andrewartha and Birch (1954) の視点
 状態空間モデルの効能：観測誤差を除去できるだけでなく，
個体群動態の「本流部分」と「派生部分」を正しく区別して議論
を行うことが可能となる。
13/13

Download Report