時系列分析における因果推定(主にConvergent Cross Mapping)

# 時系列分析における因果推定

> Convergent Cross Mapping (CCM) を中心に時系列データから因果推定を実施する際の方針のメモ(備忘録)

## 時系列データでの因果推定アプローチ

**系列Xと系列Yの関係性を考える**

上記を解く際、複数の方針が存在する。その中でも代表的なものを取り上げた。


### LiNGAM
関数形には線形性を仮定する。外生変数を非ガウス連続分布と仮定する。なので、観測変数は非ガウス分布になる。そして、回帰分析と独立性の評価を繰り返すことで、観測変数の因果的順序を推定し、得られた係数行列から残りの下三角部分を推定することで因果関係を推定していく。これをVARモデルに拡張。

### Convergent Cross Mapping (CCM)
Empirical Dynamic Modeling (EDM)の一種。EDMは時系列データを非線形力学系として捉え、変数間の作用や近未来予測を実行する手法。数式を仮定せず、観測値から系の状態空間を再構成してモデルフリーに解析をおこなう。Simplex projectionやS-mapの考え方を複数の変数の交互の予測に拡張した手法で、相互作用のある非線形力学系ならばTakens' theoremにより、お互いの時間遅れ埋め込みの情報を使って相手の時間遅れ埋め込みを再構成できる（自分の時系列のみから相手の時系列を予測できる）という特徴を利用している。

### Granger因果
Granger因果とは、「相手が存在することで、予測精度が増加するかを判断」するものである。２つの時系列X、Y (X→Y)について考えた時、２つの予測残差の平方和を比較して、Xのデータを使用することで予測残差が有意に減少したかを検定するものである。もっと言うと、XのみからXを予測した場合と比べて、XとYの情報を使ってXを予測した場合の方が予測精度が向上すれば、YはXに対してGranger因果を持つことになる。CCMとの違いは、この方法はプロセス誤差が確率過程的であれば有効だが、決定論的な力学系においては検出力が落ちる可能性がある。加えて、背後に共変量が存在し、疑似相関が発生している場合には誤った検出結果になってしまう可能性がある。

## CCMについて
### 導入
２つの時系列(X, Y)について考えてみる。
これらの関係性は主に４つに分けられる。


| 記号的関係 | 解釈的関係 | 
| -------- | -------- | 
| X -> Y     | XがYの原因の１つでYは結果     |
| Y -> X     | YがXの原因の１つでXは結果     |
| X <-> Y     | XとYとが相互的に作用     |
| X     Y     | XもYも独立     |

※ここでは疑似相関を考慮していない。


これら２つの時系列の因果関係を推定する際、その方法は大きく３つあげられる。
　・VARLiNGAM
　・Granger因果
　・CCM
生物学・生態学分野において、そのデータは決定論的なものである。その際、各時系列変動は線形性を有するものは少なく、大抵は非線形な力学系を有する。よって、既存の「モデルリング→処理・推定」を実行することは、モデルフリーな推定よりも正確性が落ちる可能性がある。CCMは数式を仮定せず、時系列の観測値のみから系の状態空間を再構成して挙動を解析する手法である。この手法は2012年にScienceで発表されたものである。今回はこのCCMについて理解を深める一助とする。

### CCMの大まかな例
２変数XとYが$ X_t=f(X_{t-1}) $、$ Y_t=g(X_{t-1}, Y_{t-1}) $に従う時、XがYに一方的に影響を与える(X→Y)力学系である。CCMで２つの事象の時系列データから因果関係を推定する時、２つの事象間には決定論的で、非線形な力学系が存在することを仮定。そして、決定論的過程を対象としているためノイズを考えない。



### もう少し詳しい解釈
e.g.) X→Yの２つの時系列をデータから類推することを考える。
関係性は以下。
　 $ X_t=f(X_{t-1}) $
　 $ Y_t=g(X_{t-1}, Y_{t-1}) $

大きく２つの方針をとる。
1. 時点Tの$ Y_T $が未観測と仮定して、これを見積もる
2. 時点Tの$ X_T $が未観測と仮定して、これを見積もる

#### 時点TのYが未観測と仮定して、これを見積もる
$ X_{T-1} $と$ X_T $に１番近い値を示す連続した時点P-1、Pを探す
↓
$ Y_P $を$ Y_T $の見積もりとする
→$ Y_t=g(X_{t-1}, Y_{t-1}) $より、$ X_{T-1} $と$ X_{P-1} $が近い値としても、$ Y_T $と$ Y_P $が近い値となることはない。あくまでもYはXや他の要因に影響されるため、Xのみで一様にYは類推することは難しい。

#### 時点TのXが未観測と仮定して、これを見積もる

$ Y_{T-1} $と$ Y_T $に１番近い値を示す連続した時点Q-1、Qを探す
↓
$ X_Q $を$ X_T $の見積もりとする
→$ Y_t=g(X_{t-1}, Y_{t-1}) $より、時点T-1からTと時点Q-1からQで$ Y_t $が近い値となった時は、$ X_{T-1} $と$ X_{Q-1} $は近い値になる。 
→$ X_{T-1} $と$ X_{Q-1} $が近い値の時は、$ X_T $と$ X_Q $も近い値となる。なぜなら、$ X_t=f(X_{t-1}) $よりXはYを関与せずにXのみで決まるため。 
→この類推は正確性が高い。

---
決定論的な力学系から生じた時系列において、「結果から原因の類推(2.)」は正確におこなえる。しかし、逆は不正確(1.)。

つまり、CCMは２つの時系列を互いに類推して正確性を測ることで因果関係を判定する。その際、この正確性は、類推した時系列と実際の時系列の平均二乗誤差・相関係数で評価する。

相関係数を考えると先程の４つの関係性は以下のようになる。


| 記号的関係|AからBの類推|BからAの類推
| --------|--------|--- 
| A -> B     |無相関|正の相関
| B -> A     |正の相関|無相関
| A <-> B|正の相関|正の相関
| A　　 B|無相関|無相関


### CCMの流れ
手順は以下。
1. 埋め込み: 時系列を空間内に反映
2. Cross Mapping: ２つの時系列の埋め込みで得た図形を参照して類推

#### 埋め込み(Embedding)
力学系のある変数の動態は、過去の他の変数の動態の情報も含まれる。そこで、ある１つの変数の時間遅れベクトルを座標にプロットする(埋め込み)ことで、全ての動態を用いたアトラクタ(時間発展が落ち着いた状態)と１対１対応し、各点の相対的な位置関係が同等のアトラクタを再構成することが可能になる(Takensの定理)。

因果関係を推定したい２つの時系列$ (X={X_1,...,X_T} $, $ Y={Y_1,...,Y_T}) $から、埋め込みにより２つのアトラクタ$ (M_X, M_Y) $を描く。このとき、アトラクタを描く空間の次元を「埋め込み次元(Embedding Dimension, Eと表記)」という。 

Xからアトラクタ$ M_X $を再構成する場合、時間遅れをτとして時間遅れベクトルを作る。
　 $ X(t)=(X_t, X_t-τ, X_t-2τ,...,X_{t-(E-1)τ}) $
このベクトルはE次元空間の中の1点として表され、これらの点の集まりを「再構成されたアトラクタ$ M_X $」という。


e.g.) 
　時刻$ t_{1-(E-1)τ} $から時刻$ t_1 $までのXの動態
　時刻$ t_{2-(E-1)τ} $から時刻$ t_2 $までのXの動態
が似ているとき、$ x(t_1) $と$ x(t_2) $の座標は近くなる。点同士の距離が近いことは、その時系列における時刻(E-1)τの時間発展の類似性を表す指標と見なせる。


#### Cross Mapping

先程と同様に時系列XとYについて(X→Yの関係)、因果関係を類推することを考える。このとき、Yは過去のYとXによって決定するものである。その上で、$ M_X $からY、$ M_Y $からXをそれぞれ類推し、その正確性を評価することで因果関係を推測する(Cross Mapping)。

例えば時系列Yにおいて、類似するパターンを示す２つの時間帯が存在したとする(m-(E-1)τからm, n-(E-1)τからn)。これはアトラクタ$ M_Y $上の近傍点で表すとy(m)およびy(n)となる。そしてy(m)およびy(n)の距離が近いことは、対応するx(m)とx(n)が近い座標であることを意味する。なぜなら、２系列はX→Yの関係であるため、Yは過去のYとXによって決定づけているからだ。

一方$ M_X $からYを類推する場合、$ M_X $上の2点が近傍だからといって対応する2点が$ M_Y $上でも近傍であるとは限らない。なぜなら、２系列はX→Yの関係であるため、XはYからの影響を受けずにXを決定づけているからだ。

したがって今回の場合であれば、$ M_Y $からXの類推は$ M_X $からYの類推よりも正確性が高くなる。よってこれをもってXが原因でYが結果と判断する。


**類推の手順($ M_X $からYを例にとる)**
1. $ x(t) $のE+1個の最近傍点(近い順に$ x(t_1), x(t_2),...,x(t_{E+1})) $を探す
→ベクトル$ x(t) $はこれらE+1個の点の位置ベクトルの加重平均
2. $ x(t_1), x(t_2),...,x(t_{E+1}) $と時間的に対応する$ M_Y $上の点$ y(t_1), y(t_2),...,y(t_{E+1}) $を用いて、$ Y_t $を類推
→$ M_x $による$ Y_t $の類推を$ \hat{Y_t} $とする
3. 全ての時刻で$ \hat{Y_t} $を求める
→類推された時系列$ \hat{Y_t} $と実際の時系列$ Y_t $との相関係数$ \rho $から、Mappingの正確性を評価


## 参照
・8/6因果フェスのプレビュー：「系列Aと系列Bの関係は？」という問いに対する4つの素敵な解法について
http://takehiko-i-hayashi.hatenablog.com/entry/2015/07/30/113328
・George Sugihara, Robert May, Hao Ye, Chih-hao Hsieh, Ethan Deyle, Michael Fogarty, Stephan Munch. (2012)Detecting Causality in Complex Ecosystems. Science 338(6106):496-500.
https://science.sciencemag.org/content/338/6106/496.full
・pyEDMを使った生態学時系列データの因果推定
https://khigashi1987.github.io/pyEDM_Maizuru/
・Convergent cross mapping の紹介：生態学における時系列間の因果関係推定法
中山 新一朗,阿部 真人,岡村 寛, 日本生態学会誌 65：241 - 253（2015）
・rEDMの使い方, 京極大助/Daisuke Kyogoku
https://sites.google.com/site/dkyogoku/redm-no-shii-fang
・rEDM を用いた Empirical Dynamic Modeling: (5) Convergent Cross Mapping (CCM)
https://ushio-ecology-blog.blogspot.com/2020/01/20200113blogger0008.html
・rEDM: an R package for Empirical Dynamic Modeling and Convergent Cross Mapping
https://mran.microsoft.com/snapshot/2018-03-30/web/packages/rEDM/vignettes/rEDM-tutorial.html
・Empirical dynamic modeling for beginners
https://link.springer.com/article/10.1007/s11284-017-1469-9