Rでサンプルサイズ設計:t検定(2群、対応あり)

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

# パラメータ設定
effect_size_dz <- 0.5 # 差分の効果量 d_z
sig_level <- 0.05 # 有意水準 alpha
power_level <- 0.80 # 検出力 1-beta

# サンプルサイズの計算
# nをNULLにすることで、必要なサンプルサイズを計算させる
sample_size_calc <- pwr.t.test(
  d = effect_size_dz,
  sig.level = sig_level,
  power = power_level,
  type = "paired", # 対応のあるt検定
  alternative = "greater" # 片側検定（スコアの向上を期待）
)

# 結果の表示
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_pre <- 50 # 事前テストの平均スコア
sd_diff <- 8 # 差分の標準偏差
mean_diff <- 4 # 期待される差分の平均 (4点向上)
mean_post <- mean_pre + mean_diff # 事後テストの平均スコア

# データ生成のための相関パラメータを設定
# (sd_diff^2 = sd_pre^2 + sd_post^2 - 2*rho*sd_pre*sd_post の関係を利用)
sd_pre <- 10 # 事前スコアの標準偏差を仮定
sd_post <- 10 # 事後スコアの標準偏差を仮定
# 相関係数 rho を逆算
rho <- (sd_pre^2 + sd_post^2 - sd_diff^2) / (2 * sd_pre * sd_post)

# 平均ベクトルと共分散行列
mu_vector <- c(mean_pre, mean_post)
sigma_matrix <- matrix(
  c(
    sd_pre^2, rho * sd_pre * sd_post,
    rho * sd_pre * sd_post, sd_post^2
  ),
  nrow = 2
)

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

# シミュレーションループ
for (i in 1:n_sim) {
  # 1. 対立仮説が真である状況の相関のあるデータを生成
  # (mvrnorm: 多変量正規分布から乱数生成)
  paired_data <- mvrnorm(n = n_required, mu = mu_vector, Sigma = sigma_matrix)

  # 2. 対応のあるt検定を実行
  t_test_result <- t.test(paired_data[, 2], paired_data[, 1],
    paired = TRUE, alternative = "greater"
  )

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

# 帰無仮説が棄却されない場合 (平均に差がない)
mu_vector_h0 <- c(mean_pre, mean_pre) # 事後も平均は同じ

p_values_h0 <- numeric(n_sim)

for (i in 1:n_sim) {
  paired_data_h0 <- mvrnorm(n = n_required, mu = mu_vector_h0, Sigma = sigma_matrix)
  t_test_result_h0 <- t.test(paired_data_h0[, 2], paired_data_h0[, 1],
    paired = TRUE, alternative = "greater"
  )
  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"))