JAGSによる状態空間モデルの実装

 JAGSを使って状態空間モデルを実装します。まず、真の個体数の時間変化を記述する状態プロセスから考えます。

状態プロセス

個体群動態を記述する式はたくさんありますが、ここでは指数増加を仮定した個体群動態モデルを用いることにします。(勝手にシジュウカラを想定し)ある年tにおけるシジュウカラの真の個体数をNt、個体群増加率をrtとすると、個体数の時間変化は以下の式で表されます。

Nt+1 = exp(rt) * Nt   (eq. 4-1)
※ exp(a)はeaという意味
※ exp(rt) = Nt+1/Nt

しかし、eq.4-1は掛け算の式になっているため、このまま計算させると計算量が膨大になってしまいます。そこで、両辺対数をとることで、eq.4-1を足し算の式に変形します(線形化)。

log(Nt+1) = rt + log(Nt)   (eq. 4-2)

Xt = log(Nt)とすると、

Xt+1 = rt + Xt   (eq.4-3)

時間tにおける個体群増加率rtは平均rmean、分散σstate2の正規分布に従うと仮定し、

rt ~ Normal (rmean, σstate2)   (eq.4-4)

となります。rtを正規分布に従う確率変数としていますが、これは個体群増加率の年変動をモデル化するためです(年変動はσstate2に制御される)。これで、指数増加モデルに基づく状態プロセスの記述は終わりです。BUGSファイルの記述には、eq.4-3、4-4を使います。

#State process
for(t in 1:(N.year-1)){
  X[t+1] <- r[t] + X[t]
  r[t] ~ dnorm(mean.r, tau.r)
}

※時間tの繰り返し構文。終わりがN.year-1となっている点に注意

観察プロセス

 次に、観察過程をモデル化し、実際の観察個体数Yt(データ)と真の個体数Xt(状態プロセス)をつなぎ合わせます。時間tにおける真の個体数Xtに誤差が加わり、観察値平均Xobs,tが得られると仮定します。

Xobs,t ~ Normal(Xt, σobs2)   (eq. 4-5)
Nobs,t = exp(Xobs,t)   (eq. 4-6)

eq.4-4のσobsは、観察誤差の程度を決めるばらつきのパラメーター(標準偏差)です。eq. 4-6でXobs,t(対数スケール)を実数スケールの個体数Nobs,tに変換し(この時点でNobs,tは正の実数;カウントデータのような正の整数ではない)、Nobs,tを平均値としたポアソン分布から観察個体数Ytが観察されるとします。

Yt ~ Pois(Nobs,t)   (eq. 4-6)

これで初めて、実際の観察データと真の個体数の関係をモデル化することができました。状態プロセスと合わせて、BUGSのモデルファイルを書き下します。

model{
# State-process
for(t in 1:(N.year-1)){
    X[t+1] <- r[t] + X[t]
    r[t] ~ dnorm(mean.r, tau.r)
  }
  
  # Observation-process
  for(t in 1:N.year){
    X.obs[t] ~ dnorm(X[t],tau.obs)
    N.obs[t] <- exp(X.obs[t])
    Y[t] ~ dpois(N.obs[t])
    
    Nest[t] <- exp(X[t]) # 真の個体数を対数スケールから実数スケールに戻す
  }

#Prior
X[1] ~ dnorm(0,0.001) #真の個体数の1番目の値は全くわからないので、事前分布を設定
tau.r ~ dgamma(0.001,0.001); sigma2.r <- 1/tau.r
tau.obs ~ dgamma(0.001,0.001); sigma2.obs <- 1/tau.obs
mean.r ~ dnorm(0,0.001)
}

このBUGSモデルは“modelSSM.R”として保存することにします。

状態空間モデルの推定

 状態空間モデルと同様のプロセスからデータYtを生成し、生成に用いたパラメーターを推定します。

###シミュレーションデータの生成###
N.year <- 20 #シミュレーションは20年分

# parameters to be estimated
mean.r <- 0.1 #20年を通じた平均個体群増加率
sigma.r <- 0.2 # 個体群増加率の年変動(標準偏差、対数スケール)σstate
sigma.obs <- 0.3 # 観察誤差の程度(標準偏差、対数スケール)σobs

# Simulated data
Y <- X <- Xobs <- N <- Nobs <- r <- NULL
X[1] <- 2 #1年目の真の個体数(対数スケール)

for(t in 1:(N.year-1)){
  #平均"mean.r"、標準偏差 "sigma.r"の正規分布から年tの個体群増加率をサンプリング
  r[t] <- rnorm(1,mean.r,sigma.r)
  X[t+1] <- r[t] + X[t]
}

for(t in 1:N.year){
  N[t] <- exp(X[t])
  
  #平均"X[t]"、標準偏差"sigma.obs"の正規分布からXobs,tをサンプリング
  Xobs[t] <- rnorm(1,X[t],sigma.obs)
  Nobs[t] <- exp(Xobs[t])
  
  #平均"Nobs[t]"のポアソン分布からYtをサンプリング
  Y[t] <- rpois(1,Nobs[t]) 
}

###modelSSM.Rを走らせるコード###
library(runjags)
Djags <- list(Y=Y, N.year=N.year)
m <- read.jagsfile("modelSSM.R")
para <- c("mean.r", "sigma2.r", "sigma2.obs", "Nest")
inits <- replicate(3,
                list(.RNG.name="base::Wichmann-Hill", .RNG.seed=round(runif(1,1,100))),
                simplify = FALSE)
post <- run.jags(m$model, monitor = para, data = Djags,
                 n.chains = 3, inits=inits, method = "parallel",
                 burnin = 500, sample = 1000, adapt = 100, thin = 3,
                 n.sims=3)
re <- summary(post)

>re
                Lower95     Median  Upper95     psrf #真の値
mean.r     -0.040826600 0.07856190 0.218723 1.068888 #(0.1)
sigma2.r    0.000213157 0.01932655 0.248413 1.124571 #(0.04)
sigma2.obs  0.000314576 0.13278200 0.316731 1.026555 #(0.09)

生成した観察値Yに基づくsigma2.r(状態変数の年変動)の推定値は、0.40となってしまいます(雑ですが、sd(Y)としたとき)。しかし、状態空間モデルでは、観察誤差を分離して評価することで、推定値は0.04付近に収まっています。こうした個体群増加率やその年変動のパラメーターは、個体群の行く末を予測するための重要なパラメーターなので、状態空間モデルのご利益は非常に大きなものとなりそうです。

推定された値を図示します。

黒の実線がNestの中央値、点線が95%CIです。

SPONSOR LINK

Spread the love
  • 4
    Shares

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です