Rでサンプルサイズ設計:t検定(2群、対応なし、同じサンプルサイズ)

cat("--- ステップ1: サンプルサイズの計算 ---\n")

# パラメータ設定
effect_size <- 0.5 # Cohen's d
sig_level <- 0.05 # 有意水準 alpha
power_level <- 0.80 # 検出力 1-beta

# サンプルサイズの計算
# nをNULLにすることで、必要なサンプルサイズ（各群）を計算させる
sample_size_calc <- pwr.t.test(
  d = effect_size,
  sig.level = sig_level,
  power = power_level,
  type = "two.sample", # 2群のt検定
  alternative = "two.sided" # 両側検定
)

# 結果の表示
print(sample_size_calc)

# 必要なサンプルサイズ（各群）を取得（切り上げ）
n_required <- ceiling(sample_size_calc$n)
cat("\n")
cat(paste0("結論: 介入群、対照群それぞれに最低 ", n_required, " 人の参加者が必要です。\n"))

cat("--- ステップ2: 検出力のシミュレーション ---\n")

# シミュレーションのパラメータ設定
set.seed(seed) # 結果を再現可能にする
n_sim <- 10000 # シミュレーション回数

# シナリオの母集団パラメータ
mean_placebo <- 140 # 対照群（プラセボ）の平均血圧 (mmHg)
sd_bp <- 12 # 両群共通の標準偏差 (mmHg)
# 介入群の平均血圧（6mmHg低下）
mean_supplement <- mean_placebo - (effect_size * sd_bp)

# 結果を格納するベクトル
p_values_h1 <- numeric(n_sim)

# シミュレーションループ
for (i in 1:n_sim) {
  # 1. 対立仮説が真である状況のデータを生成
  # 介入群データ (新サプリ)
  supplement_data <- rnorm(n = n_required, mean = mean_supplement, sd = sd_bp)
  # 対照群データ (プラセボ)
  placebo_data <- rnorm(n = n_required, mean = mean_placebo, sd = sd_bp)

  # 2. Welchのt検定を実行（不等分散を仮定する）
  t_test_result <- t.test(supplement_data, placebo_data, var.equal = FALSE)

  # 3. p値を保存
  p_values_h1[i] <- t_test_result$p.value
}

# シミュレーションによる検出力の計算
simulated_power <- mean(p_values_h1 < sig_level)

cat(paste0("\nシミュレーションによる検出力: ", round(simulated_power, 3), "\n"))

cat("--- ステップ4: 第1種の過誤のシミュレーション ---\n")

set.seed(seed)

# 結果を格納するベクトル
p_values_h0 <- numeric(n_sim)

# シミュレーションループ (帰無仮説が棄却されない)
for (i in 1:n_sim) {
  # 1. 帰無仮説が棄却されない状況のデータを生成 (平均血圧に差がない)
  group1_data_h0 <- rnorm(n = n_required, mean = mean_placebo, sd = sd_bp)
  group2_data_h0 <- rnorm(n = n_required, mean = mean_placebo, sd = sd_bp)

  # 2. Welchのt検定を実行
  t_test_result_h0 <- t.test(group1_data_h0, group2_data_h0, var.equal = FALSE)

  # 3. p値を保存
  p_values_h0[i] <- t_test_result_h0$p.value
}

# シミュレーションによる第1種の過誤率の計算
type1_error_rate <- mean(p_values_h0 < sig_level)

cat(paste("\nシミュレーションによる第1種の過誤率:", round(type1_error_rate, 3), "\n"))