２群の平均値の差の検討（t検定，ベイズ推定）

２群の平均値の差の検討（性差の検討を例に）

データの格納

x<-SNFC_1$theta[SNFC_1$sex=="2"] #第1群の指定
y<-SNFC_1$theta[SNFC_1$sex=="1"] #第2群の指定
#平均値の高い方をyにしておくと後で解釈しやすい
nx <- length(x)
ny <- length(y)
data <- list(N1=nx,N2=ny,x1=x,x2=y,c=0)

対応のないt検定

var.test(x, y) #等分散性の検定.p < 0.05 => 等分散でない．p >= 0.05 => 等分散である
t.test(x,y,var.equal=T) #t検定.等分散でない時はvar.equal=F(Welch検定)

対応のないt検定 ver.ベイズ

stanコード格納

pairedF <- '
data{
  int<lower=0> N1;
  int<lower=0> N2;
  real x1[N1];
  real x2[N2];
  real c;
}
parameters{
  real mu1;
  real mu2;
  real<lower=0> sigma1;
  real<lower=0> sigma2;
}
model{
  x1 ~ normal(mu1,sigma1);
  x2 ~ normal(mu2,sigma2);
}
generated quantities{
  real delta;
  real d_overC;
  real cohen_d;
  delta = mu2 - mu1;
  d_overC = (delta>c? 1:0);
  cohen_d = delta /((sigma1^2*(N1-1)+sigma2^2*(N2-1)) /((N1-1)+(N2-1)))^0.5;
}
'

require(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

fit<-stan(model_code=pairedF,data=data,chains=5,warmup=1000,iter=21000)
print(fit,digits_summary=3)