古典統計学とベイズ統計での2群間の平均の差の検定

今回はRの組み込みのデータセットToothGrowthを用いて
2群間の平均の差を古典統計学とベイズ統計の実装であるRstanをもちいて
検定（推定）していきたいとおもいます。
supp　サプリの種類　”VC”　”OJ”
dose　サプリの用量　0.5㎎　1.0㎎　2.0㎎
len　歯の長さ　㎜

head(ToothGrowth)
   len supp dose
1  4.2   VC  0.5
2 11.5   VC  0.5
3  7.3   VC  0.5
4  5.8   VC  0.5
5  6.4   VC  0.5
6 10.0   VC  0.5

箱ひげ図にすると次の通りです。

boxplot(len ~ supp , data = ToothGrowth)

このOJ群とVC群の平均に差があるかの検定を考えます。
ここでは等分散を仮定しないウェルチの2群間の平均の差の検定です。

 t.test(len ~ supp, data=ToothGrowth,var=F)

Two Sample t-test

data: len by supp
t = 1.9153, df = 58, p-value = 0.06039
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-0.1670064 7.5670064
sample estimates:
mean in group OJ mean in group VC
20.66333 16.96333

結果の見方ですが、95%信頼区間に0が含まれているので、
2群間の平均の差に統計学的有意差はなく、p値0.060となっています。

ここまでが古典統計学での記述統計と検定でした。
つづいてベイズ統計による推定を行います。
OJ群、VC群のサンプル数をそれぞれM、Nとします。
その後、X1にOJ群を、X2にVC群を割り付けます。

M=nrow(subset(ToothGrowth,supp=="OJ"))
N=nrow(subset(ToothGrowth,supp=="VC"))
X1=subset(ToothGrowth,supp=="OJ")$len
X2=subset(ToothGrowth,supp=="VC")$len
dat <- list(M=M, N=N, x1=X1, x2=X2)

モデル式は次の通りで、Rstanで推定しています。
マルコフ連鎖モンテカルロ法で平均の差Welch のt検定です
今回は平均、分散を一様分布に従うと仮定しています。

library(rstan)

model<-"
data{
int<lower=0> M;
int<lower=0> N;
real x1[M];
real x2[N];
}
parameters{
real<lower=0> s1;
real<lower=0> s2;
real m1;
real m2;
}
model{
for(i in 1:M)
x1[i] ~ normal(m1, s1);
for(i in 1:N)
x2[i] ~ normal(m2, s2);
}
generated quantities{
real diff;
diff=m1-m2;
}
"
fit <- stan(model_code=model, data=dat, iter=1000, chain=4)

mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
s1 6.90 0.02 0.96 5.36 6.21 6.80 7.48 9.18 1829 1
s2 8.63 0.03 1.19 6.68 7.80 8.51 9.32 11.34 2008 1
m1 20.67 0.03 1.28 18.24 19.83 20.67 21.50 23.16 2099 1
m2 17.02 0.04 1.60 13.96 15.94 17.01 18.06 20.15 1716 1
diff 3.65 0.05 2.01 -0.26 2.28 3.67 5.03 7.46 1936 1
lp__ -147.07 0.04 1.44 -150.84 -147.78 -146.73 -146.04 -145.26 1170 1

ぞれの平均と標準偏差、2群間の平均の差が推定されています。
95%ベイズ信用区間に0を含んでいるので、差はなさそうです。
まあウェルチの平均の差の検定と変わらないですね。

分散は半コーシー分布に従うことが多いので次のモデルでも検定。
平均は幅の広い正規分布、分散は半コーシー分布を指定しています。

> model0<-"
data{
int<lower=0> M;
int<lower=0> N;
real x1[M];
real x2[N];
}
parameters{
real<lower=0> s1;
real<lower=0> s2;
real m1;
real m2;
}
model{
m1~normal(0,10000);
m2~normal(0,10000);
s1~cauchy(0,5);
s2~cauchy(0,5);

for(i in 1:M)
x1[i] ~ normal(m1, s1);
for(i in 1:N)
x2[i] ~ normal(m2, s2);
}
generated quantities{
real diff;
diff=m1-m2;
}
"