rEDMの使い方

Empirical Dynamic Modeling (EDM) の使い方をエンドユーザー向けに解説します（備忘録としても）。EDMのちゃんとした解説を見ると微分同相写像とか良く分からない単語が出てきますが、以下ではそーゆー話はすっ飛ばします（すっ飛ばしても生態学データではほぼ問題ないと思われる）。私自身、エンドユーザーとして必要十分だと自分で判断した範囲内でしか勉強していません。一般性を重視しつつも、生態学っぽい実例を交えて分かりやすい解説を心がけています（その分、文字数は多くなります。また、正確性を犠牲にしている部分もあります）。間違いの指摘などは歓迎しますが、以下の内容を用いた解析等によって不利益を被ったとしても責任は負いかねますのでご了承ください。解析には R のパッケージ rEDM（CRANサイトからダウンロードする）を使います。

2018.08.31 追記

rEDM のパッケージのバージョンにより、関数の挙動が若干異なるようです。特に関数の出力オブジェクトのクラス（リストとかデータフレームとか）が異なる場合があるようで、以下のコードをコピペするだけでは上手くいかないこともあるかもしれません。出力の中身の確認やパッケージの取説の確認を適宜おこなってください。

EDMの基本的なアイデアについては George Sugihara のラボのサイトにある動画を参照してください。 http://deepeco.ucsd.edu/video-animations/

またこの解説では以下の論文を主に参考にしています。

Chang CW, Ushio M, Hsieh CH (2017) Empirical dynamic modeling for beginners. Ecological Research. DOI:10.1007/s11284-017-1469-9

中山新一朗，阿部真人，岡村寛 (2015) Convergent cross mapping の紹介：生態学における時系列間の因果関係推定法．日本生態学会誌，65，241-253．

■ 前提知識

力学系についての初歩的な知識が無いとちょっと難しいかもしれません。全く馴染みのない方には、『数理生物学入門』（巌佐庸 著 共立出版）の第4章付録Bの、離散世代個体群動態モデルの局所線形化による安定性解析の解説を読むことをお勧めします。線形代数は、とりあえず行列の加法と乗法がイメージ出来ていれば良いのではないかと思います。また、統計モデリングについての知識も無いよりかは有ったほうが理解の助けになるかもしれません。しかし EDM は既存の統計モデリングとはだいぶ違う哲学に基づくので、統計モデリングの知識があると逆にちょっと分かりにくい部分もあるかも…。

■ Simplex projection による埋込み次元の決定

埋込み次元 embedding dimension の決定はEDMを実装する上でかなり重要なステップです。Simplex projection により埋込み次元を決定する原理は以下の通りです。まずTakens’ theorem （ターケンスの埋込み定理？）により、注目している変数（例：一日にあなたが目撃するネコの数）が従う力学系の性質（どんな関数になってるか）は、その変数の有限な数の時系列データから再構成できます（i.e. 毎週日曜だけお隣のおじいちゃんがネコに餌付けしてるから日曜は近所にネコが多いとか、土曜はあまり外に出ないからそもそもネコをほとんど見ないとか、そういうことをイチイチ考慮しなくても、過去何日かのネコの観測数から明日見るネコの数は予測可能ということ。厳密に言うと”有限”の範囲が具体的にどれくらいの範囲かについての証明もあるが、力学系の真の次元数D という実証的には知りようの無さそうなパラメタが出てくるので、実証屋としては”たかだか有限個の”埋め込み次元に帰着する、ということに何となく納得しておくのが良いのではないかと思う）。この再構成作業を状態空間 State Space あるいは 多様体 Manifold の再構成と呼びます。

注目している変数の時系列 {x₁, x₂,..., x_t,...} がすでに手元にあるとします（例：目撃したネコの数を毎日手帳に書いていたとか）。Takens’ theorem から、ある自然数Eがあって、変数 x の従う力学系の状態空間はE個の時系列 {x_t, x_t-1, x_t-2,..., x_t-(E-1)} から再構成できる（実際は『時系列データを採る頻度』と『再構成する時に時系列の部分集合をどの時間間隔で作るか（通例ギリシャ文字のタウτで表されるようだ）』は別の問題ですが、ここではその問題は無視します）。

とりあえず、最適かどうかは置いておいて、任意の自然数 E での埋込み {x_t, x_t-1, x_t-2,..., x_t-(E-1)} を考えましょう。ここでは視覚的に分かりやすいように、E = 2 として2次元上の点 x = (x_t, x_t-1) が図１の曲線に沿って動くとしましょう（図１A）。この E がもし最適な埋込み次元数なら、Takens’ theorem から、この軌道は（観測されていない変数も含めて）この系を駆動してる力学系についての十分な情報を含んでいます。手持ちのデータのうち、『次の動きを予測したい点 x_P = (x_T, x_T-1)』を任意に選びます（図１A）。続いて、この x_P に近い順に、手持ちの時系列データから x₁, x₂, x₃ を探してきます（図１A）（E + 1 個の点を探す。この E + 1 個の点で x_P を囲むため。 x_P に近い順に E + 1 個の点を集めても常に x_P が囲まれる訳ではないが、囲まれないことはまれなはず）。ここで、手持ちの全ての時系列データ {x₁, x₂,..., x_t,...} から、x₁, x₂, x₃ がそれぞれ次の時刻（x = (x_t, x_t-1) で t を t + 1 とした時）にどの点に移ったかは分かります（図1B）。この x₁, x₂, x₃ それぞれが次の時間ステップに移る点を x'₁, x'₂, x'₃ としましょう（図1BC）。E 次元空間上で x'₁, x'₂, x'₃ の平均を（x_P と x₁, x₂, x₃ の距離で重みづけして）計算することで x_P:next を予測することができます（図1C）。こうして得られた予測値 x_P:next を、（手持ちの時系列データから分かる）x_P が実際に移った点 x_P' と比較することで、予測精度を評価することができます（Mean Absolute Error MAE とか相関係数 ρ とかで）。実際にはこれを色々な x_P に対して行うことで、ある E のもとでの予測精度の分布を計算します。さらにこれを様々な E に対して行うことで、最も予測精度の高い（i.e. 最適な埋込み次元数である） E を探索的に見つけることができます。

図１．Simplex projection の概念図。（A) 予測したい点 x_P を近傍の E + 1 点でかこむ。(B) 近傍の点（x₁, x₂, x₃）が次の時間ステップにどこに進んだかは、手持ちのデータから分かる。(C) x₁, x₂, x₃それぞれが次の時間ステップに進む点 x'₁, x'₂, x'₃ を使って x_P が次の時間ステップにどこに進みそうかを予測する（x_Pと x₁, x₂, x₃ それぞれとの距離による重みづけで平均を計算）。

なお、手元に十分に長い時系列データがある場合にはデータの半分を予測される側（x_P）として、のこり半分を予測する側（ライブラリと呼ぶ）として使ったりします。使える時系列が短い場合には、予測したい点を除いた全てのデータをライブラリとして使います。

E をだんだん増やしていくと（情報が増えるので）予測精度が増しますが、あまり E を増やし過ぎるとオーバーフィッティングが起きて予測精度が徐々に下がります。

Simplex projecton は R のパッケージ rEDM（CRANサイトからダウンロードする）に含まれるsimplex() で実行可能です。引数に時系列、時系列のうちどこを予測するのか・予測に使うのかと、試す埋込み次元を指定すればよい。予測する側の時系列と予測される側の時系列に同じものを指定した場合は、自動で予測点以外の点をライブラリとして用います。stats_only = “F” とすることで、統計量だけではなく各点の予測値と観測値を出力させることもできます。

以下のRコードを実行することで関数の挙動がわかるかと思います。またいろいろパラメタをいじって遊ぶこともできます（テストデータとしてはカオスを示すデータを使います。キレイな周期振動を示すデータだと、変数間の相関が強くなりすぎて、変なことが起こります）。

2020.01.02 追記

以下のコードでは simplex() 関数にデータフレームの列の形で時系列データを与えていますが、現在のバージョンではベクトルまたは2列以上の行列またはデータフレームを与えることになっているようです。行列・データフレームの場合は1列目をタイムステップのインデックスとして参照するようです。以下のコードでエラーが出る場合は適宜最新のバージョンの解説をCRANのサイトで参照してください。

# ライブラリの読み込みを忘れないように。

library(rEDM)

# 時系列データの生成：2種時間遅れ系

# 列が種、行が世代に対応するような行列 v を考える。

x.f <- function(v){

3.9*(1-(v[2,1]+0.3*v[1,1]+0.3*v[2,2])/2.4)*v[2,1]

}

y.f <- function(v){

3.8*(1-(v[2,2]+0.2*v[2,1])/2)*v[2,2]

}

g <- matrix(c(1.1,1.02,0.98,0.99),nrow=2,ncol=2)

for (i in 1:299) {

gn <- c(x.f(g[i:(i+1),]),y.f(g[i:(i+1),]))

g <- rbind(g,gn)

} # g:1列目に種x, 2列目に種yの個体数, 301 時間ステップ分

g <- cbind(g,rnorm(301,1.1,0.2)) #3列目に乱数の列であるzの動態を挿入

# 時系列のプロット

plot(200:301,g[200:301,1],type="l",ylim=c(0,2))

lines(200:301,g[200:301,2],col="red")

lines(200:301,g[200:301,3],col="blue")

# データの分割 有 vs. 無 で simplex projection: x を例として

simp.x.half <- simplex(g[,1],pred=c(251,300),lib=c(201,250),E=c(2:8))

simp.x.all <- simplex(g[,1],pred=c(201,300),lib=c(201,300),E=c(2:8))

plot(rho~E,data=simp.x.half,type="l",ylim=c(0.8,0.98))

lines(rho~E,data=simp.x.all,col="red")

# x と y の simplex projection の比較。x は時間遅れが入ってる分最適なEが大きい

simp.x <- simplex(g[,1],pred=c(201,300),lib=c(201,300),E=c(2:8))

plot(rho~E,data=simp.x,type="l",ylim=c(0.9,1))

simp.y <- simplex(g[,2],pred=c(201,300),lib=c(201,300),E=c(2:8))

lines(rho~E,data=simp.y,col="red")

# z のsimplex projection; 相関係数 rho の値がとても低い

simp.z <- simplex(g[,3],pred=c(201,300),lib=c(201,300),E=c(2:8))

plot(rho~E,data=simp.z,type="l")

■ S-map 法による非線形性の確認

S-map（エスマップ）の使い方に関して Web で検索をかけると、解散した某音楽グループについての記事ばかりが出てきて困ります。さて、S-map (Sequential locally weighted global linear map) は simplex projection 同様に、時系列データから再構成された E 次元上の軌道を使って、ライブラリからあるデータ点を予測する手法です。ただし、simplex projection が予測したい点の近傍の E + 1 点のみを使うのに対し、S-map ではすべてのライブラリデータを使います（rEDMでは予測する側の時系列と予測される側の時系列に同一のものを使う場合は、simplex projection の時のように自動で予測点を除いた点の集合がライブラリとして用いられます）。

S-map では、軌道（多様体）上でのライブラリの各点と予測したい点の距離 d を計算し、それぞれの点を exp(−θd/dm) で重みづけして回帰をすることで予測を生成します。ここで dm はライブラリから計算される距離を標準化するための定数で、θ は非線形性を規定するパラメタです。θ が大きいほど、予測したい点の近傍の点の重みが増すことになります。θ = 0 は線形回帰です。したがって、θ が大きくなるにつれて予測精度（相関係数などで評価）が 上がれば系が非線形性を持つ（i.e. 過去の時系列の効果が状態依存的である）、と判断できます。なお埋込み次元数は simplex projection の結果によって決定します。

Simplex projection のところで書いたコードに続けて以下のコードを実行することで、相互に影響を及ぼしあっている時系列（種 x: g[,1] と種 y: g[,2]）は非線形なふるまいをしていることが分かります。使う関数は s_map() で、引数として時系列、埋め込み次元E、ライブラリと予測される点の範囲、それに θ の値 (theta) をベクトルとして与えます。ここでは 0 から 20 まで 0.1 刻みでθ を試します。

E.x <- 3 # x の埋込み次元

E.y <- 2 # y の埋込み次元

smap.x <- s_map(g[,1],E=E.x,lib=c(201,300),pred=c(201,300),theta=seq(0,3,0.1))

smap.y <- s_map(g[,2],E=E.y,lib=c(201,300),pred=c(201,300),theta=seq(0,3,0.1))

plot(rho~theta,data=smap.x,type="l")

lines(rho~theta,data=smap.y,col="red")

smap.x$theta[which.max(smap.x$rho)] # 種 x の最適なθ smap.y$theta[which.max(smap.y$rho)] # 種 y の最適なθ

# 乱数である z に対して同様のことをすると滅茶苦茶な挙動を示す。

smap.z <- s_map(g[,3],E=c(2,3,4,5,6),lib=c(201,300),pred=c(201,300),theta=seq(0,30,0.1))

plot(rho~theta,data=smap.z[smap.z$E==2,],type="l",ylim=c(min(smap.z$rho),max(smap.z$rho)))

lines(rho~theta,data=smap.z[smap.z$E==3,],lty=2)

lines(rho~theta,data=smap.z[smap.z$E==4,],lty=3)

lines(rho~theta,data=smap.z[smap.z$E==5,],lty=4)

lines(rho~theta,data=smap.z[smap.z$E==6,],lty=5)

■ CCM による因果の検出

Cross Mapping は2つの時系列（１日に目撃したネコの数とイヌの数とか）の間で simplex projection による予測を行い、変数間の因果関係を推定する方法です。例えばイヌの数（の時系列）からネコの数を予測できる場合、イヌの数にはネコの数の情報が含まれていることになります。ここから、ネコの数がイヌの数に影響している（i.e. ネコの数からイヌの数に因果がある）ことが分かります。ここで注意すべき（そして分かりにくいの）は、予測が成立する方向と因果の方向が逆だということです。ネコの数がイヌの数に影響していたとしても、イヌの数は必ずしもネコの数だけで決まっているわけでは無いので、ネコの数からイヌの数を予測できるとは限りません（i.e. 因果の判定に使えない）。いっぽう、ネコの数がイヌの数の動態に影響するなら、Takens’ theorem から、イヌの数の時系列にはネコの数の動態を含む力学系の情報が含まれています。このへんの話は本当に分かりにくいです。中山ら（2015 日本生態学会誌）に分かりやすい説明があるので参考にしてみてください。また「原因とは何か」について引っかかる人は、ドーキンスの『延長された表現型』を読むと良いかもしれません（因果の操作的定義については以下で言及します）。

もう少し詳しくcross mapping を説明します。対応する２つの時系列データ {x₁, x₂,..., x_t,...}と {y₁, y₂,..., y_t,...} があるとします（イヌとネコの数とか）。１変数のSimplex projection（図１）で考えたように xt の埋込みからなる埋込み次元 E （ここでは例として E = 2 とする）のベクトル x = (x_t, x_t-1)を考えます。このとき x は時間の関数 x(t) なので、ある x(t) に対応する y_t を定義することができます（y_t-1 でも y_t-2 でも良い）。cross mapping では、予測したいある y_t の値（y* としましょう）に対応する x （x* とします）を考えます。x の再構成された多様体上において、x* の近傍にある E + 1 個の点（x₁, x₂, x₃）を選び、x* からの距離で重みづけします。この重みづけと、x₁, x₂, x₃ に対応する y の値 y₁, y₂ y₃ を使って、y* の値を予測します。この予測が cross map で、simplex projection 同様に相関係数などを使って予測精度を評価することができます。

さて、因果の判定方法についてです。予測を生成する際に、上の説明では x₁, x₂, x₃ をライブラリ全体から探してきました。因果の判定をする際には、このライブラリから一部を間引きます。疎になったライブラリを使った場合と、より密なライブラリを使った場合とで、密なライブラリを使った場合に予測精度が上がることを Convergent Cross Mapping と呼びます。EDM の枠組みではこの Convergent Cross Mapping が変数間の因果の判定基準として使われます。

この CCM は rEDM の関数 ccm() で実行することができます。順を追って関数の挙動を追ってみましょう。引数として、時系列データ、ライブラリの埋込み次元、ライブラリの列、ターゲット（ライブラリから予測するデータ）の列、間引いた後のライブラリサイズ、予測する点の時間遅れ（e.g. x(t) からy_t-1 を予測するなら tp = -1）、そして間引くことで生成するライブラリの部分集合を（各ライブラリサイズに対して）何回生成するか（num_samples）を指定します。ここでもデータは上で生成した２種の個体群動態を使います。また埋込み次元は上で得られたsimplex projection の結果を使います。

いっぺんに convergence の評価まですると大変なので、まずはccm() の挙動を以下のコードで確認してみましょう。ここでは気まぐれに、全データ点のうち100個だけを抽出して cross mapping をする、という操作を10回繰り返してみます。 lib_sizes で抽出するデータ点を指定すれば、あとは勝手に関数がうまいこと処理してくれます。

cmxy.example <- ccm(g,E=3,lib_column=1,target_column=2,

lib_sizes=100,tp=0,num_samples=10)

cmxy.example #出力の構造

10回のライブラリ生成それぞれに対してsimplex projection が行われて予測精度などが計算されているのが分かると思います。

では、ライブラリサイズを変えたときに予測精度が上がるかを試してみましょう。予測精度も高く、ライブラリが大きくなると急速に予測精度が上がる様子が見られると思います。

lib.size <- c(4:301)

cmxy.tp0 <- ccm(g, E=3, lib_column=1,target_column=2,

lib_sizes=lib.size, tp=0,num_samples=10)

plot(rho~lib_size,data=cmxy.tp0)

cmxy.tp1 <- ccm(g, E=3, lib_column=1,target_column=2, #tp=-1 でもやってみる

lib_sizes=lib.size, tp=-1,num_samples=30)

plot(rho~lib_size,data=cmxy.tp1)

cmxy.tp2 <- ccm(g, E=3, lib_column=1,target_column=2, #tp=-2 でもやってみる

lib_sizes=lib.size, tp=-2,num_samples=30)

plot(rho~lib_size,data=cmxy.tp2)

で、実際にConvergent Cross Mapping が起きているかどうか（i.e. ライブラリサイズの増大によって有意に予測精度が上がるか）は統計検定を使います。いくつか方法が考えられているようですが、基本的な流儀は2種類で、(i) ライブラリサイズと予測精度の相関が有意かを判断する方法と、(ii)予測精度とconvergence の程度が帰無モデル（surrogate data）のもとでのそれより有意に大きいかを判断する方法の2つです。各流儀内で細かな手法についてはいろいろ考えられているようです。個人的な印象としては、surrogate data 法は（2017 年現在）surrogate data の生成法についてコンセンサスが無いらしく、特殊な事情がない限り (i) の方法をとるのが良いように思います。以下では tp=0 から tp=-2 までの間で、simplex projection を行うのに最低限必要なライブラリサイズ（ここでは E = 3 なので 4 点）から手持ちのデータサイズ（301）までライブラリの大きさを変化させます（ライブラリ用のデータから 301 回の復元抽出をすることで、ライブラリサイズが最大の場合でも予測精度の“分布”を生成することができる [注]）。それぞれのライブラリサイズに対して、以下では1000回の繰返し（ライブラリの生成）を行うことで、1000個の予測精度を発生させます。この分布の両側の 2.5% がどこに来ているかを調べ、予測精度の95%信頼区間がライブラリサイズ間でかぶっているかどうか調べます。以下のコードを実行すると、いずれの時間遅れでも有意に因果が検出されるのが確認できると思います。さらに、ライブラリサイズが最大のときの 50% quantile を tp=0 から tp=-2 の間で比較することで、tp=0 でもっとも予測精度が高くなっていることが分かります。これは最初にこの時系列を発生させたときに、種 x は時間遅れ無く y から影響を受ける、とした仮定と合致しています。

[注] このようにライブラリサイズ最大の場合についても分布を生成する方法が主流のようですが、個人的には全ての手持ちのデータのもとでの予測精度（i.e. データサイズ＝ライブラリサイズで非復元抽出）が最小ライブラリサイズでの予測精度の分布から有意に外れているかを見れば良いのではないかと思う。

#tp=0 でのconvergence の確認

cmxy.tp0.min <- ccm(g, E=3, lib_column=1,target_column=2,

lib_sizes=4, tp=0,num_samples=1000)

cmxy.tp0.max <- ccm(g, E=3, lib_column=1,target_column=2,

lib_sizes=301, tp=0,num_samples=1000)

quantile(cmxy.tp0.min$rho,p=c(0.025,0.5,0.975))

quantile(cmxy.tp0.max$rho,p=c(0.025,0.5,0.975))

#tp=1 でのconvergence の確認

cmxy.tp1.min <- ccm(g, E=3, lib_column=1,target_column=2,

lib_sizes=4, tp=1,num_samples=1000)

cmxy.tp1.max <- ccm(g, E=3, lib_column=1,target_column=2,

lib_sizes=301, tp=1,num_samples=1000)

quantile(cmxy.tp1.min$rho,p=c(0.025,0.5,0.975))

quantile(cmxy.tp1.max$rho,p=c(0.025,0.5,0.975))

#tp=2 でのconvergence の確認

cmxy.tp2.min <- ccm(g, E=3, lib_column=1,target_column=2,

lib_sizes=4, tp=2,num_samples=1000)

cmxy.tp2.max <- ccm(g, E=3, lib_column=1,target_column=2,

lib_sizes=301, tp=2,num_samples=1000)

quantile(cmxy.tp2.min$rho,p=c(0.025,0.5,0.975))

quantile(cmxy.tp2.max$rho,p=c(0.025,0.5,0.975))

■ Multivariate simplex projection による予測精度の向上

さて、最初に行ったの simplex projection では x の時間遅れ埋込みだけから x の動態を予測しました。しかし CCM の結果から、 x は y からも影響を受けていることが分かりました。このようなとき、x と y を両方使って simplex projection を行うことができます（S-map も）。使う変数が一種類か複数あるかを区別するときには、univariate simplex projection とか multivariate simplex projection と呼んだりします。multivariate のときの埋込み次元数（ベクトルの次元）をどうするかについては（2017年現在）あまりコンセンサスが無いようです。univariate のときと次元数は変えずに時間遅れが大きい変数を他の変数で置き換える（以下でやります）、あるいは univariate で最適化された埋込みに因果のある変数を追加する、といった方法が（極端に偏った結果を生み出す危険が小さく）無難なようです（ただしこの点については最新の方法論に注意を払っておいたほうが良いと思う）。

ためしに、最初に生成した2種の個体数を使って試してみます。最初に行った univariate simplex projection から種 x の最適な埋込み次元は E = 3 と推定されました。また、CCM の結果から時間遅れ無く tp=0 で y から x への影響が最も強く検出されました。そこで、{x_t, x_t-1, x_t-2} で行っていた埋込みを {x_t, x_t-1, y_t}で考えてみましょう。

rEDM では multivariate の simplex projection や S-map を block_lnlp() 関数を使って実行します。simplex() や s_map() では データと、参照すべき（予測やライブラリに使う）変数を別々の引数として与えました。block_lnlp() には各行が {x_t, x_t-1, x_t-2} に対応するベクトルとなるような data.frame を与えてあげる必要があります。具体的には以下のようにしてデータを加工します。

normalized.g <- as.data.frame(apply(g,2,function(x) ((x-mean(x))/sd(x))))

multi <- cbind(normalized.g[2:301,],normalized.g[1:300,1])

ここではさらにデータを平均ゼロ、標準偏差１となるような標準化をしています。特に複数の変数を同時に扱って EDM の解析をするときには結果の解釈を容易にするためにも変数の標準化が推奨されるようです（ホントは上の解析でもそうしたほうが良いのかも）。次のコードで multivariate simplex projection が実行できます。

multi.simp.x <- block_lnlp(multi,method="simplex",lib=c(1,300),pred=c(1,300))

multi.simp.x$rho # 予測精度の表示

比較をすれば、univariate simplex projection のときよりも予測精度が上がっているのが分かると思います。これは、もともと個体数の動態を生成した式に出現する変数のみを使ってデータの解析を行っているためと解釈できます。

■ S-map 係数による変数間の関係の記述

上で説明したように、S-map は距離による重みづけを用いて局所的な回帰を繰り返します。したがって、個々の予測したい点 x_P に対して回帰係数が得られます。これを multivariate の埋込みで行うことで、ある変数から他の変数への影響を個々の時間点について計算することができます。例えば、種間相互作用の強さが時々刻々どのように変化するかを明らかにできます。ただし、ここで計算される“相互作用の強さ”は、典型的な個体群動態モデル（e.g. Lotka-Volterra model）に出てくる競争係数などとは違うものである、という点に注意が必要です（ここでの例でいえば、∂x_t+1/∂y_tとかが出てきます）。上に出てきた種 x と種 y の相互作用の強さの時間変化を試しに推定してみます。

test.theta <- block_lnlp(multi,method="s-map",num_neighbors=0, # データは上で作ったmulti

theta=seq(0.5,15,0.5),target_column=1)

plot(test.theta$rho~test.theta$theta,type="l")

test.theta[which.max(test.theta$rho),"theta"]

# とりあえず theta= 15 とする。最適かわかんないけど

smap.res <- block_lnlp(multi, method="s-map", num_neighbors=0, theta=15,

target_column=1,save_smap_coefficients=T)

smap.coef <- as.data.frame(smap.res[[1]]$smap_coefficients)

colnames(smap.coef) <- c("x0","y0","x1","intercept")

plot(smap.coef[1:100,"y0"],type="l",ylim=c(-2,1)) # 種 y の効果 # 試しにt=100までプロット

lines(smap.coef[1:100,"x0"],col="red") # 種 x の（時間遅れのない）効果

lines(smap.coef[1:100,"x1"],col="blue") # 種 x の時間遅れのある効果

■ 解析上のTips

〇 解析に必要なデータの長さ

Sugihara et al. (2012, Science) の主張によると、およそ 35-40 回のセンサスデータがあれば EDM の解析ができるそうです。ただし必要なデータサイズはモデルの複雑さ（埋込み次元 E）にも依存します。

〇 複数の時系列の結合

実際には、単一の時系列で 35~40 のデータ点が得られないことも多いと思います。そういった場合には、同じような力学系に従っていると考えられる複数の時系列データ（同じ条件で行われたレプリケーションのある種間競争の個体群動態とか）をまとめて多様体の再構成を行うことができます（2017年現在。今後このへんの方法論は新手法が出てくるかも）。rEDM であればrbind()でデータをつなぐことになるかと思います（読み込むcsvファイル上で統合してももちろん構いませんが）。単純に時系列をつなぐと、前の時系列の終わりと次の時系列の初めの間が一連の時系列であるとして扱われてしまうので、この点に注意が必要です。rEDM の simplex()、ccm()、block_lnlp() では、引数の lib（pred も同様）に2列の行列を与えることでこの問題を解決できます。各行が一つの時系列を表し、1列目が時系列の開始行（データの方の行）、2列目が時系列の終了行を指定するような行列を与えることで、あとは勝手に関数が処理してくれます。例えば30 個のデータポイントからなる観測が3反復行われている場合は、以下のような行列を指定します。

matrix(c(1,30,31,60,61,90),ncol=2,byrow=T)

■ 謝辞

長田穣さん、川津一隆さんに色々と教わりました。ただし上では自分で納得した事だけを書いたので、間違い等はすべて京極に責任があります。

■ 更新履歴

2017.Jun.23. 『rEDMの使い方』公開

2017.Jun.23. タイポの修正。Google Docs からコピペしたのがマズかったか。

2017.Jul.10. CCM判定のためのブートストラップ検定について注釈を追加。

2017.Jul.12. multivariate S-map のコードを修正。

2017.Sept.7 CCMの予測と因果の関係についての説明を修正。

2017.Dec. 1 図１Bが E = 3 の場合になっていたので本文に合わせて E = 2 に修正。

2018.Apr.15 山道真人さんから指摘をいただいた univariate simplex projection のコーディングエラーを修正しました。またS-map係数の推定をするコードがマシンによってうまく行かないそうです（どうやら rEDM のパッケージのバージョンにより関数の出力オブジェクトのデータ型〔リストとかデータフレームとか〕が異なることが原因の模様）。

2018.Jun.1 CCM のコーディングの説明が分かりにくかったので修正しました。

2021.Jul.30 旧Googleサイトから新Googleサイトへ変換。ざっと読んで変換ミス等は修正したつもりですが、見落としがあるかも。