- 被度（調査面積に対する対象植物の被覆度）は、もともとは0〜1の範囲の連続値。
- 実際には、{+, 1, 2, 3, 4, 5}などといった階級データとして記録されることが多い。
- また、測定値は目視で決められることが多い。
植生被度階級データの
モデリング
伊東宏樹
2020-03-27
はじめに
今回のコード
https://github.com/ito4303/jfs131 で公開しています。
被度階級データ
被度階級データの問題点
- 順序尺度データなので、数値的な取り扱いが簡単ではない。
- 目視による測定では、誤差が大きいと考えられる。
被度のモデリング
ベータ分布
被度（連続値）は0〜1の値をとるので、ベータ分布にあてはめるのは自然な発想。
\[ y \sim \mathrm{Beta}(\alpha, \beta) \]
平均\(\mu\)を使ったパラメータ化
\[ y \sim \mathrm{Beta}\left(\frac{\mu}{\delta}-\mu,\frac{(1-\mu)(1-\delta)}{\delta}\right) \]
\(\delta\)はもともと、pin-point法による被度測定における方形区内の分布相関 (Damgaard 2012)
しかし、被度測定における不確実性とも解釈可能。
ベータ分布(\(\mu\)=0.5のとき)
ベータ分布(\(\mu\)=0.05のとき)
被度階級
今回は以下のように定義する。
1: 0–0.01 (0を含む), 2: 0.01–0.1, 3: 0.1–0.25, 4: 0.25–0.5, 5: 0.5–0.75, 6: 0.75–1
\(\delta\)を変化させたとき、平均被度に対して、各被度階級が選ばれる確率を図示する。
\(\delta\)=0.001
\(\delta\)=0.01
\(\delta\)=0.05
\(\delta\)=0.1
\(\delta\)=0.2
\(\delta\)=0.4
被度階級のモデリング
模擬データ
- 被度 = 0.6
- \(\delta\) = 0.05 とする。
- 10回測定する(\(N = 10\))。
生成されたデータ
## [1] 5 5 5 5 4 5 5 5 5 5Stan コード
被度階級の確率分布を関数として定義
functions {
real coverclass_lpmf(int Y, vector CP, real a, real b) {
int n_cls;
real gamma;
n_cls = num_elements(CP) + 1;
if (Y <= 1) { // 0 or 1
gamma = inc_beta(a, b, CP[1]);
} else if(Y >= 2 && Y < n_cls) {
gamma = inc_beta(a, b, CP[Y])
- inc_beta(a, b, CP[Y - 1]);
} else {
gamma = 1 - inc_beta(a, b, CP[n_cls - 1]);
}
return bernoulli_lpmf(1 | gamma);
}
}モデルの定義
model {
{
real a = mu / delta - mu;
real b = (1 - mu) * (1 - delta) / delta;
for (n in 1:N)
Y[n] ~ coverclass(CP, a, b);
}
}モデルへのあてはめ
事後分布の要約
## Inference for Stan model: cover-202004211608-1-e492ef.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## mu 0.59 0 0.04 0.51 0.56 0.59 0.61 0.67 1619 1
## delta 0.04 0 0.03 0.00 0.02 0.03 0.05 0.11 1819 1
##
## Samples were drawn using NUTS(diag_e) at 火 4 21 16時08分36秒 2020.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).事後予測チェック
とくに矛盾はない
共変量のあるモデル
模擬データ
- 共変量: X (0〜1の連続値)
次式で被度が生成される。 \[ \mathrm{logit}(p) = -5 + 5 X + \epsilon \\ \epsilon \sim \mathrm{Normal}(0, \sigma) \]
- \(N = 50\) か所で測定
- 1か所につき、2回くりかえし測定(\(R = 2\))
グラフ
Stan コード
平均被度muをロジットリンクで線形予測子と結びつける
transformed parameters {
vector<lower = 0, upper = 1>[N] p; // proportion of cover
mu = inv_logit(beta[1] + beta[2] * X + sigma * err);
}モデルの定義
model {
// Observation
for (n in 1:N) {
real a = mu[n] / delta - mu[n];
real b = (1 - mu[n]) * (1 - delta) / delta;
for (r in 1:R)
Y[n, r] ~ coverclass(CP, a, b);
}
// System
err ~ std_normal();
sigma ~ normal(0, 5);
}結果の要約
## Inference for Stan model: cover2-202004211614-1-f4091b.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## delta 0.13 0.00 0.03 0.08 0.10 0.12 0.14 0.19 1232 1
## beta[1] -4.76 0.01 0.45 -5.66 -5.05 -4.75 -4.45 -3.93 1665 1
## beta[2] 4.73 0.01 0.66 3.45 4.28 4.72 5.17 6.08 2067 1
## sigma 0.59 0.01 0.23 0.10 0.45 0.61 0.75 1.00 680 1
##
## Samples were drawn using NUTS(diag_e) at 火 4 21 16時15分32秒 2020.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).事後予測チェック
1回目の測定について表示
ゼロ過剰モデル
模擬データ
- このモデルでの被度
- 0: 0, 1: 0–0.01 (0を含まない), 2: 0.01–0.1, 3: 0.1–0.25, 4: 0.25–0.5, 5: 0.5–0.75, 6: 0.75–1
- 分布している確率: \(\omega\)
- 被度1のときの発見確率: \(\psi\)
- 被度2以上のときは必ず発見できるとする(\(\psi\)=1)。
- \(N = 50\) か所で測定
- 1か所につき、2回くりかえし測定(\(R = 2\))
生成されたデータ
Stan コード
モデルの定義
model {
// Observation model
for (n in 1:N) {
real a = mu[n] / delta - mu[n];
real b = (1 - mu[n]) * (1 - delta) / delta;
if (sum(Y[n]) == 0) { // Y[n]==0 for all n
real lp[2];
lp[1] = bernoulli_lpmf(0 | omega);
lp[2] = bernoulli_lpmf(1 | omega)
+ coverclass_lpmf(1 | CP, a, b) * R
+ bernoulli_lpmf(0 | psi) * R;
target += log_sum_exp(lp);
} else {つづく
for (r in 1:R) {
if (Y[n, r] == 0) {
target += bernoulli_lpmf(1 | omega)
+ coverclass_lpmf(1 | CP, a, b)
+ bernoulli_lpmf(0 | psi);
} else if (Y[n, r] == 1) {
target += bernoulli_lpmf(1 | omega)
+ coverclass_lpmf(1 | CP, a, b)
+ bernoulli_lpmf(1 | psi);
} else {
target += bernoulli_lpmf(1 | omega)
+ coverclass_lpmf(Y[n, r] | CP, a, b);
}
}
}
}結果の要約
## Inference for Stan model: zicover-202004211616-1-5659d3.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## delta 0.14 0.00 0.03 0.08 0.11 0.13 0.16 0.21 2659 1
## omega 0.82 0.00 0.05 0.72 0.79 0.82 0.85 0.90 5357 1
## psi 0.75 0.00 0.12 0.49 0.68 0.77 0.84 0.93 4893 1
## beta[1] -4.74 0.01 0.54 -5.91 -5.09 -4.71 -4.37 -3.79 2265 1
## beta[2] 4.42 0.01 0.78 3.02 3.90 4.37 4.91 6.08 2797 1
## sigma 0.47 0.01 0.25 0.03 0.29 0.46 0.63 0.98 1071 1
##
## Samples were drawn using NUTS(diag_e) at 火 4 21 16時17分52秒 2020.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).事後予測チェック
空間データ
模擬データ
10×10=100個の方形区で、被度を生成する
↓
被度階級に変換
生成された被度（割合）
生成された被度階級
空間自己相関
隣とは近い値をとる
Stan による Conditional AutoRegressive (CAR) モデル
Exact sparse CAR models in Stan by Max Joseph
https://mc-stan.org/users/documentation/case-studies/mbjoseph-CARStan.html
Stan コード
phiの逆ロジットが平均被度mu
transformed parameters {
vector[N] mu = inv_logit(phi);
}sparse_iarで、空間自己相関の事前分布
model {
// Spatial random effects
phi ~ sparse_iar(tau, W_sparse, D_sparse, lambda, N, W_n);
// Observation model
for (n in 1:N) {
real a = mu[n] / delta - mu[n];
real b = (1 - mu[n]) * (1 - delta) / delta;
Y[n] ~ coverclass(CP, a, b);
}
// Priors
tau ~ gamma(2, 2);
}結果の要約
## Inference for Stan model: carcover-202004211619-1-ba17b0.
## 4 chains, each with iter=8000; warmup=4000; thin=1;
## post-warmup draws per chain=4000, total post-warmup draws=16000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## delta 0.02 0.00 0.01 0.01 0.01 0.02 0.02 0.04 4472 1
## tau 1.69 0.01 0.58 0.84 1.27 1.60 2.02 3.10 2703 1
##
## Samples were drawn using NUTS(diag_e) at 火 4 21 16時22分41秒 2020.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).推定された被度（割合）
真値と推定値との比較
まとめ
被度階級データ
- 扱いづらい被度階級データ
- 階層モデリング
- データが観測される過程
- データの背後にあるシステム
- 被度階級データも、適切にモデリングすることで うまく扱える