R: KFAS 1.0.4-1による一部の状態のとりだし

KFASによる状態空間モデルの解析。生態学会のときはpredict()のstates引数で状態を指定して、過去のデータについて予測をおこなうというふうにしていたのですが、KFAS 1.0.4-1ではpredict()のstatesの指定が有効でないみたい。マニュアルをよく読むと、こういうときはsignal()をつかうほうがよいのかも。

# 生態学会のときのスライドも、いまからみるといろいろ手直ししたいところがありますなあ。

KFS()で平滑化して、signal()のstatesで、状態を指定する。

###
### KFASによる、ポアソン分布の状態空間モデルの解析例
###
### 例題は、help(KFAS)より。
###
### Jouni Helske (2014). KFAS: Kalman Filter and Smoother for
### Exponential Family State Space Models.
### R package version 1.0.4-1.
### http://CRAN.R-project.org/package=KFAS
###

## KFASとggplot2を使用
library(KFAS)
library(ggplot2)

## データにはSeatbeltsを使用
data(Seatbelts)

## モデル定義
## VanKilled: バン運転手の死者数
## law: シートベルト義務化の前後をしめす変数
model.van <- SSModel(VanKilled ~ law +
                     SSMtrend(degree = 1, Q = list(matrix(NA))) +
                     SSMseasonal(period = 12, sea.type = "dummy", Q = NA),
                     data = Seatbelts, distribution = "poisson")

## 分散Qを最尤推定
fit.van <- fitSSM(model = model.van, inits = c(-4, -7, 2))

## カルマンスムーザー適用と、トレンドおよび回帰の状態のとりだし
out.van <- KFS(fit.van$model, smoothing=c("state"))
sig.van <- signal(out.van, states = c("trend", "regression"))

p <- ggplot(data.frame(Year = seq(start(Seatbelts)[1],
                                  end(Seatbelts)[1] + 11/12, 1 / 12),
                       VanKilled = c(Seatbelts[, "VanKilled"]),
                       Smooth = exp(c(sig.van$signal))))
p + geom_point(aes(x = Year, y = VanKilled)) +
    geom_line(aes(x = Year, y = VanKilled)) +
    geom_line(aes(x = Year, y = Smooth), colour = "gray35", size = 1.2) +
    ylab("No. of van drivers") +
    theme_bw(12, "Helvetica")

結果

タグ：R