どうも、スミスです。

先日、統計数理研究所公開講座、LeadingDAT養成コースの修了証が送られてきました。

めちゃくちゃいかつい感じでびっくり...！この講座をやり遂げた重みを感じました。年末年始ひたすらレポートに取り組んでいたので嬉しいです。
今回は、その復習も兼ねて、新型コロナウイルスの感染者数について状態空間モデルによるモデリングを実装していきたいと思います。今回は、理論はほどほどに、世にあるデータに対して状態空間モデルを適用させるプロセスや、その結果に基づく考察を提示することで、「実務に役立てる」ことを意識していきたいと思います。

1. 状態空間モデルとは

1-1. 状態空間モデルの構造

実装の前に、状態空間モデルの要点をまとめます。
状態空間モデルは、時系列データを扱う階層モデルの一種です。まず、階層モデルについて簡単に述べておきます。
階層モデルでは、各観測値に対して、それを生み出す確率分布のパラメータ（潜在変数）を規定します。そして、潜在変数を規定するパラメータと、各観測値を規定するパラメータとして機能する潜在変数の両者を推論します。
この構造をベースに、各潜在変数間に依存関係を仮定したものが状態空間モデルの構造になります（下図）。
このように、各観測値を前の時点での状態に基づいて規定される状態のパラメータのもとでモデリングをすることができる点で時系列データへの適用にもってこいなわけです。
さらに、観測不可能な「状態」を潜在変数として切り出すことができるため、解釈性の高いモデリングが可能であったり、欠測値に対しての扱いも簡単であったりと言うメリットもあります。

1-2. 基本的な状態空間モデル

続いて、上記で述べた構造をベースにした基本的な状態空間モデルであるローカルレベルモデルを定式化して解説します。
下図上段の通り、観測値は状態αと計測誤差ε（観測値撹乱項）でモデリングされます。状態αも下段の通り、前時点のαとその変動要素（状態撹乱項）で定義されます。各撹乱項は平均0の正規分布を仮定し、分散を未知パラメータとして、その推論によりαと観測値を求めることになります。

1-3. 状態モデルの拡張

続いて、1-2で述べたローカルレベルモデルから、今回の実装に用いる拡張モデルを紹介します。ローカルレベルモデルは、時系列データの増減トレンドについて考慮できず、長期的な予測を加味した場合に精度の悪い予測をするモデルを構築してしまう場合があります。そこで、下記に述べるように状態モデルを拡張し、予測に対して精度を保ったモデルを構築することができます。

・平滑化トレンドモデル

イメージとしては、先ほどのような階層的なモデルを状態モデルに対して適用する、というものです。前時点での状態αと、増減トレンドを加味するαの傾き成分（同様に前時点の傾き成分と、撹乱項ηで構成）に分解することで、各時点での増減の様子を調節することができるようにしたモデルです。

・ローカル線形トレンドモデル

こちらは平滑化トレンドモデルで定義した水準成分にも撹乱項を加え、より柔軟に状態の変動を表現できるようにしたモデルです。

2. KFASを用いた状態空間モデルによる新型コロナウイルス感染者数のモデリング

いよいよ、状態空間モデルを実際のデータに適用してその挙動や解釈を見ていきたいと思います。
今回は、発生から1年が経ち、医療やライフスタイルなど様々な部分で大きな影響を与え続けている新型コロナウイルス感染者数のデータを用いて、その推移をモデリングしてみます。

今回の解析は、Rの状態空間モデルパッケージ「KFAS」を使用します。

また、今回の実装に使用したソースコードについてはgithubに格納していますので、詳しくは下記ご覧ください。

モデリングの準備

今回使用するデータは、厚生労働省の公開する新型コロナウイルス陽性者数の日次データ（期間：2020.2.5〜2021.1.7）です。

データの内容は年月日と各日の陽性者数のみです。これらを週次で集計したもので今回の解析を行っています。
下のグラフは、週次でのコロナウイルス陽性者数推移です。状態空間モデルでは、モデリング前のデータの様子から、仮説を持っておくのがいいと思っています。
今回のデータでは、①5月手前付近、②8月前後付近、そして③11月下旬以降の3つの山が確認できるかと思います。時期的に明らかかと思いますが、①は1回目の緊急事態宣言下での観測値で、第1波の感染増期と見て良さそうです。同様に、②は緊急事態宣言が明けて少しした夏場に感染が増えていた時期、そして③は冬に突入し、2回目の緊急事態宣言までの感染爆発が起こっていた時期、と見られます。

モデリング結果の比較

いよいよモデリングした結果を見てみましょう。1で書いた状態空間モデルの3つのモデル（ローカルレベルモデル・平滑下トレンドモデル・ローカル線形トレンドモデル）に、月ごとの季節変動を固定・可変で調節した計6パターンのモデルを作成しました。
これらのモデルについて、「どれがいいか」を見るためにAICと最大対数尤度の比較をしてみます。
各モデルのあてはまりの良さは下表の通りです。

AICベースで見ると、「平滑化トレンドモデル・固定季節変動」が最も当てはまりのいいモデルになっていますね。つまり、過適合を避け、最適な予測値を与えてくれるモデルは、増減トレンドを考慮し、かつ月ごとに固定の季節変動を加えた形であるということになりましょうか。
一方、最大対数尤度ベースで見ると、「ローカル線形トレンドモデル・可変季節変動」が最も当てはまりのいいモデルになるようです。つまり、現在のデータに対しては、増減トレンドをより柔軟に考慮した上で、月ごとの変動を考慮した季節変動を加えた形のモデルが最適なあてはまりを示すということになりましょうか。
2つの評価指標は、前回の記事でも書かせていただきましたが、モデリングの対象となった事象の現状を微細に表現するモデルと、将来の予測を視野におき、トレンドの変化などに対応できる形で現実から多少遠い表現をするモデルという、それぞれの指標の考え方に沿ったモデルを最適と評価しているようですね。今回のモデリングでも分かりやすく現れているかなと思います。

今回実装したモデルは、まだまだ改善の余地が多くあるかと思います。
今私が思いつくだけでも、
・予測精度の検証（ちょうど2021年1月8日以降のデータが1ヶ月ほどたまっているかと思うので、そのデータで検証できますね）
・季節変動の周期を調節する（今回1ヶ月ごとの変動として定義しましたが、例えば3ヶ月ごとにして四季を周期とする変動にしてみるなどが考えられるかと思います）
・短期変動を別途切り出す（2020年11月以降はそれまでと大きく異なる増加トレンドが現れていたので、その部分を別途短期変動としてモデルに明示的に加えることもできるかと思います）
などが挙げられます。その他より良いアプローチなどありましたら、是非教えてください！！

今回は少し長くなりました。。。
仕事が少し繁忙してきましたが、継続的にアウトプットできるようやりくりしていきたいと思います。
ではまた。