続いていよいよ本題、開放個体群の個体数の推定。なぜ開放個体群から始めなかったのかというと、考えることが多いから。閉鎖個体群と異なり、開放個体群では調査期間内での個体数の変動をモデルに含める必要がある。
個体数の変動は、以下の要因によって生じる。
- 死亡
- 加入（出生）
- 移出入
開放個体群の個体数を推定するモデルは、これらをどのようにモデル化するのかによって複数の種類がある。
データの生成
開放個体群になると、データ生成だけでも長いです（汗
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
## カメラの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*Npri, 20), nrow=Ncam, ncol=Npri)
## 生存率
phi <- rep(0.8, Npri)
## 加入率
beta <- rep(1, Npri)
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)]))
}
# 調査機会ごとの個体の在不在
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[t]*z[, t-1] + nu[t]*q[,1]
} else { mu2[, t-1] <- phi[t]*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 <- 0.1
# カメラからの距離に応じた検出率の減衰関数である
# 半正規関数の分散パラメータ
sigma <- 1.5
# 各個体の活動中心と各カメラの距離を格納する箱
D2 <- matrix(0, ncol=Ncam, nrow=N)
# 検出率を格納する箱
lambda <- matrix(0, ncol=Ncam, nrow=N)
# 生成する撮影枚数データを格納する箱
y <- array(dim=c(N, Ncam, Npri))
# 撮影枚数データの生成
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:Npri) { # 調査機会
y[i,j,k] <- rpois(1,lambda[i,j]*Effort[j,k]*z[i,k])
}
}
}
# 1枚も撮影されていない個体は検出できないので、データから削除
yobs <- y[apply(y,1,sum)>0,,]
openCRによる個体数推定
捕獲履歴のデータ
開放個体群では、標本抽出の時間スケールに留意する必要がある。開放個体群、つまり個体数が死亡や加入によって変動しうるくらい長い間隔での標本抽出を1次標本抽出（Session）と呼ぶ。一方、個体数が変動しないぐらい短い期間での標本抽出を2次標本抽出（Occasion）と呼ぶ。2次標本抽出を含む標本抽出の設計は、Robust designと呼ばれる。
閉鎖個体群の場合、Sessionは1つしかなく、調査回数はそのままOccasionの回数となる。開放個体群の場合、Sessionが複数ある。そのため、捕獲履歴の用意の仕方が、閉鎖個体群より複雑になる。
# 撮影データの配列を縦方向に展開する
temp <- as.data.frame(yobs) #行は一度でも観測された個体、列はカメラ×調査機会数
colnames(temp) <- paste(rep(Cam_id, Npri), rep(1:Npri, 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 <- as.numeric(sapply(strsplit(tl$Cam_so, "\\_"), "[", 2))
tl$trapID <- sapply(strsplit(tl$Cam_so, "\\_"), "[", 1)
tl$Occasion <- 1
tl <- tl[rep(1:nrow(tl), tl$C01_1), ]
colnames(tl)[grep("individual",colnames(tl))] <- "ID"
tl <- tl[, c("Session", "ID", "Occasion", "trapID")]
# 観測位置（と努力量）の情報
library(openCR)
trapXY <- data.frame(cam_id=Cam_id, x=cam_x, y=cam_y)
trapbase <- read.traps(data=trapXY, detector="count")
rownames(trapbase) <- Cam_id
chlist <- list()
# Sessionごとにcapthistオブジェクトを作成し、
# 調査努力量もSessionごとに与える。
for (i in 1:Npri) {
usage(trapbase) <- NULL
usage(trapbase) <- as.matrix(Effort[, i], ncol=1)
chlist[[i]] <- make.capthist(tl[tl$Session==i, ], trapbase, fmt="trapID")
}
# 複数のcapthistオブジェクトを、MS.capthist()で1つのcapthistオブジェクトにする
test <- MS.capthist(chlist)
推定範囲の設定
ここは閉鎖個体群とほぼ同様。ただし、調査範囲のmaskについては複数あるので、ここでは1番目のtrapオブジェクトを使う。
# 推定範囲のポリゴンを作る
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
)
推定
使う関数は同じだが、開放個体群の個体数を推定するモデルは複数ある。ここでは、加入確率を扱うモデルとする。
secropen1 <- openCR.fit(test, type = 'JSSAsecrb',
model=list(lambda0~1, b~t, phi~1), mask = msk,
distribution="poisson", detectfn="HHN", method="Nelder-Mead"
)
結果の確認
summary()関数を使う。また、開放個体群の場合、Sessionごとの個体数をderived()関数で取り出すことができる（例えば、derived(secropen1)と打ち込んでみよう）。おおむね、設定した値を推定できているようである。
