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

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

SPONSOR LINK

状態プロセス

個体群動態を記述する式はたくさんありますが、ここでは指数増加を仮定した個体群動態モデルを用いることにします。(勝手にシジュウカラを想定し)ある年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です。

時系列データの取り扱い

時系列データの特性:「独立」仮定の破綻

Rパッケージはとても簡単に利用できるよう整備されているため、ANOVAや線形モデルを「なんとなく」利用することが多いと思います。しかし、これらの解析におけるとても重要な仮定のひとつに、「繰り返しの単位となるサンプルは“互いに独立”であること」があります(あるサンプル値iが変化しても他のサンプル値は影響を受けない)。調査研究の際、「調査地点間は~m以上離した」のように設定することが多いですが、これは地点間の距離が近すぎることで、サンプル値同士が似た値をとってしまうことを避ける意図があります(地点iのサンプル値が直近地点jのサンプル値に影響する)。“独立性”が保てないまま解析すると(あるいはそれをモデルのなかで考慮しない限り)、パラメーター推定の正確性が損なわれ、誤った結論を導きかねません。

SPONSOR LINK

それでは、生態学における時系列データはどのような特徴をもつのでしょうか。例えば、ここに個体の移出入はない閉鎖個体群があるとした場合、時間tにおける個体数Ntは、前年の個体数Nt-1に依存する繁殖や死亡の過程を経て決まると考えられます。これは、「サンプルは互いに独立」という仮定が崩れる典型といえるでしょう。この問題に対処するための様々な解析アプローチがありますが、ここでは生態学の文脈で力を発揮する「状態空間モデル(State-space model)」の基本的な考え方と実装手順を説明します。

状態空間モデル

時系列データの解析には、様々なアプローチがあります。例えば、自己回帰モデル(Auto Regressive model; AR model)や移動平均モデル(Moving Average model; MA model)、この両者を組み合わせたARMAモデルなどがありますが、いずれの手法も「誤差が年(i.e. 観察ユニット)を経るごとに蓄積する」過程を明示的に考慮する点において共通しています。しかし、これらの手法では、生態学の興味の中心である「増加率」や「増加率の年変動」を推定することは難しいという欠点を抱えています。状態空間モデルは、この意味において上記の手法より優れており、生態学の調査研究に適した手法だと思われます。なぜなら、1)生態過程(個体群の増加プロセス)を明示的に扱えるだけでなく、2)観察過程も考慮できる(調査努力量、発見効率など)という長所をもつからです。

状態空間モデルでは、このモデルの少し変わった(賢い)考え方により、こうした複雑なプロセスのモデル化を可能にしています。その考え方とは、生物の個体数Nが観察されるまでの過程を、「観察プロセス(Observation process)」と「状態プロセス(State process)」に切り分けて考える、というものです。

状態プロセスでは、生物が実際の繁殖や死亡の過程を経て増減する「真の個体数Nstate(推定対象となるパラメーター)」をモデル化します(直接観察できない“知りたい”もの)。一方、観察プロセスでは、観察個体数Nobs(データ)は、真の個体数Nstateから誤差を伴って観察されると仮定し、その過程をモデル化します。このような枠組みを取り入れることで、観察誤差(σobs)と真の揺らぎ(環境変化による増加率の年変動など;σstate)を明確に区別でき、より正確なパラメーター推定(個体群増加率など)が可能になります。