モデルの確認
本検定のためのモデルについてはこちらの姉妹記事で詳細を説明していますので、ここでは要点だけ記載しておくことにします。
$n_1$、$n_2$はそれぞれ交差点1、交差点2における車両毎の危険挙動に関するデータ（0=「危険挙動無し」 OR 1「危険挙動あり」）の標本数、$s_1$、$s_2$はそれぞれ交差点1、交差点2において危険挙動をした車両の台数です。交差点ごとに、独立の二項分布に従い危険挙動をした車両が観測されたと仮定し、発生確率を$\theta_1$、$\theta_2$とおきモデル化します。
$\theta_1$、$\theta_2$には、無情報事前分布$\mathrm{Uniform}(0,1)$を設定することとします。
このとき、我々の興味の対象となる確率変数は$\delta=\theta_1 - \theta_2$です。$\theta_1$、$\theta_2$の大小関係についての事前知識がある場合と無い場合に分け、2つの帰無仮説・対立仮説の組を設定します。交差点の比較の例では、ふたつの交差点の危険挙動発生率について大小関係が不明であると考えることもできますが、交差点2は交差点1と同じ形状であるが信号機が設置されていることから、危険挙動発生率は交差点1より少ない、すなわち$\theta_1 > \theta_2$ではないか、と考えることもできます。
これらは、従来の仮説検定における片側検定・両側検定に対応します。
Unrestricted Model
では、$\theta_1$と$\theta_2$の大小関係が明らかでないときの仮説から検定してみます。Rstanを用いて$(1)～(4)$式をモデル化し、$\delta$の周辺事後分布をMCMCによって近似推定し、Savage-Dickeyによる$(5)$式を用いて$BF_{01}$を求める方針です。またこちらの記事で導出した解析的な結果との整合性も確認してみます。
まずはサンプルデータを準備します。
# sample data
#x1:Group1
x1 <- c(rbinom(n=40, size=1, prob=0.4),rbinom(n=60, size=1, prob=0.5))
#x2:Group2
x2 <- c(rbinom(n=40, size=1, prob=0.2),rbinom(n=60, size=1, prob=0.5))
library(ggplot2)
library(gridExtra)
p1 <- ggplot(data=data.frame(X=x1)) + theme_light(base_size=11) + geom_bar(aes(x=X), fill="#004C98") +
scale_x_continuous(breaks=c(0,1), labels = c("-","危険挙動")) + xlab("Group1")
p2 <- ggplot(data=data.frame(X=x2)) + theme_light(base_size=11) + geom_bar(aes(x=X), fill="#004C98") +
scale_x_continuous(breaks=c(0,1), labels = c("-","危険挙動")) + xlab("Group2")
p <- grid.arrange(p1, p2, ncol=2)
このデータに対し、通常の比率の差の検定を関数prop.test()で行います。原理は昔紹介した分割表の一様性の検定と同じものです。
prop.test(x=c(sum(x1),sum(x2)), n=c(length(x1),length(x2)), correct=F)
## 2-sample test for equality of proportions without continuity correction
##
## data: c(sum(x1), sum(x2)) out of c(length(x1), length(x2))
## X-squared = 1.6807, df = 1, p-value = 0.1948
## alternative hypothesis: two.sided
## 95 percent confidence interval:
## -0.04549292 0.22549292
## sample estimates:
## 0.45 0.36
## prop 1 prop 2
２グループの比率が同じであるという帰無仮説を設定した検定に対し、p値$0.1948$となり、比率に差があるとは言えない結果となりました。
次にこちらの記事の$(7)$式に基づいて解析的なベイズファクターの計算による検定を行います。
data <- list(N1=length(x1), N2=length(x2), S1=sum(x1), S2=sum(x2))
# calculate BF_01 by closed form
# data[[1]]：n1
# data[[2]]：n2
# data[[3]]：s1
# data[[4]]：s2
cat(choose(data[[1]],data[[3]]) * choose(data[[2]],data[[4]]) / choose(data[[1]]+data[[2]], data[[3]]+data[[4]])*(data[[1]]+1)*(data[[2]]+1)/(data[[1]]+data[[2]]+1))
## 2.526656
対立仮説$H_1:\delta\neq0$よりも帰無仮説$H_0:\delta=0$の方がデータの説明力（周辺尤度）が約$2.53$倍良く、どちらかというと帰無仮説を支持するようです。しかしKass and Raftery(1995)の基準ではこの差は対立仮説に対する反証があまり強いとは言えない程度であると評価されました。
以上、検定結果をひととおり見た上で、次はMCMCにより$\delta$の事後分布を近似推定し、ベイズファクターを推定します。Stanによる実装例を以下に示します。
//non-restricted aalysis
//model1.stan
data {
//the number fo sample of Group1
int<lower=1> N1;
//the number fo sample of Group2
int<lower=1> N2;
//the number fo sample which occurs with a probability of theta1 of Group1
int<lower=0> S1;
//the number fo sample which occurs with a probability of theta2 of Group2
int<lower=0> S2;
}
parameters {
//the parameter of Binomial distribution of Group1
real<lower=0, upper =1> theta1;
//the parameter of Binomial distribution of Group2
real<lower=0, upper=1> theta2;
}
model {
//S1 ~ binomial(N1, theta1)
target += binomial_lpmf(S1 | N1, theta1);
//S2 ~ binomial(N2, theta2)
target += binomial_lpmf(S2 | N2, theta2);
//uninformative prior
target += beta_lpdf(theta1 | 1,1);
target += beta_lpdf(theta2 | 1,1);
}
generated quantities{
//the difference parameter between theta1 and theta2
real delta;
delta = theta1 - theta2;
}
iterationを$11000$、warmupを$1000$、chain数を$4$としMCMCを走らせ、結果を確認します。
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)
fit <- stan(file="model/model1.stan", data=data, iter=11000, warmup=1000, chain=4)
fit
## Inference for Stan model: model1.
## 4 chains, each with iter=11000; warmup=1000; thin=1;
## post-warmup draws per chain=10000, total post-warmup draws=40000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## theta1 0.45 0.00 0.05 0.36 0.42 0.45 0.48 0.55 36883 1
## theta2 0.36 0.00 0.05 0.27 0.33 0.36 0.39 0.46 36917 1
## delta 0.09 0.00 0.07 -0.05 0.04 0.09 0.13 0.22 36034 1
## lp__ -8.89 0.01 1.02 -11.60 -9.27 -8.58 -8.17 -7.90 17925 1
このシミュレーション結果を使って$\delta$の事後分布を近似推定します。推定手法は参考文献が推奨する対数スプライン推定を用いました。通常、事後分布の推定にはカーネル密度推定が用いられることが多いと思うのですが、対数スプライン推定はカーネル密度推定と比較し推定値が負の値を取らないという利点があり、この場合推定精度が比較的良いとされているようです。
library(ggplot2)
library(ggrepel)
library(logspline)
library(tidyverse)
# non_restricted analysis
d <- as.matrix(fit)[,"delta"]
fit_logspline <- logspline(d)
plot_title = ggtitle(label="equalty of proportions(unrestricted analysis)")
p <- ggplot() + theme_light(base_size=11) +
geom_path(data=data.frame(X=x <- seq(-1,1,0.01), Y=dlogspline(x,fit_logspline)), aes(x=X,y=Y), color="blue") +
geom_line(data=data.frame(x=c(-1,0,1),y=c(0,1,0)),aes(x=x,y=y)) +
geom_point(data=data.frame(x=0,y=dlogspline(0, fit_logspline)),aes(x=x,y=y), color="blue") +
geom_point(data=data.frame(x=0,y=1), aes(x=x, y=y)) +
geom_text_repel(data=data.frame(x=c(0,0),y=c(1,dlogspline(0,fit_logspline)),label=c(1,round(dlogspline(0, fit_logspline),2))),
aes(x=x,y=y,label=label)) +
annotate("text",x=-Inf,y=Inf,label=sprintf("BF01=%.2f",dlogspline(0,fit_logspline)),hjust=-.2,vjust=2) + plot_title
p
結果は、$BF_{01} \simeq 2.50$と、解析的な解とそれほど大差無い結果となりました。やはりこのデータは確信をもって比率に差があるとは言えないようですな。
同じ危険挙動発生率でも標本数が倍になると、確信をもって差があると言える程度も変わってくるはずです。ちょっと確認してみましょう。
# data with same probability but twice sample size
#new_x1:Group1
new_x1 <- rep(x1,10)
#new_x2:Group2
new_x2 <- rep(x2,10)
～（以下略）～
cat(1/dlogspline(0,fit_logspline))
## 147.1852
$BF_{10} \eqsim 147.19$と、期待通り対立仮説を強く支持する結果となりました。
Order-restricted Model
次に、$\theta_1 \geq \theta_2$が明らかな場合のOrder-restricted Modelについて考えます。モデル構成はUnrestricted Modelと同じで、あらかじめ分かっている事前情報が違うだけなので、新たに別のモデルでMCMCする必要はありません。方針としては、先ほど得たMCMCサンプルより推定した事後分布のうち、$\theta_1 \geq \theta_2$を満たす範囲のみを抽出し、事後分布とします。
# restricted analysis
fit_logspline <- logspline(d)
dens <- function(x){
dlogspline(x, fit_logspline)/(1-plogspline(0, fit_logspline))}
plot_title = ggtitle(label="equalty of proportions(order-restricted analysis)")
p <- ggplot() + theme_light(base_size=11) +
geom_path(data=data.frame(X=x <- seq(0,1,0.001), Y=dens(x)), aes(x=X,y=Y), color="blue") +
geom_line(data=data.frame(x=c(0,1),y=c(2,0)),aes(x=x,y=y)) +
geom_point(data=data.frame(x=0,y=dens(0)),aes(x=x,y=y), color="blue") +
geom_point(data=data.frame(x=0,y=2), aes(x=x, y=y)) +
geom_text_repel(data=data.frame(x=c(0,0),y=c(2,dens(0)),label=c(2,round(dens(0),2))),
aes(x=x,y=y,label=label)) +
annotate("text",x=Inf,y=Inf,label=sprintf("BF01=%.2f",dens(0)/2),hjust=1.2,vjust=2) +
xlab("δ") + ylab("density") + plot_title
p
cat(dens(0)/2)
## 1.385902
検定の結果は、Unrestricted Modelが$BF_{01} \eqsim 2.50$に対し、Order-restricred Modelが$BF_{02} \eqsim 1.39$となり、帰無仮説$H_0$を支持する傾向は同じものの、やや確信度が小さくなる結果となりました。直感的には、グループ2よりグループ1の事故発生率が大きいことから、$\theta_1 \geq \theta_2$としたOrder-Restricted Modelに対するデータの当てはまりが良く、$BF_{01} \geq BF_{02}$となりそうだと考えられ、実際には事前分布の$\delta=0$での値が$2$倍になったことにより、直感通りの結果が得られたと言えます。
またベイズファクターは相対比較ですから、$BF_{01}$と$BF_{02}$が求まれば$BF_{12}$も一意に定まります。
$$ BF_{12} =\cfrac{p(D|H_1)}{p(D|H_2)} = \cfrac{p(D|H_1)}{p(D|H_0)}\cfrac{p(D|H_0)}{p(D|H_2)} = BF_{10}BF_{02} \eqsim 0.554 $$
よって、$H_2$の方が$H_1$よりも$\cfrac{1}{0.554} \eqsim 1.80$倍データをうまく説明している、と判断することができます。
ここまで、対応の無い比率のさの検定をベイズファクターを使って実践してみました。従来通りの結果と比べて理論が簡潔で私は好みのベイズファクター、「分かりやすい」統計のために今後重要になってくるのではないでしょうか。
一方、本例のような危険挙動データの例の場合、本来、危険挙動を観測した車両の全体に占める割合は非常に小さく、通常1割にも達しません（私の経験）。 このようなデータに対し、二項分布のパラメータ$\theta$に無情報事前分布を設定すると、周辺尤度の性質から、$BF_{01}$は必ず帰無仮説$H_0:\delta=0$を強く支持する値となるでしょう。事前分布をどうするかについては、今後の研究が待たれるところです（？）。
おわりに
今回は、ベイズファクターを使って比率の差を検定する方法を実践しました。解析的に求める方法と、MCMCを使って近似的に求める方法の2通りを実践しましたが、精度の面からいうともちろん解析的に求めるのが妥当です。しかし、後半で実践したような、Order-restricted Modelで検定したいときなどは、解析解をまた別に求めてやらなばなりません。こちらの記事で導出した解析解はあくまで従来通りの検定における帰無仮説と対立仮説の（方向性をもたない）比較をしたいときの解であり、この式だけでは他の仮説の比較に応用できないのです。それならOrder-restricted Modelの解析解を求めてやれば良いのですが、目的に応じてその都度解を求めるのは大変です。MCMC法はそういった面倒な計算をすっ飛ばして近似的にベイズファクターやら事後分布やらを求めることができる、という点が魅力的なんですね（といっても今回は片側検定に対応するモデルしかMCMC法で追加検討していませんが）。
