閉鎖個体群と開放個体群に空間明示捕獲再捕獲モデルを適用するためのRのパッケージopenCRを使ってみた。登場してからあまり時間が経っておらずネットにもあまり情報が無いので、使った例を紹介する。
使うデータの生成
ここでは、自動撮影カメラを用いてある哺乳類種を撮影したデータを想定する。自動撮影カメラが格子状に設置されており、撮影画像から対象種の個体識別を行い、個体毎に撮影枚数が集計されているものとする。
以下の例では、便宜的に調査機会が4回あるが、調査期間全体を通じて閉鎖個体群であると仮定する（openCRなのにというツッコミはありそうだが^^;）。また、自動撮影カメラは数が限られているため、全ての調査機会において全ての箇所に設置できているわけではない。そのため、調査努力量は変動している。
################
# データ生成
################
set.seed(2)
## 個体数を推定したい空間的な範囲を設定
xl <- -15
xu <- 15
yl <- -15
yu <- 15
# 推定範囲に存在する真の個体数
N <- 20
# 個体ごとの活動中心の位置の設定
ac_x <- runif(N, xl, xu)
ac_y <- runif(N, yl, yu)
## 調査機会数
Nso <- 4
## カメラのID、x座標、y座標を設定
Cam_id <- paste("C", formatC(1:100, width=2, flag="0"), sep="")
cam_x <- rep(((-5):4)*2+1, 10)
cam_y <- rep(((-5):4)*2+1, each=10)
Ncam <- length(cam_x)
## 調査機会かつカメラごとの稼働状況
Effort <- matrix(rpois(Ncam*Nso, 10), nrow=Ncam, ncol=Nso)
## 個体（i）、カメラ（j）、調査機会（k）ごとの撮影枚数データの生成
# カメラとの距離が0の場合の検出率
lam0 <- 0.05
# カメラからの距離に応じた検出率の減衰関数である
# 半世紀関数の分散パラメータ
sigma <- 1.5
# 各個体の活動中心と各カメラの距離を格納する箱
D2 <- matrix(0, ncol=Ncam, nrow=N)
# 検出率を格納する箱
lambda <- matrix(0, ncol=Ncam, nrow=N)
# 生成する撮影枚数データを格納する箱
y <- array(dim=c(N, Ncam, Nso))
# 撮影枚数データの生成
for (i in 1:N) { # 個体
for (j in 1:Ncam) { # カメラ
# 各個体の活動中心と各カメラの距離の2乗を計算
D2[i,j] <- (ac_x[i]-cam_x[j])^2 + (ac_y[i]-cam_y[j])^2
# 半正規関数による、距離に応じた検出率の減衰
lambda[i,j] <- lam0*exp(-D2[i,j]/(2*sigma*sigma))
for (k in 1:Nso) { # 調査機会
y[i,j,k] <- rpois(1,lambda[i,j]*Effort[j,k])
}
}
}
# 1枚も撮影されていない個体は検出できないので、データから削除
yobs <- y[apply(y,1,sum)>0,,]
作成したデータ
生成したデータを図にしてみる。活動中心を四角で示しているが、白抜きの個体は撮影データが無いことを示している。調査範囲にある全ての個体が撮影されているわけではないことが分かる。
openCRでの個体数推定
openCRはRの別のパッケージsecrを前提に作成されている。そのため、用意するデータの形式も、secrパッケージと共通にする必要がある。基本的な手順は、以下のとおり。
- 観測点（カメラの位置など。openCRのマニュアルではdetectorと呼びます）の定義
- 推定する空間的な範囲を指定する
- 個体毎の遭遇履歴（capthistオブジェクト）の作成
- 個体密度の推定
観測点の定義
観測点（detector）は、今回は自動撮影カメラである。位置情報とともに、trapオブジェクトとして用意する。trapオブジェクトに変換してくれる関数read.traps()は、位置座標を示す列名がxとyでないと動作しないので注意すること。
調査努力量の違いは、各カメラが列、各行が調査機会となっている努力量の行列を、usage()でtrapオブジェクトに付与する。
library(openCR)
trapXY <- data.frame(cam_id=Cam_id, x=cam_x, y=cam_y)
traps <- read.traps(data=trapXY, detector="count")
rownames(traps) <- Cam_id
usage(traps) <- Effort
推定する空間的な範囲を指定する
推定する空間的な範囲は、実は勝手に設定してくれるが、自分で明示的に与えたい場合はspオブジェクトかsfオブジェクトで用意する（若干不確かです）。島などのように範囲が限定されている場合や、特定の場所は活動中心の場所として考えられない場合は、望むようなshpファイルを用意しておくといいだろう。
library(sf)
estrange <- st_polygon(list(rbind(st_point(c(-12, -12)),
st_point(c(-12, 12)),
st_point(c(12, 12)),
st_point(c(12,-12)),
st_point(c(-12, -12))
)
)
)
個体毎の遭遇履歴の作成
個体毎の遭遇履歴は、capthitsオブジェクトで用意する。各行は検出機会である。同じ個体でも、カメラや調査機会が異なれば複数行にデータがあり、同じカメラと調査機会であっても複数枚数撮影されている場合はやはり複数行にデータがある。つまり、行数がそのまま全検出数となっている必要がある。
また、すでに作成したtrapオブジェクトと合わせて、make.capthist()関数でcapthistオブジェクトを作成する。
# 配列を縦方向に展開する
temp <- as.data.frame(yobs) #行は一度でも観測された個体、列はカメラ×調査機会数
colnames(temp) <- paste(rep(Cam_id, Nso), rep(1:Nso, each=Ncam), sep="_")
tl <- reshape(temp, varying=list(1:ncol(temp)),
idvar="individual", ids=1:nrow(temp),
timevar="Cam_so", times=colnames(temp),
direction="long"
)
tl$Session <- 1
tl$trapID <- sapply(strsplit(tl$Cam_so, "\\_"), "[", 1)
tl$Occasion <- as.numeric(sapply(strsplit(tl$Cam_so, "\\_"), "[", 2))
# rep()を使って、撮影枚数だけ行を増やす
tl <- tl[rep(1:nrow(tl), tl$C01_1), ]
colnames(tl)[grep("individual",colnames(tl))] <- "ID"
tl <- tl[, c("Session", "ID", "Occasion", "trapID")]
test <- make.capthist(tl, traps, fmt="trapID")
# 推定範囲を読み込む
msk <- make.mask(traps(test), spacing = 1,
buffer = 3, type = "traprect", poly=estrange
)
個体密度の推定
ここまできたら準備完了。openCR.fit()関数でパラメータ推定を行う。今回は閉鎖個体群なので、モデルの種類を指定する引数typeはsecrDとする。
secrclose <- openCR.fit(test, type = 'secrD',
model=list(lambda0~1), mask = msk,
distribution="poisson", detectfn="HHN", method="Nelder-Mead"
)
推定結果の確認
推定値の要約量は、summary()関数で確認できる。
summary(secrclose)
$versiontime
[1] "1.5.0, run 20:21:01 25 1 2021"
$traps
Detector Number Spacing UsagePct
count 100 2 1019
$capthist
S1 S2 S3 S4 Total
Occasions 1 1 1 1 4
Detections 18 11 9 16 54
Animals 10 7 6 9 13
Detectors 100 100 100 100 100
Moves 4 3 2 3 27
$intervals
numeric(0)
$mask
Cells Spacing Area
576 1 0.0576
$modeldetails
type fixed distribution
secrD phi = 1, b = 1 poisson
$Movementmodel
[1] "static"
$AICtable
model npar rank logLik
AIC AICc
lambda0~1 superD~1 sigma~1 3 3 -196.6385 399.277 401.944
$link
lambda0 phi b superD sigma
log logit mlogit log log
$coef
beta SE.beta
lcl ucl
lambda0 -2.9943305 0.21252066 -3.4108634 -2.5777977
superD 5.5313727 0.27872276 4.9850862 6.0776593
sigma 0.3506291 0.08968483 0.1748501 0.5264082
$hessian
rankH svtol Eigen1 Eigen2 Eigen3
3 1e-05 1 0.09547 0.05959
$predicted
$predicted$lambda0
session estimate SE.estimate lcl ucl
1 1 0.05007014 0.01064094 0.03301269 0.07594106
$predicted$superD
estimate SE.estimate lcl ucl
1 252.4903 70.37479 146.2162 436.0074
$predicted$sigma
session estimate SE.estimate lcl ucl
1 1 1.419961 0.1273489 1.191068 1.692841
様々な情報が記載されているが、例えば設定した個体数を推定できているかを確認するためには、推定密度x面積で個体数を得ればよい。設定値は20なので、妥当な推定結果と言える。
summary(secrclose)$predicted$superD*summary(secrclose)$mask$Area
estimate SE.estimate lcl ucl
1 14.54344 4.053588 8.422052 25.11403
ベイズ推定する場合
上記のモデルは、ベイズ推定することも可能である。ベイズ推定する場合、データ拡大法を使う。実行コードは以下のとおり。
# データ拡大
# 観察個体の4倍の個体を加える
naug <- nrow(yobs)*5
M <- nrow(yobs) + naug
yaug <- array(0,dim=c(M,Ncam,Nso))
yaug[1:nrow(yobs),,] <- yobs
# リスト形式でデータを用意
list.data <- list(y=yaug, M=M, Ncam=Ncam, Nso=Nso,
cam_x=cam_x, cam_y=cam_y, xl=xl, xu=xu, yl=yl, yu=yu,
Effort=as.matrix(Effort), op=ifelse(Effort==0,0,1)
)
# BUGS言語によるモデルの指定
modelFilename = "secr_open_bayes.txt"
cat("
model {
# 含有率の事前分布
psi ~ dunif(0, 1)
# 個体の活動中心の位置の事前分布
for (i in 1:M) {
AC_x[i] ~ dunif(xl, xu)
AC_y[i] ~ dunif(yl, yu)
}
# 個体群が調査範囲に存在するかしないかを推定
for (i in 1:M) {
w[i] ~ dbern(psi)
}
# 検出率の事前分布
lam0 ~ dunif(0, 100)
# 半正規関数の分散パラメータの事前分布
sigma ~ dunif(0, 100)
# 個体の在不在と、カメラからの距離に応じた検出率から、
# 個体（i）かつカメラ（j）かつ調査機会（k）
# ごとの撮影枚数データを説明する
for (i in 1:M) { # 個体
for (j in 1:Ncam) { # 調査機会
# カメラと活動中心の距離の2乗を計算
D2[i,j] <- pow(AC_x[i] - cam_x[j], 2) + pow(AC_y[i] - cam_y[j], 2)
for (k in 1:Nso) {
# 半正規関数による距離に応じた検出率、カメラの稼働状況、個体の在不在
mu[i,j,k] <- lam0*exp(-D2[i,j]/(2*pow(sigma, 2)))*Effort[j,k]*w[i]
y[i,j,k] ~ dpois(mu[i,j,k]) # データはポアソン分布に従う
}
}
}
## 導出パラメータ
N <- sum(w[]) # 真の個体数は在（w=1）である個体数
D <- N/((xu-xl)*(yu-yl)) # 個体密度はN/調査範囲
} #
モデルの記述はここまで
", fill=TRUE, file=modelFilename)
# 初期値の設定
init1 <- list(sigma=2, psi=0.5,
AC_x=runif(list.data$M, list.data$xl, list.data$xu),
AC_y=runif(list.data$M, list.data$yl, list.data$yu),
w=rep(1, list.data$M), lam0=0.1
)
init2 <- list(sigma=1, psi=0.4,
AC_x=runif(list.data$M, list.data$xl, list.data$xu),
AC_y=runif(list.data$M, list.data$yl, list.data$yu),
w=rep(1, list.data$M), lam0=0.05
)
init3 <- list(sigma=4, psi=0.6,
AC_x=runif(list.data$M, list.data$xl, list.data$xu),
AC_y=runif(list.data$M, list.data$yl, list.data$yu),
w=rep(1, list.data$M), lam0=0.2
)
inits <- list(init1, init2, init3)
inits[[1]]$.RNG.name <- "base::Mersenne-Twister"
inits[[1]]$.RNG.seed <- 1
inits[[2]]$.RNG.name <- "base::Mersenne-Twister"
inits[[2]]$.RNG.seed <- 12
inits[[3]]$.RNG.name <- "base::Mersenne-Twister"
inits[[3]]$.RNG.seed <- 123
# 監視対象パラメータの指定
para <- c("sigma", "AC_x", "AC_y", "lam0", "N", "D")
# MCMCの設定
nc <- 3
ni <- 3000
nb <- 1500
nt <- 3
# jagsUIパッケージの読み込み
library(jagsUI)
# 計算開始時間を記録
start.time <- Sys.time()
# MCMC法の実行
outSECR <- jags(list.data, inits, parameters.to.save=para, parallel=TRUE,
modelFilename, n.chains=nc, n.thin=nt,
n.iter=ni, n.burnin=nb
)
#終了時間の記録と、計算時間の出力
end.time <- Sys.time()
elapsed.time <- difftime(end.time, start.time, units='hours')
cat(paste(paste('Posterior computed in ', elapsed.time, sep=''), ' hours\n', sep=''))
推定結果を確認する。全てのパラメータはRhat値が1.1以下となり収束したようだ。設定値との関係も確認する。左上の図が個体数（青の縦線が設定値、赤の縦線がデータ拡大した範囲、黒の点線が観測された個体数）。おおむね設定した個体数を推定できていると言えるだろう。
# Rhat値の確認
res <- as.data.frame(outSECR$summary)
res[res$Rhat > 1.1, ]
ベイズ推定の場合、最尤推定（openCR）では積分消去されてしまう活動中心のパラメータも取り出すことが可能である。観測値が得られている個体の活動中心は、設定した位置の周辺で推定できているようである（エラーバーは95%信用区間）。
