rEDM を用いた Empirical Dynamic Modeling: (4) S-map の周辺

S-map の周辺

前回は S-map のアルゴリズムを解説しました。S-map は様々な拡張･応用が可能で、最近でも新しい方法が提案されています。このポストでは、S-map の周辺の話題について書いていきたいと思います。

1. Simplex projection との関係

これまでに解説してきたように Simplex projection と S-map のアルゴリズムには共通する点とそうでない点があります。２つの手法では ｢時間遅れ座標系による埋め込み｣｢状態空間上でのターゲットからその他の点への距離の計算｣ は共通でした。その後は ｢数個の最近傍点の挙動を参照するのか｣ ｢全ての点に重みを付けてそれらの挙動を参照するのか｣ という違いと、予測値の計算が ｢重み付け平均か｣ ｢局所重み付け線形回帰か｣ という違いがありました。どちらも近未来の非線形予測に使われるため、使い分けが難しい部分があるかもしれません。そこで、両者の違いについて、経験的なものも含まれますがここで挙げておきます (引用をつけられないものもあるので、自分で試してみることをおすすめします)。将来的に項目を追加するかも知れません。

• EDM の枠組みでは、Simplex projection は最適埋め込み次元 $E$ の決定に用いられ、S-map は (決定した $E$ を用いて) 非線形性 $\theta$ を定量するのに用いられるのが標準的です。もちろん初めから S-map のみを用いて $E$ と $\theta$ を同時に決定する、という手順も考えられますが、この場合は試行する組み合わせが $(Eの数) \times (\thetaの数)$ となるため計算量が増大します。
• 時系列データの特性にもよりますが、S-map の方が Simplex projection より予測精度が高くなることが多いようです。
• 自分で試行して分かったこととして、Simplex projection は予測精度がしばしば S-map に劣るものの、誤差の入り方が比較的均一です。一方で S-map にはシステマティックな誤差が入ります (= 予測値が観測値に比べて小さくなりがちな領域や大きくなりがちな領域が存在する)。

2. 相互作用強度の推定

S-map の重要な利用方法として、多変数 S-map への拡張と S-map 係数 の利用があります。ここでも再び自作の２種系のデータセットを使用して説明します。

# rEDM のロード
library(rEDM)
packageVersion("rEDM") # 2019.12月時点で 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])
}

# 最初の 100 time step のデータを捨てる
d <- d0[100:1000,]

多変数 S-map

これまで Takens の埋め込み定理に基づき、一変数の時間遅れ座標系を利用して動態を再構成してきました。実はこの埋め込み定理は二変数以上の場合に一般化が可能です (Deyle & Sugihara 2011)。

上の時系列を例にすると動態 $Y(t)$ は $\{Y(t), Y(t-1)\}$ で再構成できますが、 $\{X(t), Y(t)\}$ でも再構成できる、ということです。この多変数埋め込みを利用した S-map を多変数 S-map (Multivariate S-map) と呼び、このときの局所線形回帰の係数に特別な意味をもたせることができます (Deyle et al. 2016)。どういうことか見てみましょう。

上の時系列を一変数での埋め込みを用いて S-map を実行したときの局所重み付け線形回帰の式は以下の通りでした。

$$\hat{Y}(t_{target}+1) = C_0 + C_1 \times Y(t_{target}) + C_2 \times Y(t_{target}-1)$$

そして、多変数 S-map を用いると次のようになります。
$$\hat{Y}(t_{target}+1) = C_0 + C_1 \times Y(t_{target}) + C_2 \times X(t_{target})$$

$C_i$ は S-map 係数 (S-map coefficients) と呼ばれますが、$C_i$ が何を表しているか考えてみましょう。例えば、$C_2$ は $X(t_{target})$ が少し増加したときに $Y(t_{target}+1)$ に及ぼす影響、と考えることができます。別の表し方をすれば、$C_2$ は $\partial Y(t_{target}+1)/{\partial X(t)}$ と書けるでしょう。つまり、S-map 係数 $C_2$ は $X(t)$ が $Y(t+1)$ に及ぼす影響 と解釈できます。$X$ や $Y$ がどちらも生物であれば S-map 係数が 種間相互作用強度 を表していると解釈できますし、$X$ が環境変数、$Y$ が生物であれば、環境が生物に与える影響を表していると解釈することが可能です。さらに S-map では $C_i$ を各時間ごとに計算するため、時間に伴って変数間の相互作用強度が変わる場合にもそれを柔軟に追跡することが可能です (Deyle et al. 2016)。

rEDM による多変数 S-map と S-map 係数の出力

では、rEDM を用いて多変数 S-map を実行してみましょう。まず、埋め込みに使用する変数を指定して、その後 block_lnlp() という関数を実行します (多変数 S-map は s_map() 関数では実行できません)。

# 多変数埋め込みを作成
smap_block <- cbind(d$x, d$y)

# 多変数 S-map を実行し最適な theta を決定
m_smap_res <- block_lnlp(smap_block, method = "s-map",
                         theta = c(0, 1e-04, 3e-04, 0.001, 0.003, 0.01, 0.03, 0.1, 0.3, 0.5, 0.75, 1, 1.5, 2, 3, 4, 6, 8),
                         silent = TRUE)

上では、通常の S-map と同様にいくつかの $\theta$ を試して最も最適な $\theta$ を探しています。$\theta$ と RMSE の関係は以下のとおりです。

plot(m_smap_res$theta, m_smap_res$rmse, type = "b")

$\theta=8$ が最適なようなので、それを使って多変数 S-map を実行し、予測結果や S-map 係数を出力します。

# 多変数 S-map を実行
m_smap_res2 <- block_lnlp(smap_block, method = "s-map", theta = 8, silent = TRUE, save_smap_coefficients = TRUE)

m_smap_res2 の中身は以下のとおりです (出力の後半は省いています)。

> m_smap_res2
  embedding tp nn theta num_pred       rho         mae        rmse
1      1, 2  1  0     8      900 0.9999472 0.003973874 0.004498086

ちなみに一変数で埋め込んだときの結果は以下の通りです。

> s_map(d$y, E = 2, theta = 8, silent = TRUE)
  E tau tp nn theta num_pred       rho         mae        rmse
1 2   1  1  0     8      899 0.9989587 0.006919698 0.008887697

比べると多変数 S-map では予測精度が向上していることが分かります (一変数での RMSE = 0.0088…, 二変数での RMSE = 0.0044…)。埋め込みに使用する変数 (この例では $X$) が予測する変数 ($Y$) に因果があるとき、予測精度は多変数 S-map によって向上する傾向にあります。

S-map 係数は m_smap_res2$smap_coefficients 内に格納されています。

> head(m_smap_res2$smap_coefficients[[1]])
         c_1          c_2       c_0
1  0.5747497  0.024614030 0.6633858
2 -3.0510028 -0.009597196 3.0771957
3  1.7387918 -0.029947627 0.2782008
4 -1.7066457 -0.088475705 2.0510882
5 -1.9924075  0.042627392 2.1762527
6 -1.4671749  0.005263682 1.8050193

また、S-map 係数の分布を箱ひげ図で図示すると以下のとおりです。

boxplot(m_smap_res2$smap_coefficients[[1]][,1:2], las = 1, ylab = "S-map coefficients")
abline(h = 0, lty = 2)

$C_1$ は $Y(t)$ から $Y(t+1)$ への効果、$C_2$ は $X(t)$ から $Y(t+1)$ への効果です。$C_2$ は少しわかりにくいですが、平均的には負のようです。元の方程式から考えると、$C_1$ も $C_2$ も負で $C_1 \ll C_2$ ですから、S-map 係数の結果と合っているようです (ただし、元の方程式の相互作用強度は個体あたりの相互作用強度で、S-map 係数は個体群あたりの影響なので、単純に両者が対応するわけではないので注意が必要です)。

rEDM::block_lnlp() のオプション

ここで、block_lnlp() のオプションについて解説しておきます。

block_lnlp(block,                             # 時系列データ
           lib = c(1, NROW(time_series)),     # 埋め込みに使用するデータの範囲
           pred = lib,                        # 予測値を出すデータの範囲
           norm = 2,                          # 重み付けのための距離の計算方法
           method = c("simplex", "s-map"),    # Simplex projection を利用するか、S-map を利用するか
           tp = 1,                            # 何ステップ先を予測するか
           num_neighbors = switch(match.arg(method), simplex = "e+1", `s-map` = 0),
                                              # 予測に利用する近傍点の個数 (デフォルトは方法によって自動設定)
           columns = NULL,                    # 状態空間の再構成に使用する列
           target_column = 1,                 # 予測のターゲットになる変数の列
           stats_only = TRUE,                 # 予測の統計情報のみを出力するのか、それとも予測値も実際に出力するのか
           first_column_time = FALSE,         # block の第一列目が時間かどうか
           exclusion_radius = NULL,           # 最近傍点を探す条件として、予測したい点とどの程度時間的に近い点を除くか (数値を指定する)
           epsilon = NULL,                    # 最近傍点を探す条件として、予測したい点とどの程度距離的に遠い点を除くか (数値を指定する)
           theta = NULL,                      # 試行する theta の値
           silent = FALSE,                    # Warning message を出力するかどうか
           save_smap_coefficients = FALSE)    # S-map の局所重み付け線形回帰の係数を出力に含めるかどうか

ほとんど simplex() や s_map() と同じですが、いくつか違う点があります。まず、埋め込み次元を指定する E がありません。block_lnlp() では block をこれまでの一変数のベクトルではなく、データフレームもしくは行列として与えることで E を同時に指定していることになるためです。つまり、block が２列のデータフレームであれば自動的に E = 2 と指定したことになります (columns のオプションを指定することで与えたデータフレームのどの列を埋め込みに用いるかの指定も可能です)。また、method で Simplex projection を実行するか、S-map を実行するか指定できます。多変数 S-map を実行したいときは method = "s-map" と指定する必要があります。

EDM における最近の発展手法

• S-map は局所的に見ると単なる (重み付け) 重回帰であり、そのため様々な統計手法を取り入れることで容易に拡張できます。例えば、LASSO や Ridge など正則化を入れると、正則化 S-map (Regularized S-map; Cenci et al 2019) となり、こちらでいくつか試したところ予測精度の向上が見られました。Cenci et al. (2019) に R の解析コードが付属していますので、興味ある方は試してみて下さい。また、自分でも異なる実装をして Github で公開していますので、良ければこちらもお試し下さい (https://github.com/ong8181/random-scripts)。

• EDM に明示的に階層性や空間構造を取り込むために ガウス過程 (Gausiaan Process) を組み込んだ手法も提案されています。Munch et al. (2017) においてそれらが解説されていますし、rEDM にもすでに block_gp() という関数が実装されておりガウス過程 EDM が実行できます (空間構造を明示的に指定はできないようですが)。少し試してみたのですが予測精度がかなり良く、内挿に強い、という印象です。

• 多変数埋め込みをさらに拡張し予測精度を向上させた Multi-view embedding (多視点埋め込み?とでも訳したらよいでしょうか…) という手法も考案されています (Ye & Sugihara 2016)。様々な埋め込みにより少しずつ違う視点で状態空間を再構成し、予測精度の良い埋め込みを選抜した後、それらの予測を統合する、という方法により、多変数 S-map 法よりも良い予測を行います。こちらも rEDM に multiview() 関数として実装されています。アイデアがすごく好きなのですが、計算時間がかかることが多く、今の所自分はあまり利用していません。

• Multi-view embedding と似ていますが、更に大規模かつランダムにたくさんの埋め込みを作り、それらの予測を統合することで少ない時間点･多変数の時系列データを用いて正確かつ長期の予測が可能になる、という Randomly Distributed Embedding (RDE) という方法も提案されています (Ma et al. 2018) (まだ自分で試したことはありません)。

参考文献

• Cenci S, Sugihara G, Saavedra S (2019) Regularized S-map for inference and forecasting with noisy ecological time series. Methods in Ecology and Evolution 10, 650–660
• Deyle ER, May RM, Munch SB, Sugihara G. (2016) Tracking and forecasting ecosystem interactions in real time. Proc. R. Soc. Lond. B 283, 20152258
• Deyle ER, Sugihara G (2011) Generalized theorems for nonlinear state space reconstruction. PLoS ONE 6, e18295.
• Ma H, Leng S, Aihara K, Lin W, Chen L (2018) Randomly distributed embedding making short-term high-dimensional data predictable. PNAS, 115 E9994–10002
• Munch SB, Poynor V, Arriaza JL (2017) Circumventing structural uncertainty: A Bayesian perspective on nonlinear forecasting for ecology. Ecological Complexity 32, 134–143
• Sugihara G. (1994) Nonlinear forecasting for the classification of natural time series. Philos. Trans. R. Soc. A 348, 477–495
• Takens, F. (1981) Detecting strange attractors in turbulence. In D. Rand & L.-S. Young (Eds.), Lecture Notes in Mathematics (Vol. 898, pp. 366–381).
• Ye H, Sugihara G. (2016). Information leverage in interconnected ecosystems: Overcoming the curse of dimensionality. Science, 353, 922–925
• Ye H, et al. (2018) rEDM: Applications of Empirical Dynamic Modeling from Time Series. doi:10.5281/zenodo.1935847

Written with StackEdit.