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

# パラメータ設定
effect_size <- 0.5 # Cohen's d
sig_level <- 0.05 # 有意水準 alpha
power_level <- 0.80 # 検出力 1-beta
n_ratio <- 2 # サンプルサイズの比 n2/n1

# n1を小さい値から1つずつ増やし、検出力が目標を超える最初の値を見つける
n1_current <- 2 # 探索を開始するn1の初期値
calculated_power <- 0

# 検出力が目標値(power_level)に達するまでループを続ける
while (calculated_power < power_level) {
  n1_current <- n1_current + 1
  n2_current <- n1_current * n_ratio

  # 現在のn1, n2における検出力を計算
  # power = NULLとすることで、検出力を計算させる
  power_result <- pwr.t2n.test(
    d = effect_size,
    sig.level = sig_level,
    power = NULL,
    n1 = n1_current,
    n2 = n2_current
  )

  # 計算された検出力を更新
  calculated_power <- power_result$power
}

# ループを抜けた時点のn1, n2が求めるサンプルサイズ
n1_required <- n1_current
n2_required <- n1_required * n_ratio

# 結果の表示
print(paste("介入群(n1)に必要なサンプルサイズ:", n1_required))
print(paste("対照群(n2)に必要なサンプルサイズ:", n2_required))
print(paste("その際の検出力:", round(calculated_power, 3)))

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

# シナリオの母集団パラメータ
mean_g1 <- 7.5 # 介入群の平均点の上昇分
mean_g2 <- 0 # 対照群の平均点の上昇分
sd_val <- 15 # 両群共通の標準偏差
# (μ1 - μ2) / σ = (7.5 - 0) / 15 = 0.5 となり、効果量と一致

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

# シミュレーションループ
for (i in 1:n_sim) {
  # 1. 対立仮説が真である状況のデータを生成
  # 介入群データ (新アプリ)
  group1_data <- rnorm(n = n1_required, mean = mean_g1, sd = sd_val)
  # 対照群データ (従来学習)
  group2_data <- rnorm(n = n2_required, mean = mean_g2, sd = sd_val)

  # 2. Welchのt検定を実行 (不等分散を仮定)
  t_test_result <- t.test(group1_data, group2_data, var.equal = FALSE)

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

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

print(paste("シミュレーションによる検出力:", round(simulated_power, 3)))

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

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

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

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