dlmによる時変係数モデル

最終更新:2017年6月1日

dlmパッケージを使って、「傾きが時間によって変化する回帰モデル」を推定します。

Stanを使って「傾きが時間によって変化する回帰モデル」を推定する解説記事を以前書きました。
しかし、Stanを使うのはちょっと敷居が高いという方も多いと思います。

そこで、このサイトでもよく取り上げている、dlmパッケージを使って、”慣れ”の統計モデリングをしてみます。
dlmパッケージは本当によくできているので、短いコードですますことができますよ。

なお、この記事は、状態空間モデルをある程度知っている方を対象としています。
「状態空間モデルってなに？」という方は、リンクを張った以下の記事を参照してください。上から順番に、徐々に難易度が上がっていきます。ローカルレベルモデルまでの知識があればOKです。

●状態空間モデル関連のページ
なぜ状態空間モデルを使うのか
状態空間モデル：状態空間モデルのことはじめ
ｄｌｍの使い方 ：Ｒで正規線形状態空間モデルを当てはめる
ローカルレベルモデル：ｄｌｍパッケージを使ってローカルレベルモデルを当てはめる
季節とトレンド：ｄｌｍパッケージを使って季節成分とトレンドの入ったモデルを作る
dlmによる時変形数モデル：ｄｌｍによる「時間によって係数が変化する回帰モデル」の作成

Pythonによる状態空間モデル：R言語ではなくPythonを使いたい方はこちらをどうぞ

スポンサードリンク

 

目次

1．時変係数モデルのご利益
2．今回扱うモデルの説明
3．シミュレーションデータの作成
4．dlmによる時変係数モデルの推定
おまけ:状態方程式、観測方程式を読み解く
5．結果の確認

 

1．時変係数モデルのご利益

今回は、状態空間モデルを使って「傾きが時間によって変化する回帰モデル」を推定します。

傾きが時間によって変化するだろうと考えられることはよくあります。
例えば、広告効果の測定が挙げられます。広告をだすと、最初は商品が良く売れるようになるけれども、徐々に慣れてきて、広告の効果が薄れていくと考えられます。その時、以下のような回帰モデルを作ったとします。

売り上げ ～ 切片 ＋ 広告(広告有りが１、無しが０) × 広告の効果

このとき、広告を出している間に、徐々に「広告の効果がなくなっていく」ということを、普通の回帰モデルでは表現することができません。

あるいは、薬の効果が徐々に薄れていくという場合など「傾きが時間によって変化する」と考えられることはよくあります。
そこで、状態空間モデルの出番です。

dlmパッケージを使って、お手軽に「傾きが時間によって変化する回帰モデル」を推定しましょう。
Rさえあれば、すぐに推定ができます。

 

2．今回扱うモデルの説明

今回は、以下のようなモデルを組みます。
英語で書かれた部分は、後で書くRのコードと対応しています。

response　=　intercept　＋　explanatory × slope + 観測誤差

ただし、観測誤差は、平均0で、標準偏差が「observationErrorSD」の正規分布に従います。

intercept(切片)は、ランダムウォークします。切片が時間によって変わるとお考えください。
また、データは1日ごとに取られているとします。

今日のintercept = 昨日のintercept ＋ 過程誤差その1

ただし、過程誤差その1は、平均0で、標準偏差が「processErrorSD」の正規分布に従います。

また、slope(傾き)もランダムウォークしていると仮定します。

今日のslope = 昨日のslope ＋ 過程誤差その2

ただし、過程誤差その2は、平均0で、標準偏差が「coefErrorSD」の正規分布に従います。

dlmパッケージの「dlm」とは「Dynamic Linear Model：動的線形モデル」の略です。
これは「誤差分布が正規分布、状態方程式と観測方程式が線形」である正規線形状態空間モデルを指します。

まずは、このモデルをそのまま使って、シミュレーションデータを作りましょう。

 

3．シミュレーションデータの作成

まずは、モデルのパラメタの設定をしておきます。
ここでは、観測誤差と過程誤差の標準偏差を指定します。

# パラメタの設定

# サンプルサイズ
N <- 500

# 観測誤差の標準偏差
observationErrorSD <- 20

# 過程誤差の標準偏差
processErrorSD <- 10

# 係数の過程誤差の標準偏差
coefErrorSD <- 0.5

 

次は、説明変数をシミュレーションで作ってみます。
これは適当に、標準偏差10の正規分布に従う乱数としておきます。

# 各種シミュレーションデータの生成

# 説明変数をシミュレーションで作る
set.seed(1)
explanatory <- rnorm(n=N, mean=10, sd=10)

 

次はslope(傾き)です。こいつは毎日過程誤差が加わって、変化していくようにしています。

# 各種シミュレーションデータの作成

# slopeのシミュレーション
set.seed(12)
slope <- cumsum(rnorm(n=N, mean=0, sd=coefErrorSD)) + 10

plot(slope, main="時間によるslopeの変化", xlab="day")

slopeは、「今日のslope = 昨日のslope ＋ 過程誤差その2」というように毎日変化していきます。
過程誤差2は「rnorm(n=N, mean=0, sd=coefErrorSD)」で表されます。
過程誤差が毎日つみあがっていくので、累積和をとるcumsum関数を使っています。

cumsum関数の動きは以下の通りです。
> cumsum(c(1,2,3,4,5))
[1] 1 3 6 10 15
数値を前から順番に足し合わせていることがわかります。これを利用して過程誤差を積み上げて行きます。

結果を図示しておきます。

 

同じようにして、intercept(切片)も作ります。

# 各種シミュレーションデータの作成

# interceptのシミュレーション
set.seed(3)
intercept <- cumsum(rnorm(n=N, mean=0, sd=processErrorSD))

plot(intercept, main="時間によるinterceptの変化", xlab="day")

 

最後は、responseです。

# 各種シミュレーションデータの作成

# responseのシミュレーション
set.seed(4)
response <- intercept + explanatory*slope + rnorm(n=N, mean=0, sd=observationErrorSD)

plot(response, main="responseのシミュレーション結果", xlab="day")

忘れないように、モデル式をもう一度書いておきます。
このモデル式とRのコードが対応していることを確認してください。
response　=　intercept ＋ explanatory × slope + 観測誤差
今日のintercept = 昨日のintercept ＋ 過程誤差その1
今日のslope = 昨日のslope ＋ 過程誤差その2
ちなみに、観測誤差の標準偏差は20、過程誤差その1は10、過程誤差その2は0.5となっています。

 

4．dlmによる時変係数モデルの推定

次は、dlmパッケージを使ってモデルを推定ましょう。
dlmパッケージの詳しい使い方はこちらをご覧ください。

まずは解析の準備からです。
パッケージを読み込みます。パッケージのダウンロードがまだの方は、Rを管理者として実行したうえで「install.packages(“dlm”)」を実行してください。そうやってパッケージをインストールできたうえで、以下のコードを実行します。

# 解析の準備


# ==================================================
# dlmパッケージによる、状態空間モデルの推定
# ==================================================

# パッケージの読み込み
library(dlm)

 

次は、モデルの型を決めます。
状態空間モデルと一口にいっても、かなりのバリエーションがあります。
今回は、説明変数を使って回帰モデルを組むので「dlmModReg()」という関数を使います。

# モデルの推定

#　Step1 モデルの型を決める
buildDlmReg <- function(theta){
  dlmModReg(
    X=explanatory , 
    dV=exp(theta[1]), 
    dW=c(exp(theta[2]), exp(theta[3]))
  )
}

dlmModReg()の引数「X」は説明変数を指定します。
dVは観測誤差の分散を指定します。分散の値は負にならないので「exp」をとっています。
dWは過程誤差の分散を指定します。過程誤差は2つあるので、パラメタも2つあります。

Step1で作成された関数を使って、尤度を最大にする3つのパラメタを推定しましょう。

# モデルの推定

#　Step2　パラメタ推定
fitDlmReg <- dlmMLE(
  response, 
  parm=c(2, 1, 1), 
  buildDlmReg,
  method = "SANN"
)

dlmMLE関数を使うことで、最尤法によるパラメタの推定ができます。
method = “SANN”を指定するとちょっと時間がかかるんですが、そのままだと警告文が出たので指定しました。

 
結果を確認します。

# 結果の確認

> # 推定されたパラメタの確認
> fitDlmReg
$par
[1]  6.064804  4.663469 -1.923347

$value
[1] 1945.618

$counts
function gradient 
   10000       NA 

$convergence
[1] 0

$message
NULL

> sqrt(exp(fitDlmReg$par))
[1] 20.7470061 10.2957837  0.3822527

推定されたパラメタにexpをとると、分散の値になります。
分散の平方根をとると、シミュレーションの時に設定した標準偏差になりますね。
『[1] 20.7470061 10.2957837 0.3822527』は左から順番に、観測誤差、過程誤差1、過程誤差2の標準偏差を表します。
正解は「20、10、0.5」だったので、まずまずというところです。

推定されたパラメタを使って、モデルを組みなおします。

# モデルの推定

# 推定された分散を使って、モデルを組みなおす
modDlmReg <- buildDlmReg(fitDlmReg$par)

 

あとは、フィルタリングしたり、スムージングしたりして、完了です。

# モデルの推定

# Step3
# カルマンフィルター
filterDlmReg <- dlmFilter(response, modDlmReg)

# Step4
# スムージング
smoothDlmReg <- dlmSmooth(filterDlmReg)

 

おまけ:状態方程式、観測方程式を読み解く

※数式が出てきます。苦手な方は飛ばしてもらってOKです。

動的線形モデルの解説文に、必ずと言ってよいほど出てくる数式があります。

$$\begin{array}{l}
Y_t = F_t θ_t + \mathit{v}_t, \hspace{ 27pt }\mathit{v}_t \sim N_m(0,V_t) \hspace{ 20pt }(1)\\
θ_t = G_t θ_{t-1} + \mathit{w}_t, \hspace{ 15pt }\mathit{w}_t \sim N_m(0,W_t) \hspace{ 13pt }(2)
\end{array}$$

式(1)を観測方程式、式(2)を状態方程式と呼びます。
添え字の「t」は時間を表しています。tが今日なら、t-1は昨日です。
Yが観測値で、Gは遷移行列と呼ばれます。またθは状態のパラメタ、ｖやｗは各々、観測誤差と過程誤差を表しています。
観測方程式は「今日の状態から今日の観測値を求める方程式」です。
状態方程式は「昨日の状態から今日の状態を求める方程式」です。

この方程式により、以下のモデル式を表します。
response　=　intercept ＋ explanatory × slope + 観測誤差
今日のintercept = 昨日のintercept ＋ 過程誤差その1
今日のslope = 昨日のslope ＋ 過程誤差その2

ここでは、GとかFとか出てきている方程式と、私が日本語で書いた方程式との対応を説明します。
特に回帰モデル(dlmModReg)の場合は、dlmの出力がちょっとややこしいので気を付けてください。

 

推定されたモデル式は、以下のようになっています。

# dlmの出力

> modDlmReg
$FF
     [,1] [,2]
[1,]    1    1

$V
         [,1]
[1,] 430.4383

$GG
     [,1] [,2]
[1,]    1    0
[2,]    0    1

$W
         [,1]      [,2]
[1,] 106.0032 0.0000000
[2,]   0.0000 0.1461172

$JFF
     [,1] [,2]
[1,]    0    1

$X
     [,1] 
[1,] 3.735
[2,] 11.84
[3,] ...  

$m0
[1] 0 0

$C0
      [,1]  [,2]
[1,] 1e+07 0e+00
[2,] 0e+00 1e+07

まず目につくのが「$FF」です。
$FFは、このようになっています。
$FF
　　　[,1]　[,2]
[1,]　　1　　1
これが、観測方程式のFtと一致しているように見えますし、一致することもあります。
しかし、dlmModRegの場合は一致しません。

真ん中のほうにある$JFFに注目してください。
$JFF
　　　[,1]　[,2]
[1,]　　0　　1
$JFFは「説明変数がどこに入るか」を表した行列です。
$JFFの1行1列目は「0」になっています。よって、ここには説明変数は入りません。
$JFFの1行2列目が「1」になっています。すると、$FFの1行2列目には、1番目の説明変数が入ります。

説明変数は$Xであらわされています。
$X
　　　[,1]
[1,] 3.735
[2,] 11.84
[3,] …
この$Xの1列目が、$FFに入ることになります。
時間がtであれば、$Xのt行目の値が入ります。

よって、観測方程式は以下のように書き直されます。
$$Y_t = \begin{bmatrix} 1 & X_t \end{bmatrix} θ_t + \mathit{v}_t, \hspace{ 27pt }\mathit{v}_t \sim N(0,\mathit{v}_t)$$

分散は、observationErrorSDの2乗ですね。dlmの出力の$Vに対応します。
$$Y_t = \begin{bmatrix} 1 & X_t \end{bmatrix} θ_t + \mathit{v}_t, \hspace{ 27pt }\mathit{v}_t \sim N(0,observationErrorSD^2)$$

 

次はθと状態方程式に移ります。
まず、θがinterceptとslopeの2つに分かれます。各々をμとβとおいておきます。
$$Y_t = \begin{bmatrix} 1 & X_t \end{bmatrix}
\begin{bmatrix} \mu_t \\ \beta_t \end{bmatrix} +
\mathit{v}_t,
\hspace{ 27pt }\mathit{v}_t \sim N(0,observationErrorSD^2)$$

状態方程式の遷移行列Gは、dlmの$GGと、今回はまったく同じになります。
分散は$Wの行列と対応しています。
$$\begin{bmatrix} \mu_t \\ \beta_t \end{bmatrix} =
\begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix}
\begin{bmatrix} \mu_{t-1} \\ \beta_{t-1} \end{bmatrix} +
\begin{bmatrix} \mathit{w}_{t,1} & 0 \\ 0 & \mathit{w}_{t,2} \end{bmatrix},
\begin{array}{l}
\hspace{ 27pt }\mathit{w}_{t,1} \sim N(0,processErrorSD^2) \\
\hspace{ 27pt }\mathit{w}_{t,2} \sim N(0,coefErrorSD^2)
\end{array}$$

 
 

まとめます。
最初の観測・状態方程式はこちら。
$$\begin{array}{l}
Y_t = F_t θ_t + \mathit{v}_t, \hspace{ 27pt }\mathit{v}_t \sim N_m(0,V_t) \hspace{ 20pt }(1)\\
θ_t = G_t θ_{t-1} + \mathit{w}_t, \hspace{ 15pt }\mathit{w}_t \sim N_m(0,W_t) \hspace{ 13pt }(2)
\end{array}$$

回帰モデルの場合はこうなります。
$$\begin{array}{l}
Y_t = \begin{bmatrix} 1 & X_t \end{bmatrix}
\begin{bmatrix} \mu_t \\ \beta_t \end{bmatrix} +
\mathit{v}_t,
\hspace{ 100pt }\mathit{v}_t \sim N(0,observationErrorSD^2) \\[5pt]
\begin{bmatrix} \mu_t \\ \beta_t \end{bmatrix} =
\begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix}
\begin{bmatrix} \mu_{t-1} \\ \beta_{t-1} \end{bmatrix} +
\begin{bmatrix} \mathit{w}_{t,1} & 0 \\ 0 & \mathit{w}_{t,2} \end{bmatrix},
\begin{array}{l}
\hspace{ 27pt }\mathit{w}_{t,1} \sim N(0,processErrorSD^2) \\
\hspace{ 27pt }\mathit{w}_{t,2} \sim N(0,coefErrorSD^2)
\end{array}
\end{array}$$

 
最後に、行列であらわされた式をばらしてみましょう。
$$\begin{array}{l}
Y_t = \mu_t + X_t \times{\beta_t} + \mathit{v}_t,
\hspace{ 20pt } \mathit{v}_t \sim N(0,observationErrorSD^2) \\
\mu_t = \mu_{t-1} + \mathit{w}_{t,1},
\hspace{ 47pt }\mathit{w}_{t,1} \sim N(0,processErrorSD^2) \\
\beta_t = \beta_{t-1} + \mathit{w}_{t,2},
\hspace{ 47pt }\mathit{w}_{t,2} \sim N(0,coefErrorSD^2)
\end{array}$$

あとは、Yにresponse、Xにexplanatory、μにintercept、βにslopeを当てはめれば、日本語の式と一致します。
動的線形モデルは、モデル式を行列であらわす上に、dlmパッケージの出力がちょっと読みづらいです。
しかし、この行列が読めると、モデルの解釈が大分としやすくなると思います。

 

5．結果の確認

まずはresponseの状態の推定結果からみていきます。
今回はスムージングされた結果を評価します。

スムージングされた結果は「smoothDlmReg$s」に格納されています。
こちらは2列あるんですが、1列目がinterceptの状態の推定結果、2列目がslopeの状態の推定結果となります。
というわけで、以下のように計算してやれば、responseが求まります。

# ==================================================
# 推定結果の図示と確認
# ==================================================

#　推定されたresponseの状態
estimatedLevel <- 
  dropFirst(smoothDlmReg$s)[,1] + 
  explanatory*dropFirst(smoothDlmReg$s)[,2]

推定結果の頭には余計な値が入っているので「dropFirst」しています。詳しくはこちらの記事を参照してください。
ここでは、intercept + explanatory*slopeが計算されていることを確認してください。

図示します。

# 図示

#　元データ
plot(response, col=1, main="response", pch="。", xlab="day")

# 推定された状態
lines(estimatedLevel, col=4)

# 観測誤差の入っていない、正しい値
lines(intercept + explanatory*slope, col=2)

# 凡例
legend("topleft", legend=c("スムージングの結果","正しい値"), lty=1, col=c(4,2))

大体一致しているようです。

次は、slopeの時間変化を見ていきます。
slopeは、先ほど説明したように「smoothDlmReg$s」の2列目の値です。
また、今回はフィルタリングされた結果「filterDlmReg$m」も併せて表示させます。

# 図示

#　推定された傾き
plot(dropFirst(filterDlmReg$m)[,2], type="l", ylim=c(0, 15), col=2, xlab="day", ylab="slope")
lines(dropFirst(smoothDlmReg$s)[,2], type="l", ylim=c(0, 15), col=4)

# 正しい傾き
lines(slope, col=1)

# 凡例
legend("topleft", legend=c("フィルタリングの結果","スムージングの結果","正しい値"), lty=1, col=c(2,4,1))

回帰係数の時間的変化をうまくとらえることができました。

 

参考文献

Dynamic Linear Models with R (Use R!)   私がdlm関連の記事を書くときは、いつもこの本を参考にしています。 dlmパッケージを使いこなしたい方はぜひどうぞ。

Rによるベイジアン動的線形モデル (統計ライブラリー)   Dynamic Linear Models with Rの日本語訳です。日本語版は読んだことないですが、日本語が良い方はこちらをどうぞ。

状態空間時系列分析入門   このサイトでいつも紹介している、状態空間モデルの入門書です。 状態空間モデルについて知りたければ、まずはこの本から始めるのが良いと思います。最も読みやすい状態空間モデルの入門書です。

カルマンフィルタの基礎   今回の記事では端折ってしまいましたが、カルマンフィルタについて詳しく知りたい方はこちらをどうぞ。多々あるカルマンフィルタの本の中では最も読みやすいです。

書籍以外の参考文献

・Giovanni Petris (2010), An R Package for Dynamic Linear Models. Journal of Statistical Software, 36(12), 1-16.
英語ですが、無料で読めるdlmの解説文です。ただし、dlmModRegの説明は少ないです。。。

・Commandeur & Koopman「状態空間時系列分析入門」をRで再現する
参考文献としても上げた状態空間時系列分析入門をRで再現した記事です。

 

うちのサイトのほかの記事を参照したい方はこちらをどうぞ。
状態空間モデルのカテゴリだけをまとめたリンクです。

スポンサードリンク