平均の差を線形回帰で表現

引き続いて線形回帰で確認します。
t検定と線形回帰は数学的に同値です。

fit.lm=lm(len~supp,data=ToothGrowth)

結果を確認します。

summary(fit.lm)

Call:
lm(formula = len ~ supp, data = ToothGrowth)
Residuals:
Min 1Q Median 3Q Max
-12.7633 -5.7633 0.4367 5.5867 16.9367

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 20.663 1.366 15.127 <2e-16 ***
suppVC -3.700 1.932 -1.915 0.0604 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 7.482 on 58 degrees of freedom
Multiple R-squared: 0.05948, Adjusted R-squared: 0.04327
F-statistic: 3.668 on 1 and 58 DF, p-value: 0.06039

OJ群の平均は20.6、VC群の平均は16.3、平均の差(傾き)は3.70で、
95%信頼区間に0を含んでいますし、p値0.060です。
線形回帰での傾きに統計学的有意差はない。

引き続いてマルコフ連鎖モンテカルロ法を用いてベイズ推定します。
線形回帰を指定します。

まず、処置（supp）を0と1の数値に変換し、Stanに渡します。Stanでは線形回帰モデルの構造をシンプルに指定できます。

ToothGrowth$supp_numeric <- ifelse(ToothGrowth$supp == "VC", 1, 0) 

Stanに渡すデータリストを指定します。

stan_data <- list( N = nrow(ToothGrowth), len = ToothGrowth$len, supp = ToothGrowth$supp_numeric )

モデルを指定します

model.lm =
"
data {
  int<lower=0> N;             // サンプルサイズ
  vector[N] len;              // lenデータ
  int<lower=0, upper=1> supp[N]; // 処置グループ（0: OJ, 1: VC）
}

parameters {
  real alpha;                 // 切片（OJ群の平均）
  real beta;                  // suppの効果（OJとVCの平均の差）
  real<lower=0> sigma;        // 観測ノイズの標準偏差
}

model {
  for (i in 1:N) {
    len[i] ~ normal(alpha + beta * supp[i], sigma);
  }
}

generated quantities {
  real mean_diff = beta; // VCとOJの平均の差
}
"

モデルを実行します。

fit <- stan( file = model.lm, data = stan_data,
  iter = 2000, chains = 4 )

結果を確認します。

print(fit, pars = c("alpha", "beta", "sigma", "mean_diff"))

mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
alpha 20.65 0.03 1.45 17.88 19.68 20.64 21.60 23.53 2108 1
beta -3.68 0.04 2.04 -7.71 -5.06 -3.71 -2.31 0.34 2118 1
sigma 7.67 0.02 0.75 6.37 7.14 7.62 8.13 9.34 2369 1
mean_diff -3.68 0.04 2.04 -7.71 -5.06 -3.71 -2.31 0.34 2118 1

OJの平均値は20.７でVCとの差は傾きbetaで3.68となっており、
ベイズ95％信用区間は-7.71ー0.34で０を含んでおり、
平均の差はないと考えて良さそうです。