NIMBLEによる条件付き自己回帰 (CAR) モデル

R
NIMBLE
作者

伊東宏樹

公開

2023年6月20日

久保さんの『データ解析のための統計モデリング入門』11章「空間構造のある階層ベイズモデル」にある、条件付き自己回帰 (CAR) モデルのあてはめを、NIMBLEで再現してみます。

モデルの詳細は、同書をご覧ください。

データ解析のための統計モデリング入門

準備

まず、使用するパッケージの読み込みと疑似乱数系列の固定です。

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:stats' からマスクされています:
## 
##     simulate
library(coda)
library(posterior)
## This is posterior version 1.4.1
## 
##  次のパッケージを付け加えます: 'posterior'
##  以下のオブジェクトは 'package:stats' からマスクされています:
## 
##     mad, sd, var
##  以下のオブジェクトは 'package:base' からマスクされています:
## 
##     %in%, match
library(stringr)
library(ggplot2)
set.seed(123)

データ

同書のサポートページからデータをダウンロードします。

# data
data_url <- "https://kuboweb.github.io/-kubo/stat/iwanamibook/fig/spatial/Y.RData"
data_file <- "Y.RData"

if (!file.exists(data_file)) {
  download.file(data_url, data_file, "auto")
}

load(data_file)

データの確認

データを確認します。1〜50の番号がついた場所が1次元で並んでいて、各場所で、ある生物の個体数が得られているというデータです。

N_site <- length(Y)
p <- ggplot(data.frame(location = seq_len(N_site), y = Y)) +
  geom_point(aes(x = location, y = y)) +
  ylim(0, 25) +
  labs(x = expression(paste("Location ", italic(j))),
       y = expression(paste("Population ", italic(y)[italic(j)]))) +
  theme_bw(base_family = "Helvetica")
print(p)

次に、空間自己相関の計算につかう数値を用意します。各場所(1〜N_site)について順番に、adjに隣接する場所の番号、weightsに重み（ここではすべて1）、numに隣接する場所の数を入れておきます。

adj <- c(2,
         sapply(2:(N_site - 1), function(i) c(i - 1, i + 1)) |> c(),
         N_site - 1)
weights <- rep(1, length(adj))
num <- c(1, rep(2, N_site - 2), 1)

モデル

NIMBLEのモデルです。CARモデルにはいくつかアルゴリズムがありますが、ここではintrinsic Gaussian CARモデルを使用しています。WinBUGSではcar.normalという名前でしたが、NIMBLEではdcar_normalという名前になっています。

car_code <- nimbleCode({
  for (j in 1:N_site) {
    Y[j] ~ dpois(mean[j])
    log(mean[j]) <- beta + r[j]
  }

  r[1:N_site] ~ dcar_normal(adj[1:L], weights[1:L], num[1:N_site],
                            tau, zero_mean = 1)
  beta ~ dnorm(0, 1.0e-4)
  tau <- 1 / (s * s)
  s ~ dunif(0, 1.0e+4)
})

初期値設定のための関数を定義します。

init_fun <- function() {
  list(beta = runif(1, -2, 2),
       s = runif(1, 0, 2),
       r = runif(N_site, -2, 2))
}

MCMCの実行です。

car_mcmc <- nimbleMCMC(car_code,
                       constants = list(N_site = N_site, L = length(adj),
                                        adj = adj, weights = weights,
                                        num = num),
                       data = list(Y = Y),
                       inits = init_fun,
                       monitors = c("beta", "s", "mean"),
                       niter = 3000, nburnin = 1000,
                       nchains = 3,
                       samplesAsCodaMCMC = 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
## Checking model calculations
## 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...
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|

結果

結果をposteriorパッケージのdraws_df形式に変換して、要約をcar_summaryに格納します。betaとsについて値を表示します。

car_summary <- car_mcmc |>
  as_draws() |>
  summarise_draws()
car_summary |>
  dplyr::filter(variable %in% c("beta", "s"))
## # A tibble: 2 × 10
##   variable  mean median     sd    mad    q5   q95  rhat ess_bulk ess_tail
##   <chr>    <num>  <num>  <num>  <num> <num> <num> <num>    <num>    <num>
## 1 beta     2.27   2.27  0.0492 0.0523 2.19  2.35   1.00     437.     786.
## 2 s        0.229  0.225 0.0483 0.0476 0.155 0.314  1.02     153.     399.

Yの期待値(mean)の事後分布の中央値および5%, 95%分位点を取り出したものに、場所の番号の列を加えます。

mean_summary <- car_summary |>
  dplyr::filter(str_starts(variable, "mean")) |>
  dplyr::select(median, q5, q95) |>
  cbind(location = 1:N_site)

プロット

元の図に、期待値の事後分布の中央値および90%区間を重ねます。

p +
  geom_ribbon(data = mean_summary,
              mapping = aes(x = location,
                            ymin = q5,
                            ymax = q95),
              alpha = 0.5) +
  geom_line(data = mean_summary,
            mapping = aes(x = location, y = median))