TMBによる、変量効果ありのN混合モデル

前回作成したTMBのN混合モデルに変量効果を加えてみました。

準備

前回と同様、TMBとggplot2を読み込みます。擬似乱数も固定します。

library(TMB)
library(ggplot2)
set.seed(123)

データ

模擬データを生成します。

R <- 80           # Number of replications (sites)
T <- 3            # Number of survey (times)
G <- 10           # Number of groups
log_lambda <- 1.1 # Log mean of N individuals (lambda = 3.0)
p <- 0.6          # Detection probability
G_ind <- rep(0:(G - 1), each = R / G)  # Group index
log_sigma <- -0.7 # Log sd of group random effect (sigma = 0.50)
ranef <- rnorm(G, 0, exp(log_sigma)) # Group random effect

N <- rpois(R, exp(log_lambda + ranef[G_ind])) # True N
Y <- sapply(1:R, function(i) rbinom(T, N[i], p)) |>
  t()             # Observed N

データのプロット

作成した模擬データを確認します。横軸はサイト番号、黒丸は真の個体数、赤丸は観測個体数です。

data.frame(Site = 1:R, N = N) |>
  ggplot(mapping = aes(x = Site, y = N)) +
  geom_point(size = 2.5) +
  geom_point(data = data.frame(Site = rep(1:R, T),
                               Y = c(Y)),
             mapping = aes(x = Site, y = Y),
             color = "red", alpha = 0.6) +
  scale_x_continuous(breaks = seq(0, R, R / G),
                     minor_breaks = NULL) +
  scale_y_continuous(breaks = seq(0, max(N), 2))

モデル

TMBのC++のコードです。"models/nmix_ranef.cpp"に保存しておきます。

群変量効果のepsilonが加わっています。

観測値データのYは、前回のモデルではDATA_MATRIX型にしていましたが、整数値なので今回はDATA_IMATRIX型にしました。その影響で、各サイトでの観測個体数の最大値を求めるのに前回使っていたTMBのmax関数が使えなくなりました。なので、自前で計算するようにしています（23〜27行）。

また、dbinom関数中ではType型に変換しています（34行）。

nmix_ranef.cpp

// N-mixture model
#include <TMB.hpp>

template<class Type>
Type objective_function<Type>::operator() ()
{
  DATA_IMATRIX(Y);           // Counted numbers of individuals
  DATA_IVECTOR(G_i);         // Group index
  DATA_INTEGER(K);           // Maximum number of individuals
  int R = Y.rows();          // Number of replications (sites)
  int T = Y.cols();          // Times of survey occasions
  PARAMETER(log_lambda);     // Log mean number of individuals
  PARAMETER(p);              // Detection probability
  PARAMETER_VECTOR(epsilon); // Random effect
  PARAMETER(log_sigma);      // log s.d. of random effect
  Type lp = 0;               // negative log likelihood

  lp += -sum(dnorm(epsilon, 0, exp(log_sigma), true));
  for (int i = 0; i < R; i++) {
    vector<Type> pr(K);
    int max_y = 0;

    for (int j = 0; j < Y.cols(); j++) {
      if (max_y < Y(i, j)) {
        max_y = Y(i, j);
      }
    }
    pr.setZero();
    for (int k = 0; k < K; k++) {
      if (k >= max_y) {
        Type theta = exp(log_lambda + epsilon(G_i(i)));
        pr(k) = dpois(Type(k), theta, false);
        for (int t = 0; t < T; t++) {
          pr(k) *= dbinom(Type(Y(i, t)), Type(k), p, false);
        }
      }
    }
    lp += -log(sum(pr));
  }
  return lp;
}

コンパイルと最適化

モデルをコンパイルして、できたライブラリをロードします。

model_name <- "nmix_ranef"
file.path("models", paste(model_name, "cpp", sep = ".")) |>
  compile()
file.path("models", dynlib(model_name)) |>
  dyn.load()

最適化を実行します。Kの値には200を与えました。MakeADFun関数のrandom引数で変量効果のepsilonを指定します。

data <- list(Y = Y, G_i = G_ind, K = 200)
parameters <- list(log_lambda = 0, p = 0.5,
                   epsilon = rep(0, G),
                   log_sigma = 0)
obj_nmix <- MakeADFun(data, parameters, DLL = model_name,
                      random = "epsilon")
opt_nmix <- nlminb(obj_nmix$par, obj_nmix$fn, obj_nmix$gr)

結果

結果です。推定されたパラメータの値は、だいたいデータを生成した元の値を再現できていました。

print(opt_nmix)
## $par
## log_lambda          p  log_sigma 
##  1.0174550  0.6564568 -0.8788686 
## 
## $objective
## [1] 349.2008
## 
## $convergence
## [1] 0
## 
## $iterations
## [1] 15
## 
## $evaluations
## function gradient 
##       21       16 
## 
## $message
## [1] "relative convergence (4)"