隠れマルコフモデルによる多状態モデルのモデル化

R
Stan
NIMBLE
作者

伊東宏樹

公開

2023年6月18日

更新日

2023年8月27日

『BUGSで学ぶ階層モデリング入門』9章「多状態モデルを用いた捕獲再捕獲による生存と移動の推定」のうち、もっとも簡単なモデルをStanのhmm_marginalを使ってあてはめてみます。比較のため、NIMBLEでもやってみます。

モデル

詳細は、上記書籍をご覧いただきたいのですが、捕獲再捕獲により動物の生存確率と2サイト間の移動確率を推定するというものです。モデルは、以下のような潜在状態と観測値をもつ隠れマルコフモデルになっています。

潜在状態

  1. サイトAで生存
  2. サイトBで生存
  3. 死亡

観測値

  1. サイトAで発見
  2. サイトBで発見
  3. 不発見

データ

使用するライブラリの読み込みなどです。

library(extraDistr)
library(cmdstanr)
## This is cmdstanr version 0.5.3
## - CmdStanR documentation and vignettes: mc-stan.org/cmdstanr
## - CmdStan path: /usr/local/cmdstan
## - CmdStan version: 2.32.2
knitr::knit_engines$set(stan = cmdstanr::eng_cmdstan)
library(nimble)
## nimble version 1.0.1 is loaded.
## For more information on NIMBLE and a User Manual,
## please visit https://R-nimble.org.
## 
## Note for advanced users who have written their own MCMC samplers:
##   As of version 0.13.0, NIMBLE's protocol for handling posterior
##   predictive nodes has changed in a way that could affect user-defined
##   samplers in some situations. Please see Section 15.5.1 of the User Manual.
## 
##  次のパッケージを付け加えます: 'nimble'
##  以下のオブジェクトは 'package:extraDistr' からマスクされています:
## 
##     dcat, dinvgamma, pinvgamma, qinvgamma, rcat, rinvgamma
##  以下のオブジェクトは 'package:stats' からマスクされています:
## 
##     simulate

サイトごとの生存率・移動確率・発見確率と、個体数と観測回数を設定します。

phiA <- 0.8   # survival probability in A
phiB <- 0.7   # survival probability in B
psiAB <- 0.3  # movement probability from A to B
psiBA <- 0.5  # movement probability from B to A
pA <- 0.7     # detection probability in A
pB <- 0.4     # detection probability in A

N <- 160   # number of individuals
T <- 6     # number of occasions

すると、推移行列は以下のようになります。

transition <- matrix(c(phiA * (1 - psiAB), phiA * psiAB,       1 - phiA,
                       phiB * psiBA,       phiB * (1 - psiBA), 1 - phiB,
                       0,                  0,                  1),
                     ncol = 3, nrow = 3, byrow = TRUE)
print(transition)
##      [,1] [,2] [,3]
## [1,] 0.56 0.24  0.2
## [2,] 0.35 0.35  0.3
## [3,] 0.00 0.00  1.0

出力確率は以下のようになります。

emission <- matrix(c(pA,  0,  1 - pA,
                     0,   pB, 1 - pB,
                     0,   0,  1),
                   ncol = 3, nrow = 3, byrow = TRUE)
print(emission)
##      [,1] [,2] [,3]
## [1,]  0.7  0.0  0.3
## [2,]  0.0  0.4  0.6
## [3,]  0.0  0.0  1.0

模擬データを生成します。まずは潜在状態です。

set.seed(1234)

x <- matrix(0, nrow = N, ncol = T)
x[1:N, 1] <- rcat(N, c(0.5, 0.5))

for (n in 1:N) {
  for (t in 2:T) {
    x[n, t] <- rcat(1, transition[x[n, t - 1], ])
  }
}

観測値を生成します。ただし、1回も観測されなかった個体は除外します。

y <- matrix(0, nrow = N, ncol = T)

for (n in 1:N) {
  for (t in 1:T) {
    y[n, t] <- rcat(1, emission[x[n, t], ])
  }
}

detected  <- (apply(y, 1, sum) < T * 3)
y <- y[detected, ]

モデル中で使用するので、最初に観測された調査回を求めておきます。

first <- sapply(1:nrow(y),
                function(i) {
                  which(y[i, ] %in% c(1, 2)) |>
                    min()})

Stan

モデル

Stanのモデルです。hmm_marginalを使用します。このモデルではpAとpBは既知としています。

// Multistate model using hmm_marginal

data {
  int<lower=0> N; // number of individuals 
  int<lower=0> T; // number of occasions
  array[N, T] int<lower=1, upper=3> Y; // observations
  array[N] int<lower=1, upper=T> f;    // first detection
  real<lower=0, upper=1> pA; // detection probability in A
  real<lower=0, upper=1> pB; // detection probability in B
}

transformed data {
  matrix[3, 3] p = [[pA, 0,  1 - pA], // output probability
                    [0,  pB, 1 - pB],
                    [0,  0,  1     ]];
  array[N] matrix[3, T] log_omega;   // log density for each output
  array[N] vector[3] rho;            // initial state

  for (n in 1:N)
    for (t in 1:T)
      log_omega[n, 1:3, t] = log(p[, Y[n, t]]);

  for (n in 1:N)
    if (Y[n, f[n]] == 1)
      rho[n] = [1, 0, 0]';
    else
      rho[n] = [0, 1, 0]';
}

parameters {
  real<lower=0, upper=1> phiA;  // survival probability in A
  real<lower=0, upper=1> phiB;  // survival probability in B
  real<lower=0, upper=1> psiAB; // movement probability from A to B
  real<lower=0, upper=1> psiBA; // movement probability from B to A
}

transformed parameters {
  matrix[3, 3] Gamma =  // transition matrix
    [[phiA * (1 - psiAB), phiA * psiAB,       1 - phiA],
     [phiB * psiBA,       phiB * (1 - psiBA), 1 - phiB],
     [0,                  0,                  1       ]];
}

model {
  for (n in 1:N)
    target += hmm_marginal(log_omega[n, 1:3, f[n]:T], Gamma, rho[n]);
}

あてはめ

fit1 <- model$sample(data = list(N = nrow(y), T = T, Y = y, 
                                  f = first, pA = pA, pB = pB),
                      iter_warmup = 1000, iter_sampling = 1000,
                      chains = 4, parallel_chains = 4,
                      refresh = 0)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 finished in 3.6 seconds.
## Chain 4 finished in 3.5 seconds.
## Chain 2 finished in 3.7 seconds.
## Chain 3 finished in 3.7 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 3.6 seconds.
## Total execution time: 3.8 seconds.

結果

結果です。

fit1$print(c("phiA", "phiB", "psiAB", "psiBA"))
##  variable mean median   sd  mad   q5  q95 rhat ess_bulk ess_tail
##     phiA  0.80   0.80 0.04 0.04 0.74 0.86 1.00     3529     2613
##     phiB  0.64   0.64 0.06 0.05 0.55 0.73 1.00     3436     2898
##     psiAB 0.32   0.32 0.05 0.05 0.25 0.40 1.00     2911     2821
##     psiBA 0.54   0.54 0.07 0.07 0.42 0.66 1.00     4186     2774

だいたい元の値を再現できました。

pAとpBも推定するモデルも試したのですが、識別可能性の問題が発生するようで、うまくいきませんでした。

NIMBLE

モデル

こちらは、pAとpBも推定します。

multistate_code <- nimbleCode({
  phiA ~ dunif(0, 1)
  phiB ~ dunif(0, 1)
  psiAB ~ dunif(0, 1)
  psiBA ~ dunif(0, 1)
  pA ~ dunif(0, 1)
  pB ~ dunif(0, 1)

  Gamma[1, 1] <- phiA * (1 - psiAB)
  Gamma[1, 2] <- phiA * psiAB
  Gamma[1, 3] <- 1 - phiA
  Gamma[2, 1] <- phiB * psiBA
  Gamma[2, 2] <- phiB * (1 - psiBA)
  Gamma[2, 3] <- 1 - phiB
  Gamma[3, 1] <- 0
  Gamma[3 ,2] <- 0
  Gamma[3, 3] <- 1
  p[1, 1] <- pA
  p[1, 2] <- 0
  p[1, 3] <- 1 - pA
  p[2, 1] <- 0
  p[2, 2] <- pB
  p[2, 3] <- 1 - pB
  p[3, 1] <- 0
  p[3, 2] <- 0
  p[3, 3] <- 1

  for (n in 1:N) {
    z[n, f[n]] <- y[n, f[n]]
    for (t in (f[n] + 1):T) {
      z[n, t] ~ dcat(Gamma[z[n, t - 1], 1:3])
      y[n, t] ~ dcat(p[z[n, t], 1:3])
    }
  }
})

あてはめ

zの初期値を設定して、あてはめます。

z_init <- y
for (n in 1:nrow(y)) {
  z_init[n, y[n, ] == 3] <- rcat(1, c(0.5, 0.5))
}
fit2 <- nimbleMCMC(multistate_code,
                  constants = list(N = nrow(y), T = T, f = first),
                  data = list(y = y),
                  inits = list(z = z_init),
                  niter = 11000,
                  nburnin = 1000,
                  nchains = 3,
                  setSeed = 1,
                  summary = TRUE)
## Defining model
## Building model
## Setting data and initial values
## Running calculate on model
##   [Note] Any error reports that follow may simply reflect missing values in model variables.
## Checking model sizes and dimensions
##   [Note] This model is not fully initialized. This is not an error.
##          To see which variables are not initialized, use model$initializeInfo().
##          For more information on model initialization, see help(modelInitialization).
## Checking model calculations
## [Note] NAs were detected in model variables: phiA, logProb_phiA, phiB, logProb_phiB, psiAB, logProb_psiAB, psiBA, logProb_psiBA, pA, logProb_pA, pB, logProb_pB, Gamma, p, logProb_z, logProb_y.
## Compiling
##   [Note] This may take a minute.
##   [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.
## running chain 1...
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## running chain 2...
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## running chain 3...
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|

結果

結果です。

fit2$summary$all.chains
##            Mean    Median    St.Dev. 95%CI_low 95%CI_upp
## pA    0.6570738 0.6326177 0.10338352 0.5165141 0.9269316
## pB    0.6296234 0.6361255 0.20874289 0.2680916 0.9795913
## phiA  0.7826865 0.7797642 0.03870145 0.7151124 0.8691169
## phiB  0.6248509 0.6262631 0.06243422 0.4980178 0.7427421
## psiAB 0.2470975 0.2134588 0.10639427 0.1191031 0.5084560
## psiBA 0.6062178 0.6137067 0.10121906 0.4021126 0.7841810

pAとpBも含めて、うまく推定できたようです。