■

平均の差の推定シリーズ
古典統計学とベイズ統計学の比較してみました。

今回は手持ちのデータで薬剤sと薬剤Tを使用したときに

どの程度parameterであるPRdeltaが変化したかのデータです。

2群間の平均の差をいろんな方法で比較していきます。

まずデータをBBとして読み込みます。

その後、X1に薬剤Tを使用した群、X2に薬剤sを使用した群を割り付けます。

> X1=subset(BB,Medicine=="T")$PRdelta　　
> X2=subset(BB,Medicine=="s")$PRdelta
> M=nrow(subset(BB,Medicine=="T"))　　　　　
> N=nrow(subset(BB,Medicine=="s"))

まずは通常通りWelch のt検定です。

> t.test(X1,X2)

Welch Two Sample t-test

data: X1 and X2
t = 0.68164, df = 45.735, p-value = 0.4989
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-3.652168 7.391298
sample estimates:
mean of x mean of y
13.00000 11.13043

もちろん古典統計学では統計学的有意差なし、

非劣性マージンを”2”としたときに非劣性は証明できません。

 

次に、ベイズ統計学です。
マルコフ連鎖モンテカルロ法で平均の差Welch のt検定です。

データをlist形式にしてdatに格納します。

> dat <- list(M=M, N=N, x1=X1, x2=X2)

モデルを指定します。

今回は平均、分散は一様分布に従うと仮定しています。

> 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;
+ }
+ "

> fit1 <- stan(model_code=model, data=dat, iter=1000, chain=4)
> print(fit1)

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
s1 12.10 0.04 1.76 9.25 10.86 11.92 13.11 16.15 1884 1
s2 8.23 0.03 1.35 6.09 7.26 8.04 9.02 11.29 1936 1
m1 12.96 0.06 2.36 8.27 11.41 13.00 14.49 17.59 1773 1
m2 11.11 0.04 1.71 7.71 9.99 11.13 12.19 14.47 1660 1
diff 1.85 0.07 2.88 -3.82 -0.06 1.76 3.76 7.69 1800 1
lp__ -134.65 0.04 1.43 -138.30 -135.35 -134.33 -133.59 -132.83 1033 1

 

このモデルでは95%ベイズ信用区間からは有意差は認めません。

2群間の平均の差は1.85程度と推定されています。

さらに別のモデルを指定します。
平均は幅の広い正規分布、分散は半コーシー分布を指定しています。

> 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;
+ }
+ "

> fit0 <- stan(model_code=model0, data=dat, iter=1000, chain=4)
> print(fit0)

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
s1 11.76 0.04 1.73 9.05 10.48 11.56 12.82 15.64 1605 1
s2 7.91 0.03 1.23 5.97 7.06 7.75 8.64 10.82 1797 1
m1 13.05 0.05 2.21 8.64 11.58 13.02 14.54 17.46 1827 1
m2 11.13 0.04 1.67 7.83 10.00 11.14 12.23 14.41 1714 1
diff 1.92 0.06 2.76 -3.51 0.23 2.01 3.68 7.36 2189 1
lp__ -137.69 0.05 1.45 -141.29 -138.42 -137.37 -136.63 -135.87 850 1

このモデルでは95%ベイズ信用区間からは有意差はないと判断されます。

また2群間の平均の差は1.92値程度と推定できます。

 

次回は線形回帰で確認していきたいと思います。