Summary

状態空間モデル（State Space Modeling）とは、モデル構造として状態過程と観察過程をもつ統計モデルの総称。JAGSとRスクリプトの備忘録。

State Space Modeling

生態学の時系列データ解析では、状態空間モデルがとても役立つ。なぜなら、「状態過程（State process）」と「観察過程（Observation process）」という二つプロセスを明示的に切り分けたモデル構造を持つからである。状態過程は、状態の時間変化をモデル化する（状態とは、集団サイズの増減のような実際には直接観察できないものを指す）。一方、観察過程では、状態に何らかの誤差が加わり、その結果としてデータが得られるという過程をモデル化する。例えば、状態過程では集団動態を表すモデルを当てはめ、そこに何らかの観察誤差が加わった結果として手元のデータが得られているという仮定を置く。これにより、集団増加率やその年変動など、将来予測に大きな影響を及ぼすパラメーターをより適切に推定できる（バイアスが少ない）。

JAGS script

次のような15年分の時系列データ（魚の個体数）があるとする。これに対し、状態空間モデルを当てはめてみる。

plot(dat$X, type = "o", ylab = "Fish abundance", xlab = "Year")

State process：状態過程では、集団動態を指数増加関数としてモデル化することにする（対象に応じて変える）。各年の集団増加率\(r_t\)は、正規分布に従って年変動すると仮定し、その平均は\(\mu_r\)、分散は\(\sigma^2_r\)とする（対数スケールであることに注意）。
\[
ln~n_{t+1} = r_{t} + ln~n_{t}\\
r_t \sim Normal(\mu_{r},\sigma^2_r)
\]

Observation process：この集団動態モデルから予測される値に対し、観察誤差が加わった結果として観察値\(X_t\)が得られるとする。ここでは、観察誤差は過分散を考慮したポアソン分布（Overdispersed-Poisson distribution）に従うと仮定する。
\[
X_t \sim Poisson(n_{obs,t})\\
ln~n_{obs,t} \sim Normal(ln~n_t, \sigma^2_{obs})\\
\]

JAGSのスクリプトで書き下す。

model{
  # Prior
  mu.r ~ dnorm(0, 0.001)

  tau.r <- pow(sigma.r, -2)
  sigma.r ~ dunif(0, 100)

  tau.obs <- pow(sigma.obs, -2)
  sigma.obs ~ dunif(0, 100)

  ln.n[1] ~ dnorm(0, 0.001)

  # Likelihood
  ## State process
  for(t in 1:(Nt-1) ){
   ln.n[t+1] <- r[t] + ln.n[t]
   r[t] ~ dnorm(mu.r, tau.r)
  }

  ## Observation process
  for(t in 1:Nt){
    ln.n.obs[t] ~ dnorm(ln.n[t], tau.obs)
    n.obs[t] <- exp(ln.n.obs[t])
    X[t] ~ dpois(n.obs[t])

    ### backtransformation to a real-scale
    n[t] <- exp(ln.n[t])
  }

}

R script

上記のJAGSスクリプトをRから走らせる。

# R script
# Read data
X <- dat$X # time series data
Nt <- length(dat$X)

# Dataframe for JAGS
Djags <- list( X = X, Nt = Nt )

# Claim parameters to be monitored
para <- c("mu.r", "sigma.r", "sigma.obs", "n")

# MCMC setting
n.ad <- 100 # number of adaptation
n.iter <- 1E+4 # number of total iteration
n.thin <- max(3, ceiling(n.iter/500)) # number of thinning
burn <- ceiling(max(10, n.iter/2)) # number of burn-in
Sample <- ceiling(n.iter/n.thin) # number of sample per chain

# Intial values
inits <- replicate(3, list(.RNG.name = "base::Wichmann-Hill", .RNG.seed = NA), simplify = FALSE)
for(i in 1:3){ inits[[i]]$.RNG.seed <- i }

# Run JAGS
library(runjags)
m <- read.jagsfile("BayesModel/model_SSM.R") # read model file
post <- run.jags(model = m$model, monitor = para, data = Djags,
                 n.chains = 3, inits =inits, method = "parallel",
                 burnin = burn, sample = Sample, adapt = n.ad, thin = n.thin,
                 n.sims = 3)

## Calling 3 simulations using the parallel method...
## Following the progress of chain 1 (the program will wait for all
## chains to finish before continuing):
## Welcome to JAGS 4.3.0 on Sun Jun 09 18:19:26 2019
## JAGS is free software and comes with ABSOLUTELY NO WARRANTY
## Loading module: basemod: ok
## Loading module: bugs: ok
## . . Reading data file data.txt
## . Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 15
##    Unobserved stochastic nodes: 33
##    Total graph size: 100
## . Reading parameter file inits1.txt
## . Initializing model
## . Adapting 100
## -------------------------------------------------| 100
## ++++++++++++++++++++++++++++++++++++++++++++++++++ 100%
## Adaptation successful
## . Updating 5000
## -------------------------------------------------| 5000
## ************************************************** 100%
## . . . . . Updating 10000
## -------------------------------------------------| 10000
## ************************************************** 100%
## . . . . Updating 0
## . Deleting model
## . 
## All chains have finished
## Simulation complete.  Reading coda files...
## Coda files loaded successfully
## Calculating summary statistics...
## Calculating the Gelman-Rubin statistic for 18 variables....
## Finished running the simulation

re <- summary(post)
round(re[1:3,1:3], 3)

##           Lower95 Median Upper95
## mu.r       -0.007  0.079   0.184
## sigma.r     0.000  0.100   0.307
## sigma.obs   0.001  0.095   0.274

Comparison with “truth”

先に示したデータは、パラメータが既知のシミュレーションデータである。シミュレーションに用いたパラメータの値は以下の通りである。推定結果と比べると、95%CIに真の値が含まれていることがわかる。
\[
\mu_r = ln~1.1 = 0.095…\\
\sigma_r = 0.1\\
\sigma_{obs} = 0.1\\
\]

実際の集団動態と推定結果を比べてみる。ポイントが推定に用いたデータ、点線が「真」の個体数（観察誤差が加わる前の個体数\(n_t\)）、実線がモデルの推定値である。実線が真の値をよく反映していることが見て取れる。

ID <- which(rownames(re) == "n[1]"):which(rownames(re) == "n[15]")
n.est <- re[ID, "Median"]
plot(dat$X, ylab = "Fish abundance", xlab = "Year")
lines(dat$x, lty = 2)
lines(n.est, lty = 1)

一方、こちらは観察値をそのまま各種パラメータの推定に用いた場合である。

mu.r <- mean(diff(log(dat$X) ) )
round(mu.r, 3)

## [1] 0.097

sigma.r <- sd(diff(log(dat$X) ) )
round(sigma.r, 3)

## [1] 0.405

（このシンプルな例では）集団増加率の平均の推定に大きな問題がなさそうだが、年変動はかなり大きめに推定されている。これは、集団増加率の年変動と観察誤差が一緒になっているためだろう。