Liboschik, T, Fokianos, K, Fried, R. (2017) tscount: An R Package for Analysis of Count Time Series Following Generalized Linear Models. Journal of Statistical Software, 82(5), 1–51.
カウント時系列データ分析のRパッケージtscountの解説。このパッケージを実践投入する予定はないんだけど、カウント時系列の分析方法について知りたくて手に取った。
本文だけで34ページある…
1. イントロダクション
カウント時系列をどうモデル化するか。いくつかアプローチがある。
- GLM。カウントのための適切な分布とリンク関数を選ぶ。
- 整数ARMA(INARMA)モデル。thinning演算子に基づき、ARMAモデルの構造を模倣する。Waib(2008 AStA Advances in Statistical Analysis)というレビューがある。
- 状態空間モデル。Fokianos(2011 Statistics), Jung & Tremayne(2011 AStA Advances in Statistical Analysis), Fokianos(2012 Annals.ISM), Tjostheim(2012 TEST), Fokianos(2015 Chap.)をみよ。
GLMはINARMAに比べ、共変量の効果と負の相関を簡単に表現できるし、ツールキットが沢山ある。状態空間はモデルは柔軟だけどモデル指定が複雑。
というわけで、tscountパッケージはカウント時系列をGLMの枠組みで扱い、尤度ベースの方法を提供します。[…]
2. モデル
カウント時系列を\(Y_t\)、\(r\)次元共変量ベクトルを\(\mathbf{X}_t = (X_{t1}, \ldots, X_{tr})^\top\)とする。
条件付き平均\(E(Y_t|\mathcal{F}_{t-1})\)を過程\(\lambda_t\)としてモデル化する。\(\mathcal{F}_{t}\)というのは、同時過程\({Y_t, \lambda_t, X_{t+1}}\)の時点\(t\)までの履歴のことである(共変量については\(X_{t+1}\)を含んでいることに注意)。\(\mathcal{F}_{t-1}\)のもとで\(Y_t\)の条件付き分布がどうなるか、というのを仮定するわけです。
モデルを一般的に書こう。リンク関数を\(g: \mathbb{R}^+ \rightarrow \mathbb{R}\)、なんらかの変換関数を\(\tilde{g}: \mathbb{N}_0 \rightarrow \mathbb{R}\)として、$$ g(\lambda_t) = \beta_0 + \sum_{k=1}^p \beta_k \tilde{g}(Y_{t-i_k}) + \sum_{l=1}^q \alpha_l g(\lambda_{t-j_l}) + \eta^\top \mathbf{X}_t $$ GLMの用語に従い、\(\nu_t = g(\lambda_t)\)を線形予測子と呼ぶ。説明変数としてどの過去の観察を使うかを一般的に表すため、整数\(0 \lt i_1 \lt i_2 \lt \ldots \lt i_p \lt \infty\)を決めている。\(j\)についても同様。以下、\(P = \{i_1, \ldots, i_p\}, Q = \{j_1, \ldots, j_q\}\)とする。
[いきなり一般的に書きすぎだよう… 左辺\(\nu_t\)は観察変数\(Y_t\)の期待値 \( \lambda_t = E(Y_t|\mathcal{F}_{t-1}) \)をリンク関数で変換した連続潜在変数、ロジスティック回帰でいうところのlogit(p)だと思う。第二項は、季節性などを表現するために自由度の高い書き方になっているが、単純に\(\sum_k^p \beta_k \tilde{g}(Y_{t-k})\)と読み替えると、過去カウント(を実数に変換した変数)を\(p\)次ラグまで投入するということ。第三項も同様で、過去の\(\lambda_t\)を\(q\)次ラグまで投入している。第二項と第三項を両方入れているのは、たとえば左辺をポアソン分布のパラメータ(つまり期待値かつ分散)としてGARCHモデルをカバーしたいという発想であろう。最後のは時間変動共変量で、これはラグがない。えーっと… 書いてないけど、\(Y_t|\mathcal{F}_{t-1}\)の撹乱項はiidってことでいいんですよね？]
例を挙げよう。[こっちを先に書いてくれないと困るよね]
- \(g(x) = \tilde{g}(x) = x, P=\{1, \ldots, p\}, Q = \{1, \ldots, q\}, \eta = 0\)とすれば、$$ \lambda_t = \beta_0 + \sum_k^p \beta_k Y_{t-k} + \sum_l^q \alpha_l \lambda_{t-l}$$ 過去のもとでの\(Y_t\)の条件付き分布がポアソンだとすれば(後述)、これは整数値GARCHモデル INGARCH(p, q)である。自己回帰条件付きポアソン(ACP)モデルとも呼ばれている。
- \(g(x) = \log(x), \tilde{g}(x) = \log(x+1), P=\{1, \ldots, p\}, Q = \{1, \ldots, q\}, \eta = 0\)だとすれば対数線形モデルになる。\(v_t = \log(\lambda_t)\)として$$ v_t = \beta_0 + \sum_{k=1}^p \log(Y_{t-k} + 1) + \sum_{l=1}^q \alpha_l v_{t-l} $$ \(\tilde{g}(x) = \log(x+1)\)とするのは、線形予測子\(v_t\)とスケールを揃えるためである。時間依存共変量を入れたときに、整数値GARCHモデルでは加法的になるがこのモデルでは乗法的になる。詳細はFokianos & Tjostheim (2011 J.MultivariateAnal.) をみよ。
過去のもとでの\(Y_t\)の条件付き分布がポアソンだ、つまり\(Y_t | \mathcal{F}_{t-1} \sim Poisson(\lambda_t)\)だということは、つまり$$ P(Y_t = y | \mathcal{F}_{t-1}) = \frac{\lambda^y_t \exp(-\lambda_t)}{y!}, \ y = 0,1,\ldots $$ということである。この場合、\(Var(Y_t | \mathcal{F}_{t-1}) = E(Y_t|\mathcal{R}_{t-1}) = \lambda_t\)である。
いっぽう、条件付き分散は条件付き平均より大きくなることもある(overdispersion)。これを許すため、負の二項分布をつかって\( Y_t | \mathcal{F}_{t-1} \sim NegBin(\lambda_t, \phi)\)とする方法もある。すなわち $$ P(Y_t = y | \mathcal{F}_{t-1}) = \frac{\Gamma (\phi+y)}{\Gamma(y+1)\Gamma(\phi)} \left( \frac{\phi}{\phi+\lambda_t} \right)^\phi \left( \frac{\lambda_t}{\phi + \lambda_t} \right)^y$$ このとき\(Var(Y_t | \mathcal{F}_{t-1}) = \lambda_t + \lambda^2_t / \phi \)である。
実は、負の二項分布は混合ポアソン過程というクラスに属している。混合ポアソン過程とは、\(\{N_t\}\)を強度1のiidなポアソン過程、\(\{Z_t\}\)を\(\{Y_t\}\)と独立なiidな確率変数(平均1, 分散\(\sigma^2\)として、\( Y_T = N_t (0, Z_t \lambda_t] \)という過程である。
[← なにこの表記！？ 意味がわかんないよ！ ここで躓いてしまい、この段落の内容自体が理解できなくなってしまったので、段落ごとスキップする。読んでないので推測だけど、負の二項分布をポアソン分布のガンマ混合とみるという話じゃないかと思う]
上記のモデルでは、共変量は過去の観察への回帰を通じて将来の観察に影響するし、過去の条件付き平均への回帰を通じても将来の観察に影響する。これを内的共変量効果と呼ぶことにする。
それに対して、過去の観察への回帰のみを通じて将来の観察に影響する、という定式化もありうる。これを外的共変量効果と呼ぶことにする。
両方を含めて定式化するとこうなる。共変量が効くかどうかを表す二値ベクトルを\(\mathbf{e}\)として、$$ g(\lambda_t) = \beta_0 + \sum_{k=1}^@ \beta_k \tilde{g}(Y_{t-i_k}) + \sum_{l=1}^q \left( g(\lambda_{t-j_l}) – \eta^\top diag(\mathbf{e}) \mathbf{X}_{t-j_l}\right) + \eta^\top \mathbf{X}_t $$ 内的効果と外的効果の比較については Liboschik et al.(2016 Int.J.ComputerMath.)をみよ。経験的には、どっちにしても大差ない。[なあんだ]
3. 推定と推論
[パス。ややこしくて付き合いきれない…]
4. 予測
MSEの観点からいえば、最適な一期先予測\(\hat{Y}_{n+1}\)は\(\lambda_{n+1}\)の条件付き期待値である。条件付き分布はモデル次第で、ポアソンだったりNBDだったりする。
\(n\)期先予測はその再帰予測である。条件付き分布は解析的には未知なので、パラメトリック・ブートストラップで近似する。
実際には\(\lambda_{n+1}\)の代わりに、回帰パラメータを\(\hat{\theta}\)として\(\hat{\lambda}_{n+1} = \lambda_{n+1}(\hat{\theta})\)を使うわけだし、NBDだったら\(\phi\)の代わりに\(\hat{\phi}\)を使うわけで、そのせいで予測分布に不確実性が加わる。以下では予測区間の構築について説明するけど、このパラメータ推定の不確実性は考慮していない。
予測区間の構築には２通りのやりかたがある。(1)近似された予測分布の\(\alpha/2\)分位点と\((1-\alpha/2)\)分位点をとる。(2)予測分布からカバレッジ率が少なくとも\(1-\alpha\)になるような最短の区間を選ぶ。どっちを使っても結果はそうそう変わらない。
5. モデル評価
観察にあてはめた値\(\hat{\lambda} = \lambda_t(\hat{\theta})\)は分布に依存しない。残差にはいろんな定義がありうる。tscountパッケージは以下に対応しています。
- response: \(r_t = y_t – \hat{\lambda}_t\)
- pearson: \(r^P_t = (y_t – \hat{\lambda}_t) / \sqrt{ \hat{\lambda}_t – \hat{\lambda}^2_t \hat{\sigma}^2} \)
- anscombe. 標準化Anscombe残差というやつで、ピアソン残差よりもより対称的になる。$$ r^A_t = \frac{ \frac{3}{\hat{\sigma}^2} \left( (1+y_t \hat{\sigma}^2)^{2/3} – (1+ \hat{\lambda}_t \hat{\sigma}^2)^{2/3} \right) + 3(y^{2/3}_t – \hat{\lambda}^{2/3}_t) }{2 (\hat{\lambda}_t + \hat{\lambda}^2_t \hat{\sigma}^2)^{1/6}}$$ [ここまで汚い式というのもそうそうお目にかかれないと思って、つい写経してしまった。なんなのこれ]
まずは経験ACFをみるのをお勧めする。それから時間-残差プロットで、どっかでDGPが変わってないかどうかチェック。\(r^2_t\)と\(\hat{\lambda}_t\)の散布図もチェック。45度線の周りにあったらポアソン、二次曲線っぽかったらNBDが示唆される。[あっそうか、なるほど。\(\hat{\lambda}_t\)は分布の仮定に依存してないからね]
カウント時系列の予測性能を評価するため、次のようなツールを御用意しています。以下、CDFを\(P_t(y) = P(Y_t \leq y | \mathcal{F}_{t-1}) \)、PDFを\(p_t(y) = P(Y_t = y | \mathcal{F}_{t-1})\)、予測分布(ポアソン分布だかNBDだか)のSDを\(v_t = \sqrt{Var(Y_t|\mathcal{F}_{t-1}} \)。
- pit. prequential principle という考え方があって、それによれば、予測分布が正しいとき、probability untegral transform(PIT)というのが一様分布に従う。カウントデータの非ランダム化PITは次のように求められる: $$ F_t(u|y) = \begin{cases} 0 & u \leq P_t(y-1) \\ \frac{u-P_t(y-1)}{P_t(y) – P_t(y-1)} &
P_t(y-1) \lt u \lt P_t(y) \\ 1 & u \geq P_t(y) \end{cases} $$ これを\(t\)を通じて平均して\(\bar{F}(u)\)とする。で、\(u=0, 0.1, \ldots, 1\)について求めて差をとって、10ビンのヒストグラムを書く。それがU字型になったら予測分布のunderdispersion, 逆U字型になったらoverdispersionである。[ひいいいい。なにいってんだかぜんぜんわかんない…!] - marcal. marginal calibrationの略。[このパート、ほぼ1パラグラフしかないけど、なんだかめんどくさそうなのでスキップ]
- scoring. [この話、面白いので細かくメモする]
Gneiting et al.(2007 JRSS-B)によれば、PITヒストグラムやmarginal calibration plotは理想的予測の必要条件ではあるが十分条件ではない。彼らはかわりに、十分にカリブレートされたすべてのモデルのうち、シャープネスが最大であるモデルを採用することを勧めている。シャープネスとは、予測分布が中心に集まっている程度のことで、予測区間の幅で測れる。カリブレーションとシャープネスを同時に評価した単一の数値的スコアを、プロパー・スコアリング・ルールによって実現できる。
予測分布\(P_t\)、観察\(y_t\)に対するスコアを\(s(P_t, y_t)\)とする。プロパー・スコアリング・ルールの定義はいろいろあり、表1に7種類を示す。たとえば二乗誤差スコアは\((y_t – \lambda_t)^2\)、基準化二乗誤差スコアは\((y_t – \lambda_t)^2 / v^2_t\)、対数スコアは\(-\log(p_t(y_t))\)である。二乗誤差スコア以外の6つは分布に依存する。
これを\(t\)を通じて平均し、低いモデルを選ぶわけである(基準化二乗誤差スコアの場合は1に近いモデルを選ぶ)。
[なるほど… こういうのもproper scoring ruleっていうんだ…] - AIC, BIC, QIC。[パス]
6. 介入分析
[ここ、もっとも関心あるパートなので細かくメモする]
Folianos & Fried(2010 J.TimeSeriesAnal., 2012 Stat.Modelling)は、位置に影響する介入をモデル化するために決定論的な共変量を導入している。その形式は、発生時点を\(\tau\)、既知の減衰率を\(\delta\)としてて$$ \delta^{t-\tau} \mathbf{1}(t \geq \tau_m) $$ である。\(\delta = 0\)ならスパイク、1までなら指数減衰、1なら永続的なレベルシフトになる。tscountパッケージではこの形式を表す関数 interv_covariateを用意している。
共変量の場合と同じで、介入の効果というのは線形モデルなら加法的だし対数線形モデルなら乗法的である。しかし、介入は過程のダイナミクスに入ってくるので、線形予測子に対するその効果は純粋には線形でない[←???]。こうした介入効果について検定するためのツールもご用意している。
介入が\(s\)タイプあり、それぞれのサイズを\(\omega_1, \ldots, \omega_s\)として、$$ g(\lambda_t) = \beta_0 + \sum_{k=1}^p \beta_k \tilde{g}(Y_{t-i_k}) + \sum_{l=1}^q \alpha_l g(\lambda_{t-j_l}) + \eta^top \mathbf{X}_t + \sum_{m=1}^s \delta^{t-\tau_m}_m \mathbf{1}(t \geq \tau_m)$$ とモデル化する。線形モデルならレベルが\(\omega_m\)だけ増え、対数線形モデルならレベルが\(\exp(\omega_m)\)倍になる。
- iterv_test関数。\(H_0: \omega_1 = \ldots = \omega_s = 0\), \(H_1:\) なんらかの\(l \in {1, \ldots, s}\)について\(w_l \leq 0\) という検定を行う。いくつかの条件を仮定すると、\(H_0\)のもとで、スコア検定統計量\(T_n(\tau_1, \ldots, \tau_s\)は近似的に自由度\(s\)のカイ二乗分布に従う。
- interv_detect関数。未知の時点\(\tau\)の下で、あるタイプの単一の介入が効果を持ったかどうかを検定する。\(T_n(\tau)\)の最大値を求め、パラメトリック・ブートストラップで\(p\)値を得る。介入が起こりうる時点の集合を指定できる。もっとも、こういう推定量は変動性が大きい。なお、本来はそれぞれのブートストラップ標本からモデル・パラメータを推定するのだが、先にモデルパラメータを決めてからブートストラップ表ホウンを推定してブートストラップ検定統計量を求めるという方法もある。計算は速いけど保守的になる。
- interv_multiple関数。複数の介入が疑われるが、時点もタイプもわからない、という場合に使える。介入がありうる時点の集合を\(D\)、ありうる介入タイプの集合(たとえば{0, 0.8, 1})を\(\Delta\)とする。まず、iterv_detect関数でそれぞれのタイプの介入について検定しボンフェローニ調整する。介入が検出されたら、\(p\)が最小である介入タイプを採用して、その効果を除外する。これを繰り返す。最終的には複数タイプの介入が複数時点で行われているというモデルが手に入るわけだが、注意して扱うこと。[そうねえ… ものすごいバイアスがありそうだわね]
現実場面では、減衰率\(\delta\)も推定したいことが多い。しかし\(\omega\)がゼロなら\(\delta\)は識別できない。この場合、このパラメータを通じて尤度をプロファイリングすることで推定できるだろう。たとえば介入効果が単一だったとする。\(\delta\)をなにかの値に決め、他のすべてのパラメータについて(擬似)最尤推定量を得る。これをすべての\(\delta\)について繰り返し、対数尤度が最大となる値を得る。もっともこの方法は、他のパラメータについての推論の妥当性にも影響してしまう。
7. パッケージの使用方法
[2つのデータ例を使って例示している。読まずに飛ばしたけど、二番目の例は、Seatbelts … おまえはSeatbeltsデータではないか！ 状態空間モデルについて勉強しているときに散々使ったやつだ。あのときはカウントの対数を目的変数にしていたけれど、確かに、本来はカウントでしたね…]
8. 他のパッケージとの比較
8.1 独立なデータのためのパッケージ
時系列データの場合でも、\(Y_{t-i}\)を共変量に加えて独立データのパッケージを使うということもできなくはない。ただし、定常性を保証するための制約とか、確率的季節性とかは扱えない。
- stat::glm(), MASS::glm.nb()。[メモ省略]
- gamlssパッケージ。位置・スケール・形状についての一般化加法モデル(GAMLSS)。
- ZIMパッケージ。ゼロ過剰モデル。
- pscl::zeroinfl(), pscl::hurdle()。
- VGAMパッケージ。ベクトルGLM。vglm()は一変量ならglm()と同じ。vgamという関数もある(ベクトルGAM)。Yee (2015, 書籍)をみよ。
まとめると、基本的にカウント時系列はこういうパッケージではなくてtscountパッケージでモデル化したほうがよいが、ゼロ過剰なときはZIMパッケージなどを検討してもよい。overdispersionが時変する場合はgamlss、多変量カウント時系列はVGAMが選択肢になる。
8.2 時系列データのためのパッケージ
- acpパッケージ。自己回帰条件付きポアソン回帰モデル(つまりINGARCHモデル)の最尤推定ができる。tscountと比べると、tscountのほうがパラメータ制約ができるし(3章で説明)、NBDモデルとか対数線形モデルとかが使えるし、確率的季節性も扱えるし(7章で説明)、外的共変量・内的共変量を扱える。
- glarmaパッケージ。一般化線形ARMA(GLARMA)モデル。tscountと比べると、glarmaは対数リンクのみで…[ほかにいろいろ細かいことが書いてあるけど、読み飛ばした]
- gamlss.utilパッケージ。一般化ARMA(GARMA)モデル。GLARMAとどうちがうかというと…[疲れちゃったのでパス]
- VGAM::garma()もGARMAモデルを推定できるが、まだ開発途上でいろいろ制限がある。
ここまでのモデルは、ランダムネスのソースがひとつだけだった(たとえばポアソン分布からの観察のドローに伴うランダムネス)。ここからは、イノベーション過程が複数あるモデルについて紹介する。モデルは柔軟になるが予測においては不利になる。詳細はKFASパッケージの解説をみよ。
- surveillance::hhh()。カウント時系列の変化点のオンライン検出ができる。モデルは$$ Y_t \sim NegBin(\lambda_t, \phi)$$ $$ \lambda_t = \exp(\beta_0 + \delta t + \gamma_t) + \beta_1 Y_{t-1}$$ ただし、\(\gamma_t\)は季節効果。\(\exp\)がトレンド・季節にのみかかって自己回帰要素にはかかってない。GLMではなくてGAMとみるべきである。
- surveillance::hts()。階層時系列モデル(hts)。モデルは$$ \lambda_t = \exp(\beta_{0,t} + \delta t + \gamma_t + \eta^\top \mathbf{X}_t) $$ $$ \Delta_d \beta_{0,t} | \beta_{0, t-1}, \ldots, \beta_{0, t-d} \sim N(0, \kappa^{-1}_{\beta_0})$$ [えええ？ 切片項\(\beta_{0,t}\)の\(d\)次差分がWNってこと？ なぜそんな発想を？] 推論はINRAパッケージでなされる。ベイジアンアプローチとtscountの尤度アプローチを比べると、予測区間が推定の不確実性を表現できるという点ではベイジアンアプローチが優れているが、計算時間がかかるのが難点。一般的に言うと、一期先予測ならGLMベースのほうがよくて、階層モデル予測ならinlaのほうがよい。この点についてはさらなる研究が必要。
- KFASパッケージ。指数分布族の多変量時系列の状態空間モデル。INLAと違って最尤推定。[… きちんと読んでないけど、tscountに近い模様]
- pompパッケージ。部分観察マルコフ過程(POMP)という状態空間モデルを使う。たとえば、$$ Y_t = Poisson(\phi N_t) $$ $$ N_{t+1} + r N_t \exp(-N_t + \epsilon_t)$$ $$ \epsilon_t \sim N(0, \sigma^2_\epsilon$$ とする。パッケージがやたらと使いにくいのが難点。
- gcmrパッケージ。Gaussian copula marginal regressionという枠組みによるパッケージで、カウント時系列もモデル化できる。共変量ベクトル\(\mathbf{X}_t\)のもとでのカウント時系列\(Y_t\)の周辺分布を、平均\(\lambda_t\)のポアソン分布なりNBDなりだとして、\(g(\lambda_t) = \beta_0 + \eta^\top \mathbf{X}_t\)とする(ここには時間依存性がはいらない)。で、未観察の誤差過程\(\{\epsilon_t\}\)を考え、ポアソン分布なりNBDなりのCDFを\(F^{-1}_t\), 標準正規分布のCDFを\(\Psi\)として$$ Y_t = F^{-1}_t(\Psi(\epsilon))$$とする。つまり、\(Y_t\)はポアソン分布なりNBDなりの分位点なのである。時間依存性は\(\{\epsilon_t\}\)のARMAモデルで実現する。[ひえええええ。なんでこんなわけわかんないことを考えるの？]
- ZIMパッケージにも、いわゆる動的ゼロ過剰モデルというのがあるね。
- tsintermittentパッケージ。intermittent demand time seriesという手法を提供している。[いろいろ書いてあるけど、疲れちゃったのでパス]
9. 今後の課題
- NBDモデルの代替案としてQuasi-Poissonモデルが考えられる。NBDでは条件付き分散を\(\lambda_t + \phi \lambda^2_t\)とするが、かわりに\(\phi \lambda_t\)とするモデルである。
- INGARCH時系列モデルでゼロ過剰を扱うという試みもある。
- イベントの回数ではなく、単位当たりイベント数に関心があることもある。そういうときは率について対数リンクを使うことができる。分母の対数がいわゆるオフセットになる。tscountパッケージでも今後はオフセットを実装してもいいかも。
- 閾値モデルというのもある。
- 外れ値への対処。
- 長期的には、二値・カテゴリカル時系列や多変量カウント時系列への対応も考えられる。
- 全然系統がちがうけど、thinningベースの時系列モデルとか、INLA, KFAS, pompへのカウント時系列用インタフェイスとか。
——————
途中で疲れ切ってしまった。正直、舐めてました。単に「ARMAモデルのリンク関数をポアソンにするんでしょ。かんたんじゃん」と思っていたのである。この話、滅茶苦茶難しいじゃんか… どうすんだよ…