rEDM を用いた Empirical Dynamic Modeling: (5) Convergent Cross Mapping (CCM)

Convergent Cross Mapping (CCM)

ここでは変数間の因果関係を検出する Convergent Cross Mapping (CCM; Sugihara et al. 2012) について解説します。

CCM で使用されるアルゴリズムは基本的に Simplex projection に非常に似ており、従って Simplex projection を理解できていればアルゴリズムを理解することは比較的容易です。

1. CCM の概念的な説明

｢因果がある｣ とは

｢相関は必ずしも因果を意味しない｣ というのはあまりにも有名です (例えば、Wikipedia による解説も御覧ください)。ここではそれについては改めて詳しく説明することはしませんが、それではそもそも ｢因果がある｣ とはどういうことでしょうか？この疑問について深く考えると非常に難しく、哲学や宗教も絡んできそうですが、ここではひとまずそのようなことは考えず ｢同じ力学系に属していれば因果がある｣ と思うことにしましょう。例えば ｢時系列データを用いた動態の再構成｣ で用いたローレンツアトラクターは以下の微分方程式から現れる動態でした。

$$\frac{dx}{dt} = -pX + pY\\ \frac{dy}{dt} = -XZ + rX - Y\\ \frac{dz}{dt} = XY - bZ$$

この系は一つの力学系ですが、$X$, $Y$, $Z$ が同じ力学系に属して (お互いに因果を持ちながら) 動態を形成しています。このような状況のときには、例えば $X$ の時間変化に $Y$ が影響を与えていますし、直感的な意味で $Y$ は $X$ に因果がある、影響があると考えて良さそうです。もっと簡単に考えると、$Y$ の値が変化するときに、$X$ も変化する、ということができます。

同じ力学系に属しているかどうかを調べる

さて、CCM は２つの変数が同じ力学系に属しているかどうかを調べようとします (他の解釈も可能ですが、図的に説明しやすいのでここではこの解釈を説明します)。では、どうやって同じ力学系に属しているかどうか調べるのでしょうか？

下の図を覚えているでしょうか？左は $X(t)$ の時間遅れを用いて再構成した動態 ($M_x$ とします)、右は $Y(t)$ の時間遅れを用いて再構成した動態 ($M_y$) です。

以前説明したように、$M_x$ はもともとの $X$, $Y$, $Z$ を用いて再構成した動態 (便宜的に真の動態、と呼びます) と同じかたちとみなせます。同様に $M_y$ も真の動態と同じかたちとみなせます。従って、$M_x = [真の動態] = M_y$ となり、もし $X$ と $Y$ が同じ力学系に属していれば $M_x$ と $M_y$ は同じかたちになるはず です (下図)。

CCM では、Simplex projection を応用して、再構成された２つの動態が同じかたちかどうかを調べることで変数間の因果関係を検出しようとします。

注意: 最近、とある研究会で ｢動態の再構成についてローレンツアトラクタ−を例として説明するのはあまりよくないのではないか｣ とのご指摘をいただきました。その方はカオス時系列などの研究をされていた方でいくつかの根拠を持ってそのような指摘をしてくれました。この点について現在勉強中で、ひょっとすると将来説明の仕方が変わる可能性があります。

2. CCM のアルゴリズム

では、Simplex projection を応用して、どのようにして ｢同じかたちかどうか｣ をチェックするのか、R を利用して解説していきます。

マニュアルコードによる解説

例によって自作の２種系のデータセットを使用して説明します (四度目)。

# rEDM のロード
library(rEDM)
packageVersion("rEDM") # 2020年1月時点で v0.7.5

# 2種系のデータセットを作成
d0 <- as.data.frame(matrix(rep(NA,3*1000),ncol=3))
colnames(d0) <- c("time","x","y")
d0[,1] <- 1:1000
d0[1,2:3] <- c(0.1, 0.1)
for(i in 1:999){
        d0[i+1,2] <- d0[i,2]*(3.8 - 3.80*d0[i,2] - 0.02*d0[i,3])
        d0[i+1,3] <- d0[i,3]*(3.5 - 0.10*d0[i,2] - 3.50*d0[i,3])
}

# 最初の 99 点を捨てて、時間 100 から 1000 のデータを使用
d <- d0[100:1000,]

最適埋め込み次元の決定

まずは Simplex projection を利用して最適埋め込み次元 (E) を決定します。ここでは rEDM の simplex 関数を利用します。

simp_res <- simplex(y)
plot(simp_res$E, simp_res$rmse, type = "b", las = 1, xlab = "E", ylab = "RMSE")

RMSE に基づいて、最適埋め込み次元は 2 を使用することにします。

Ey <- simp_res$E[which.min(simp_res$rmse)]

最近傍点の特定

最適埋め込み次元を決定した後は Simplex projection と同様に $Y(t)$ の埋め込みを作り、ターゲットの点を決めて最近傍点を特定します。

# Y(t) の埋め込みを作成して、ターゲットの点を指定
y_embed <- make_block(y, max_lag = Ey)
t_target <- 240　# ターゲットの点を指定

# t = 240 の点から各点への距離を計算
distance <- sqrt(rowSums((y_embed[,2:3] - c(y_embed[t_target, 2:3]))^2))
nn_index <- order(distance)[2:4] # 自分自身を除いて 3 点近い点を選ぶ

# 最短距離を取得
min_distance <- distance[nn_index][1]
# 最短距離で補正した重みを計算
weights <- exp(-distance[nn_index]/min_distance)
# 合計の重みを計算
total_weight <- sum(weights)

Cross mapping

通常の Simplex projection ではこの後、以下のように $Y(t_{target})$ の近傍点の情報 (nn_index) から、$Y(t+1)$ を予測します。

# Y(t+1) を予測
pred <- (weights %*% as.matrix(y_embed[nn_index + 1,3:2])) / total_weight

一方、CCM では以下のようにして $Y(t_{target})$ の近傍点の情報 (nn_index) を元にして $X(t)$ の値を予測します。この

# Y(t) の近傍点の情報を利用して X(t) を予測
pred <- (weights %*% as.matrix(x[nn_index])) / total_weight

予測値を見てみましょう。かなり近い値を予測できていることが分かります。

> pred # 予測値
[1,] -1.21572
> x[t_target] # 実測値
[1] -1.254002

この操作は一体何をしているのでしょうか？これは、元の時系列に戻ってみると分かりやすいでしょう。下の図を使って説明します。

図の左上のパネルは $Y(t)$ のターゲット点周辺の様子です ($t_{target} = 240$ ですが、生成したデータの最初の 99 点を捨てているので実際の time は 339 となっています)。２次元で動態を再構成しているので、$Y(t)$ と $Y(t-1)$ の点の ｢動き｣ に注目していることになります。最近傍点の特定とは、ターゲットの点の ｢動き｣ と似た点を $Y(t)$ の時系列から探していることになります。これは図中では (1) Identify nearest neighbors に相当します。その後、特定した $Y(t)$ の時間と同じ時間の $X(t)$ を特定します (図中 (2) Cross-mapping)。$Y(t)$ の情報を元に特定された $X(t)$ の値を重み付け平均して $X(t)$ の値を予測しています (図中 (3) Prediction)。(図中の (2) と (3) を合わせて Cross-mapping と呼ぶべきかも知れませんが、ここでは説明の都合上２つを分けています)

Convergence のチェック

CCM では上記の手順を $X(t)$ の全てのデータ点で行い、$Y(t)$ の情報のみから $X(t)$ を予測できるかを調べます。さらに CCM では Simplex projection では通常行わない計算が加わります。それは、状態空間の再構成に利用するデータ点の数 (rEDM パッケージ内では library size と呼ばれます) を変化させ、データ点数の増加に伴って予測精度が上昇するかどうかをチェックします。この現象は Convergence と呼ばれていて、CCM の名前の由来にもなっています。

rEDM の ccm 関数を利用してこの Convergence をチェックしてみましょう。ccm 関数で求めた Cross-mapping による予測精度を ccm_means 関数で見やすい形にまとめています。

# Convergence のチェック
ccm_all_res <- ccm(cbind(y, x), E = Ey, lib_sizes = seq(Ey+1, nrow(d), by = 30))
ccm_res <- ccm_means(ccm_all_res)

# 結果の図示
plot(ccm_res$lib_size, ccm_res$rho, type = "l",
     xlab = "Library size", ylab = expression(rho), las = 1)
plot(ccm_res$lib_size, ccm_res$rmse, type = "l",
     xlab = "Library size", ylab = "RMSE", las = 1)

Library size の増加に伴って、相関係数 ($\rho$) の増加や RMSE の減少が見られます。これは $X$ から $Y$ への因果があることを示唆しています。

注意点: CCM では $Y$ の情報を使って $X$ を予測できたとき、それは $X$ から $Y$ への因果を意味しています。$X → Y$ のときは、$X$ の情報が $Y$ に流れており、そのため $Y$ の情報を用いた $X$ の予測が可能になります。

また、ここで ccm 関数の引数の説明もしておきます。

ccm(time_series,                             # 時系列データ
    lib = c(1, NROW(time_series)),           # 埋め込みに使用するデータの範囲
    pred = lib,                              # 予測値を出すデータの範囲
    norm = 2,                                # 重み付けのための距離の計算方法
    E = 1,                                   # 試行する埋め込み次元の範囲
    tau = 1,                                 # 時間遅れの単位
    tp = 0,                                  # 何ステップ先を予測するか
    num_neighbors = "e+1",                   # 使用する最近傍点の数
    lib_sizes = seq(10, 100, by = 10),       # library size の指定
    random_libs = TRUE,                      # library からベクトルを選ぶときにランダムに選ぶかどうか
    num_samples = 100,                       # ある library size に対して何回 library を作るか
    replace = TRUE,                          # library を選ぶ際に同じベクトルを複数回選択することを許すかどうか
    lib_column = 1,                          # cross-map 元 (= library) の変数 (減員側の変数)
    target_column = 2,                       # cross-map 先の変数 (因果側の変数)
    first_column_time = FALSE,               # 時系列のデータフレームの最初の列が時間かどうか
    RNGseed = NULL,                          # library の再構成に使用するデータを選ぶ際の random seed の設定
    exclusion_radius = NULL,                 # 最近傍点を探す条件として、予測したい点とどの程度時間的に近い点を除くか (数値を指定する)
    epsilon = NULL,                          # 最近傍点を探す条件として、予測したい点とどの程度距離的に遠い点を除くか (数値を指定する)
    stats_only = TRUE,                       # 予測の統計情報のみを出力するのか、それとも予測値も実際に出力するのか
    silent = FALSE)                          # Warning message を出力するかどうか

たくさんの引数がありますが、解析上重要なのは、time_series, E, lib_sizes, tp などです。特に tp を変化させることで時間遅れの効果による因果の判定もできます (Ye et al. 2015)。

なぜ Convergence が重要か？

では、なぜ Convergence をチェックすることで因果のあるなしが、すなわち２つのかたちが同じかどうか判定できるのでしょうか？ Cross-mapping を利用して片方の動態からもう片方の値を予測するということは、ローレンツアトラクターを例にすると以下のように２つの動態の中のある点とある点の対応をチェックしていることになります。

Library size が小さいときは、疎な状態のアトラクターを利用して最近傍点を特定し、Cross-mapping を行います。このときはたとえ二変数の間に因果があっても比較的遠くの点しか最近傍点として利用できず、当然 Cross-mapping による予測精度は低くなります。一方、だんだんと Library size を大きくすると、アトラクターが密な状態になっていき、より近い点を最近傍点として選べるようになります。仮に無限のデータ点が利用でき、そのデータに誤差が全く含まれていない場合は、無限に近い点を最近傍点として利用できることになり、最終的には完全な予測 (RMSE = 0) が達成できることになります。このとき、二つのアトラクター上の点と点は完全に 1:1 に対応していると判断でき、そのため二つのかたちは同じということができます。実際の時系列データは有限で誤差も含まれると思いますが、予測精度の上昇傾向を見て、最終的に二つのアトラクター上の点が 1:1 に対応しそうかどうかを調べている、ということです。

Sugihara et al. (2012) の補足資料である以下のビデオも参考になると思います。https://www.youtube.com/watch?v=NrFdIz-D2yM&t=10s

典型的な解析手順

CCM について解説しましたが、典型的な解析手順を再度まとめておきます。

1. 時系列データの標準化
2. Simplex projection による最適埋め込み次元 ($E$) の推定
3. 決定した $E$ を用いて ccm 関数で Cross-mapping を実行。この際、library_sizes 引数を指定して Convergence を同時にチェックします。典型的な library size の範囲は E + 1 から全てのデータ点 (length(y) や nrow(d)) までです。
4. ccm_means 関数で結果をまとめ、Convergence が見られるかチェックする。少なくとも最小のライブラリ数による予測精度を最大のライブラリ数による予測精度が上回っていなければ因果なしという判定になります。
5. Convergence が見られた場合は、サロゲートデータ (ある一定のルールに従って元の時系列データをシャッフルしたもの) を作成し、その Convergence が有意かどうか検定する。この際、rEDM の make_surrogate_data 関数が使用できます。どのサロゲートデータを使用するかも重要な問題ですが、初学者はまずは make_surrogate_ebisuzaki を試してみるのが良いかと思います。このサロゲートや検定に関してはまた時間があるときに解説してみようと思います。

解析上の注意点

CCM はうまく使用すれば非常に有用な手法ですが、比較的若い手法でまだまだ未解決問題やそれらに関連した注意点がいくつかあります。開発も続いているので、近い将来状況が変わる可能性は大いにありますが、現段階での注意点･未解決問題を列挙しておきます。自分のデータが以下の問題に当てはまるときは、EDM･CCM の原理を理解した上で慎重に解析を実行する必要があります。

• データに強い季節性がある (特に変数間でその季節性が同調している): このとき、単純な CCM では因果を誤検出する可能性があり、季節性を考慮した解析が必要です。現段階ではしばしば make_surrogate_seasonal 関数によって季節性を保存したサロゲートデータを作ることで対応しますが、これでも十分ではないかも知れません。Deyle et al. (2016), Baskerville & Cobey (2017), Sugihara et al. (2017), Cobey & Baskerville (2016) などを御覧ください。
• 欠損値が含まれる: 欠損値が含まれると時間遅れ埋め込みの際に状態空間中で位置を確定できないベクトルが増加します。そのため、CCM を行う際にライブラリとして使用するための十分なデータ点数を確保できなくなる可能性があります。多少の欠損値であれば何らかの方法で値を内挿し値を埋めることもありますが、一番はやはり欠損値がないようなデータセットを準備することです。
• 時系列データの時間間隔が一定ではない: 時系列データの時間間隔が一定でない場合は時間遅れ埋め込みによる状態空間の再構成が使えません。従って、複数のデータ点の平均値を取るなどして時間間隔を調整するか、EDM の適用を諦めることになります。
• データが相対優占度 (%) である: 微生物群集の DNA データなどにありがちですが、データが絶対量ではなく、相対優占度でしか手に入らない場合があります。相対優占度のデータは計算上、どうしても二変数間に因果が出てしまいます。つまり２種類の生物 A と B しかいない場合、A と B の間に因果関係がなくても、相対優占度は常に A = 1 - B となるように計算されており、明らかに A の相対優占度は B の相対優占度の影響を受けます。従って、CCM の安易な適用はできません。絶対量のデータを取得する、というのがもっとも直接的な解決方法で、Ushio et al. (2018) や Ushio (2019) はそれを目的として行われた手法開発です。
• 時系列データは短いが、時系列データに反復がある: 一般的には EDM は最低でも 35–40 点程度の時系列データが必要と言われていますが (もちろん動態によります)、これは生態学の分野ではなかなか厳しい制約です。たとえば、データが 1 年に 1 回しか取得できなければ 35–40 年分のデータが最低でも必要ということになります。ただし、データ点がそれよりも少ないが、調査地点が複数あるときは別の選択肢が考えられます。例えば、時系列データは 20 点しかないが、そのような調査が行われた地点が 10 地点ある、などの場合は Spatial convergent cross mapping (Clark et al. 2015) のような手法を適用することも考えられます。

終わりに

５回に渡って EDM について書いてきましたが、Simplex projection, S-map, CCM について理解できれば一通り EDM についての入門は終わりと言っても良いでしょう。あとは自分の時系列データに適用しながら理解を深めていって下さい。

EDM はまだまだ開発途上でこれから先数年のうちにまた新しい手法が出てくることが期待されます。そうなった場合にはまた解説記事を書くかも知れません。また今回の記事ではカバーできなかった部分も要望などがあれば再び追加で記事を書くかもしません (サロゲートデータや検定についてなど)。そうなったときはまた目を通していただければ幸いです。

参考文献

• Baskerville EB & Cobey S (2017) Does influenza drive absolute humidity? PNAS 114:E2270-E2271
• Cobey S & Baskerville EB (2016) Limits to causal inference with state-space reconstruction for infectious disease. PLoS ONE 11: e0169050
• Clark AT, et al. (2015) Spatial convergent cross mapping to detect causal relationships from short time series. Ecology 96: 1174-1181
• Deyle ER, Maher MC, Hernandez RD, Basu S, Sugihara G (2016) Global environmental drivers of influenza. PNAS 113:13081–13086
• Sugihara G, Deyle ER, Ye H (2017) Reply to Baskerville and Cobey: Misconceptions about causation with synchrony and seasonal drivers. PNAS 114: E2272-E2274
• Sugihara G, et al. (2012) Detecting causality in complex ecosystems. Science 338: 496–500
• Ushio M, et al. (2018) Quantitative monitoring of multispecies fish environmental DNA using high-throughput sequencing. Metabarcoding & Metagenomics 2: e23297
• Ushio M (2019) Use of a filter cartridge combined with intra-cartridge bead-beating improved detection of microbial DNA from water samples. Methods in Ecology and Evolution 10:1142-1156
• Ye H, Deyle ER, Gilarranz LJ, Sugihara G (2015) Distinguishing time-delayed causal interactions using convergent cross mapping. Scientific Reports 5:14750
• Ye H, et al. (2018) rEDM: Applications of Empirical Dynamic Modeling from Time Series. doi:10.5281/zenodo.1935847

Written with StackEdit.