頻度論的な統計手法のオルタナティブ

相関分析

データの格納

score_1<-SNFC_1$btheta
score_2<-SNFC_1$kakusan
score<-scale(cbind(score_1,score_2))
ni<-1317
nt<-2
data<-list(y=score,ni=ni,nt=nt)
par<-c("phi")
war<-1000
ite<-21000
see<-1234
dig<-3
cha<-5

ピアソンの積率相関

cor(score_1,score_2)

ピアソンの積率相関 ver.ベイズ

stanコード格納

corstan <- '
data{
	int ni;
	int nt;
	matrix[ni,nt] y;
}
parameters{
	vector[2] mu;
	corr_matrix[2] phi;
}
model{
	for(i in 1:ni){
		y[i]~multi_normal(mu,phi);
	}
}
'

MCMC実行

fit_valid<-stan(model_code=corstan,data=data,pars=par,verbose=F,seed=see,chains=cha,warmup=war,iter=ite)
print(fit_valid,pars=par,digits_summary=dig)

このstanコードは豊田(2017)を引用

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

データの格納

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)