統計数理（2018）
第 66 巻 第 2 号 319–337©2018 統計数理研究所［研究ノート］

トータルパワー寄与率を用いた
海洋生態システムにおける因果性推測
ソルヴァン加藤 比呂子1 ・Subbey Sam2,3（受付 2017 年 3 月 16 日；改訂 2018 年 7 月 30 日；採択 8 月 2 日）

要

旨

多変量時系列データの変量間相互関係を推測する方法として，Ozaki（2012）は，Granger（1969）と Geweke（1982）のペアワイズ因果性推測（Partial pairwise causality）と赤池（1968）のを統合した因果性推測を提案した．我々は赤池情提案した相対パワー寄与率（Total causality）
報量規準の枠組みで，因果関係の有意性に関する規準を加え，シミュレーションデータによりその手法の検証をおこなった．また，実データ分析として，バーレンツ海域の生態系におけで重要な 4 種類の海洋生物（シシャモ，タラ，オキアミ，ニシン）のバイる食物連鎖（food web）
オマス時系列データを用い，シシャモの年齢毎に生物種間のフィードバック関係を考察した．
2–3 年齢のシシャモを含むフィードバックシステムは，1 年齢もしくは 4 年齢のシシャモを含むフィードバックシステムよりもより多くの生物種間の相互関係を示した．それらはバーレン
ツ海の food web に関するこれまでの先行研究を裏付け，シシャモがバーレンツ海域生態シス
テムの food web に関連する生物種間の重要な駆動源になることを明らかにした．本稿が提案する手法は，海洋学研究で対象となる複雑な生態系において，生物種間，環境因数間の因果関係を推測する一手段として有用であると考えられる．
キーワード：多変量自己回帰モデル，多変量時系列データ，フィードバックシステム，
Granger の因果性，海洋生態システム，バーレンツ海．

はじめに

1.

近年の地球温暖化に伴う海洋生態系と海洋生物種の資源量の変化，生物種と環境要因との因果関係を予測することは，学術的・商業的に重要である．バーレンツ海域では，その漁業管理のために，ロシアとノルウェーの共同評価（WGIBAR, 2016）がおこなわれている．その調査をと海洋生物種（Biotic）に関する高次元の時系列データが観測され通じて，環境的要因（Abiotic）ている．しかし，それらのデータ分析手法としては，時系列モデルに基づくものは殆どみられ（ICES）でスタンダードとず，その多くは，International Council for the Exploration of the Sea
1Marine Mammals Research Group, Institute of Marine Research, P.O.Box 1870, Nordness N-5817 Bergen,

Norway2

Research Group on Fisheries Dynamics, Institute of Marine Research, P.O.Box 1870, Nordness N-5817,

Bergen, Norway
3

USA

Department of Natural Resources, Cornell University, Room120 Fernow Hall, Ithaca, New York 14853,

320

統計数理

第 66 巻 第 2 号 2018

推奨されている主成分分析や対応分析が用いられている（WGIBAR, 2016）．時系列データにおける主成分分析は，時間領域か周波数領域に適用するかで，観測データを独立同一分布として扱えるかどうかが変わり，それらの結果の解釈は全く異なることも気をつけなければならない（川崎, 2001）．ICES で推奨している通常の主成分分析は，サンプル間を独立と仮定しているため，理論的には時間軸上のデータにそのまま適用するにはふさわしくなく，まして変量間の因果関係を推察するには適切ではない．Planque and Arneberg（2017）は，主成分分析による
Abiotic と Biotic を含む高次元時系列データの分析結果の信憑性が極めて薄いことを指摘している．
多次元時系列間の因果性推測については，1956 年の Wiener の論文にさかのぼる（Wiener,．2 次元の時系列データについて，1 つの時系列に対しもう一方の時系列からの関係を考1956）慮した方が最小二乗値が小さい場合，2 変量間に因果関係が存在するということを提唱した．
その後 Granger がその考えを経済時系列データ間の因果性分析に適用した（Granger, 1969）．
Granger の因果性は基本的に時間軸領域における分析で，2 変量自己回帰モデルを適用しその予測誤差の比較により推測をおこなっている．一方で，1968 年に，赤池は秩父セメントのロー
タリーキルンの駆動制御のために，多変量自己回帰モデルを用いたフィードバックシステム解析を考案した（Akaike, 1968）．この手法はシステムを撹乱すると考えられるノイズに注目している．実データ分析では有色ノイズであるが，モデルの導出過程で白色化を施し，最小二乗法
でモデルの係数の推定を可能にした．係数と予測誤差の共分散行列を用い，時間領域では閉
ループと開ループのインパルス応答を，周波数領域では相対パワー寄与率を計算することで，双方の領域から変数間の因果関係を推察できる分析手法であった．周波数領域で提案された相対パワー寄与率は，長/中/短周期において，ある変量のパワーに対し各変量からのノイズ（パワー）の寄与がどの程度の割合か算出することができる．この手法は分析対象となるシステムに関連する多次元時系列データをまとめて取り扱えるため，医学，工学，経済学，生物学等
様々な分野での適用例が多数報告されている（Akaike and Kitagawa, 1994）．
多次元時系列データ間の相互関係を分析するには，このような多変量時系列モデルを用いた因果性分析の方が，主成分分析による変動の類似性に基づく手法よりも，変量から変量への方向性を踏まえた関係を推測するのに有用である．時系列間の相互関係に対し多変量自己回帰モデルを適用して，多変量間のコヒーレンスを用いた手法（Baccalá and Sameshima, 2001）も提案されているが，それらは多変量自己回帰モデルの係数のみ，すなわち伝達特性のみに注目を対象としたものである．しかし，実際にはシステムを構成するサブシステムを撹乱するノイズ項を考慮することは大切で，システムのダイナミクスを左右するイノベーションと言われている（Ozaki, 2012）．多変量間のフィードバックを想定した場合，伝達特性とイノベーションの寄与双方から考察することが時系列間の因果性を推測する上で重要である（Bosch-Bayard et．
al., 2012）
Ozaki（2012）の 14 章では，赤池の相対パワー寄与率を全周波数領域で積分し（Ozaki（2012）と呼んでいる），関連した全変量を対象としたパワーの寄与からあで「赤池の Total causality」
る変数の関係を外した寄与について定式化している．全変量を対象としたパワーの線形和に対し，周波数領域において総和をとると，変動の分散そのものになり，その対数表示はまさにモデルの対数尤度を示す．さらに，Granger-Geweketype causality の，2 変量間因果性の有無に関して予測誤差の対数を比較する手法（Ozaki（2012）で「Partial pairwise causality」と呼んでいる）
を取り上げ，赤池の Total causality の枠組みで計算し，ある変量と全変量とのペアワイズに比較する手法を提案した．これは，Granger が 2 変量のみの時系列モデルによる予測誤差を頼りにするだけでは多変量全ての間の因果推測をするには不足であるということと，赤池が全変量を対象としているが，複雑なフィードバックシステムにおける因果性推測には，ペアワイズでトータルパワー寄与率を用いた海洋生態システムにおける因果性推測

321の比較が予備的に因果関係を探す手段として役立つという，双方の困難な点を補える利点がある．データの背後に潜む本当の因果関係は，限られたサンプル内のデータで予備的に探られた結果と，データに関する知見や背景の情報と合わせて推測されるだろう．この Ozaki（2012）の提案手法は，実データ分析においてシステマティックで実用的な方法と期待されるが，未だ具体的な適用例が少なく，複雑な海洋生態システムのダイナミクスの解明においては皆無である．
多次元時系列データ間に内在する因果性を分析する手法が取り上げる問題は，因果関係の有無，因果関係の強度，因果関係の可視化，因果関係による機能の分析等が考えられる．先にあ，Geweke（1982）
はペアワイズのモデル化に基づいて，また，Akaike は関げた，Granger（1969）
連した全ての変量の関係を含むモデルで因果関係の有無を調べる手法と考えられる．Granger（1969）のコヒーレンス関数や Akaike の相対パワー寄与率やインパルス応答関数のグラフ化は因果関係を可視化と言えよう．因果関係による機能の分析では，赤池・中川（1972）のフィード
バックシステムの動的制御や，加藤・石黒（1997）によるフィードバック経路切断による経済シ
ステムの動的解析があげられる．
の手法に，さらに因果関係の有無を査定する規準を情報量規準の本論文では，Ozaki（2012）
枠組みから導入し，海洋生態システムの Abiotic や Biotic の因果関係の有無や強度，可視化の可能性を提案する．まずシミュレーションデータによる検証を行い，次に実データ分析として，バーレンツ海生態システムの food web において重要な 4 種類の海洋生物のバイオマスを用い，生物種間の因果関係を推測する．
ノイズ寄与率は，赤池の他に inverse スペクトルを適用したパワー寄与率を考慮した．実デー
タ分析として，バーレンツ海生態システムにおける food web に関する生物資源のデータに適用し，そのシステムにおける生物種間の因果関係を推測した．
2.

統計的手法

観測された k 次元時系列データ xt = (x1 (t), x2 (t), . . . , xk (t)) , t = 1, . . . , N とする（ここで，(·) は転置の記号）．これらのデータは，以下に示すような多変量自己回帰（Multivariate
過程の実現値と仮定する：
auto-regressive，MAR）
(2.1)

xt =

M


Am xt−m + εt ,m=1

ここで M は自己回帰の次数，Am は自己回帰モデルの係数，そして εt は平均ゼロベクトル，
分散共分散行列 Σ に従う多変量正規分布に従うとする．自己回帰係数の推定には多数の方法があり，最小二乗法や Yule-Walker 法，また，Yule-Walker 法を効率よく計算できる Levinson
アルゴリズム等が提案されている（Ozaki, 2012, 4 章）．自己回帰の最適な次数 M は Akaike（AIC）（Akaike, 1974）
等の情報量規準による統計的モデル選択で同定でき，Information Criteria
予測誤差系列により共分散行列が計算できる．本稿におけるモデル選択では，AIC（AIC = −2 ×モデルの最大対数尤度 + 2 × モデルの自由パラメータ数）を用いる．
2.1 赤池の相対パワー寄与率
自己回帰係数とフーリエ変換から，周波数 f に対する周波数応答関数 Ff が求まり，以下のようなパワースペクトルが求まる：

統計数理

322

⎛
(2.2)

p11f
⎜ p21f⎜Pf = Ff ΣFf∗ = ⎜

⎝

pk1f

第 66 巻 第 2 号 2018

p12f
p22f
..
.
pk2f

⎞

···
···

p1kf
p2kf ⎟
⎟

···

pkkf

⎟,
⎠

0 ≤ f ≤ 0.5Δ,

ここで Ff∗ は Ff の共役複素数，Δ は観測値のサンプリング間隔とする．Pf の非対角成分はク
ロスパワースペクトルである．もし，Σ の非対角成分が非常に小さい，すなわち各変量の予測
誤差は独立であることが仮定できると，i 番目の変量 xi のパワースペクトルは他の変量 xj か2の影響を含む項の和
らの周波数応答関数 Fijf と予測誤差の分散 σjj
222
piif = |Fi1f |2 σ11+ · · · + |Fiif |2 σii+ · · · + |Fikf |2 σkk(2.3)の形で示すことができ，j 番目の変量 xj からのノイズの影響 rijf を以下のように示すことが
できる：
(2.4)

rijf =

2|Fijf |2 σjj∈ [0, 1] .
|piif |

（2.4）式で示された比を赤池の相対パワー寄与率と呼ぶ．予測誤差を独立としない場合の相対パワー寄与率は Tanokura and Kitagawa（2014）により提案されている．
Ozaki（2012）
では，多次元時系列データに MAR モデルを適用しその推定値と予測誤差の情報を用いて（2.4）を算出するので，全時系列データの全ての因果関係を含んだ情報（total causality）
を用いて（2.4）式を算出したものを，partial innovation contribution と呼んでいる．それにより（partial causality）と呼んでいる．
推測される結果を Partial 因果性
2.2 Granger と Geweke 型の因果性観測された 2 つの時系列 xt と yt に対して時系列モデルを考えることにする．モデル A を適用した際の時系列の予測誤差の分散を Var(∗|A) と表すことにすると，Granger と Geweke の因：果性は以下のように定義する（Ozaki, 2012）
xt− = (xt−1 , xt−2 , . . .) ，yt− = (yt−1 , yt−2 , . . .) と す る と ，
定 義 1．（Granger の 因 果 性 ）Var(xt |xt− ) − Var(xt |xt− , yt− ) > 0 ならば，時系列 xt の変動はもう一方の時系列 yt の変動に起因している．

定義 2．（Geweke の因果性）
log |Var(xt |xt− )| − log |Var(xt |xt− , yt− )| > 0 ならば xt の変動は
yt の変動に起因している．
ここで，Var(xt |xt− ) は xt− によるモデルで予測された xt の予測誤差の分散を，Var(xt |xt− , yt− )は x1,t− と yt− を含むモデルで予測された yt の予測誤差の分散を示す．定義 2 で示された式，
左辺の第一項で示した予測誤差を算出する予測モデルを Model(0) とし，同じく左辺第二項に示した予測誤差を算出する予測モデルを Model(1) とすると，Ozaki（2012）では，Granger とGeweke の定義は，本質的に以下のような予測モデルの尤度の比較に相当すると示している：
定義 3．ある観測された時系列 xt の変動が時系列 yt に起因されているということは，モデ(0)(1)
ルの対数尤度の観点から，−2 log lModel (xt ) − (−2) log lModel (xt ) > 0 となる．
(0)

ここで，lModel(1)(xt ) と lModel(xt ) は Model(0) と Model(1) で予測した xt の尤度を示す．

トータルパワー寄与率を用いた海洋生態システムにおける因果性推測

323

2.3 Ozaki（2012）の因果性推測は，予測誤差 σ 2 とパワースペクトル p(f ) との関係をさて，Kolmogorov（1941）
1/2

σ 2 = exp

(2.5)

log p(f )df−1/2と示した．その周波数領域において積分したものをトータルパワーの全周波数領域の総和と言（2.3）式でパワースうことにしよう．多次元時系列データに MAR モデルを適用する話に戻り，ペクトルが求められたとすると，今 i 番目の変数のトータルパワーの全周波数領域の総和は
1/2

(2.6)

log σi2 =1/2

log pii (f )df =
−1/2

log
−1/2

K


|αik (f )|2 σk2 df

k=1

1/2

=
−1/2

2
log(|Fi1 (f )|2 σ12 + |Fi2 (f )|2 σ22 + · · · + |FiK (f )|2 σK)dfと示すことができる．また，j 番目の変数からの影響を除いた場合のトータルパワーの全周波
式及び
（2.6）
式に示されている全変数の周波数応答とノイズの分散の積
数領域の総和は，（2.3）の線形和から，|Fij (f )|2 σj2 を除いた式を計算することになる．例えば j = 2 とすると，
1/2

(2.7)

2
log σi∧2
=

−1/2

(j)

log pii (f )df

1/2

=

log
−1/2

j−1


|αik (f )|2 σk2 +

k=1

k


|αik (f )|2 σk2

df

k=j+1

1/2

=
−1/2

2
log(|αi1 (f )|2 σ12 + |αi3 (f )|2 σ32 + · · · + |αiK (f )|2 σK)dfとなる．これらを以下のように比較すると，
1/2

(2.8)

2
2
log σi∧j− log σii
=

−1/2

(j)

log pii (f )df −

1/2

(j)

log

=
−1/2



1/2

log

=
−1/2

pii (f )df =pii (f )1/2

log pii (f )df
−1/2
1/2

log
−1/2

2
|αij (f )| σjj
1−
pii (f )

2pii (f ) − |αij (f )| σjj
df
pii (f )dfとなり，定義 3 で示したような，ある変数からの影響を含むか含まないかの 2 つのモデルの比式の最終項にあるように，全変数からある 1 つの変数の関較と等価になる．また，それは（2.7）（全周波数領域の総和をとった）をみることになる．本稿
係を外した場合の Partial パワー寄与
では，（2.7）式で得られる値を対数尤度の差の値ということで，本稿では ΔLL 値として示すことにする．
ΔLL 値がどの値になると変数 j から変数 i に関係していると判断するか．本稿では，モデルの比較として AIC の観点からみることにする．変数 j から変数 i の関係を外した（2.6）のように示した擬似（対数）
尤度とトータル寄与の擬似（対数）尤度をもつモデル，すなわち 1 変数だけ外したモデルとの比較になるので，

統計数理

324

(2.9)

AICσ2

i∧j

第 66 巻 第 2 号 2018

2
2
− AICσ2 = −2 × log σi∧j
+ 2 × (k − 1) − −2 × log σii
+2×k



ii

= ΔLL − 2となり，（2.7）式の値が 2 より小さい場合は差があるとは言い難くなる．
モデルの対数尤度の比較という点からは，ΔLL を尤度比検定の枠組みに拡張することもできるかもしれない．その場合帰無仮説は変数 j から変数 i に関係はない，すなわち |Fij (f )|2 σj2 = 0ということになる．尤度比検定量を調べるには，リサンプリングしたデータに基づいてモデルを適用した後に ΔLL を算出し，それを何回か繰り返すことによって ΔLL の帰無分布を作成しなければならない．しかし帰無仮説が |Fij (f )|2 σj2 = 0 となるには，MAR モデルに含まれる全変量の自己回帰係数が関係していて，予測誤差の共分散については，実データに適用したときに独立と仮定できてもわずかの非対角要素も存在するので，帰無分布を求めるための厳密な性質を踏まえたデータをリサンプリングすることは難しい．その点に関しては，理論的な背景を踏まえ別稿で議論することにする．従って，本稿では，因果性を調べたい対象となるある変量（2.9）をもとにこの Total パワーに対し，関係を調べたい変量とのペアごとに ΔLL を計算し，の値が 2 から 3 以上の場合，有意な因果関係を推測できるものとする．
2.4 Bosch-Bayard の pNCR
赤池のパワー寄与率と似たアプローチで，Bosch-Bayard により partial noise contribution（pNCR）が提案されている（Bosch-Bayard et al., 2012）．この手法は（2.2）式で示したパ
ratio
（2.4）
式と同様な比率を算出する．ある対象とする変数
ワースペクトルの逆行列 Pf−1 を用いてのパワースペクトルの逆行列に対し，他の変数からの寄与の割合を計算する．すなわち
b
rijf
=

(2.10)∗
2
Fijf )/σjj(Fijf|piif |−1と表される．赤池のパワー寄与率と同様，多次元時系列データに対して MAR モデルから推定された係数と予測誤差の共分散を用い，対象とする変数とそれへの影響を調べるための変数のbをプロットすることで因果関係の可視化が可能となる．本稿でペアで周波数領域において rijf
は，赤池の手法によるトータルパワーに加え，Bosch-Bayard のトータルパワーも検討することにする．
3.

シミュレーションデータ分析

ここでは，2 変量のみ対象とする Granger（1969）と Geweke（1982）の方法と，全変量の情報の ΔLL について，3，5 変量シミュレーションデータにより比較する．
を用いる Ozaki（2012）
また，ΔLL について，サンプル数に関する Sensitivity と，因果関係の強度との関係について検討する．
3.1 Granger と Geweke の因果性分析と ΔLL の実装の比較—3 変量データの場合
データは以下に示される多変量自己回帰モデルから生成された 200 時点の 3 変量（xt ，yt ，zt ，t = 1, . . . , 200）
時系列，sim1 と sim2 を用いる：
sim1: xt = 0.3xt−1 − 0.45yt−1 + εx,tyt = 0.4yt−1 + 0.35zt−1 + εy,tzt = 0.25zt−1 + εz,tsim2: xt = 0.3xt−1 − 0.45zt−1 + εx,tyt = 0.4yt−1 + 0.35zt−1 + εy,tzt = 0.25zt−1 + εz,tここで εx,t ，εy,t ，εz,t は平均 0 分散 1 の独立同一分布に従うものとする．生成されたデータのトータルパワー寄与率を用いた海洋生態システムにおける因果性推測325

図 1．MAR（1）model によって生成された 3 変量シミュレーションデータ sim1 と sim2 とその変量間の関係．

プロットと変量間に仮定しているの関係を示したダイアグラムを図 1 に示す．
Granger と Geweke による方法では，各変量に対して以下のようなモデルを適用する：
xt に対し，xt−1 からのみ説明されるモデル
xt−1 と yt−1 から説明されるモデル
xt−1 と zt−1 から説明されるモデル

xt = axx xt−1 + εxx,txt = axx xt−1 + axy yt−1 + εxy,txt = axx xt−1 + axz zt−1 + εxz,tyt に対し，yt−1 からのみ説明されるモデル
yt−1 と xt−1 から説明されるモデル
yt−1 と zt−1 から説明されるモデル

yt = ayy yt−1 + εyy,tyt = ayx xt−1 + ayy yt−1 + εyx,tyt = ayy yt−1 + ayz zt−1 + εyz,tzt に対し，zt−1 からのみ説明されるモデル
zt−1 と xt−1 から説明されるモデル
zt−1 と yt−1 から説明されるモデル

zt = azz zt−1 + εzz,tzt = azx xt−1 + azz zt−1 + εzx,tzt = azy yt−1 + azz zt−1 + εzy,ta∗∗ はモデルの係数で，最小二乗法から推定される．それぞれのモデルに対する予測誤差の分散を計算し，定義 1 に示したように分散間の差と，定義 2 に示したように分散の対数間の差を算出する．その結果を表 1 にまとめた．
まず予測誤差の分散の差はどの場合も関係が入らない変量との差が 0 より大きくなる．sim1
については，各変量ごとでそれぞれの誤差間の分散の差を比較すると，xt に対して yt−1 を含む

326

統計数理

第 66 巻 第 2 号 2018

表 1．sim1 と sim2 データにおける Granger と Geweke の方法を適用した結果．

モデルの方が zt−1 を含むモデルより誤差間の差が大きく，yt に対しては zt−1 から説明されるモデルの方が xt−1 を含むモデルより誤差間の差が大きくなるため，xt と yt の変動には，yt−1
からと zt−1 からの関わりがあると推察される．しかし，zt に関しては，他の変量からの関わりがあると言えるのか言えないのか不明である．誤差間の差の大小関係からだと xt−1 からの関わりが若干ありそうである．sim2 については，各変量毎にそれぞれの誤差間を比較するとかなり大きな差がみられるので，xt と yt に zt−1 から関わりがあると推察されるが，zt に xt−1と yt−1 からの関係はないはずだが，これらの比較では yt−1 からの関係があるかどうか判断が難しい．xt に yt−1 からの関係もありそうである．定義 1 と 2 に沿った誤差分散の差の比較を明解にするために，2 つのモデルの残差和（誤差の総和）を比較する F 検定が有用である．本稿（vars）にある関数 causality を実装してみる．各検定から得られた p 値を表 2
では R の library（pxy は変量 yt−1 から xt への関係に関する値）にまとめる．
5% 未満の p 値を有意な差であるとすると，表 2 では，概ね正解に近い関係を示せているとみられるが，sim1 に関しては，zt−1 から xt と xt−1 から zt への関係が正解と異なっていた．
さて，ΔLL による因果推測については，Akaike のパワー寄与率（RPC）と Bosch-Bayard の
pNCR を用いた場合の ΔLL を算出する．
（2.9）式で言及したことを考慮し，ΔLL > 2.5 を満たす変量からの関係は有意であるとする．図 2 に ΔLL の値を変量毎にプロットした図と，その値から推察されるダイアグラムを示す．それぞれのグラフのタイトルに示した変量に対し，x軸に示した変量を外した場合の ΔLL 値を × でプロットし，ΔLL > 2.5 の場合は○をつけていトータルパワー寄与率を用いた海洋生態システムにおける因果性推測327

表 2．sim1 と sim2 データにおける Granger 因果性の F 検定の結果．

る．その場合は，外された変量はタイトルの変量に寄与していることになる．どの図も，自分自身からの影響に関しては非常に大きいため，それが外されると y 軸で示している範囲を大幅に超えているのでプロットされていない．RPC も pNCR による場合も，ΔLL はシミュレー
ションデータを生成するときに用いた関係と同じ関係を導き出している．
3.2 Granger と Geweke の因果性分析と ΔLL の実装の比較—5 変量データの場合
3.1 節の例と比べ，もう少し複雑な関係を含む 5 変量のシミュレーションについても実験をおこなった．以下がデータを生成するモデルである：
x1,t = 0.7x1,t−1 + 0.5x2,t−1 + 0.6x4,t−1 + ε1,tx2,t = 0.7x2,t−1 + ε2,tx3,t = 0.7x3,t−1 + 0.6x5,t−1 + ε3,tx4,t = 0.7x4,t−1 + 0.6x5,t−1 + ε4,tx5,t = 0.7x5,t−1 + ε5,t .

生成されたデータ 200 点について，上記と同様の Granger 因果性 F 検定をおこなうと，以下のような p 値が得られる：

⎛

p12
⎜ p21
⎜
⎜ p31
⎜
⎝ p41
p51

p13
p23
p32
p42
p52

p14
p24
p34
p43
p53

⎞

⎛

p15
2.0 × 10−4⎜
p25 ⎟⎟ ⎜ 0.91⎜=
p35 ⎟⎟ ⎜ 0.35⎠⎝ 0.0053
p45
p54
0.58

1.93 × 10−6
0.11
0.87
0.64
0.96

< 2.2 × 10−16
0.84
0.040
0.28
0.72

⎞

1.1 × 10−5⎟
0.97⎟
< 2.2 × 10−6 ⎟⎟< 2.2 × 10−16 ⎠0.69

変量数の増加に伴い検定回数も増えることによる多重検定問題を回避するために，5% の False（FDR）（Benjamini and Hockberg, 1995）を適用し有意となった関係を図 3 中央
Discovery Rateに示す．図 3 左に描かれた正解のダイアグラムと比較すると，3 箇所余計な関係が推定されていることがわかる．さて，このデータに対して ΔLL を算出してみると，図 3 右のような関係が描ける．5 → 4 → 1 を 5 → 1 と両方推定してしまう．このような indirect/direct な関係を最終的にどちらが正しいか見極めるには，それらの関係を含んだモデルを再度適用し，AIC で比較することによりふさわしい関係を判定することが可能である．例えばこの例の場合，5 → 4 → 1
を 5 → 1 を含んだモデルの AIC は 579.0（対数尤度 −264.5）
5 → 1 を外すと AIC は 528.8（対数尤度 −265.4）となるため，5 → 1 は冗長であると判断される．3.1，3.2 の結果から，ペアワイズに基づく Granger と Geweke の方法よりも，全変量の情報に基づいた ΔLL の方が，シミュレーションデータに仮定した関係をより正しく推測していた．
3.3 サンプル数の違いによる ΔLL の Sensitivity
ところで，これまでのシミュレーションデータは MAR モデルを適用するには十分なサンプ

328

統計数理

第 66 巻 第 2 号 2018

図 2．sim1，sim2 データにおける，Akaike の RPC と Bosch-Bayard の pNCA により算出された ΔLL 値と，それらの結果を元に作成した変数間のダイアグラム．矢印のわきには RPC と pNCA から得られた ΔLL 値を記してある．

ル数 200 点を有していた．そこで，サンプル数の変化に伴う ΔLL の Sensitivity を調べることにする．具体的には，3.2 節で生成した 5 変量データで，サンプル数を 50，100，150，200 点変
え，10000 回生成することにより，ΔLL がどの程度正解を得ることができるか調べてみた．変量間の関係の正解数と誤正解数の平均 ± 標準偏差を表 3 に示す．正解数に関しては，サンプル数が 50 点の段階から正解数に近く 150 点以上ではほぼ正解数を示しているが，50 点では誤正解である関係も推定しやすくなる．150 点以上でとった誤正解として推定した関係は，3.2 節にもみられた indirect/direct の関係であった．実データの観点から考えると，例えば地球温暖化に関連する海洋調査データでは，その殆どが，温暖化を意識され始めた年前後含めて 30 年から
多くて 50 年分程度の時系列データしか観測されていない場合が多い．また，生物種間によってトータルパワー寄与率を用いた海洋生態システムにおける因果性推測329

図 3．5 変量シミュレーションデータに仮定した関係（左），Granger 因果性 F 検定より推定
（右）．
された関係
（中），ΔLL により推定された関係
表 3．5 変量データ（正解数は 4）に関して，サンプル数を変化させた場合における，関係の正解数と誤正解数．

は，indirect/direct の関係のどちらかを選ぶ必要がないものもある．従って限られたデータ数の関係を推定するには，ΔLL によりシステマティックに推測された変量間の関係に加え，さらに生物学的知見とその生態系の専門知識，また関連した先行研究による裏づけが必要となる．
3.4 ΔLL による因果性の強度
先に述べたように，多次元時系列データ間に内在する因果性を分析する手法が取り上げる問
題として，因果関係の有無，因果関係の強度，因果関係の可視化等が考えられる．ここでは，
ΔLL と変量間の寄与の強度との関係について検討する．以下の数式を用いてデータを作成し，
xt = 0.3xt−1 − 0.15zt−1 + εx,tyt = 0.3yt−1 + ayz zt−1 + εy,tzt = 0.3zt−1 + εz,tayz を 0.1 から 0.9 に変化させた場合のΔLL 値を計算し，そのプロットを図 4 にまとめる．9 つ

図 4．z から y への MAR モデルの係数を 0.1 から 0.9 まで変化させた場合の，ΔLL 値の変化．x 軸は変化させた係数の値，y 軸は ΔLL 値．

統計数理

330

第 66 巻 第 2 号 2018の図のタイトルがどの変量からどの変量への寄与かを示す．例えば xx，yy ，zz は変量自身からの寄与，xy は y から x への寄与を示す．z から x への影響が xz の図に示されており，zz と比べると，xx は自身からの影響が低めに出ている．yz の図に注目すると，z から y への影響
は係数の値が大きくなると ΔLL 値も大きくなるのがわかる．それに伴い，y 自身は ΔLL 値が下がってきて，y 以外の変量からの係わり合いが含まれていく変化がみられる．ΔLL は対数スケールなので，yy と yz を線形に足し合わせて一定の値にはならない．他の変量の係数については，MAR で全極モデルとして係数を推定するために，若干の変動はみられるが，y 軸のスケールを合わせると，それらの変動は yz の変化よりも微量であることがわかる．これらの結果から，多次元データが，線形で定常であることが仮定でき MAR モデルを適用した場合，
ΔLL は変量からの寄与の強さを示すことができると示唆される．
4.

実データ分析

バーレンツ海における生態系では，シシャモの資源量の維持が他の生物種が関連する food
（食物連鎖）において重要な役割を示している．シシャモは 3 から 4 年程度の寿命の遠洋魚
webで，バーレンツ海における北東北極タラとニシンの稚魚の餌食として知られている（Gjøsæter
．シシャモの消費に伴う資源量の衰退はバーレン
et al., 2009; Hallfredsson and Pedersen, 2009）
ツ海の food web のダイナミクスに大きく影響を与えると考えられている．例えば，Hjermann
et al.（2004）はバーレンツ海の food web を，シシャモを中心に栄養段階の高低レベル双方に．
連鎖しているという仮説を示した（図 5）
本稿では，バーレンツ海におけるシシャモ，タラ，ニシン，オキアミ間の関係に対し，フィー
ドバックシステムを仮定し，ΔLL によりその因果関係を分析する．時系列データは Working（WGIBAR）（ICES, 2016）で取り上げGroup on the Integrated Assessments of the Barensts Sea
1 から 4 年齢の年ている，1972 年から 2014 年まで収録された，シシャモ（Capelin または Cap）とニシン（Herring）の全年間バイオマス，オキアミ（Krill）はバーレン
間バイオマス，タラ（Cod）
ツ海のバイオマス密度を用いる（図 6）．モデルを適用する際，データは標準化されシシャモの各年齢におけるデータとタラ，オキアミ，ニシンを合わせて 4 次元の時系列データとして取り扱うことにする．
まず MAR モデルを適用する前に，これらのデータを標準化し相互相関を調べる．図 7 にシ
シャモ 1 年齢と他の 3 種類の生物種との場合の相互相関のプロットを示す．シシャモの変動に対する他の変量との相互相関は，ラグが低いところにピークがみられるが，タラやオキアミに関する相互相関ではラグの高いところにもピークがみられる．本稿の図には示していないが，他の年齢においても似たような傾向がみられた．図 6 にみられるように，タラやニシン，オキ
アミはシシャモの変動と比べると若干非定常性や非線形性がみられるため，長いラグに対して

図 5．シシャモを中心としたバーレンツ海の食物連鎖に対する仮説（Hjermann et al.（2004）
から図を再描画した）．

トータルパワー寄与率を用いた海洋生態システムにおける因果性推測

331

図 6．シシャモ 1–4 年齢のバイオマス（y 軸 ×106 kg），タラ，オキアミのバイオマス（y 軸kg），ニシンのバイオマス密度（y 軸 × 103 g/m2 ），x 軸は年を示す．

図 7．シシャモ年齢 1，タラ，ニシン，オキアミ 4 系列に関する相互相関（同じ種同士は自己相関）関数のプロット．各プロットのタイトルはどの変量間の相互相関かを示すもので，
例えば c12 はシシャモとタラに関する相互相関を示す．下付の数字は図全体のタイトルに示した通り．x 軸に −10 から 10 までのラグを，y 軸に相互相関関数値を示す．

相互関係のピークが現れているかもしれない．
MAR モデルを適用した結果，AIC により，シシャモ 1 と 4 年齢の場合は MAR（1）モデル，モデルが同定された（Solvang et al., 2017 の Table2 をシシャモ 2 と 3 年齢の場合は MAR（2）

332

統計数理

第 66 巻 第 2 号 2018

図 8．ΔLL から推察される 4 生物種（シシャモ Capelin，タラ Cod，オキアミ Krill，ニシン
Herring を表す）間の関係について（右ダイアグラム，矢印横数字は ΔLL 値を示す）．

参照）．推定値より計算した ΔLL を用い，ΔLL > 2.5 となる関係を実線で，ΔLL が 2.4 から
2.5 となる関係を点線で示したダイアグラムを図 8 にまとめる．シシャモ 1 年齢よりも 2，3 年
齢の方が，4 生物種間の関係が多様に現れ，シシャモが 4 年齢になるとそれらの相互関係が希薄になる．これは，4 年齢のシシャモが，実データに見られるように，バーレンツ海において
ドミナントではないことによると考えられる．また，シシャモ 1 から 3 年齢の変動はタラの変動に高く影響を及ぼすことを示していて，これは図 5 の仮説で示されているシシャモがタラの餌食となることに対応すると考えられる．シシャモ 2 年齢は最も高くオキアミで代表される動物性プランクトンの変動に影響を与えることを示している．これは，シシャモとオキアミが強（例えば Baum and Worm, 2009）を支持している．シシャモ 3 年
い top-down の関係をもつこと
齢の場合には，ニシンの変動は動物性プランクトンの変動への影響がみられるが，2 年齢の場
合は小さい．ニシンの変動はタラへ影響を及ぼす方がその逆よりも若干高いが数値的にはあまり大きな差はない．図 5 に示された food web の仮説においても，動物性プランクトンとニシン，ニシンとタラは点線なので，実線で示された関係より弱い関係を示している．シシャモとニシンの関係については，先行研究で，バーレンツ海における成魚でないニシンは，シシャモの卵を餌食するので，ニシンの幼魚の資源量が 1 年齢未満のシシャモの生存を左右することが報告されている（Harme, 2000）．シシャモ 2 年齢と 3 年齢の変動がニシンの変動に関係している点については，その年齢のシシャモが動物性プランクトンやオキアミをよく餌食するので，ニシンと食糧の上で競合する可能性が影響していると考えられる．そしてニシンはシシャ
モ 3 年齢を餌食することはないので，さらにオキアミに関して競合関係になると考えられる．
（Dalpadado and Skjoldal, 1996）

トータルパワー寄与率を用いた海洋生態システムにおける因果性推測

333

図 9．シシャモとオキアミ，ニシンとオキアミのフィードバック経路を切断した場合のオキア
ミのシミュレーションデータ．上図がシシャモ 2 年齢，下図がシシャモ 3 年齢の場合の推定された MAR の係数を用いてデータを生成した．

この点について，加藤・石黒（1997）におけるシステム解析で，フィードバック経路を切断（マスキング）
した場合の各変量のシミュレーション値を比較する方法で検証してみる．
シシャモ 2 年齢と 3 年齢の場合で推定された MAR モデルの係数について，シシャモとオキアミのフィードバック経路（a13 (m), a31 (m), m = 1, 2），ニシンとオキアミのフィードバック
を 0 に置き換え，シミュレーションデータを生成し，図 9 にプ
経路（a34 (m), a43 (m), m = 1, 2）
ロットした．上図はシシャモ 2 年齢，下図は 3 年齢の場合のオキアミのシミュレーション値である．若干 3 年齢の値が高く出る傾向がみられる．すなわち，シシャモ 3 年齢の場合，シシャ
モとニシンのフィードバック関係がないと，オキアミの資源量はシシャモ 2 年齢のときよりも上がる傾向があるということである．それは，オキアミに対する競合関係がシシャモ 3 年齢のほうがより強いという傾向を示唆している．
近年，バーレンツ海におけるシシャモの回遊や資源量が変化し，シシャモを餌食とする生物種，例えば鯨類等がオキアミを主に餌食し始めているという報告がある（Haug et al., 2002;
Solvang et al., 2017）．そこで，仮にシシャモがバーレンツ海の生態システムから存在しなくなった場合を想定した分析として，シシャモを含まない他 3 生物種で構成されるフィードバッ
クシステムを仮定し，同様の分析をおこないその結果を図 9 にまとめた．タラとオキアミのフィードバック関係が高く，タラとニシン，オキアミとニシンとの直接的な関係はそう高い値によると，1–2 年齢のタラが 3–6 年齢のタラよりもよりを示さなかった．Bogstad et al.（2015）
多くオキアミを食し，3–6 年齢のタラは，オキアミよりより多くのシシャモを食することを述べていた．また，北大西洋タラとニシンの稚魚がシシャモを好んで餌食する（Gjøsæter et al.,と報告されているが，シシャモがいなくなった場合，タ2009; Hallfredsson and Pedersen, 2009）
ラやニシンの餌食としてオキアミは重要であることも指摘されている（Gjøsæter et al., 2009）．
図 10 での分析では，データとしてニシンの稚魚を扱っていないが，成魚のニシンとまったく無関係という結果ではなく，また，タラとオキアミはこれらの先行研究を裏付けていたといえる．従って，タラやニシンの餌食となるシシャモの変動が，次の餌食と考えられるオキアミの変動に関連すると推察される．

統計数理

334

第 66 巻 第 2 号 2018

図 10．シシャモを含まない生態システムを仮定した場合における他 3 生物種のフィードバッ
ク関係．

5.

おわりに

多次元時系列データ間の相互関係を調べる手法として，Ozaki（2012）は，赤池のトータルパ
ワー寄与を基にした，Partial pairwise contribution による因果推測を提案した．我々のシミュ
レーション実験では，Granger と Geweke の 2 変量間だけのペアワイズ比較よりも，トータルパワーに基づく比較の方がより正確な関係をとらえやすいことを示せた．ΔLL と情報量規準に基づく有意性の規準から，因果性の有無と度合いが判断でき，可視化も可能であることが示せた．MAR の全極を対象とするモデルのため，直接や非直接的な関係や，特にサンプル数が少ない場合，本来考慮しない関係をとる場合があることも示した．サンプル数が 100 点以上の場合，得られた関係をマスキングしたモデルか，もしくはその関係に注目し係数を数値的最適化で求めたモデル間の AIC を比較することで適した関係を示すモデルを再考でき，因果推測が可能となるだろう．実データ解析では，限られたサンプル数であったが，バーレンツ海の生態系における food web で重要な 4 生物種，43 年間のバイオマス時系列データを用い，food
web を駆動するシシャモの資源の重要性を裏付けた．提案手法は，複雑なフィードバックシス
テムの関係を予備的に予測することが可能で，あらかじめ検討をつけた関係について，さらに詳しく関係を踏まえたモデルを再考し，正確な関係をとらえるのに有効であると考えられる．
また，赤池の相対パワー寄与率は，本稿で扱ったようなあまり長くない時系列データの場合，
長中短期周波数領域における結果の解釈が困難な場合がある．提案手法は周波数領域に関する総和をとるため，そのような細かい解釈は必要なくなる．これまで ICES が海洋学関連の時系
列データ分析に主成分分析をスタンダードとしていたが，本稿でとりあげた手法がより適切であり，より具体的に変量間の関係をとらえることができるので，海洋生態系の変化に伴う様々な要因を探ることに貢献できると考えられる．
謝辞

本稿執筆に際し，突発的な質問にも拒まれず，常に丁寧に議論にお付き合いいただき，多くの助言と激励を頂いた尾崎統博士に，また多くの有益なコメントを頂いた査読者の方々に感謝の意を表します．

トータルパワー寄与率を用いた海洋生態システムにおける因果性推測

参

考

文

335

献

Akaike, H. (1968). On the use of a linear model for the identiﬁcation of feedback systems, Annals of
the Institute of Stattistical Mathematics, 20, 425–439.
Akaike, H. (1974). A new look at the statistical model identiﬁcation, IEEE Transactions on AutomaticControl, 19, 716–723.
Akaike, H. and Kitagawa, G. (eds.) (1994). The Practice of Time Series Analysis, Springer-Verlag,
New York.
赤池弘次，中川東一郎 (1972)．『ダイナミックシステムの統計的解析と制御』，サイエンスライブラリ（情，サイエンス社，東京．
報電算機 9）Baccalá, L.A. and Sameshima, K. (2001). Partial directed coherence: a new concept in neural structure
determination, Biological Cybernetics, 84, 463–474.
Baum, J.K. and Worm, B. (2009). Cascading top-down eﬀects of changing oceanic predator abundances,
Journal of Animal Ecology, 78, 699–714.
Benjamini, Y. and Hochberg, Y. (1995). Controlling the false discovery rate: A practical and powerfulapproach to multiple testing, Journal of the Royal Statistical Society. Series B (Methodological), 57(1), 289–300.
Bogstad, B., Gjøsæter Haug, H., T. and Lindstrøm, U. (2015). A review of the battle for food in
the Barents Sea: cod vs. marine mammals, Frontiers in Ecology and Evolution, 25, doi:
10.3389/fevo.2015.00029.
Bosch-Bayard, J., Wong, K.F.K., Okazaki, S., Oshio, R., Galka, A., Ozaki, T. and Sadato, N. (2012).
Directed causality for non-stationary time series based on Akaike’s noise contribution ratio,
FORMATH, 11, 121–131.
Dalpadato, P. and Skjoldal, H.R. (1996). Abundance, maturity and growth of the krill species Thysanoessa inermis and T. longicaudata in the Barents Sea, Marine Ecology Progress Series, 144,
175–183.
Geweke, J. (1982). Measurement of linear dependence and feedback between multiple time series,
Journal of the American Statistical Association, 77, 304–313.
Gjøsæter, H., Bogstad, B. and Tjelmeland, S. (2009). Ecosystem eﬀects of the three capelin stock
collapses in the Barents Sea, Marine Biology Research, 5, 40–53.
Granger, C.W.J. (1969). Investigating causal relations by economietric models and cross-spectral methods, Econometrika, 37, 424–438.
Hallfredsson, E.H. and Pedersen, T. (2009). Eﬀects of predation from juvenile herring (clupea harengus) on mortality rates of capelin (mallotusvillosus) larvae, Canadian journal of Fisheries andAquatic Sciences, 66, 1693–1706.
Harme, J. (2000). Capelin and herring as key species for the yield of north-east Arctic cod. Resultsfrom multispecies model runs, Scientia Marina, 67(Suppl 1) 315–323.
Haug, T., Lindstrøm, U. and Nilssen K.T. (2002). Variations in minke whale (Balaenoptera acutorostrata) diet and body condition in reponse to ecosystem changes in the Barents Sea, Sarsia, 87,
409–422.
Hjermann, D.ø., Strenseth, N.C. and Ottersen, G. (2004). Indirect climate forcing of the Barents Sea
capelin: A cohort eﬀect, Marine Ecology Progress Series, 273, 229–238.
加藤比呂子，石黒真木夫 (1997). 多変量時系列モデルによる経済システムの動的解析，統計数理, 45(2),301–318.
川崎能典 (2001). 多変量時系列に対する主成分・因子分析, 統計数理, 49(1), 109–131.
Kolmogorov, A.N. (1941). Stationary sequences in Hilbert space, Moscow University Mathematics
Bulletin, 5, 3–14.
Ozaki, T. (2012). Time series modeling of neuroscience data, Chapman & Hall/CRC, Boca Raton,

統計数理

336

第 66 巻 第 2 号 2018

Florida.
Planque, B. and Arneberg, P. (2017). Principal component analyzes for integrated ecosystem assessments may primarily reﬂect methodological artefact, ICES Journal of Marine Science,
doi:10.1093/icesjms/fsx223.
Solvang, H., Subbey, S. and Frank, A.S.J. (2017). Causal drivers of Barents Sea capelin (Mallotus villosus) population dynamics on diﬀerent time scales, ICES Journal of Marine Science,
doi:10.1093/icesjms/fsx179.
Tanokura, Y. and Kitagawa, G. (2004). Power contribution analysis for multivariate time series withcorrelated noise sources, Advances and Applications in Statistics, 4, 65–95.
WGIBAR. (2016). Final report of the Working Group on the Integrated Assessments of the Barents
Sea.
Wiener, N. (1956). The theory of prediction, Modern Mathematics for Engineers (ed., E.F.
Beckenback), Chapter2, 165–190, McGraw-Hill, New York.

Proceedings of the Institute of Statistical Mathematics Vol. 66, No. 2, 319–337 (2018)337

Causal Inference for Marine Ecosystems Based on Total Power
Contribution
Hiroko Kato Solvang1 and Subbey Sam2,3
1 Marine Mammals Research Group, Institute of Marine Research
2 Research Group on Fisheries Dynamics, Institute of Marine Research
3 Department of Natural Resources, Cornell University

We introduce a statistical methodology that integrates Granger’s pair-wise causal
analysis and its expansion to causality basedon the log-likelihood (Partial pairwise causality), and Akaike’s power contribution approach whole frequency domain (Total causality).
Although the initial idea was proposed by Ozaki (2012), it has hitherto not been applied
to complex marine ecosystem dynamics. In this article, we implement the approach byadding a criterion to assess signiﬁcance to detect causal relationship. We perform a simulation study to verify the eﬃcacy and sensitivity of the method, using data generated bythree autoregressive models with three and ﬁve dimensions. We also applied the method toreal observations to investigate causal drivers of Barents Sea capelin population dynamics.
The goal of this analysis was to explore inter-species relationships, which are importantfood web drivers in the Barents Sea ecosystem. We present results demonstrating that the
proposed methodology is a useful tool in early-stage causal analysis of complex feedback
systems.

Key words: Multivariate auto-regressive model, multivariate time series data, feedback system, Granger’s
causality, marine ecosystem, Barents Sea.