TMBによるN混合モデル – 伊東生態統計研究室

TMBでN混合モデル(N-mixture model)を実装してみました。TMBのexamplesにも実装例はあるのですが、R風の密度関数をつかって、より簡潔に書いてみました。

準備

例によって、TMBとggplot2を読み込みます。擬似乱数も固定します。

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

データ

R <- 40      # Number of replications (sites)
T <- 3       # Times of surveys
lambda <- 3  # Poisson mean of individuals
p <- 0.6     # Detection probability

# Number of individuals for each site
N <- rpois(R, lambda)

# Counted number of individuals for each site and survey
Y <- sapply(1:R, function(i) rbinom(T, N[i], p)) |>
  t()

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

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)

モデル

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

Royle (2004)の式(3)をだいたいそのままTMBのコードにしています。真の個体数の上限を無限にするのは不可能なので、上限をKとしています。

TMBの型の扱いがまだよくわかっていないので、型変換でちょっと苦労しました。

nmix.cpp

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

template<class Type>
Type objective_function<Type>::operator() ()
{
  DATA_MATRIX(Y);    // Counted numbers of individuals
  DATA_INTEGER(K);   // Maximum number of individuals
  int R = Y.rows();  // Number of replications (sites)
  int T = Y.cols();  // Times of survey occasions
  PARAMETER(lambda); // Mean number of individuals
  PARAMETER(p);      // Detection probability
  Type lp = 0;       // negative log likelihood

  for (int i = 0; i < R; i++) {
    vector<Type> pr(K);
    Type max_y = max((vector<Type>)Y.row(i));

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

コンパイルと最適化

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

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

最適化を実行します。

Kの値には200を与えました。

data <- list(Y = Y, K = 200)
parameters <- list(lambda = 1, p = 0.5)
obj_nmix <- MakeADFun(data, parameters, DLL = model_name)
opt_nmix <- nlminb(obj_nmix$par, obj_nmix$fn, obj_nmix$gr)

結果

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

このアルゴリズムをStanで実装してMCMCでパラメータ推定すると、それなりに時間がかかるのですが、TMBによる最尤推定は非常に高速でした。

print(opt_nmix)
## $par
##    lambda         p 
## 3.1408005 0.6314746 
## 
## $objective
## [1] 179.004
## 
## $convergence
## [1] 0
## 
## $iterations
## [1] 9
## 
## $evaluations
## function gradient 
##       12       10 
## 
## $message
## [1] "relative convergence (4)"

参考文献

Royle, J.A (2004) N-mixture models for estimating population size from spatially replicated counts. Biometrics, 60:108-115. doi:10.1111/j.0006-341X.2004.00142.x