2020年12月30日水曜日

Shiny と rEDM による Empirical Dynamic Modeling のウェブアプリ

20201230_Blogger0014

以前紹介した非線形時系列解析 Empirical Dynamic Modeling (EDM) (https://ushio-ecology-blog.blogspot.com/2019/12/20191205blogger0004.html など) ですが、R の Shiny を用いてコーディングなしで EDM を体験できるウェブアプリケーションを作成しました。このアプリでは EDM の基本的なツールである Simplex projection (Sugihara & May 1990) ・ S-map (Sugihara 1994) と、発展的な手法で変数間の因果関係が推定できる Convergent Cross-Mapping (CCM; Sugihara et al. 2012) を実行できます。計算パッケージとしては rEDM (Ye et al. 2018) を使用しています。

ウェブアプリケーションは以下からアクセスできます。
Empirical Dynamic Modeling with Shiny App

*注意点: このアプリは EDM を体験することを主目的に作りました。そのため、計算をスピーディに進めることを重視しており、データ解析としては甘い部分があります。特に CCM の結果を確定的なものと捉えないでください。

*注意点: 公開には現在のところ shinyapp.io の無料アカウントを使用しており、月当たり一定のユーザー数・使用時間数を超えるとアクセスできなくなります。アクセス不可が頻発するようでしたら他の公開方法に切り替える可能性もあります。

データの形式とアップロード

このアプリで受け付けるのは以下のような形式の CSV ファイルです。1列目に “time_index” という列名で時間のインデックスを入力してください。この列名は他のものは受け付けません。

2列目・3列目に CCM で因果関係を推定したい変数の時系列データを入力します。これらの列名は何でも構いません。

enter image description here

手持ちの時系列データがない場合はアプリ中の「Two-species demo data」をクリックしてデモデータをダウンロードし、それをそのままアップロードしてください。

enter image description here

このデモデータは以下の式に従って作成した X が Y に影響を与えている時系列データです。Y から X への影響はありません。

X[t+1]=X[t]×(3.83.8×X[t])Y[t+1]=Y[t]×(3.53.5×Y[t]0.1×X[t])X[t+1] = X[t] \times (3.8 - 3.8 \times X[t]) \\ Y[t+1] = Y[t] \times (3.5 - 3.5 \times Y[t] - 0.1 \times X[t])

デモを走らせるには time_index と2変数のデータが格納されたデータセットが必要です。1変数のみのデータでも Simplex projection や S-map は可能ですが、それを行いたい場合はダミー変数でもいいので (0, 1, 0, 1, … とか) 1列分余計にデータを追加してください (なお、ダミー変数を入れた状態で CCM の実行ボタンを押すとアプリが停止することがあります)。

1変数のみのデータをアップロードすると以下のように警告が出ます。

enter image description here

正しい形式のデータをアップロードしたら、「Show standardized timse series」のボタンを押すと時系列データを図示します。表示範囲はスライダーで 0 – 300 の範囲で適当に変更できます。300 時間点以上のデータセットもアップロード & 解析可能ですが、CCM の計算時間が長くなるかもしれません。また、短い時系列データもアップロード & 解析可能ですが、概ね 30 – 40 時間点以上の時系列データが望ましいです。短い時系列 & 多数の反復を持つデータセットは反復間に NaN を挟むことで Simplex projection と S-map が可能ですが、CCM には今の所対応していません (この場合は time_index は NaN にせず、時系列の値だけ NaN にしてください)。

enter image description here

Simplex projection の実行

次に、Simplex projection と S-map で解析する変数を指定します。この変数は CCM の際に「結果」の候補変数としても使用され、指定しなかった変数が「原因」の候補変数となります。また、Simplex projection では埋め込み次元 (過去の EDM に関する記事を参照してください) の推定を行いますが、この際に使用する予測精度の指標も指定します。特に好みがなければデフォルトの RMSE を指定してください。

その後、「1. Perform simplex projection: Check embedding dimension」をクリックします。

enter image description here

ボタンを押すと E = 1 – 20 までを使用したときの予測結果が出力されます。予測値は Leave one-out cross-variation (LOOCV) によるものです。下の図の赤丸が最適な予測を示す E を表しています。ここも詳細については過去の記事を参照してください (https://ushio-ecology-blog.blogspot.com/2019/12/20191211blogger0005.html)。

図のタイトル部分に変数名・Simplex projection による E の推定結果・使用した予測指標を表示しています。

enter image description here

enter image description here

S-map の実行

Simplex projection で推定した最適埋め込み次元を用いて S-map を
実行します (S-map の詳細)。「2. Perform S-map: Check nonlinearity」のボタンをクリックすると S-map を実行し、Simplex projection と同様に予測結果を出力します。最適な非線形パラメーター (theta) を決定しますが、この theta は次の CCM では特に使用していません。アップロードしたデータに非線形があるかどうかをチェックしています。最適な theta > 0 であれば EDM の枠組みでは「非線形な動態」とされます。

ここでも図のタイトル部分に変数名・S-map による theta の推定結果・使用した予測指標を表示しています。

enter image description here

enter image description here

CCM の実行

最後に最初に指定した変数に、もう一つの変数が因果を及ぼしているかどうかを CCM で推定します (CCM の詳細)。CCM の前に結果判定に使用する予測指標を指定します。オリジナルの論文 (Sugihara et al. 2012) では rho (相関係数) が使用されていますが、MAE や RMSE も指定できます。特に好みがなければ rho のままにしておいて構いません。なお、rEDM::ccm()には tpという引数がありますが、ここでは -1 を指定しています。

enter image description here

CCM の結果は以下のように出力されます。実線と点がオリジナルデータを元にした cross-mapping です。グレーの領域は rEDM::make_surrogate_ebisuzaki()で生成したサロゲートデータ 100 本に基づいた 95% 信頼区間です。本来はもう少し多くのサロゲートデータに基づいて計算すべき、かつ、データによって使用するサロゲートを検討すべきですが、計算時間を短縮するためこのようにしています。

因果の判定は「最大ライブラリサイズでの予測精度 >> 最大ライブラリサイズでの 95% 信頼区間」&「最大ライブラリサイズでの予測精度 >> 最小ライブラリでの予測精度」に基づいています。この指標もちゃんと解析するときはもう少し検討すべきですが、簡略化しています。

ここでも図のタイトル部分に CCM の結果が表示されます。モデルデータでは、「X が Y に影響を与えている」というモデル式通りの結果が得られます。

enter image description here

最後に

「自分のデータに EDM を適用してみたいけどハードルが高そう」という方に、まずこのアプリで遊んでもらえば、と思って作成してみました。適用した結果、Simplex projection や S-map での予測精度がよく、CCM で因果がありそうな結果が出れば、EDM の適用を本格的に考えてもいいと思います。

そのような場合に一人では解析に自信がない、ということであれば連絡いただければ何らかの助けができるかもしれません。

また、このアプリに何らかのバグがあれば教えていただけるとありがたいです。時間が取れるときがあれば修正したいと思います。

Shiny でのアプリ作成、2作目ですが、前回より少し分かってきました。面白いです。まだもやもやしている部分はありますが、またアイデアが出てきたら何か作ってみたいです。

参考情報

アプリ内部で使用したパッケージのバージョン

  • R 0.4.3
  • shiny: v1.5.0
  • rEDM: v0.7.5
  • tidyverse: v1.3.0
  • cowplot: v1.1.0
  • ggsci: v2.9

関連ウェブサイト

参考文献

  • Sugihara, G. et al. (2012) Detecting causality in complex ecosystems. Science 338, 496–500
  • Sugihara G, May RM. (1990) Nonlinear forecasting as a way of distinguishing chaos from measurement error in time series. Nature 344, 734–741
  • Sugihara G. (1994) Nonlinear forecasting for the classification of natural time series. Philos. Trans. R. Soc. A 348, 477–495
  • Ye H, et al. (2018) rEDM: Applications of Empirical Dynamic Modeling from Time Series. doi:10.5281/zenodo.1935847

Written with StackEdit.

Empirical Dynamic Modeling with rEDM: (2) Near-future forecasting (Simplex projection)

20231015_Blogger0021 *This is an English version of my previous post (first translated by ChatGPT, checked and edited ...