Rでベイズファクター

# 必要パッケージの読み込み
library(BayesFactor)

# 再現性のための乱数シードの設定
seed <- 20260819
set.seed(seed)

# サンプルサイズの範囲（10から1000まで10刻み）
n_sizes <- seq(10, 1000, by = 10)

# 結果を格納するベクトルの準備
p_h0 <- numeric(length(n_sizes))
bf_h0 <- numeric(length(n_sizes))

p_h1 <- numeric(length(n_sizes))
bf_h1 <- numeric(length(n_sizes))

# 対立仮説が正しい場合の「真の効果量（平均値）」
true_effect <- 0.3

for (i in seq_along(n_sizes)) {
  n <- n_sizes[i]

  # --- パターンA: 帰無仮説 (H0) が正しい場合 (真の平均 = 0) ---
  data_h0 <- rnorm(n, mean = 0, sd = 1)

  # 従来のt検定 (p値)
  p_h0[i] <- t.test(data_h0)$p.value
  # ベイズファクター (BF10) の抽出
  bf_h0[i] <- extractBF(ttestBF(data_h0), onlybf = TRUE)

  # --- パターンB: 対立仮説 (H1) が正しい場合 (真の平均 = 0.3) ---
  data_h1 <- rnorm(n, mean = true_effect, sd = 1)

  # 従来のt検定 (p値)
  p_h1[i] <- t.test(data_h1)$p.value
  # ベイズファクター (BF10) の抽出
  bf_h1[i] <- extractBF(ttestBF(data_h1), onlybf = TRUE)
}

# グラフの描画設定（2行2列の表示）
par(mfrow = c(2, 2), mar = c(4.5, 4.5, 2.5, 1))

# 1. H0が真のときのp値
plot(n_sizes, p_h0,
  type = "b", pch = 19, col = "royalblue",
  xlab = "サンプルサイズ (N)", ylab = "p値",
  main = "H0が真のとき: p値の推移", ylim = c(0, 1)
)
abline(h = 0.05, col = "red", lty = 2) # 有意水準 5% の線

# 2. H0が真のときのBF10 (y軸は対数スケール)
# BF10 < 1 は帰無仮説を支持することを意味します
plot(n_sizes, bf_h0,
  type = "b", pch = 19, col = "forestgreen",
  xlab = "サンプルサイズ (N)", ylab = "BF10 (対数軸)", log = "y",
  main = "H0が真のとき: BF10の推移"
)
abline(h = 1, col = "red", lty = 2) # BF10 = 1 (証拠なし) の基準線

# 3. H1が真のときのp値
plot(n_sizes, p_h1,
  type = "b", pch = 19, col = "royalblue",
  xlab = "サンプルサイズ (N)", ylab = "p値",
  main = "H1が真のとき: p値の推移", ylim = c(0, 1)
)
abline(h = 0.05, col = "red", lty = 2)

# 4. H1が真のときのBF10 (y軸は対数スケール)
plot(n_sizes, bf_h1,
  type = "b", pch = 19, col = "forestgreen",
  xlab = "サンプルサイズ (N)", ylab = "BF10 (対数軸)", log = "y",
  main = "H1が真のとき: BF10の推移"
)
abline(h = 1, col = "red", lty = 2)

# 描画設定の初期化
par(mfrow = c(1, 1))