多次元時系列分析ツールを公開しました

先日我々の論文がPhilosophical Transactions of the Royal Society Aから出版されました。

T. Ezaki, T. Watanabe, M. Ohzeki, N. Masuda,
"Energy landscape analysis of neuroimaging data", Phil. Trans. R. Soc. A 375, 20160287 (2017)

 

この論文はEnergy Landscape Analysis [1-3] という分析手法についてまとめた論文です。それに合わせて今回、誰でも(MATLABさえあれば)簡単にこの手法が使えるようにプログラムを公開しましたのでこの機会に詳しく書いてみたいと思います。

 

何ができる分析なのか?

多次元の時系列(例えば脳の様々な領域からとってきた活動データなど)を(単に個々の活動の総和としてではなく)全体として分析したほうが本質的なことがわかる場合があります。しかし、多次元の時系列は状態空間の次元が高く、何がどうなっているか理解するためには何らかの方法で低い次元での動きに読み直してやる必要があります。

 

これを解決する方法の一つがEnergy Landscape Analysis [1-3] です。この手法ではデータから尤もらしいモデル分布を推定し、それを元にして状態空間をエネルギー地形に見立てて粗視化することができます。

 

尚、脳のfMRIデータを例としていますが、(基本的には)多次元の時系列であればなんのデータにでも使えます。

 

手法の詳細

 

f:id:tEzk:20170710125035p:plain

 

この分析では、まず2値化した(活動が高い=+1,低い=–1とした)時系列を考えます。

 

一つの領域(変数)ごとに+1/–1で2パターンあるので、N個の領域を分析すると2^N個の活動パターンが発生しえます。データの中でそれぞれのパターンが何回現れたかをカウントすることで出現の確率分布を得ることができます。

 

次に、この確率分布をボルツマン分布でfittingします。Maximum Entropy Distributionとも言います。また、統計物理の文脈ではIsing(スピングラス)模型の定常確率分布です。このfittingですが、時系列の一次と二次のモーメント(=平均と相関)はデータとモデルでぴったり一致させることができます。ただ、Nが大きくなると厳密にこれを求めるのは現実的ではなくなり、近似手法に頼ることになります。ちなみに、三次以降の項は無視して大丈夫なのかというと、少なくともfMRIのデータでは(今のところ)大体よさそうです。

 

このモデル分布は式の定義の中に「エネルギー」が入っていて、フィッティングができた時点でそれぞれの活動パターンにひとつずつエネルギーの値が割り振られたことになります。

 

エネルギーが低い状態は出現頻度が高く、したがって重要な状態である可能性が高いです。そこで、それらの状態がエネルギー的にどのような位置関係になっているかを分析するとdisconnectivity graphと呼ばれるデンドログラムが書けます。

また、エネルギーの谷底(local minimum)になっているところの周り(basin)にどのような状態が集まっているのかも決定することができます。

 

こうして、すべての2^N個の状態にエネルギーの谷底(local minimum)かその周囲(basin)のラベルが付きます。すると、元の時系列が今そのEnergy Landscapeの中でどこにいるのかを定義することができ、このエネルギー地形を状態が動いていくという解釈をすることが可能となります。

 

 

公開したプログラム(Energy Landscape Analysis Toolkit: ELAT)

ダウンロードはこちらから

 

今回、この分析を簡単に利用してもらえるように時系列データを読み込むだけで全部自動的に出してくれるMATLABのプログラムコードを公開しました。ここではその使い方とできることについて紹介したいと思います。

 

準備

1.分析する時系列データを用意してください。現状、変数の数は20以下くらいでないとうまく分析できません。ファイルは数字だけからなるデータで各行が各変数の時系列になるようにします。つまり横(列)方向に時間が進んでいく向きで準備してください。連続値のデータの場合はELATが指定された閾値で二値化しますが、それが嫌な場合(データによってはもっと良い方法があるはずなので)はご自分で二値化した後のデータを準備します。

2.変数名リストのファイルを用意していただけると出力結果に反映されます。なくても動きます。

3.私のサイトからELATをダウンロードしていただき、適当なところに解凍してください。

 

ELATの使い方

1.MATLABを起動し、”StartProgram.m”を実行します。すると次のようなウィンドウが出ます。

f:id:tEzk:20170710134530p:plain

 

 

2.データファイルを読み込みます。テストデータとして"testdata.data"が入っているのでとりあえず動くか試してみたい方はそれを読み込んでください。

f:id:tEzk:20170710134800p:plain

 

3.テストデータは既に二値化されているのでData typeは"Binarized data"を選択します。

f:id:tEzk:20170710134934p:plain

連続値のデータを読み込んだ場合は、"Continuous data"の方を選択します。この時Thresholdが選択できますが、特段の理由がなければデフォルトの0のままで良いと思います。尚、fMRIデータでgroup内の被験者をひとまとめにpoolしたデータの場合、二値化は被験者ごとにやったほうが良い(base lineの値が違ったらこまる)ので、先に別途二値化したものを用意することをお勧めします。

 

 

4.変数名(ここでは、脳の領域の名前)のリストが(あれば)読み込みます。ここではテストデータ"roiname.dat"を読み込みます。

f:id:tEzk:20170710135457p:plain

 

 

5.最後に結果が保存されるフォルダを指定します。後の分析のためにbasinの定義(どの状態がどのエネルギーの谷に属するか)も保存したい場合は"Save Basin List"にチェックを入れます。

f:id:tEzk:20170710135724p:plain

 

6.Executeを押します。

 

 

結果

f:id:tEzk:20170710135926p:plain

上段の二枚は2^N個の状態をlocal minimumとそのbasinに分類した絵です。

数字はそれを二進数表示にして、1 → +1, 0 → –1と読み替えたときに出てくるパターンに対応します。色と高さはエネルギーの値を表します。

 

下段の右はdisconnectivity graphでlocal minimumになっている状態の間をエネルギー地形の上で行き来するときにどれだけエネルギー差があるかをまとめたものです。横軸は状態(定義は下段左)で縦軸はエネルギーの値を表します。

 

 

 

この後は何をすればいいのか?

この分析では、状態空間を粗視化したものとエネルギー地形的な解釈を与えますが、それだけで分析が終了ということではなく、それを足掛かりにしてより詳細な分析を行います。問題の対象や設定によって異なると思いますが、例えば

1.disconnectivity graphの形が条件AとBで異なっているかをみる

2.あるlocal minimumの周りでの状態の滞在時間や遷移レートなどの指標に着目してデータを分析する

などが考えられます(具体例については論文 [1] のイントロで紹介した文献をご参照ください)。

 

 

その他注意事項

1.このコードを実行するためにはMATLABの2016bよりも新しいバージョンおよび Statistics and Machine Learning Toolboxが必要です。

 

2.時系列の変数の数ですが20より小さいくらいでないと動きません。ただ、ボルツマン分布を求める部分については数百まで増やしても動きます(が、データによってパフォーマンスはまちまちなようです)。

 

3.同梱されているmain.mで一通りの分析が行えますので、適宜ご利用ください。

 

4.理論的なことや技術の詳細に関しては論文 [1](OAなので誰でも読めます)に書いてありますのでご参照ください。

 

5.バグや質問等ありましたら江崎でご連絡ください。

 

 

[1] T. Ezaki, T. Watanabe, M. Ohzeki, and N. Masuda, "Energy landscape analysis of neuroimaging data", Phil. Trans. R. Soc. A 375, 20160287 (2017).

[2] T. Watanabe, N. Masuda, F. Megumi, R. Kanai, G. Rees, "Energy Landscape and dynamics of brain activity during human bistable perception", Nat. Commun. 5, 4765 (2014).

[3] T. Watanabe, S. Hirose, H. Wada, Y. Imai, T. Machida, I Shirouzu, Y. Miyashita, N. Masuda, "Energy Landscapes of resting-state brain networks", Front. Neuroinform. 8, 12 (2014).

 

 

さきがけ

この度、「ボルツマンマシンを利用した脳の機能障害ダイナミクスの理解」という課題でJSTのさきがけに採用されました。最近更新が滞っていましたが、本プロジェクトに関する内容など含めポストしていきたいと思います。

 

新しい環境

四月からとあるプロジェクトの研究員として働くことになり、都内のとある国立研究所に異動となった。ようやく研究室の引っ越しや様々な手続きがひと段落して研究体制が整った。研究環境としては非常によく、かなり集中して研究に取り組めると思う。

 

以前書いていた論文は少々の紆余曲折ののち無事レビューに回り現在査読中となっている。良い知らせを待ちたい。

 

3月くらいからデータ分析に関する新しい研究を始めている。これはうまくはまれば相当なインパクトが出るような研究で、非常に楽しみである。一通りの(新しい)既存手法でできることに関しては何とかコードを実装でき、実際のデータに適用してみて結果がどうか、ということを見はじめたところである。とはいえ、まだ実装が正しくできているか不安な面もあり、現在は正しさをチェックしている段階である。

本日からかなり本腰で研究の作業が再開できた。共著者に質問を投げたり、一通りのコードを整理して足固めがだいぶ進んだ気がする。明日もこの調子で進めていきたい。

 

Closing

 前回の日記から2か月弱が経過したが、件の論文はほぼ投稿間近となった。記録によると12月7日くらいから原稿を書く作業を開始したらしいので、年末年始がなかったことを考えると二か月くらいかかっていることになる。その間原稿は20回程度アップデートされているので(無論、共著者の方のおかげだが)いままで自分が書いてきた原稿と比較すると比べ物にならないほど細かいところまで詰められているのではないかと思う。

 こういうふうにひとつのプロジェクトにどれだけ時間を使ったかを記録していくことは結構重要なのではないか、ということに最近思い至った。日単位ではなく、時間単位で簡単に記録したいと思い、togglというウェブサービスを使い始めてみた。これで論文が書きあがるまでにどれくらい時間がかかったのか後でみればわかり、今後の計画などに利用できるような気がする。

 やってみて気づいたのが、研究室にいる時間に比べて(思ったよりは)少ない時間しか使えていないということである。単純に時間の効率化という面でもいいのかもしれない。最近やり始めたばかりなので、1プロジェクトまるまるこれを使って記録したものが上がってくるのはもう少し先になるかもしれないが、辛抱強くつづけてみたい。

 

Manuscript

件の論文については原稿を前回の日記のあと12月7日あたりから書き始めた。そこから大幅な結果の増量などもあり、現在は数回やり取りしたところである。自分にとっては新しい分野のため参考文献の知識が乏しいのがつらい。ただ共著者の方が毎度正しい方向性を示してくださるのでその意味ではかなり効率よくレベルアップできている気がする。今年度中にもう一度アップデートして年末を迎えたい。

Cooperation

前回の記事から一週間、人間の協調行動についての研究がスタートした。モデルを使った数値計算の研究である。

 とりあえずあたりを付けていたモデル自体は初日に実装できたのだがいまいち実験結果と整合する結果が出ない。二つの量の変化を見ているのだが、概ね期待通りの振る舞いをしている一方で大小関係が逆転するところがどうしても出てしまうのだ。

 これを除去するのが一苦労であった。モデルを修正したりパラメータ範囲を変化させたりして逆転が出ないようにするとその他の(実験と整合するという意味で)望ましい性質も失われてしまった。結論から言うと元のモデルを形の上では変更せず、あるパラメータの解釈を変えて変化可能な領域を変えた。これで概ねこれまで報告されている(多少のバリエーションがある)実験結果をあるパラメータの範囲の中で再現できるところまでこぎつけた。

 明日以降は、論文の後半部分(内容としては追加的な部分、だがまあまあ作業量はある)をやる。

時間スケール

 先週、知り合いの研究者の方にいろいろ教えてもらい、(仮に狙い通りの結果が出れば)手っ取り早く新しいことができそうなテーマが見つかった(というかほぼ提供してもらった)のでそれをやっていた。結論からいうと仮説は誤りで、この方針で進めても大したcontributionにならないということであった。ただ、モデルのコーディングを始めてから2日でこの結論に至ることができたので(このあたりのことを勉強できたという意味では)時間当たりのコスパは良かったように思う。

 同じ方から別のテーマ(こちらの研究の狙いはそこまでリスキーではない、がインパクトは比較的高そう)を一緒にやらないかというお誘いを頂いたので今日はこのあたりの(教えてもらった)論文を幾つか読み、大体のイメージを掴んだ。わからないことや方針を即座に教えてもらえているので、ここのところ非常にテンポよく研究が進められている。この調子でじゃんじゃん論文を書いていきたい。