p値を計算したくなる検定の数々を試しにStanによるベイジアンモデリングで代替してみた

この記事は、やたらはてブを稼いでしまった前回の記事の続きです。

ASAのプレスリリース及び声明の中には、確かに「p値に依拠しない新たなアプローチの例」として予測値を重視するアプローチ*5、ベイジアンモデリング、決定理論的アプローチ*6およびfalse discovery rate*7といったものを用いるべき、という趣旨のコメントが入っています。とは言え、重回帰分析とか機械学習のような多変量モデリング（なおかつサンプルサイズも大きい）を伴うテーマならともかく、統計学的仮説検定のようなサンプルサイズも小さい（データも少ない）シチュエーションでどうやるんだよ的な疑問を持つ人も多いのではないかと。

そんなわけで、実際にそれっぽい各種検定の数々をStanによるベイジアンモデリングで代替してみたので、この記事ではその結果をつらつら紹介してみようと思います。テーマは前々回のこちらの記事の1節で取り上げた2種。

とりあえず今回はt検定とカイ二乗検定だけ取り上げてみます。ANOVAはそもそも線形回帰モデルで代替可能なので、ここでは割愛しますよということで。

（対応のある）t検定をベイジアン モデリングでやってみる

@TJO_datasci https://t.co/JTkxUd3J4q
の記事における最初のStanコードは、t検定ではなく対応のあるt検定の形になっているように思うのですが如何でしょう？
— Dr. おやすみくん (@SuperOyasumi) November 27, 2016
というご指摘の通り、下記のStanコードだとpaired t-testになってしまいます。これをunpairedに改める方法はこの節の最後に追記しました。

そもそもt検定ってどうやってモデルで表現するんだっけ？と思ったんですが、よくよく考えたら「平均値の差」のモデリングなのでこれで十分なんですね。

$x_2 \sim x_1 + m + N(0, \sigma)$

要は2つのデータの間の差 $m$ の事後分布が0から十分に離れているというサンプリング結果になれば、 $x_1$ と $x_2$ とは十分に差があると推論できるわけです。これをStanコードで表現するとこうなります。

data {
	int<lower=0> N;
	real<lower=0> x1[N];
	real<lower=0> x2[N];
}

parameters {
	real s;
	real m;
}

model {
	for (i in 1:N)
		x2[i]~normal(x1[i]+m,s);
}

このStanコードをttest.stanという名前で保存した上で、以下のRコードでkickします。

> d<-read.csv('https://raw.githubusercontent.com/ozt-ca/tjo.hatenablog.samples/master/r_samples/public_lib/DM_sampledata/ch3_2_2.txt',header=T,sep=' ')
> dat1<-list(N=nrow(d),x1=d$DB1,x2=d$DB2)
> fit1<-stan(file='ttest.stan',data=dat1,iter=1000,chains=4)

SAMPLING FOR MODEL 'ttest' NOW (CHAIN 1).

Chain 1, Iteration:   1 / 1000 [  0%]  (Warmup)
# ... #
Chain 1, Iteration: 1000 / 1000 [100%]  (Sampling)# 
#  Elapsed Time: 0.008053 seconds (Warm-up)
#                0.007259 seconds (Sampling)
#                0.015312 seconds (Total)
# 
# ... omitted ... #

SAMPLING FOR MODEL 'ttest' NOW (CHAIN 4).

Chain 4, Iteration:   1 / 1000 [  0%]  (Warmup)
# ... #
Chain 4, Iteration: 1000 / 1000 [100%]  (Sampling)# 
#  Elapsed Time: 0.007944 seconds (Warm-up)
#                0.008107 seconds (Sampling)
#                0.016051 seconds (Total)
# 
# ... omitted ... #

> fit1
Inference for Stan model: ttest.
4 chains, each with iter=1000; warmup=500; thin=1; 
post-warmup draws per chain=500, total post-warmup draws=2000.

      mean se_mean   sd  2.5%   25%   50%   75% 97.5% n_eff Rhat
s     0.73    0.01 0.15  0.50  0.62  0.70  0.81  1.08   577 1.01
m     0.72    0.01 0.19  0.35  0.59  0.72  0.84  1.09   824 1.00 # ←ココ！
lp__ -1.78    0.05 1.08 -4.57 -2.29 -1.43 -0.97 -0.67   509 1.00

Samples were drawn using NUTS(diag_e) at Wed Mar  9 22:06:55 2016.
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).
> fit1.coda<-mcmc.list(lapply(1:ncol(fit1),function(x) mcmc(as.array(fit1)[,x,])))
> plot(fit1.coda)

f:id:TJO:20160309221455p:plain

ということで、一応 $m$ の事後分布の2.5%点が0より大きいこと、そして綺麗な単峰性分布でサンプリングがきちんと収束しているとみなせることから、2つのデータの平均値の差は0より十分に大きいと言って良さそうです。

追記 (Mar 15, 2017)

以下のStanスクリプトに替えれば、unpairedな検定になります。

data {
	int<lower=0> N;
	real<lower=0> x1[N];
	real<lower=0> x2[N];
}

parameters {
	real s;
	real m;
}

model {
	real q[N];
	for (i in 1:N)
		q[i] <- x2[i] - x1[i];
	q~normal(m,s);
}

> d<-read.csv('https://raw.githubusercontent.com/ozt-ca/tjo.hatenablog.samples/master/r_samples/public_lib/DM_sampledata/ch3_2_2.txt',header=T,sep=' ')
> dat1<-list(N=nrow(d),x1=d$DB1,x2=d$DB2)
> fit1<-stan(file='ttest.stan',data=dat1,iter=1000,chains=4)
> fit1
Inference for Stan model: ttest.
4 chains, each with iter=1000; warmup=500; thin=1; 
post-warmup draws per chain=500, total post-warmup draws=2000.

      mean se_mean   sd  2.5%   25%   50%   75% 97.5% n_eff Rhat
s     0.73    0.01 0.18  0.49  0.61  0.70  0.80  1.18   224 1.02
m     0.74    0.01 0.20  0.34  0.61  0.74  0.85  1.16   424 1.01 # ←ココ！
lp__ -1.87    0.11 1.44 -6.02 -2.24 -1.41 -0.96 -0.67   181 1.03

Samples were drawn using NUTS(diag_e) at Wed Mar 15 16:45:13 2017.
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).
> fit1.coda<-mcmc.list(lapply(1:ncol(fit1),function(x) mcmc(as.array(fit1)[,x,])))
> plot(fit1.coda)

f:id:TJO:20170315164849p:plain

pairedの場合に比べて若干事後分布が広くなっていますが、それでも2.5%点が0を超えているのでOKと見て良さそうです。

カイ二乗検定をベイジアン モデリングでやってみる

こちらは平均値の差のモデリングで済んだt検定に比べるとやや複雑で、二項ロジットの枠組みの中でモデリングすることになります。式に書き下すとこんな感じでしょうか。

$cv \sim Bi(cv + ncv, p)$

$p = \displaystyle \frac{exp(a + bx + q)}{1 + exp(a + bx + q)}$

（ただし $cv$ はCV数、 $ncv$ は非CV数、 $p$ は二項分布の確率パラメータ、 $a$ は定数項、 $b$ は介入効果パラメータ、 $x$ は介入の有無を表す二値パラメータ、 $q$ は介入試験ごとの個体差、2行目はいわゆるinverse logit関数）

この場合は介入効果パラメータである $b$ の事後分布が十分に0から離れていれば、施策介入によるCV引き上げ効果ありと推論できることになります。で、Stanでの演算ですが実はカイ二乗検定のメタアナリシスで使ったStanコードがそのまま転用できます。

ただし、以下のコード中個体差を表すcl（上のモデル式中の $q$ ）はそもそも0になるはず*1で不要なので、計算負荷が無駄に増えるのが気に入らないという人は削除しちゃってください。もしかしたらわずかにサンプリング結果が変わるかもですが、最終的な結論には影響しないはずです。

data {
	int<lower=0> N;
	int<lower=0> ncv[N];
	int<lower=0> cv[N];
	real<lower=0,upper=1> iv[N];
}

parameters {
	real cl[N/2];
	real s;
	real a;
	real b;
}

model {
	real p[N];
	s~uniform(0,1e4);
	cl~normal(0,s);
	for (i in 1:N/2){
		p[2*i-1]<-inv_logit(a+b*iv[2*i-1]+cl[i]);
		p[2*i]<-inv_logit(a+b*iv[2*i]+cl[i]);
	}
	for (i in 1:N)
		cv[i]~binomial(ncv[i]+cv[i],p[i]);
}

以前の記事同様にbinom_hb_gen.stanという名前で保存した上で、以下のRコードでkickします。

> x<-data.frame(ncv=c(117,32),cv=c(25,16),iv=c(0,1))
> dat2<-list(N=2,ncv=x$ncv,cv=x$cv,iv=x$iv)
> fit2<-stan(file='binom_hb_gen.stan',data=dat2,iter=1000,chains=4)

SAMPLING FOR MODEL 'binom_hb_gen' NOW (CHAIN 1).

Chain 1, Iteration:   1 / 1000 [  0%]  (Warmup)
# ... #
Chain 4, Iteration: 1000 / 1000 [100%]  (Sampling)# 
#  Elapsed Time: 2.40675 seconds (Warm-up)
#                2.91676 seconds (Sampling)
#                5.32351 seconds (Total)
# 
# ... omitted ... #
 
> fit2
Inference for Stan model: binom_hb_gen.
4 chains, each with iter=1000; warmup=500; thin=1; 
post-warmup draws per chain=500, total post-warmup draws=2000.

         mean se_mean     sd     2.5%     25%     50%     75%   97.5% n_eff Rhat
cl[1] -283.14  275.47 560.17 -1371.49 -688.51 -158.22   63.54  768.06     4 1.42
s      902.56  325.75 808.28    41.42  279.46  671.98 1297.48 3132.93     6 1.79
a      281.58  275.48 560.18  -769.40  -65.01  156.69  686.87 1370.13     4 1.42
b        0.87    0.04   0.40     0.13    0.58    0.88    1.15    1.68    83 1.03 # ←ココ！
lp__  -104.30    0.28   1.37  -107.43 -105.07 -104.29 -103.44 -101.66    24 1.18

Samples were drawn using NUTS(diag_e) at Wed Mar  9 22:08:31 2016.
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).
> fit2.coda<-mcmc.list(lapply(1:ncol(fit2),function(x) mcmc(as.array(fit2)[,x,])))
> plot(fit2.coda)

f:id:TJO:20160309221515p:plain

f:id:TJO:20160309233539p:plain

ということで、介入効果パラメータ $b$ の事後分布の2.5%点が0より大きいこと、そして概ねまともな単峰性分布でサンプリングが収束しているとみなせることから*2、施策による介入効果は十分に（0より）大きいと言って良さそうです。

やってみて思ったこと

要はこれってconfidence intervalを明示的に出すというやり方のさらに細かい方法論なので、普通にCI出せばいいのかなと思いました。ただしCIが出しにくいケースでどうするかというのは無論あって、それは課題かなと。

そして、この程度の検定にわざわざStan引っ張り出してきてベイジアンモデリングするのって、ぶっちゃけ完全にオーバーキルだと思うんですよね（汗）。。。とは言え、場合によってはただ漫然と検定するだけでなく、このようにきちんとモデリングした結果を見せる必要があるのかもと思えば、覚えておいて損はないかもです。

*1:個体というかstudyが1回のみなので、個体差は出ない

*2:他のパラメータがきちんと最頻値と最大値が一致してなおかつ0が点推定値になるものの、事後分布自体は微妙な二峰性みたいな感じだったりするのが何とも。。。