Stanによる標準距離標本抽出法データの解析(2)

R
Stan
作者

伊東宏樹

公開

2025年3月2日

更新日

2025年3月4日

まえにも1回やりましたが、『生態学のための階層モデリング』第8章8.3.1にある、標準距離標本抽出法(conventional distance sampling)で採取したデータの解析をふたたびStanでやってみました。本にある実データを使用して、ビン分割データの解析もやってみます。距離データで解析するほうはまえとほぼ同じです。

今回はRStanを使用します。

library(rstan)
## 要求されたパッケージ StanHeaders をロード中です
## 
## rstan version 2.32.6 (Stan version 2.32.2)
## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)
## For within-chain threading using `reduce_sum()` or `map_rect()` Stan functions,
## change `threads_per_chain` option:
## rstan_options(threads_per_chain = 1)
library(ggplot2)
set.seed(123)

データ

長さ60kmのトランセクトの中央から両側で観測されたインパラまでの距離(m)がデータとなります。このデータからトランセクト内のインパラの個体数（と密度）を推定します。

X <- c(71.93, 26.05, 58.47, 92.35, 163.83, 84.52, 163.83, 157.33,
       22.27, 72.11, 86.99, 50.8, 0, 73.14, 0, 128.56, 163.83, 71.85,
       30.47, 71.07, 150.96, 68.83, 90, 64.98, 165.69, 38.01, 378.21,
       78.15, 42.13, 0, 400, 175.39, 30.47, 35.07, 86.04, 31.69, 200,
       271.89, 26.05, 76.6, 41.04, 200, 86.04, 0, 93.97, 55.13, 10.46,
       84.52, 0, 77.65, 0, 96.42, 0, 64.28, 187.94, 0, 160.7, 150.45,
       63.6, 193.19, 106.07, 114.91, 143.39, 128.56, 245.75, 123.13,
       123.13, 153.21, 143.39, 34.2, 96.42, 259.81, 8.72)

ヒストグラムを表示します。

ggplot(data.frame(X = X), aes(x = X)) +
  geom_histogram(binwidth = 50, center = 24.999) +
  labs(x = "Distance (m)", y = "Frequencty")

トランセクト幅の半分(m)（観測距離の最大より大きく設定）は500mとします。

B <- 500

観測された個体数です。

(N_obs <- length(X))
## [1] 73

データ拡大の設定

200個のゼロデータでデータ拡大をおこないます。

N_aug <- 200

距離データをそのまま使う場合

モデル

距離データをそのままデータとして使用する場合は以下のようなモデルになります。

発見確率pは、半正規関数にしたがって減少していくと仮定してモデリングしています。

コンパイルして、model1というオブジェクトに格納しておきます。

data {
  int<lower=1> N_obs;
  int<lower=1> N_aug;
  vector<lower=0>[N_obs] X;
  real<lower=0> B;
}

transformed data {
  real<lower=0> Area = B * 2 / 1000 * 60;  // 調査面積(km^2; Length=60km)
}

parameters {
  real<lower=0> sigma;                     // 半正規関数の尺度
  real<lower=0, upper=1> psi;              // データ拡大パラメータ
  vector<lower=0, upper=B>[N_aug] x_aug;   // 拡大個体の距離
}

transformed parameters {
  vector<lower=0, upper=1>[N_obs + N_aug] p; // 発見確率

  for (i in 1:N_obs) {
    p[i] = exp(-X[i]^2 / (2 * sigma^2));
  }
  for (i in (N_obs + 1):(N_obs + N_aug)) {
    p[i] = exp(-x_aug[i - N_obs]^2 / (2 * sigma^2));
  }
}

model {
  x_aug ~ uniform(0, B);
  sigma ~ uniform(0, 1000);

  // 観測データ
  target += bernoulli_lpmf(1 | psi * p[1:N_obs]);

  // 拡大データ
  for (i in (N_obs + 1):(N_obs + N_aug)) {
    target += log_sum_exp(bernoulli_lpmf(0 | psi),
                          bernoulli_lpmf(1 | psi)
                          + bernoulli_lpmf(0 | p[i]));
  }
}

generated quantities {
  int<lower=0, upper=N_obs + N_aug> N = N_obs;
  real D; // 密度(個体/km^2)

  for (i in (N_obs + 1):(N_obs + N_aug)) {
    real lpnz =  bernoulli_lpmf(1 | psi) + bernoulli_lpmf(0 | p[i]);
    real lp = lpnz - log_sum_exp(lpnz, bernoulli_lpmf(0 | psi));
    N = N + bernoulli_rng(exp(lp));
  }
  D = N / Area;
}

あてはめと結果

データを用意してrstan::sampling関数であてはめます。N, sigma, Dの事後分布の要約を表示します。

stan_data <- list(N_obs = N_obs, N_aug = N_aug, X = X, B = B)
inits <- function() {
  list(psi = runif(1, 0, 1),
       sigma = runif(1, 40, 200))
}
fit1 <- sampling(model1, data = stan_data, init = inits,
                 chains = 4, iter = 4000, warmup = 2000,
                 refresh = 0)
print(fit1, pars = c("N", "sigma", "D"))
## Inference for Stan model: anon_model.
## 4 chains, each with iter=4000; warmup=2000; thin=1; 
## post-warmup draws per chain=2000, total post-warmup draws=8000.
## 
##         mean se_mean    sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
## N     222.92    0.51 24.38 175.00 206.00 223.00 241.00 267.00  2256    1
## sigma 131.31    0.16 10.69 112.72 123.89 130.36 137.94 154.29  4389    1
## D       3.72    0.01  0.41   2.92   3.43   3.72   4.02   4.45  2256    1
## 
## Samples were drawn using NUTS(diag_e) at Sun Mar  2 12:01:28 2025.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

だいたい本にあるWinBUGSと同様の結果が得られました。

ビン分割されたデータの場合

次に、距離そのままではなく、距離が区間ごとの階級データになっている場合です。

データ

区間の幅を50mとします。

Delta <- 50

区間の分割点です。

xg <- seq(0, B, Delta)

距離データを距離カテゴリーに変換します。

D_class <- X %/% Delta + 1

区間数です。

(N_D <- length(xg) - 1)
## [1] 10

図示するとこのようなデータになりました。

factor(D_class, levels = 1:N_D) |>
  table() |> as.data.frame() |>
ggplot(aes(x = Var1, y = Freq)) +
  geom_col() +
  labs(x = "Distance class", y = "Frequency")

距離区間の中点を計算しておきます。

Mid_pt <- xg[-1] - Delta / 2

モデル

ビン分割されたデータ用のモデルです。コンパイルして、model2というオブジェクトに格納しておきます。

data {
  int<lower=1> N_obs;
  int<lower=1> N_aug;
  real<lower=0> B;
  int<lower=1> N_D;
  array[N_obs] int<lower=1, upper=N_D> D_class;
  real<lower=0, upper=B> Delta;
  vector<lower=0, upper=B>[N_D] Mid_pt;
}

transformed data {
  real<lower=0> Area = B * 2 / 1000 * 60;  // 調査面積(km^2; Length=60km)
  simplex[N_D] pi = rep_vector(Delta / B, N_D); // 各区間に含まれる確率
}

parameters {
  real<lower=0> sigma;          // 半正規関数の尺度
  real<lower=0, upper=1> psi;   // データ拡大パラメータ
}

transformed parameters {
  vector<lower=0, upper=1>[N_D] p = exp(-(Mid_pt .* Mid_pt) // 発見確率
                                    / (2 * sigma^2));
}

model {
  sigma ~ uniform(0, 1000);

  // 観測データ
  for (i in 1:N_obs) {
    target += bernoulli_lpmf(1 | psi)
              + bernoulli_lpmf(1 | p[D_class[i]])
              + categorical_lpmf(D_class[i] | pi);
  }

  // 拡大データ
  for (i in (N_obs + 1):(N_obs + N_aug)) {
    real lp[N_D + 1];
    lp[1] = bernoulli_lpmf(0 | psi);
    for (g in 1:N_D) {
      lp[g + 1] = bernoulli_lpmf(1 | psi)
                  + bernoulli_lpmf(0 | p[g])
                  + categorical_lpmf(g | pi);
    }
    target += log_sum_exp(lp);
  }
}

generated quantities {
  int<lower=0, upper=N_obs + N_aug> N = N_obs;
  real D;

  for (i in 1:N_aug) {
    real lp;
    real lp2[N_D];

    for (n in 1:N_D) {
      lp2[n] = bernoulli_lpmf(0 | p[n]) + categorical_lpmf(n | pi);
    }
    lp = bernoulli_lpmf(1 | psi) + log_sum_exp(lp2)
         - log_sum_exp(bernoulli_lpmf(1 | psi) + log_sum_exp(lp2),
                       bernoulli_lpmf(0 | psi));
    N = N + bernoulli_rng(exp(lp));
  }
  D = N / Area;
}

あてはめと結果

モデルをあてはめて、結果を表示します。

stan_data <- list (N_obs = N_obs, N_aug = N_aug, D_class = D_class, B = B,
                   Delta = Delta, N_D = N_D, Mid_pt = Mid_pt)
inits <- function() {
  list(psi = runif(1, 0, 1),
       sigma = runif(1, 40, 200))
}
fit2 <- sampling(model2, data = stan_data, init = inits,
             chains = 4, iter = 4000, warmup = 2000,
             refresh = 0)
print(fit2, par = c("N", "sigma", "D"))
## Inference for Stan model: anon_model.
## 4 chains, each with iter=4000; warmup=2000; thin=1; 
## post-warmup draws per chain=2000, total post-warmup draws=8000.
## 
##         mean se_mean    sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
## N     219.23    0.56 25.51 170.00 201.00 219.00 238.00 267.00  2080    1
## sigma 134.31    0.29 11.85 114.75 126.11 133.03 141.17 161.94  1725    1
## D       3.65    0.01  0.43   2.83   3.35   3.65   3.97   4.45  2080    1
## 
## Samples were drawn using NUTS(diag_e) at Sun Mar  2 12:02:17 2025.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

こちらも、だいたい本にあるWinBUGSと同様の結果が得られました。