Rでサンプルサイズ設計:t検定(1群の平均値の検定)

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

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

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

# 結果の表示
print(sample_size_calc)

# 必要なサンプルサイズ（ペットボトルの本数）を取得（切り上げ）
n_required <- ceiling(sample_size_calc$n)
cat(paste0("\n結論: この検査には最低 ", n_required, " 本のペットボトルが必要です。\n"))

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

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

# シナリオの母集団パラメータ
mu_0 <- 500 # 基準値 (ml)
sd_val <- 4 # 標準偏差 (ml)
# 対立仮説が真である場合の平均値（2mlずれている）
mu_1 <- mu_0 - (effect_size * sd_val) # 498mlを想定

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

# シミュレーションループ
for (i in 1:n_sim) {
  # 1. 対立仮説が真である状況のデータを生成
  sample_data <- rnorm(n = n_required, mean = mu_1, sd = sd_val)

  # 2. 1群のt検定を実行
  t_test_result <- t.test(sample_data, mu = mu_0)

  # 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")

p_values_h0 <- numeric(n_sim)

for (i in 1:n_sim) {
  # 1. 帰無仮説が棄却されない状況のデータを生成 (平均は基準値と同じ500ml)
  sample_data_h0 <- rnorm(n = n_required, mean = mu_0, sd = sd_val)

  # 2. 1群のt検定を実行
  t_test_result_h0 <- t.test(sample_data_h0, mu = mu_0)

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

type1_error_rate <- mean(p_values_h0 < sig_level)

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