対応のある分散分析（１要因）forベイズ

対応のある分散分析（１要因）

データの読み込み

#作業ディレクトリの設定（適宜変更）
setwd("C:/Users/turid/OneDrive/ドキュメント")
#データ読み込み
library(readr)
SNFC_1 <- read_csv("SNFC.csv")
library(rstan)                 #パッケージrstanの呼び出し
dat1=na.omit(SNFC_1)
x <- subset(dat1,select = c("kakusan","tokushu","tankyugouri"))
n <- nrow(x)
a=3
data <-list(N = n, X = x)

stanコード格納

stancode <- "
data{
int<lower=0> N;
vector[3] X[N];
}

parameters{
vector[3] mu;
real<lower=0> sig[3];
real<lower=-1,upper=1> rho;
}

transformed parameters{
matrix[3,3] Sigma;
Sigma[1,1] = sig[1]*sig[1];
Sigma[2,2] = sig[2]*sig[2];
Sigma[3,3] = sig[3]*sig[3];
Sigma[1,2] = sig[1]*sig[2]*rho;
Sigma[1,3] = sig[1]*sig[3]*rho;
Sigma[2,3] = sig[2]*sig[3]*rho;
Sigma[2,1] = Sigma[1,2];
Sigma[3,1] = Sigma[1,3];
Sigma[3,2] = Sigma[2,3];
}

model{
X ~ multi_normal(mu, Sigma);
rho ~ uniform(-1,1);
sig ~ cauchy(0,2.5);
}

generated quantities{
real delta21;
real delta31;
real delta32;
real d_overC1;
real d_overC2;
real d_overC3;

delta21 = mu[2] - mu[1];
delta31 = mu[3] - mu[1];
delta32 = mu[3] - mu[2];
d_overC1 = (delta21>0? 1:0);
d_overC2 = (delta31>0? 1:0);
d_overC3 = (delta32>0? 1:0);
}
"

rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
fit <- stan(model_code = stancode, iter=400, warmup=100, chains=4, data = data)
print(fit, digits_summary = 2) 

連言命題が正しい確率

df<-data.frame(extract(fit))
t<-sum((df$mu.3 -df$mu.1 > 0) & (df$mu.3 -df$mu.2 > 0) =="TRUE")
round(t/ nrow(df),3)
t<-sum((df$mu.1 -df$mu.3 > 0) & (df$mu.1 -df$mu.2 > 0)  =="TRUE")
round(t/ nrow(df),3)

特定のパラメーターだけプロット

library(ggmcmc)
df_param <- ggs(fit)
df_param %>% {
  print(class(.))
  dplyr::tbl_df(.)
}
fit_param_mu <- df_param %>% dplyr::filter(grepl("^mu", Parameter))
ggs_caterpillar(fit_param_mu)