R: KFASによる状態空間モデル (6)

二項分布のばあい。

ライブラリのよびだしとか、逆ロジット関数の定義など。

library(ggplot2)
library(KFAS)

inv.logit <- function(x) 1 / (1 + exp(-x))

模擬データの生成。一定数の調査地があって、そのうちのいくつで対象とする生物が発見されるか、というイメージです。発見率はロジットスケールでランダムウォークして、さらに変数xの値によって発見率が変化するとします。

set.seed(5)

N.t <- 40                # 測定回数
N.u <- 20                # 調査地の数
x <- rbinom(N.t, 1, 0.5) # 変数
alpha <- rep(NA, N.t)    # xの効果をふくまない出現率のロジット
alpha[1] <- -2
for (t in 2:N.t) {
  alpha[t] <- rnorm(1, alpha[t - 1], 0.3)
}
N <- rbinom(N.t, N.u, inv.logit(alpha + 1.0 * x))

グラフに表示してみます。

df <- data.frame(Time = 1:N.t, x = x, N = N)
p11 <- ggplot(df, aes(x = Time, y = N)) +
      geom_line() +
      theme_bw()
print(p11)

KFASで解析してみます。二項分布は、SSModel関数でdistribution引数に "binomial"を指定します。全体の数はu引数で指定します。

model6 <- SSModel(N ~ x + SSMtrend(degree = 1, Q = list(NA)),
                  data = df,
                  u = N.u,
                  distribution = "binomial")
fit6 <- fitSSM(model = model6, inits = c(1, 1), method = "BFGS")
out6 <- KFS(fit6$model, smoothing = c("state", "mean"))

発見数の期待値と、変数xの影響がないばあい(x=0のばあい)の発見数の期待値をかさねてプロットしました。

df$Smooth <- fitted(out6) * N.u
df$Level <- inv.logit(c(signal(out6, states = "trend")$signal)) * N.u
p12 <- ggplot(df) +
       geom_line(aes(x = Time, y = N), color = "black") +
       geom_line(aes(x = Time, y = Smooth), color = "red") +
       geom_line(aes(x = Time, y = Level), color = "blue", alpha = 0.5) +
       theme_bw()
print(p12)

fitted関数を適用してえられる平滑化された値は0〜1の発見率になっていますので、全体の数N.uをかけて発見数の期待値にしました。また、変数xの影響がないばあいの発見数の期待値は、ローカルレベルモデル部分のみの信号成分をとりだし、逆ロジット関数を適用して発見率にしたうえで、全体の数をかけて求めています。

赤が期待値、青がxの影響がない場合の期待値です。

タグ：R KFAS