発見率モデル(除去法)

除去法による発見率モデル

発見率を推定する方法として、繰り返し調査の他に、除去法による推定方法があります。除去法では、調査区間から個体を一時的に除去していき(例えば、バケツに入れておくなど)、個体数の減少パターンから発見率を推定します。個体は閉鎖区間からいずれいなくなるのですが(全数調査:この場合推定はいらない)、それを実行するのは現実的ではありません。そのため、2回~3回の調査を行い、その限られた情報から真の個体数を推定することになります。
 除去法で得られたデータは、次のような過程を経て観察されていると考えることができます(簡単のため、3回除去を行った場合を想定)。まず、3回調査したときにとれた総個体数nは、真の個体数Nに採捕効率π(0–1の値をとる:真の個体数に対する採捕個体数の割合)をかけたものなので、

SPONSOR LINK

n = πN

と表すことができます。次に、1回目の調査で採捕された個体数n1は、真の個体数Nに対し、「調査一回あたりの採捕効率p」を掛け合わせたものになります。

n1 = Np = nπ−1 p

2回目の調査で採捕される個体数n2は、真の個体数Nから1回目の調査で採捕された個体数を差し引いたもの(N− n1)に対し、採捕効率pを掛け合わせたものとして表します。

n2 = (N− n1)p = (nπ−1 − nπ−1 p)p = nπ−1p(1−p)

同様に考えて、3回目の調査で採捕される個体数n3は以下の式で表せます。

n3 = (N− n1− n2)p = nπ−1p(1−p) (1−p)

つまり、k回目の採捕個体数を一般化した形で表すと、

nk = nπ−1p(1−p)k−1

となります。

このように得られるデータは、多項分布(発見プロセス)とポアソン分布(真の個体数Nの空間的ばらつき)の混合モデルとして表現することができます。

[ni1, ni2, ni3|Ni, π] ~ Multinomial(Ni, π)

Ni ~ Pois(λ)

SPONSOR LINK

BUGSによる実装

実際にスクリプトを書く場合、BUGSの多項分布は限られた観察数n(総採捕数)に対する各要素の比率(nk/n)を推定することしかできないので、真の個体数Nと総採捕個体数nの間の関係を記述する際に、二項分布を仮定することになります。

model{
###Priors###------------------------------------------------
p0 ~ dbeta(1,1)
log(lambda) <- alpha
alpha ~ dnorm(0,ninfo)
  
###Likelihood###--------------------------------------------
  for(i in 1:Nsite){
    #p[i]: probablity of being captured with a single pass
    p[i] <- p0
   
    mu[i,1] <- p[i] # proportion of fish removed at 1st pass
    mu[i,2] <- p[i]*(1-p[i])
    mu[i,3] <- p[i]*(1-p[i])*(1-p[i])
   
    #total catchability
    pi[i] <- sum(mu[i,])
    
    for(j in 1:J){# J = 3
      muc[i,j]	<- mu[i,j]/pi[i]
    }
    y[i,1:J] ~ dmulti(muc[i,1:J], ncap[i])
    ncap[i] ~ dbin(pi[i], N[i]) #ncap: total number of fish captured
  }
  
  for(i in 1:Nsite){ N[i] ~ dpois(lambda) }
}

発見率について、地点間でランダムにばらつく場合を想定する場合、以下のように書き換えます。

model{
###Priors###------------------------------------------------
logit.p0 <- logit(p0)
p0 ~ dbeta(1,1)
log(lambda) <- alpha
alpha ~ dnorm(0,ninfo)
tau.eps <- pow(sigma.eps,-2)
sigma.eps ~ dunif(0,100)

###Likelihood###--------------------------------------------
  for(i in 1:Nsite){
    #p[i]: probablity of being captured with a single pass
    logit(p[i]) <- logit.p0 + eps[i]
    eps[i] ~ dnorm(0,tau.eps)
   
    mu[i,1] <- p[i] # proportion of fish removed at 1st pass
    mu[i,2] <- p[i]*(1-p[i])
    mu[i,3] <- p[i]*(1-p[i])*(1-p[i])
   
    #total catchability
    pi[i] <- sum(mu[i,])
    
    for(j in 1:J){# J = 3
      muc[i,j]	<- mu[i,j]/pi[i]
    }
    y[i,1:J] ~ dmulti(muc[i,1:J], ncap[i])
    ncap[i] ~ dbin(pi[i], N[i]) #ncap: total number of fish captured
  }
  
  for(i in 1:Nsite){ N[i] ~ dpois(lambda) }
}

そのほか、説明変数に関する情報があれば、それらの発見率に対する効果なども見ることが可能です。

発見率モデル(繰り返し法)

ロバストデザインによる調査

生物の分布や個体数を調べ、対象生物のある地点における存在確率を推定する際、「いない」ことを証明するのはとても難しい問題です。なぜなら、生物調査で「本当はいる個体」を100%発見することは不可能だからです。この問題をそのままにすると存在確率を過小評価し、パラメーター推定に偏りをもたらすことが知られています。

この問題に対応するため、最近ではロバストデザインによる調査が一般的になってきています。ロバストデザインとは、同じ調査地点で複数回調査する調査デザインのことで、これにより発見率に起因する「在・不在」や「個体数」データのばらつきを統計的に補正することができます。基本的な原理としては、繰り返し調査の間は「真の在・不在」の状態は変化しないという仮定を置くことで、「同じ地点における調査回ごとの在・不在データのばらつきは発見率に由来する」というモデリングを行います。

BUGSによるモデリング(応答変数が在・不在)

調査回tにおける調査地点iでの在・不在をYti、真の在不在をzi(調査回tの間ではzは変化しないと仮定するので、tは省く)、発見率をpとしたとき、在・不在データが得られる確率は下記のように整理されます。

真の状態は在、データも在:zi*p
真の状態は在、データは不在:zi*(1-p)
真の状態は不在、データも不在:zi*p or zi*(1-p)

これらをモデル化するためには、まずは真の状態(zi)を記述し、次に真の状態ziに依存して観察値Ytiが得られるプロセスをモデル化します。

model{
#Likelihood
for(i in 1:Nsite){
 Y[i] ~ dbin(mu[i],T.survey)  #Y[i]はT.survey回調査したうち、在となった回数
 mu[i] <- p*z[i] #データが在となる確率は真の状態と発見率の積
 z[i] ~ dbern(psi) #真の存在確率psi
}
# Prior
p ~ dunif(0,1)
psi ~ dunif(0,1)

}

BUGSモデルでは、Ytiは発見回数として表されます(例えば、3回調査して2回見つかった場合はY=2)。前半の二項モデルの部分(Y[i] ~ dbin(mu[i],T.survey))では、発見率に応じて3回中何回が「在」となるそうかをモデル化しています。

真の在・不在の状態に影響する要因を調べたい場合は、下記のようにスクリプトを書き替えることで容易にパラメーター推定を行えます。

model{
#Likelihood
for(i in 1:Nsite){
 Y[i] ~ dbin(mu[i],T.survey)
 mu[i] <- p*z[i]

 z[i] ~ dbern(psi[i])
 logit(psi[i]) <- b[1] + b[2]*x[i]
}

# Prior
p ~ dunif(0,1)
for(i in 1:2){ b[i] ~ dnorm(0, 0.001) }

}

なお、実際にモデルを回す際には、zの初期値の与え方に注意しなければなりません。zがありえない値(例えばY=3なのにz=0となるなど)になるとBUGSはエラーコードを返してきます。そのため、zの初期値はz = rep(1, N.sample)などとしておくほうが無難です。

BUGSによるモデリング(応答変数が個体数)

応答変数が個体数などの離散値の場合でも、少しスクリプトを書き替えるだけで対応することが可能です。調査回tにおける調査地点iでの発見個体数をYti、真の在不在をNiとすると:

model{
#Likelihood
for(i in 1:Nsite){
 for(t in 1:Nsurvey){
  Y[t,i] ~ dbin(p, N[i])
 }
 N[i] ~ dpois(lambda) #真の個体数の平均値lambda
}

# Prior
p ~ dunif(0,1)
log(lambda) <- log.lambda
log.lambda ~ dnorm(0, 0.001)
}

なお、ここでの発見率pは、発見された個体数を真の個体数で除したもの(見つかった個体の割合)として解釈することができます。上と同様に、lambdaを線形予測子と関連付けることで、真の個体数に影響する要因を調べることができます。

発見率のばらつき

上記のモデルでは、場所によらず発見率は一定としていましたが、様々な理由で発見率が地点によってばらつくことが考えられます。そのような場合には、発見率をなんらかの確率分布に従うランダム変数とすることで、そのばらつきを考慮することが可能です。上記の在・不在モデルにおいて、発見率がランダムにばらつくと仮定したモデルは次のように書くことができます。

model{
#Likelihood
for(i in 1:Nsite){
 Y[i] ~ dbin(mu[i],T.survey)
 mu[i] <- p[i]*z[i]
 z[i] ~ dbern(psi)

 logit(p[i]) <- logit.mean.p + eps[i]
 eps[i] ~ dnorm(0, tau.p)
}
# Prior
logit(mean.p) <- logit.mean.p
mean.p ~ dunif(0,1)
psi ~ dunif(0,1)
tau.p ~ dgamma(0.001,0.001)

}

上のモデルではロジットスケールで正規分布を仮定した誤差を考えていますが、ベータ分布によるモデリングなども可能です。

事前分布によく用いる確率分布

事前分布に用いられる確率分布をリストします。パラメーターが実際にどの範囲の値をとりうるのか、という視点から整理します。

とる値に制限のないパラメーター

該当するパラメーター:偏回帰係数など
正規分布 dnorm
一様分布 dunif

#Prior
b ~ dnorm(0,0.0001)
b ~ dunif(-1000,1000)

正の値をとるパラメーター

該当するパラメーター:分散
(逆)ガンマ分布 dgamma
一様分布 dunif
正規分布(truncationつき)dnorm

#likelihood
y[i] ~ dnorm(mu, tau) # tau is a precision parameter

#Prior
tau ~ dgamma(0.0001,0.0001) # 精度(分散の逆数)に事前分布にガンマ分布を指定することで、分散の事前分布は逆ガンマ分布となる
tau <- pow(sigma,-2); sigma ~ dunif(0,100) # SDに対して0-100の範囲をとる一様分布を事前分布に指定
tau <- pow(sigma,-2); sigma ~ dnorm(0,0.0001)T(0,100) # SDに対して0-100の範囲をとるフラットな正規分布を事前分布に指定

備考:逆ガンマ分布について、Gelmanさんらにより無情報ではないという指摘がされていますが、この問題は分散の推定値が小さいときに顕著になるようです(Gelman and Hill 2007)。一方、一様分布は無情報すぎて、複雑なモデルでは収束しにくい傾向にあります。経験上truncated Normalが一番使い勝手がいいように思います。なお、ランダム効果のグループ数が少ない場合、ばらつきの推定が非常に難しくなります。そのような場合はhalf Cauchy分布を用いるとよいらしいです(Gelman and Hill 2007)。

Gelman and Hill 2007 Data Analysis Using Regression and Multilevel/Hierachical Models

0~1の値をとるパラメーター

該当するパラメーター:生存率、発見率などの確率パラメーター
ベータ分布 dbeta
一様分布 dunif

#Prior
s ~ dbeta(1,1) #無情報事前分布に該当
s ~ dunif(0,1)

尤度関数によく用いる確率分布

JAGSで利用できる一般的な確率分布をリスト化します。尤度関数部分のみを書いていますので、実装する場合は適宜事前分布を設定してください。

正規分布

パラメーター:平均、分散
使用条件:正規性の仮定が求められる
備考:推定が早い

y[i] ~ dnorm(mu[i], tau)
mu[i] <- b[1] + b[2]*x[i]

ポアソン分布

パラメーター:平均(=分散)
使用条件:非負の離散値、平均と分散が同程度であることが必要
備考:正規乱数の誤差項を追加することでOverdispersed Poissonに拡張可能

# Poisson model
y[i] ~ dpois(lambda[i])
log(lambda[i]) <- b[1] + b[2]*x[i]

# Overdispersed Poisson model
y[i] ~ dpois(lambda[i])
log(lambda[i]) <- b[1] + b[2]*x[i] + eps[i]
eps[i] ~ dnorm(0,tau)

負の二項分布

パラメーター:サイズr、確率p(平均=r(1-p)p, 分散=r(1-p)/p2
使用条件:非負の離散値
備考:ポアソン分布では過分散する場合によく用いられる。確率pで生じる現象がr回生じるまでの試行回数(yに相当)を表現。サイズパラメーターは正の実数。幾何分布(r = 1)の一般化。

y[i] ~ dnegbin(p[i], r)
p[i] <- r/(r+mu[i])
log(mu[i]) <- b[1] + b[2]*x[i]

二項分布

パラメーター:試行回数n、確率p
使用条件:非負の離散値で上限値が既知
備考:過分散がある場合、正規乱数の誤差項を追加することでNormal-binomialに拡張可能。ベルヌーイ分布(n = 1)の一般化。

# Binomial model
y[i] ~ dbin(p[i], n)
logit(p[i]) <- b[1] + b[2]*x[i]

# Normal-binomial model
y[i] ~ dbin(p[i], n)
logit(p[i]) <- b[1] + b[2]*x[i] + eps[i]
eps[i] ~ dnorm(0,tau)

打ち切り関数

打ち切りのあるデータ

野外のデータ、あるいは実験データを扱う際、「ある特定の値以上であることはわかるけれども、本当の値はわからない」という状況に出くわすことがあります。このようなデータタイプは「打ち切り」と呼ばれ、統計モデルを組むうえで厄介な問題です。該当する例はとして、以下のようなものが挙げられます。

  • ある一定量のエサを与え、一定時間後に無くなった量を調べる(全部消費された場合、もともと与えたエサ量「以上」を消費したかったことは確か)
  • 標識再捕獲を行い、ある時点で生きていたかどうかを調べる(最後の再捕獲で生きていた場合、その時点より”長く”生きていたことは確か)

他にも例は挙げればたくさんありますが、こういう性質のデータを扱うのは頭を抱えます。

センサリング(Censoring)とは

センサリングとは、そのまま「打ち切り」を意味する言葉なのですが、このテクニックにより上記の問題に対応することが可能です。センサリングを適用する場合、上限値に達した観測値は一時的に欠測、NAとして扱われます。BUGSソフトウェアでは、(通常の線形モデルの場合)応答変数の欠測は自動的に推定されるべきパラメーター(つまり、ほかの偏回帰係数などの推定には用いられない情報)として扱われ、ランダムな値があてがわれます。しかし、打ち切りの場合は部分的な情報(y >= 60など)を持っていますので、上限値に達した値に対しては、こちらで定義する上限値以上の値が(ランダムに?)割り振られます。そうした値が割り振られたうえで、線形モデルの各種パラメーターを推定することになります。

JAGSによる実装

基本的な実装の手順は、普通のBUGSモデルのコードに1行追加するだけで済みます。こちらのモデルファイルを参考にBUGSコードを書いてみます。

model{
#Likelihood
 for(i in 1:N.sample){ # loop for i
   isCensored[i] ~ dinterval(y[i],C[i])
   y[i] ~ dpois(lambda[i]) # poisson model
   log(lambda[i]) <- alpha + beta * X[i] # linear model structure
 }
 
#Prior
 alpha ~ dnorm(0, 0.001)
 beta ~ dnorm(0, 0.001)
}

dintervalが打ち切り関数になっており、isCensored[i](こちらでRから与えるデータ)はデータiが打ち切りかどうかを示す変数(0:打ち切り無し、1:打ち切りあり)、C[i]は打ち切った場合の上限値(こちらデータを与える;任意の変数)になります。なお、打ち切りありの場合、yは”NA”として与えることになります。例えば、観測の上限がy=7だった場合、データは下記のようになります。

> head(D)
   y isCensored C         x
1 NA          1 7 12.056546
2 NA          1 7 11.205303
3  4          0 7  9.865988
4 NA          1 7 12.633539
5 NA          1 7  9.116709
6  6          0 7 12.067588

図にするとこんなふうになります。

打ち切りのあるデータプロット。白丸は真の観察値、グレイは打ち切りされた観察値。

センサリングモデルと普通のポアソンモデルで推定値を比較してみます。

#true values; alpha=0.3, beta=0.15

#with censoring model
      Lower95 Median Upper95 Mean    SD Mode  MCerr MC%ofSD SSeff AC.200 psrf
alpha   0.135   0.59    1.08 0.59 0.241   NA 0.0167     6.9   208  0.100    1
beta    0.083   0.13    0.18 0.13 0.024   NA 0.0017     7.0   202  0.097    1

#with no censoring model
      Lower95 Median Upper95  Mean    SD Mode MCerr MC%ofSD SSeff  AC.200 psrf
alpha   0.587  0.970    1.41 0.972 0.212   NA 0.010       5   407 -0.0027    1
beta    0.026  0.071    0.11 0.071 0.021   NA 0.001       5   403 -0.0037    1

センサリング関数を使わない場合、y>=7の値はy=7とされてしまうので、推定値がかなり過小評価されています。それに対し、センサリング関数を用いたモデルでは、推定値の偏りが解消され、その95%CIも真の値を含んでいます。

注意点:センサリング関数を用いる場合、打ち切りされたデータについては、「打ち切り上限値以上の値」を初期値として与える必要があります。上の例で挙げると、y=NAとなっているデータについては、y=8などの値を初期値として与えます。そうしないとエラーメッセージが帰ってきます。

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です。

時系列データの取り扱い

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

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

それでは、生態学における時系列データはどのような特徴をもつのでしょうか。例えば、ここに個体の移出入はない閉鎖個体群があるとした場合、時間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)を明確に区別でき、より正確なパラメーター推定(個体群増加率など)が可能になります。

Varying slope model: 回帰係数に線形予測子を加える

Varying intercept modelでは切片のみがグループ間で変化することをモデル化していますが、ある変数xに対するyの応答が、グループによって異なることもあります。例えば、魚の種数yは、局所スケールの生息地面積xに従って増えますが、その応答の強さは流域スケールの変数である気温xgで変わるかもしれません。

例:生息地面積xに対する魚の種数yの応答。傾きは流域(色の違い)によって異なっている。

このような「空間スケールをまたがった交互作用」を扱いたい場合は、上の例ので考えると、生息地面積xの効果b(回帰係数)は、流域の平均気温xgによって変わる、という階層構造をもったモデルを構築することになります。

model{
  ninfo <- 0.0001
  
  # Prior
  for(j in 1:2){
    for(k in 1:2){
      B[j,k] ~ dnorm(0,ninfo)
    }
  }
  tau ~ dgamma(ninfo, ninfo); sigma <- sqrt(1/tau)
  TAU[1:2,1:2] ~ dwish(W[,], 3)
  SIGMA[1:2,1:2] <- inverse(TAU[,])
  
  # Likelihood
  for(i in 1:N){
    y[i] ~ dnorm(mu[i], tau)
    mu[i] <- b[group[i],1] + b[group[i],2]*x[i]
  }
  
  for(j in 1:Ng){
    b[j,1:2] ~ dmnorm(b.mu[j,], TAU[,])
    for(k in 1:2){
      b.mu[j,k] <- B[1,k] + B[2,k]*xg[j]
    }
  } 
}

上記のモデルでは、切片と偏回帰係数bがグループ(流域)で変化すると仮定しています。ただし、切片と偏回帰係数において、同じグループ間のばらつきを考えることになるので、それらを独立と考えるのは無理があります。そのため、b[,1]とb[,2]の誤差(グループレベルの説明変数で説明できないばらつき)は、共分散を考慮する必要があります。このスクリプトでは、bの誤差構造には多変量正規分布を用いており(ハイライト箇所21行目)、その分散共分散行列の逆数(TAU[,])の事前分布にウィシャート分布を指定しています(Wはdiag(2)としてRから与えた)。

Varying intercept model: 切片に線形予測子を加える

ランダム効果入りのベイズ推定モデルでは、グループ間の「ランダムなばらつき」を表現するモデルを作りましたが、そのグループ間のばらつきが何によって変わるのかを知りたい場合もあります。例えば、川にすむ魚の個体数をモデリングするとき、小さい空間スケールの個体数は隠れ場所の量によってよく説明できたとしても、河川ごとの平均的な個体数は河川単位で変化する他の要因、例えば降水量で決まるかもしれません(降水量が多い河川は氾濫が多いため、個体数は全体的に少ないなど)。

色の違いが川の違い;X軸は生息場所の量、Y軸は魚の密度。

このような問題には、ランダム効果に線形予測子を加えることで対処できます。下記スクリプトのハイライト箇所で、グループレベルの線形予測子がモデル化されています。なお、このグループレベルの説明変数だけで説明できないばらつきは、パラメーターsigma.Rで考慮されることになり、これによりランダムなグループ間のばらつきも表現されます。

model{
  ninfo <- 0.0001
  
  # Prior
  for(i in 1:2){ b[i] ~ dnorm(0, ninfo) }
  for(i in 1:2){ B[i] ~ dnorm(0, ninfo) }
  tau ~ dgamma(ninfo, ninfo); sigma <- sqrt(1/tau)
  tau.R ~ dgamma(ninfo, ninfo); sigma.R <- sqrt(1/tau.R)
  
  # Likelihood
  for(i in 1:N){
    y[i] ~ dnorm(mu[i], tau)
    mu[i] <- alpha[group[i]] + b[1]*x[i]
  }
  
  for(j in 1:Ng){
    alpha[j] ~ dnorm(mu.R[j], tau.R)
    mu.R[j] <- B[1] + B[2]*xg[j]
  } 
}

なお、注意点として、グループレベルの説明変数は、繰り返しの単位がグループIDと対応することになります。例えば、ここに10の河川から5地点ずつ集められたデータがあるとします(N=5*10=50)。グループレベルでは「河川」が繰り返しとのなるので、添え字jは1から10(Ng: Number of groups)の繰り返し構文となります。したがって、RからJAGSに渡すデータセットは二つのタイプに分かれます(y、x、groupは長さNのベクター、xgは長さNgのベクター)。

ランダム効果入りのベイズ推定モデル

ここでは一般的なGLMMをJAGSで実装するためのスクリプトを紹介します。GLMM、すなわちグループ(ランダム効果)によって切片の値が異なるとするモデルですが、その中身は以下のようなものになります。

Yi ~ Normal(μi, σ2)
μi = α + β*Xi + Rj(i)

上の式ではRj(i)がランダム効果ですが、これは平均0、分散σR2の正規分布に従う確率変数になります。j(グループID)ごとに正規分布に従う乱数が生成され、その値に応じて切片(グループごとのYの平均値)の値が変化します。平均は0に固定されており、問題となるσR2はデータから推定することになります(グループ間のばらつき)。このようなモデルをJAGSで実装する場合、ランダム効果も数字表記(Numeric)でなければならないので、その点に注意が必要です(lmerなどRの関数では因子型として指定;BUGSでは因子型変数は利用不可)。

データセットを用意します。

> head(D,10)
             x group         y
1   1.89336207     1 3.1029784
2   1.96143183     1 3.0784519
3   1.46765994     1 2.8607315
4  -1.61391038     1 2.1712377
5  -0.57316145     1 2.3753531
6   0.70379683     2 0.6382498
7   2.72372835     2 1.0363066
8   0.97064235     2 0.7026342
9   0.02245384     2 0.6256619
10  1.63627604     2 0.8633311

モデルスクリプトを書きます(model_random.Rとして保存)。ランダム効果は入れ子として表記(R[group[i]])されます。データセットはiをユニットとした繰り返しになるのに対し、ランダム効果はjをユニットとした繰り返しになります。

model{
  ninfo <- 0.0001
  
  # Prior
  for(i in 1:2){ b[i] ~ dnorm(0, ninfo) }
  tau ~ dgamma(ninfo, ninfo); sigma <- sqrt(1/tau)
  tau.R ~ dgamma(ninfo, ninfo); sigma.R <- sqrt(1/tau.R) #ランダム効果の精度の(超)事前分布
  
  # Likelihood
  for(i in 1:N){
    y[i] ~ dnorm(mu[i], tau)
    mu[i] <- b[1] + b[2]*x[i] + R[group[i]]
  }
  
  for(j in 1:Ng){ R[j] ~ dnorm(0, tau.R) } 
}

上記のモデルは正規分布になってますが、該当箇所をdpoisやdnegbinとすることで、ポアソンや負の二項分布のモデルでもランダム効果を加えることができます。とてもフレキシブルです。下にいくつかほかの確率分布の場合の例を挙げてみます。モデルを走らせるRコードも下につけますが、詳細はこちらを参照してください。

model{
  ninfo <- 0.0001
  
  # Prior
  for(i in 1:2){ b[i] ~ dnorm(0, ninfo) }
  tau ~ dgamma(ninfo, ninfo); sigma <- sqrt(1/tau)
  tau.R ~ dgamma(ninfo, ninfo); sigma.R <- sqrt(1/tau.R) #ランダム効果の精度の(超)事前分布
  
  # Likelihood
  for(i in 1:N){
    y[i] ~ dpois(mu[i])
    log(mu[i]) <- b[1] + b[2]*x[i] + R[group[i]]
  }
  
  for(j in 1:Ng){ R[j] ~ dnorm(0, tau.R) } 
}
model{
  ninfo <- 0.0001
  
  # Prior
  for(i in 1:2){ b[i] ~ dnorm(0, ninfo) }
  tau ~ dgamma(ninfo, ninfo); sigma <- sqrt(1/tau)
  tau.R ~ dgamma(ninfo, ninfo); sigma.R <- sqrt(1/tau.R) #ランダム効果の精度の(超)事前分布
  size ~ dnorm(0,ninfo)T(0,1000)

  # Likelihood
  for(i in 1:N){
    y[i] ~ dnegbin(p[i],size)
    p[i] <- size/(mu[i]+size)
    log(mu[i]) <- b[1] + b[2]*x[i] + R[group[i]]
  }
  
  for(j in 1:Ng){ R[j] ~ dnorm(0, tau.R) } 
}
###MCMC setting
n.iter <- 3E+4
n.thin <- max(3, ceiling(n.iter/500))
burn <- ceiling(max(10, n.iter/2))
n.ad <- 100
Sample <- ceiling(n.iter/n.thin)

# data for JAGS
Djags <- list(N = nrow(D), Ng = length(unique(D$group)), y = D$y, x = c(scale(D$x, scale=F)), group=D$group)

# parameters to be monitored
para <- c("b", "sigma", "sigma.R")
  
# set initial values
inits <- replicate(3, list(.RNG.name="base::Wichmann-Hill", .RNG.seed=NA), simplify = FALSE)
for(i in 1:3){ inits[[i]]$.RNG.seed <- 2*i+2 }

# Run JAGS
library(runjags)
m <- read.jagsfile("model_random.R")
post <- run.jags(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)