「曲がりくねった共変量関係」をStanを使ってガウス過程でモデル化

R
Stan
作者

伊東宏樹

公開

2025年3月5日

『生態学のための階層モデリング』第10章10.4節にある「曲がりくねった共変量関係」をStanを使ってモデル化してみました。本ではBUGSを使って罰則つきスプラインでモデル化していますが、ここではガウス過程でモデル化しました。

データ

データはAHMbookパッケージを利用して生成します。その他、ggplot2とrstanパッケージを読み込みます。

library(AHMbook)
library(ggplot2)
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)

データの生成

wigglyOcc関数でデータを生成します。サイト占有モデルで、占有確率と検出確率がそれぞれ共変量に対して曲がりくねった関係となっているというものです。

data <- wigglyOcc(seed = 1, show.plot = FALSE, verbose = TRUE)
##    True number of occupied sites: 141 
##    Observed number of occupied sites: 114 
##    Proportional underestimation of distribution: 0.19

データには以下のような要素が含まれています。

  • M: サイトの数
  • J: サイトあたりの反復調査の数
  • Xsite: サイト共変量
  • Xsurvey: 調査共変量
  • psi: サイトの占有確率
  • z: 在・不在
  • p: 検出確率
  • y: 観測値（調査ごとの検出・不検出）
summary(data)
##           Length Class  Mode   
## M           1    -none- numeric
## J           1    -none- numeric
## Xsite     240    -none- numeric
## Xsurvey   720    -none- numeric
## psi       240    -none- numeric
## z         240    -none- numeric
## p         720    -none- numeric
## y         720    -none- numeric
## x.index   720    -none- numeric
## p.ordered 720    -none- numeric

データの確認

サイト共変量Xsiteと占有確率psi（赤線）、在・不在z（点）の関係を図示します。

このように占有確率は共変量に対して曲がりくねっています。

site_df <- data.frame(x = data$Xsite, psi = data$psi, z = data$z)
plot_site <- ggplot(site_df) +
  geom_point(aes(x, z), color = "black", size = 3, alpha = 0.5) +
  geom_line(aes(x, psi),
            linewidth = 1, colour = "red") +
  labs(x = "Site covariate (Xsite)", y = "Occupancy prob.") +
  theme_bw()
plot(plot_site)

調査共変量Xsurveyと検出確率p（赤線）、観測値y（点）の関係を図示します。

このように検出確率も共変量に対して曲がりくねっています。

survey_df <- data.frame(x = c(data$Xsurvey), y = c(data$y), p = c(data$p))
plot_survey <- ggplot(survey_df) +
  geom_point(aes(x, y),
              colour = "black", size = 3, alpha = 0.5) +
  geom_line(aes(x, p),
            linewidth = 1, colour = "red") +
  labs(x = "Survey covariate (Xsurvey)", y = "Detection prob.") +
  theme_bw()
plot(plot_survey)

モデル

このような関係をガウス過程でモデル化します。Stanのコードは、User’s guideのGaussian Processesを参考にしています。

このモデルをコンパイルして、stan_modelというオブジェクトに格納しておきます。

data {
  int<lower=1> N_site1;         // Number of sites
  int<lower=1> N_rep;           // Number of replications
  int<lower=1> N_obs1;          // Number of obs. (N_site1 x N_rep)
  vector[N_site1] X_site1;      // Site covariate
  vector[N_obs1] X_survey1;     // Survey covariate
  array[N_site1, N_rep] int<lower=0, upper=1> Y; // Observations
  int<lower=1> N_site2;         // Number of predictions (site)
  int<lower=1> N_survey2;       // Number of predictions (survey)
  vector[N_site2] X_site2;      // Covarites for predictions (site)
  vector[N_survey2] X_survey2;  // Covarites for predictions (survey)
}

transformed data {
  real delta = 1e-9;
  int N_site = N_site1 + N_site2;
  int N_obs = N_obs1 + N_survey2;
  array[N_site] real X_site;
  array[N_obs] real X_survey;
  array[N_site1] int Y_sum;

  for (n in 1:N_site1)
    X_site[n] = X_site1[n];
  for (n in 1:N_site2)
    X_site[N_site1 + n] = X_site2[n];
  for (n in 1:N_obs1)
    X_survey[n] = X_survey1[n];
  for (n in 1:N_survey2)
    X_survey[N_obs1 + n] = X_survey2[n];
  for (n in 1:N_site1)
    Y_sum[n] = sum(Y[n]); 
}

parameters {
  // for occupancy parameter (psi)
  real<lower=0> rho_psi;
  real<lower=0> alpha_psi;
  real a_psi;
  vector[N_site] eta_psi;
  // for detection parameter (p)
  real<lower=0> rho_p;
  real<lower=0> alpha_p;
  real a_p;
  vector[N_obs] eta_p;
}

transformed parameters {
  vector[N_site1] logit_psi;       // Logit occ. prob.
  vector[N_site2] logit_psi2;
  matrix[N_site1, N_rep] logit_p;  // Logit detection prob.
  vector[N_survey2] logit_p2;

  { // for occupancy parameter (psi)
    vector[N_site] f;
    matrix[N_site, N_site] L_K;
    matrix[N_site, N_site] K = gp_exp_quad_cov(X_site, alpha_psi, rho_psi);

    // diagonal elements
    for (n in 1:N_site)
      K[n, n] = K[n, n] + delta;
    L_K = cholesky_decompose(K);
    f = L_K * eta_psi;
    logit_psi = a_psi + f[1:N_site1];
    logit_psi2 = a_psi + f[(N_site1 + 1):N_site];
  }

  { // for detection parameter (p)
    vector[N_obs] f;
    matrix[N_obs, N_obs] L_K;
    matrix[N_obs, N_obs] K = gp_exp_quad_cov(X_survey, alpha_p, rho_p);

    // diagonal elements
    for (n in 1:N_obs)
      K[n, n] = K[n, n] + delta;
    L_K = cholesky_decompose(K);
    f = L_K * eta_p;
    logit_p = to_matrix(a_p + f[1:N_obs1], N_site1, N_rep);
    logit_p2 = a_p + f[(N_obs1 + 1):N_obs];
  }
}

model {
  rho_psi ~ inv_gamma(5, 5);
  alpha_psi ~ std_normal();
  a_psi ~ std_normal();
  eta_psi ~ std_normal();
  rho_p ~ inv_gamma(5, 5);
  alpha_p ~ std_normal();
  a_p ~ std_normal();
  eta_p ~ std_normal();
  for (n in 1:N_site1) {
    if (Y_sum[n] > 0) {  // detected
      target += bernoulli_logit_lpmf(Y[n, 1:N_rep] | logit_p[n, 1:N_rep])
                + bernoulli_logit_lpmf(1 | logit_psi[n]);
    } else { // not detected
      target += log_sum_exp(bernoulli_logit_lpmf(0 | logit_psi[n]),
                            bernoulli_logit_lpmf(Y[n, 1:N_rep] |
                                                 logit_p[n, 1:N_rep])
                            + bernoulli_logit_lpmf(1 | logit_psi[n]));
    } 
  }
}

generated quantities {
  vector[N_site2] psi = inv_logit(logit_psi2);
  vector[N_survey2] p = inv_logit(logit_p2);
}

あてはめ

Stanにわたすデータを用意して、モデルをあてはめます。実行にはそれなりに時間がかかります。

new_Xsite <- seq(-2, 2, length = 41)
new_Xsurvey <- seq(-2, 2, length = 41)

stan_data <- list(N_site1 = data$M, N_rep = data$J,
                  N_obs1 = data$M * data$J,
                  X_site1 = data$Xsite,
                  X_survey1 = c(data$Xsurvey),
                  Y = data$y,
                  N_site2 = length(new_Xsite),
                  N_survey2 = length(new_Xsurvey),
                  X_site2 = new_Xsite,
                  X_survey2 = new_Xsurvey)
fit <- rstan::sampling(stan_model, data = stan_data,
                       iter = 2000, warmup = 1000,
                       chains = 4, cores = 4,
                       refresh = 0)

結果

占有確率のガウス過程に関係するパラメータの事後分布の要約です。

print(fit, par = c("rho_psi", "alpha_psi", "a_psi"))
## Inference for Stan model: anon_model.
## 4 chains, each with iter=2000; warmup=1000; thin=1; 
## post-warmup draws per chain=1000, total post-warmup draws=4000.
## 
##           mean se_mean   sd  2.5%   25%  50%  75% 97.5% n_eff Rhat
## rho_psi   0.94    0.01 0.25  0.53  0.76 0.91 1.09  1.51  1072    1
## alpha_psi 1.65    0.01 0.48  0.87  1.30 1.59 1.94  2.71  2807    1
## a_psi     0.09    0.01 0.74 -1.38 -0.39 0.10 0.58  1.51  2506    1
## 
## Samples were drawn using NUTS(diag_e) at Tue Mar  4 20:46:01 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).

検出確率のガウス過程に関係するパラメータの事後分布の要約です。

print(fit, par = c("rho_p", "alpha_p", "a_p"))
## Inference for Stan model: anon_model.
## 4 chains, each with iter=2000; warmup=1000; thin=1; 
## post-warmup draws per chain=1000, total post-warmup draws=4000.
## 
##          mean se_mean   sd  2.5%   25%   50%   75% 97.5% n_eff Rhat
## rho_p    0.54    0.00 0.10  0.37  0.47  0.53  0.60  0.77  1390    1
## alpha_p  1.66    0.01 0.46  0.94  1.33  1.60  1.94  2.68  2179    1
## a_p     -0.63    0.02 0.67 -1.93 -1.09 -0.66 -0.20  0.71  1568    1
## 
## Samples were drawn using NUTS(diag_e) at Tue Mar  4 20:46:01 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).

ガウス過程による占有確率の予測値を青色（線: 中央値、領域: 95%信用区間）で図示します。

psi_summary <- summary(fit, pars = "psi")$summary
pred_df <- data.frame(x = new_Xsite,
                      median = psi_summary[, "50%"],
                      lower = psi_summary[, "2.5%"],
                      upper = psi_summary[, "97.5%"])
plot_site +
  geom_ribbon(data = pred_df,
              mapping = aes(x = x,
                            ymin = lower,
                            ymax = upper),
              fill = "blue", alpha = 0.3) +
  geom_line(data = pred_df,
            mapping = aes(x = x, y = median),
            linewidth = 1, colour = "blue")

ガウス過程による検出確率の予測値を青色（線: 中央値、領域: 95%信用区間）で図示します。

p_summary <- summary(fit, pars = "p")$summary
pred_df <- data.frame(x = new_Xsurvey,
                      median = p_summary[, "50%"],
                      lower = p_summary[, "2.5%"],
                      upper = p_summary[, "97.5%"])
plot_survey +
  geom_ribbon(data = pred_df,
              mapping = aes(x = x,
                            ymin = lower,
                            ymax = upper),
              fill = "blue", alpha = 0.3) +
  geom_line(data = pred_df,
            mapping = aes(x = x, y = median),
            linewidth = 1, colour = "blue")

わりといい感じに推定ができたのではないかと思います。