unmarkedでTMBを使う – 伊東生態統計研究室

N-mixtureモデルや占有モデルなど、生物の個体数推定の階層モデルをあつかうunmarkedパッケージで、エンジンにTMBが使えるようになっていましたので、試してみました。

準備

unmarkedパッケージを読み込みます。

library(unmarked)

データ

データとして、“Bayesian Population Analysis using WinBUGS”（日本語訳「BUGSで学ぶ階層モデリング入門」）の12.2.2節のデータを用います。共変量のあるN-mixtureモデル（二項混合モデル）です。

データを読み込んだら、unmarkedFramePCount関数で、N-mixtureモデルへのあてはめをおこなうpcount関数用のデータを作成します。

url <- "https://raw.githubusercontent.com/stan-dev/example-models/refs/heads/master/BPA/Ch.12/binmix_cov.data.R"
source(url)
d <- unmarkedFramePCount(y, siteCovs = data.frame(X = X))

データを確認します。

plot(d)

あてはめ

pcount関数で、N-mixtureモデルにあてはめます。engine = "TMB"とすることで、エンジンにTMBが使用されます。

結果をパイプでsummary関数に渡しています。

pcount(~ X ~ X, data = d, K = 100, mixture = "P", engine = "TMB") |>
summary()
## 
## Call:
## pcount(formula = ~X ~ X, data = d, K = 100, mixture = "P", engine = "TMB")
## 
## Abundance (log-scale):
##             Estimate     SE    z  P(>|z|)
## (Intercept)     1.09 0.0996 10.9 1.08e-27
## X               3.35 0.1883 17.8 8.38e-71
## 
## Detection (logit-scale):
##             Estimate    SE      z  P(>|z|)
## (Intercept)    -0.33 0.145  -2.27 2.29e-02
## X              -5.00 0.265 -18.87 1.94e-79
## 
## AIC: 1101.482 
## Number of sites: 200
## optim convergence code: 0
## optim iterations: 57 
## Bootstrap iterations: 0

推定値は、データ生成に使われた値（alpha0 = 1, alpha1 = 3, beta0 = 0, beta1 = -5）をだいたい再現できています。

TMBを使う利点ですが、ほかのエンジン(“C”および”R”)とくらべて、TMBの方がかなり高速のように思われます。