Rstanで2×2クロス集計表の検定をしてみる

さて、統計でRstanを使ってみましょう。
薬のA症状に対する効果とB症状に対する効果でどうだったかなんかを検証するときは2×2のクロス集計表に対する統計検定で行われると思います。

通常の検定ではχ二乗検定が使われますが数の小さなマスがあるのでFisherの直接確立検定が用いられるところではあるかなと思います。この際はp=0.0016というなかなか素敵な値が出たわけですが次の場合はどうでしょう

この場合はχ二乗検定でp=0.077となりp>0.05なので棄却されます。ただ、サンプル数が27というのが多いのか少ないのか。同じ割合でサンプル数が倍の(C症状、D症状）に対し効果あり(20,10)と効果なし(6,18)になるとp=0.0056となります。一般的にサンプル数が増えると微妙な差で判定が変わることが多いです。じゃあどの程度のサンプル数が適切かなどを議論するのですが後ろ向き研究だとサンプル数が決まっていて変えられないなんてことがあるわけです。サンプル数の問題はめちゃくちゃめんどくさいです。

SO・KO・DE

　Rstanを使ったハミルトニアンモンテカルロ法を使ってみましょう。C症状に効果があるなし問いのは、ある期間の間に起こった事象の数というポアソンの分布に従いそうです。D症状も同様です。この場合は事前確率をガンマ分布にすることが一般的そうです。そうすると得られたデータがガンマ分布に従っているものとして計算します。

最低限のプログラムは必要ですので泣く泣く.stanのファイルを作ります。。しくしく。

data {
 int<lower=0> x[2];
 int<lower=0> 2[2];
 int<lower=0> x2[N];
}

parameteres {
 real<lower=0, upper=1> p[2];
}

transformed parameteres{
}

model {
for (i in 1:2){                            ←2つのカテゴリーだからiは1から２まで代入
 x[i] ~ binomial(n[i],p[i]);        ←乱数を発生する
 }
}

generates quantities{
 int <lower=0> xaste[2];
 real p_sa;
 real p_hi;
 real Odds_hi;
 real Odds[2];
 p_sa = p[1]-p[2];
 p_hi = p[1]/p[2];
 Odds_hi[1] = p[1]/(1-p[1]);
 Odds_hi[2] = p[2]/(1-p[2]);
 Odds = Odd[1]/Odds[2]
for (i in 1:2){
 xaste[i] = binomial_rng(n[i],p[i])
 }
}

最初のdataはデータ処理です。parameters〜transformed parameteresはxが変化するときに変化する事後確率pを求めています。次に出力をしていきます・・・（白目）generates quantitiesで2群にどんくらい差があるのか計算します。χ二乗検定やFisherの直接確率検定ではなかなか2群の差が体感できないのでこの点は大きな利点があります。stanのファイルは型のある統計ごとに作っておけば再利用できます。

Rを起動して打ち込んでいきます。。データ入力以外はほぼ流れ作業。。
library(rstan)   #Rstanの起動
rstan_options(auto_write=T)  #高速化 
options(mc.cores=parallel::detectCores())  #マルチコア処理
x<-c(10,5); n<- c(13,14)   #データ入力：x=効いた数、n=総数
dat <- list(x=x,n=n)
scr <- "（ファイル名）.stan"  #作ったstanスクリプトファイルを読み込む
prob <- c(0.025,0.05,0.5,0.95,0.975)     #2.5%,5%,50%,95%,97.5%で打ち出す
cha <- 4         #4回シミュレーションをして分布重なりを比較（一致するか）
war <- 1000       #収束しそうにない最初の1000回の計算を切り捨て
ite <- 21000      #乱数の数（21000-1000)で20000個残る予定
fi <- NA
fit <-stan(file=scr, data=dat, chains=cha, warmup=war, iter=ite, fit=fi)
ext <- extract(fit, permuted=FALSE);
print(fit, probs = prob)

そうすると・・・
Inference for Stan model: stan.
5 chains, each with iter=21000; warmup=1000; thin=1;
post-warmup draws per chain=20000, total post-warmup draws=1e+05.

mean se_mean sd 2.5% 5% 50% 95% 97.5% n_eff Rhat
p[1] 0.73 0.00 0.11 0.49 0.53 0.74 0.90 0.92 85743 1
p[2] 0.38 0.00 0.12 0.16 0.19 0.37 0.58 0.62 87074 1
xaste[1] 9.52 0.01 2.11 5.00 6.00 10.00 13.00 13.00 92478 1
xaste[2] 5.25 0.01 2.41 1.00 2.00 5.00 9.00 10.00 93070 1
p_sa 0.36 0.00 0.16 0.02 0.08 0.36 0.61 0.65 84724 1
p_hi 2.20 0.00 0.97 1.04 1.15 1.98 3.95 4.63 72258 1
Odds_hi 7.30 0.03 8.04 1.10 1.40 5.02 20.53 27.33 62326 1
Odds[1] 3.66 0.01 2.88 0.97 1.15 2.90 8.60 10.90 59775 1
Odds[2] 0.67 0.00 0.37 0.20 0.24 0.59 1.37 1.61 80629 1
lp__ -20.32 0.00 1.04 -23.13 -22.40 -20.00 -19.34 -19.31 43683 1

Samples were drawn using NUTS(diag_e) at Thu Dec 15 16:48:25 2022.
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).

こんな感じのが出ます。p1 73%、p2 38%、p1-p2 36%
p1-p2>0となれば検定としては有意差ありなので5%点をみると0.08と、95%の範囲で有意差ありとなります。95%確信区間(2.5%〜97.5%)での差が2〜65%で平均36%というのをどのように捉えるかですがちょっと目にみえる差が実感できるかと思います。

他にはbrmsパッケージを使う方法なんかも・・・。

えっと・・・みなさん帰ってきてますか・・・？

ちなみにJASPという統計ソフト（無料！）を使うとBayesA/Bテストはできますが詳しい差や比などは出すことができないのかめんどくさいかどちらかのようです。Rstanも面倒ではありますが。これからの更新に期待です。t検定は普通にできそうです（試してない）。