今度は、1度の捕獲機会において、1トラップで最大でも1個体しか捕獲されず、かつ捕獲された個体が同じ調査機会において他のトラップに捕獲されない状況での個体数推定を考える。ネズミ類などの生体捕獲で標識再捕獲をする場合などが、このような標本抽出になるだろう。
カメラ調査とは異なり、生体捕獲調査では1次調査期間（個体群が開放）が厳密に等間隔にできないことが多い。そして、1次調査期間の中に複数回の反復調査（2次調査期間、個体群は閉鎖）があることが多い。このような、2段階の標本抽出で、かつ1次調査期間の間隔が一定でない場合を考える。
この場合、生存率の扱いに注意が必要。調査期間が違えば生存率もそれに従って変動するので。openCR.fit()では、調査期間の感覚に応じて生存率（phi）は規格化されると記述している。元論文をあたると、おそらく日あたりの生存率を推定し、それを期間長分べき乗している。
データの生成
コードがさらに長くなります（汗
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)
## 1次調査機会数
Npri <- 6
## 2次調査機会数
Nsec <- 2
## トラップのID、x座標、y座標を設定
Trap_id <- paste("C", formatC(1:100, width=2, flag="0"), sep="")
trap_x <- rep(((-5):4)*2+1, 10)
trap_y <- rep(((-5):4)*2+1, each=10)
Ntrap <- length(trap_x)
## 調査機会かつトラップごとの稼働状況
Effort <- matrix(1, nrow=Ntrap, ncol=Npri)
Effort <- Effort[, rep(1:ncol(Effort), each=Nsec)]
## 上記の通り、1次調査期間の間隔が一定でないので、
## 1日の生存率を設定する
phi_day <- 0.999
## 1次調査機会間の間隔（ここで、少し異なる期間を設定します）
Intpri <- rpois(Npri-1, 100)
## 加入率
beta <- c(5, rpois(Npri-1, 3))
b <- numeric()
for (t in 1:Npri) {
b[t] <- beta[t]/sum(beta[1:Npri])
}
# 加入確率に変換
nu <- numeric()
nu[1] <- b[1]
for (t in 2:Npri) {
nu[t] <- b[t]/(1 - sum(b[1:(t-1)]))
}
# 1次調査機会ごとの個体の在不在
z <- matrix(0, ncol=Npri, nrow=N)
z[,1] <- rbinom(N, 1, nu[1])
q <- matrix(0, ncol=Npri-1, nrow=N)
mu2 <- matrix(0, ncol=Npri-1, nrow=N)
for (t in 2:Npri) {
q[, t-1] <- 1 - z[, t-1]
if (t==2) { mu2[, t-1] <- phi_day^Intpri[t-1]*z[, t-1] + nu[t]*q[,1]
} else { mu2[, t-1] <- phi_day^Intpri[t-1]*z[, t-1] + nu[t]*apply(q[,1:(t-1)],1, prod) }
z[,t] <- rbinom(N, 1, ifelse(mu2[,t-1]>1,1,mu2[,t-1]))
}
## 個体（i）、トラップ（j）、調査機会（k）ごとの撮影枚数データの生成
# トラップとの距離が0の場合の検出率
lam0 <- 1
# トラップからの距離に応じた検出率の減衰関数である
# 半世紀関数の分散パラメータ
sigma <- 1.5
# 検出率
p <- 0.9
# 各個体の活動中心と各トラップの距離を格納する箱
D2 <- matrix(0, ncol=Ntrap, nrow=N)
# 検出率を格納する箱
mu <- matrix(0, ncol=Ntrap, nrow=N)
# 生成する撮影枚数データを格納する箱
y <- array(0, dim=c(N, Ntrap, Npri, Nsec))
# 撮影枚数データの生成
for (i in 1:N) { # 個体
for (j in 1:Ntrap) { # トラップ
# 各個体の活動中心と各トラップの距離の2乗を計算
D2[i,j] <- (ac_x[i]-trap_x[j])^2 + (ac_y[i]-trap_y[j])^2
# ハザード半正規関数による、距離に応じた検出率の減衰
mu[i,j] <- 1-exp(-(lam0*exp(-D2[i,j]/(2*sigma*sigma))))
for (k in 1:Npri) {
for (m in 1:Nsec) {
y[i,j,k,m] <- rbinom(1, 1, mu[i,j]*Effort[j,Nsec*(k-1)+m]*z[i,k]*p)
}
}
}
}
# ある個体を捕獲したtrapは、同一occasionにおいて他の個体を捕獲できない
for (k in 1:Npri) {
for (m in 1:Nsec) {
for (j in 1:Ntrap) {
for (i in 2:N) {
if(sum(y[1:(i-1),j,k,m])>0) {y[i,j,k,m] <- 0}
}
}
}
}
# あるtrapで捕獲された個体は、同一occasionにおいて他のtrapで捕獲できない
for (k in 1:Npri) {
for (m in 1:Nsec) {
for (i in 1:N) {
for (j in 2:Ntrap) {
if(sum(y[i,1:(j-1),k,m])>0) {y[i,j,k,m] <- 0}
}
}
}
}
# 1枚も撮影されていない個体は検出できないので、データから削除
yobs <- y[apply(y,1,sum)>0,,,]
openCRによる推定
# 撮影データの配列を縦方向に展開する
temp <- as.data.frame(yobs) #行は一度でも観測された個体、列はトラップ×1次調査機会数x2次調査機会数
colnames(temp) <- paste(rep(Trap_id, Npri*Nsec), rep(rep(1:Npri, each=Ntrap), Nsec), rep(1:Nsec, each=Ntrap*Npri),sep="_")
tl <- reshape(temp, varying=list(1:ncol(temp)),
idvar="individual", ids=1:nrow(temp),
timevar="Trap_sc", times=colnames(temp),
direction="long"
)
tl$Session <- as.numeric(sapply(strsplit(tl$Trap_sc, "\\_"), "[", 2))
tl$trapID <- sapply(strsplit(tl$Trap_sc, "\\_"), "[", 1)
tl$Occasion <- as.numeric(sapply(strsplit(tl$Trap_sc, "\\_"), "[", 3))
tl <- tl[tl$C01_1_1==1, ]
colnames(tl)[grep("individual",colnames(tl))] <- "ID"
tl <- tl[, c("Session", "ID", "Occasion", "trapID")]
# 観測位置（と努力量）の情報
library(openCR)
trapXY <- data.frame(trap_id=Trap_id, x=trap_x, y=trap_y)
# これまでと異なり、detectorの種類はsingleにします
trapbase <- read.traps(data=trapXY, detector="single")
rownames(trapbase) <- Trap_id
chlist <- list()
for (i in 1:Npri) {
usage(trapbase) <- Effort[, (Nsec*(i-1)+1):(Nsec*(i-1)+Nsec)]
chlist[[i]] <- make.capthist(tl[tl$Session==i,], trapbase, fmt="trapID")
}
test <- MS.capthist(chlist)
# 推定範囲のポリゴンを作る
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))
)
)
)
msk <- make.mask(traps(test)[[1]], spacing = 1,
buffer = 3, type = "traprect", poly=estrange
)
# 引数intervalで、1次調査期間の間隔を指定します
secropen2 <- openCR.fit(test, type = 'JSSAsecrb',
model=list(lambda0~1, b~t, phi~1), mask = msk,
distribution="binomial", detectfn="HHN", method="Nelder-Mead",
interval=Intpri, details=list(control=list(maxit=10000))
)
推定結果
設定した値と推定値の関係は以下のとおり。今回の条件では、比較的設定値を推定できているようである。ただ、調査期間や乱数の種によっては、うまくいかないことも。
