こちらはQiitaのStanアドベントカレンダー2020 第19日目用の記事になります。
近くの山にカブトムシ調査に行った状況を想定します。
知りたいことは、
・どのくらいカブトムシがいるのか
・カブトムシの見つけやすさはどのくらいなのか
・見つけやすい時間帯はあるのか
です。
調査デザインとしては、山を適当なサイズのプロットに区切り、それぞれのプロットでカブトムシの在不在状況を記録します。
夜の山は危険なので、カブトムシ取りは朝〜夕方の間に行う予定で、事前調査もその時間帯で行います。
また、カブトムシの存在量を数量化することも可能ですが、ここでは簡単のため在・不在情報として記録することとします。
ここで、この調査において最も重要な考え方について、いくつか記述します。
1. カブトムシの検出は「プロットiにカブトムシが存在する確率psi(i)とカブトムシの検出率p(i)」の積から決定される
つまり、「カブトムシが検出された」という事実は、疑いもなくそこにカブトムシが存在することを意味するが、「カブトムシが検出されなかった」ことは「カブトムシが存在しなかった」状態と「存在したが見つけられなかった」状態の2つのパターンが考えられる、ということです。
この対応として、同じプロットで複数回の調査を繰り返します。一回の「不検出」からは「不在」と「在→不検出」の区別ができませんが、例えば調査を3回繰り返し「不検出・不検出・検出」となれば、「存在はしたが1回目と2回目は見つけられなかった」という事実が浮かび上がり、存在率と検出率の分離が可能になります。
2. 各プロット内での閉鎖個体群を想定する
これは1の考え方を可能にする仮定になります。閉鎖個体群とは、「個体の誕生・死亡・移出・移入」が発生しない、ということです。
現実にこの仮定が満たされることはないので、実際上は「閉鎖個体群として近似できるほどの期間の間に調査する」ことが必要となります。幸い、甥っ子とのカブトムシ取りと今回の事前調査は日中に行うため、夜行性のカブトムシは１日の昼間は閉鎖個体群として近似できるでしょう。
まずは、もっともシンプルな状況でシミュレーションしましょう（といってもずっとシミュレーションですが）。
これは調査時間は考慮せず、一定存在率・一定検出率を仮定します。
擬似データを生成します。
library(rstan)
set.seed(123)
R <- 100 # 調査サイト数
T <- 5 # 時間的調査反復数
psi <- 0.6 # 占有率（存在率）
p <- 0.4 # 検出率
z <- rbinom(R, 1, psi) # R地点の真の在不在情報
y <- array(NA, dim=c(R,T)) # 調査記録枠組み
for (i in 1:R){ # 調査記録（検出/不検出）格納
for (t in 1:T){
y[i,t] <- rbinom(1, 1, z[i] * p) # 検出成功率は z(在不在) * p(検出率)
}
}
調査から得られるデータは以下のようになります。
行がそれぞれのサイト、列が反復です（6行目まで表示しています）。
100プロットそれぞれについて5反復の調査を実施し、検出ならば1、不検出ならば0となります。
さて、存在率をpsi[i]、検出率をp[i,t]、プロットiの真の存在状態をz[i] (二値)、調査結果の観測値をy[i,t] (二値)としてモデルを整理すると、以下のようになります。
ただし、ここではpsiとpは一定なので、添字は必要ありません。
z[i] ~ Bernoulli(psi[i])
y[i,t] | z[i] ~ Bernoulli(z[i] * p[i,t])
Stanでは離散パラメータは扱えないので、以下のようなコードとなります。
data{
int R; // プロット数
int T; // 反復数
int<lower=0, upper=1> y[R,T]; // 観測データ
}
transformed data{
int sum_y[R]; // プロットiの合計検出回数
for (i in 1:R)
sum_y[i] = sum(y[i]);
}
parameters{
real<lower=0, upper=1> psi; // 存在率
real<lower=0, upper=1> p; // 検出率
}
model{
for (i in 1:R){
if (sum_y[i] > 0) // 一度でも検出された場合（存在）
target += bernoulli_lpmf(1 | psi)
+ bernoulli_lpmf(y[i,] | p);
else // 一度も検出されなかった（不在 or 存在するが不検出）
target += log_sum_exp(bernoulli_lpmf(0 | psi),
bernoulli_lpmf(1 | psi)
+ bernoulli_lpmf(y[i,] | p)
);
}
}
キックするコードは以下で、結果を図示します。
data <- list(R=nrow(y), T=ncol(y), y=y)
fit <- stan("beetle0.stan", data=data, seed=123)
真値はpsi=0.6, p=0.4ですので、まぁいい感じと思われます。
先の例はシンプルなシミュレーションなので、もう少し複雑にします。
具体的には、
・カブトムシの見つけやすい時間帯を推定したい
・複数人に調査を手伝ったもらったので、調査者ごとの検出率差も考慮したい
ということにします。
プロット数100、反復5、調査者4人で手分け、として生成された擬似データは以下のようになります。
なお、データ生成コードや図示コードが長くなってしまったので、ページの最後に載せます。
time_y1列とPerson_y1列が、y1列に対応する観測時刻と観察者、という見方になります。
timeについては、午前7時からの経過時間を分単位で記録しています（0 < time < 600）
また、観測時間及び4人の観察者による真の検出率pは以下のようになっています
12時前後が最も見つけづらい、という風になっています。
基本的には、先のシンプルな例をアレンジするだけです。
Stanコードは以下になります。
data{
int R; // プロット数
int T; // 反復数
matrix[R,T] Time; // 観測時刻
int<lower=1, upper=4> Person[R,T]; // 観測者ID(4人)
int<lower=0, upper=1> y[R,T]; // 観測データ
int T_new; // 時刻と検出率図示用
vector[T_new] Time_new; // 同上
}
transformed data{
int sum_y[R]; // 検出回数合計
matrix[R,T] Time2; // Timeの2乗
vector[T_new] Time2_new; // Time_newの2乗(図示用)
for (i in 1:R)
sum_y[i] = sum(y[i]);
Time2 = Time .* Time;
Time2_new = Time_new .* Time_new;
}
parameters{
real<lower=0, upper=1> psi; // 存在率
real beta0;
real beta1;
real beta2;
vector[4] eps; // 観測者ごと検出率誤差
real<lower=0> sigma; // その標準偏差
}
transformed parameters{
matrix[R,T] logit_p; // ロジット軸での検出率p
for(i in 1:R)
for(t in 1:T)
logit_p[i,t] = beta0 + beta1 * Time[i,t]
+ beta2 * Time2[i,t] + eps[Person[i,t]];
}
model{
sigma ~ normal(0, 1); // 事前分布
eps ~ normal(0, sigma); // 観測者ごと誤差
for (i in 1:R){
if (sum_y[i] > 0)
target += bernoulli_lpmf(1 | psi)
+ bernoulli_logit_lpmf(y[i,] | logit_p[i,]);
else
target += log_sum_exp(bernoulli_lpmf(0 | psi),
bernoulli_lpmf(1 | psi)
+ bernoulli_logit_lpmf(y[i,] | logit_p[i,])
);
}
}
generated quantities{
vector[T_new] p_pred; // 検出率と時間の効果
p_pred = inv_logit(beta0 + beta1 * Time_new + beta2 * Time2_new);
}
調査者4人には、事前にカブトムシの見つけ方講座を行っていたものとし、調査者ごとの検出率誤差に関わるsigmaにNormal(0,1)と幅の小さめの事前分布を適用しています。
見つけやすい、あるいは見つけにくい時刻を探すために、Timeの2乗を用いて回帰しています。
また、観測時刻と検出率の関係を知りたいので、generated quantitiesブロックで計算します。
推定結果は以下のようになります。
観察しやすい、あるいは観察しづらい時間については、Timeの2乗にかかる係数beta2がほぼ正の値を取ることから、ある時間で最も観察しづらい時間が存在することが分かります。これも図示してみましょう。
真値は
存在率psi = 0.6, beta0 = -1, beta1 = -0.3, beta2 = 0.5, eps[1] = 0.6, eps[2] = 0.3, eps[3] = -0.3, eps[4] = -0.6
と設定しており、悪くはなさそうです。
時間と検出率の関係も良さそうです。
ちゃんと12時前後でもっとも見つけづらくなっています。
時刻と検出率の図示については、下に凸であることを強調するため、MCMCサンプル300個を図示する方法を取りました。
このモデルはいろいろな拡張が可能です。
今回は閉鎖個体群を想定していますが、現実的には開放個体群を想定したいところでしょう。
その場合、開放個体群を想定する調査期間と、閉鎖個体群を想定する調査期間を組み合わせるロバストデザインが用いられるようです。
これにより、個体の死亡率や移入率なども推定できる場合があります。
今回の記事に当たっては以下の書籍を参考にさせていただきました。
- 作者:K´ery,Marc,Schaub,Michael
- 発売日: 2016/03/25
- メディア: 単行本
今回紹介した内容とその拡張も含めて、生態学（特に個体数や確率論的個体群動態）での様々な事例を、主に状態空間モデルから解析した内容になっています。
Winbugsを用いた解析ですが、全ての例についてStanでの実装がGithub上に挙げられています。
何より、日本語訳が素晴らしいです。考え方が分かりやすく、的確な統計学的用語で解説してくれていると思いました。
不満があるとすれば、Stanでは離散パラメータが簡単に扱えないことから、Winbugsだとこんなに簡単なコードなのに、この複雑なStanコードは一体...となることが度々あったことでしょうか。
まぁそれは仕方ないですね...一度理解できればなんとかなります。
基本的にはStanの方が高速なので、その引き換えと思う他ありません。
では最後に、擬似データ生成、Stanのキック、図示に関わるRコードをバーっと載せときます。
もっと簡単にできそうな気がしますが、技術力不足ゆえ長くなっている気がします。
library(rstan)
library(tidyverse)
set.seed(123)
#### 擬似データ生成 ####
R <- 100 # 調査サイト数
T <- 5 # 時間的調査反復数
psi <- 0.6 # 占有率（存在率）
Person_logit_p <- c(0.6, 0.3, -0.3, -0.6) # 調査者による検出率誤差（ロジット軸）
time <- array(NA, dim=c(R, T)) # 観測時間
Time <- array(NA, dim=c(R, T)) # 標準化された観測時間
Person <- array(NA, dim=c(R,T)) # 調査者ID
p <- array(NA, dim=c(R,T)) # 検出率
y <- array(NA, dim=c(R,T)) # 観測値
for(i in 1:R){
for(t in 1:T){
time[i,t] <- sample(0:600, 1) # 午前7時からの経過時間（分）
Person[i,t] <- sample(1:4, 1) # 調査者ID (1~4)
}
}
for(i in 1:R){
for(t in 1:T){
Time[i,t] <- (time[i,t] - mean(time)) / sd(time) # 観測時間を標準化
p[i,t] <- plogis( -1 - 0.3 * Time[i,t] + 0.5 * Time[i,t] * Time[i,t]
+ Person_logit_p[Person[i,t]]) # 観測時間と観察者による検出率
}
}
z <- rbinom(R, 1, psi) # 真の在不在情報
for(i in 1:R){ # 観測データ生成
for(t in 1:T){
y[i,t] <- rbinom(1, 1, z[i] * p[i,t]) # 観測成功率 = 在不在(1/0) * 検出率p
}
}
##### 検出率と時刻、調査者の関係図示 ####
P <- array(NA, dim=c(3*R, 2))
for(i in 1:R){
P[i,1] <- time[i,1]
P[i,2] <- p[i,1]
P[i+R, 1] <- time[i,2]
P[i+R, 2] <- p[i, 2]
P[i+2*R, 1] <- time[i,3]
P[i+2*R, 2] <- p[i,3]
}
P <- as.data.frame(P)
colnames(P) <- c("time", "P")
ggplot(P, aes(x=time, y=P))+
geom_point()+
ylim(0,1)
#### 観測データをまとめる ####
d <- cbind(y, time, Person) # 調査で得られるデータ
colnames(d) <- c("y1", "y2", "y3", "y4", "y5",
"time_y1", "time_y2", "time_y3", "time_y4", "time_y5",
"Person_y1", "Person_y2", "Person_y3", "Person_y4", "Person_y5")
#### 推定 ####
time_new <- seq(0, 600, 1) # generated quantitiesブロック用
Time_new <- (time_new - mean(time)) / sd(time) # 標準化
T_new <- length(Time_new) # generated quantitiesブロック用
data <- list(R=nrow(y), T=ncol(y), K=4, Time=Time, Person=Person, y=y,
Time_new=Time_new, T_new=T_new)
fit <- stan("beetle.stan", data=data, seed=123,
control=list(adapt_delta=0.995))
fit
#### 検出率と時刻の関係を図示 ####
ms <- rstan::extract(fit)
p_pred <- ms$p_pred
# 時刻ごと中央値を算出
p50 <- apply(p_pred, 2, quantile, probs=0.5)
P50 <- data.frame(Time=time_new, p50)
# MCMCサンプル300個を取り出す
p_p <- t(p_pred)
P_p <- as.data.frame(p_p[,1:300])
P_p$Time <- time_new
P_t <- pivot_longer(P_p, cols=-Time, names_to="MCMC", names_prefix="V", values_to="p")
# 図示
ggplot()+
theme_classic()+
geom_line(data=P50, aes(x=Time, y=p50), colour="blue", size=1.6)+
geom_path(data=P_t, aes(x=Time, y=p, group=MCMC), alpha=0.3, size=0.2)+
ylab("Detection rate")
