« 読了: Petris & Petrone (2011), Petris (2010) dlmパッケージとそのライバルたち | メイン | 読了：「インドでキャバクラ始めました(笑)」「かごめかごめ」「甘々と稲妻」「イノサン」「36歳無職さん」「孤独のグルメ」「木曜日のフルット」「それでも街は廻っている」 »
2014年10月24日 (金)
Holmes, E.E., Ward, E.J., Scheuerell, M.D. (2014) Analysis of multivariate time-series using the MARSS package. version 3.9. Northwest Fisheries Science Center, Seattle, WA.
RのMARSSパッケージのユーザーズ・ガイドに相当する文書で、３部構成、全16章、200頁以上に及ぶ。MARSSパッケージとは要するに、多変量時系列の背後に少数の自己回帰系列を考える状態空間モデルのソフトで、いわゆる動的回帰やベクトル自己回帰を扱うことができる。パラメータ推定を著者らが提案している一種のEMアルゴリズムで行うのが特色。
第一部。
1章はイントロ。えーっと、MARSSのモデルは下記の通り。
状態方程式: $x_t = B_t x_{t-1} + u_t + C_t c_t + w_t$
状態方程式の撹乱項: $w_t \sim MVN(0, Q_t)$
観察方程式: $y_t = Z_t x_t + a_t + D_t d_t + v_t$
観察方程式の撹乱項: $v_t \sim MVN(0, R_t)$
初期状態: $x_1 \sim MVN(\pi, \Lambda)$ ないし $x_0 \sim MVN(\pi, \Lambda)$
時系列の長さを$T$, 状態変数を$m$本, 観察変数を$n$本とする。$c_t$, $d_t$は外生変数で、それぞれ$p$本、$q$本。観察変数$y_t$には欠損を許すが外生変数$c_t$, $d_t$には許さない。
章末に他のソフトの紹介が載っている。著者ら曰く、MARSSパッケージは速度を最適化してない、特に時点数が多いときには堪忍な、とのこと。
2章、関数紹介。主役はMARSS()関数である。3章は著者らが開発したEMアルゴリズムの説明(パス！)。
第二部。
4章がたぶんいちばん大事な章で、MARSS()に対するデータとモデルの渡し方。このパッケージの特徴として、渡すものはすべてユーザの責任で正しくつくらないといけない。
えーっと... データはn行T列の横長な行列で渡す。モデルは、dlmパッケージのdlmMod*()やKFASパッケージのSSM*()のようなヘルパー関数はなくて、自分でシステム行列を書きまくりリストにして渡す必要がある。モデルの要素は行列なのだが、なかに数値と文字列を混在させたいので、list()をmatrix()で並べるのが望ましい。つまり、たとえばmatrix(list(1, "a", "b"), 3, 1)ってな風に書けってことである。
MARSS()に渡すモデルにいれられる要素は以下の通り。初期状態 x0, V0を除き、すべて時間変動可能であって、そのときはarray()で渡すのだが(5.3節)、ややこしいので省略。ここでは時間不変の場合についてメモする。
- U: 状態に乗せる傾き$u_t$。m行1列の行列で渡す。たとえばU=matrix(list(0.01, "u1", "u2"), 3, 1)とすると、$x_t$の一本目には毎時点で0.01が, 二本目には毎時点で定数u1が、三本目には毎時点で定数u2が加算されることになる。ほかに、U="unequal", U="equal", U="zero"というショートカット表現がある。デフォルトはU="unequal"。
- A: 観察に乗せる傾き$a_t$。Uと同様に指定。さらに、A="scaling"というショートカットがある。こうすると、$a_t$の要素は通常は0となり、同一の状態変数に対して係数を持つ観察変数が複数個存在するときに限り、2本目以降の観察変数での要素が定数a1, a2, ... になる。デフォルトはA="scaling"。
- x0: 初期状態の平均$\pi$。Uと同様に指定。面白いことに、ここに変数名をつけてパラメータ扱いすることができるのである(つまり、初期状態の事前分布次第で結果が変わってきちゃうよね、という議論を回避できるわけだ)。デフォルトはよくわかんないんだけど(4.5節にはx0じゃなくてpiって書いてある...)、たぶんx0="unequal"だと思う。
- Q: 状態撹乱項の分散行列$Q_t$。m行m列の行列を渡す。気をつけないといけないのは、Q=diag(c("sa", "sa", "sb"))というような書き方をしてはいけないということ。Q=matrix(list(0), 3, 3)とやって、あとでdiag(Q)=c("sa", "sa", "sb")とするのはオッケー。ショートカットは"diagonal and equal", "diagonal and unequal", "unconstrained", "equalvarconv", "identity", "zero"。デフォルトは"diagonal and unequal"。
- R: 観察撹乱項の分散行列$R_t$。Qと同様。
- V0: 初期状態の分散行列$\Lambda$。Qと同様。V0の指定には細心の注意を払え、さもないと破壊的結果を招くぞ、とのこと。デフォルトは"zero", これは未知を表す由(えええ... それは教えてくれないとわかんないね...)
- B: 状態遷移行列$B$。Qと同様。数値をいれる場合、すべての固有値の絶対値が1以下になるようにしとかないといけない。デフォルトは"identity"。
- Z: 観察方程式における状態変数への係数行列$Z$。要するにこれはn行m列の行列で、頑張って作るしかないのだが、以下の場合は手が抜ける。(1)状態変数の数$m$と観察変数の数$n$が同じ場合。Qと同様の書式が使える。(2)$Z$はデザイン行列だという場合(つまり、要素は0か1で縦に合計すると必ず1だという場合)。このときは、factor()を使って謎めいた指定ができるのだが、さっぱり理解できないので省略。
- C, c, D, d: 4章には記載がないが、まあ使う段になればどうにかなるだろう。
- tinitx: 初期状態をt=0とするかt=1とするか。デフォルトはtinitx=1。
5章は、簡単なコード例(5.1節)、状態変数の数が観察変数の数と異なるコード例(5.2節)、時間変動パラメータのコード例(5.3節)、共変量$C, c, D, d$についての短い説明(5.4節)、結果のみかた(5.5節)、信頼区間の求め方(5.6節)、推定結果の取り出し方(5.7-5.9節)、パラメータのブートストラップ推定(5.10節)、初期状態をランダムに与えてモンテカルロ法で初期化する方法(5.11節)、シミュレーション用のデータ作成(5.12節)、ブートストラップAIC(5.13節)、収束条件について(5.14節)。
第三部。ここからは事例紹介の連続爆撃である。どの例も生態学(?)の話なので、ちょっと覚悟が必要だ。ざっとめくって、タイトルをメモしておくと、
- 6章: curruptedデータを使ったcountベースのpopulation viability分析。最初からナンノコッチャという感じだが、なにかの生き物の個体数を毎年数えているようなデータがあって、絶滅しないかどうか調べる、というような話らしい。モデルをみたら、単変量のローカル線形トレンドモデルで傾きが定数の奴だった。ちゃんと読んだらすごく勉強になりそう。読まないけど。
- 7章: 複数のsiteのデータを結合し、地域の個体数トレンドを推定する(populationの訳ってここでは個体数でいいのかしらん？)。観察変数は5本。状態変数はひとつからはじめて、だんだん複雑にしていく。これも面白そう... 読まないけど。
- 8章: 空間のpopulation structureと共分散を同定する。なにいってんだかわかんないけど、地域別にアザラシの個体数の時系列があって、地域がどう固まっているか、というようなことを調べたいらしい。これはめんどくさそうだ。ぜったい読まない。
- 9章: 動的因子分析(DFA)。正直なところこれが目当てで読み始めたのだが、きちんと読むには時間がかかるので(ここだけで18頁ある)、先に全体のメモを取っている次第である。
- 10章: ノイズの多い動物トラッキングデータの分析。えーっと、ビッグ・ママという愛称のウミガメの位置の時系列の分析らしい(観察変数は2本、状態変数も2本)。あー、そりゃ面白いな、ランダム・ウォークというか、ランダム・スイミングというか。読まないけどさ。
- 11章: 外れ値と構造的ブレークの検出。データはご存知ナイル川データ。これは仕事とも関係するので、こんどきちんと読もう。
- 12章: 共変量のとりこみ。季節効果の推定もこの枠組みで扱うわけで、いかん、ここは読まざるをえない...
- 13章: 種の相互作用の強さの推定。オオカミとヘラジカの個体数の時系列を分析する。ああそうか、多変量自己回帰(MAR)で扱えるわけだ、なるほどねー。こういうの、いつか仕事で使えるような気がするが、とりあえず後回し。
- 14章: 複数の時系列の統合。よくわかんないんだけど、ある種の個体数を陸から数えた時系列と空から数えた時系列があって...というような話らしい。一因子のDFAみたいなもんだろか。
- 15章: 単変量動的回帰モデル(DLM)。係数が時間変動する話なんでしょうね、読んでないけど。
- 16章: MARSSによるラグ-pモデル。AR(2)とかMAR(2)とかをどう扱うかという話らしい。
というわけで、これからどの章をきちんと読まねばならんか見当がついたので、とりあえず良しとしよう。疲れたし。
論文：データ解析(-2014) - 読了：Holmes, Ward, & Scheuerell (2014) MARSSパッケージで多変量時系列分析