Rでサンプルサイズ設計:効果量･検定力･第一種過誤･第二種過誤

# パラメータ設定
mu_h0 <- 0 # 帰無仮説の母平均の差
mu_h1 <- 0.8 # 対立仮説の母平均の差（効果量大）
sd <- 1 # 母標準偏差（両群共通と仮定）
n <- 25 # 各群のサンプルサイズ
alpha <- 0.05 # 有意水準

# 標準誤差の計算
se <- sd * sqrt(1 / n + 1 / n)

# 棄却域の境界値（臨界値）を計算
critical_value <- qnorm(1 - alpha / 2, mean = mu_h0, sd = se)

# 分布のデータを作成
dist_df <- tibble(
  x = seq(-3 * se, mu_h1 + 3 * se, length.out = 400),
  h0_density = dnorm(x, mean = mu_h0, sd = se),
  h1_density = dnorm(x, mean = mu_h1, sd = se)
)

# プロット
ggplot(dist_df, aes(x = x)) +
  # 帰無仮説の分布
  geom_line(aes(y = h0_density, color = "帰無仮説 (H0)"), linewidth = 1) +
  geom_area(
    data = . %>% filter(x > critical_value),
    aes(y = h0_density), fill = "red", alpha = 0.3
  ) +
  geom_area(
    data = . %>% filter(x < -critical_value),
    aes(y = h0_density), fill = "red", alpha = 0.3
  ) +
  # 対立仮説の分布
  geom_line(aes(y = h1_density, color = "対立仮説 (H1)"), linewidth = 1) +
  geom_area(
    data = . %>% filter(x < critical_value & x > -critical_value),
    aes(y = h1_density), fill = "skyblue", alpha = 0.5
  ) +
  # 凡例と注釈
  geom_vline(xintercept = c(critical_value, -critical_value), linetype = "dashed", color = "black") +
  annotate("text", x = critical_value + 0.1, y = 0.6, label = "α/2", parse = TRUE, size = 5) +
  annotate("text", x = -critical_value - 0.1, y = 0.6, label = "α/2", parse = TRUE, size = 5) +
  annotate("text", x = 0.4, y = 0.1, label = "β", parse = TRUE, size = 5, color = "blue") +
  annotate("text", x = 1.0, y = 0.2, label = "検定力 (1-β)", parse = TRUE, size = 5, color = "navy") +
  scale_color_manual(name = "分布", values = c("帰無仮説 (H0)" = "red", "対立仮説 (H1)" = "blue")) +
  labs(
    title = "第一種過誤(α)と第二種過誤(β)の概念図",
    subtitle = paste0("サンプルサイズ(n)=", n, ", 効果量(d)=", mu_h1, ", 有意水準(α)=", alpha),
    x = "標本平均の差",
    y = "確率密度"
  ) +
  theme_minimal()

# シミュレーションを実行する関数を定義
run_power_simulation <- function(n, effect_size, n_sim = 500, alpha = 0.05) {
  mu1 <- 100
  sd <- 15
  mu2 <- mu1 + effect_size * sd

  p_values <- replicate(n_sim, {
    sample1 <- rnorm(n, mean = mu1, sd = sd)
    sample2 <- rnorm(n, mean = mu2, sd = sd)
    t.test(sample1, sample2)$p.value
  })

  power <- sum(p_values < alpha) / n_sim
  return(power)
}

# シミュレーションのパラメータを設定
sample_sizes <- seq(10, 200, by = 10)
effect_sizes <- c(0.2, 0.5, 0.8) # 小、中、大の効果量

# expand.gridで全ての組み合わせを作成
sim_params <- expand.grid(
  n = sample_sizes,
  d = effect_sizes
)

# シミュレーションを実行
# mapplyを使って各パラメータの組み合わせでシミュレーションを実行
set.seed(seed)
sim_results <- mapply(
  run_power_simulation,
  n = sim_params$n,
  effect_size = sim_params$d
)

# 結果をデータフレームに格納
results_df <- sim_params %>%
  mutate(
    power = sim_results,
    effect_size_label = factor(paste0("d = ", d))
  )

# 結果をプロット
ggplot(results_df, aes(x = n, y = power, color = effect_size_label)) +
  geom_line(linewidth = 1.0) +
  geom_point(size = 2.5) +
  geom_hline(yintercept = 0.8, linetype = "dashed", color = "red") +
  annotate("text", x = 150, y = 0.85, label = "目標検定力 (80%)", color = "red") +
  scale_y_continuous(labels = scales::percent, limits = c(0, 1)) +
  labs(
    title = "サンプルサイズと効果量が検定力に与える影響",
    x = "各群のサンプルサイズ (n)",
    y = "検定力 (Power)",
    color = "効果量 (Cohen's d)"
  ) +
  theme_bw(base_size = 14)

# 必要サンプルサイズの計算
# ------------------------------------------------------------------
# パラメータ設定
# ------------------------------------------------------------------
target_power <- 0.8 # 目標とする検定力 (80%)
sig_level <- 0.05 # 有意水準 (5%)

# ------------------------------------------------------------------
# 1. 効果量 大 (d = 0.8) の場合
# ------------------------------------------------------------------
power_calc_large <- pwr.t.test(
  d = 0.8,
  power = target_power,
  sig.level = sig_level,
  type = "two.sample", # 2つの独立したサンプルの検定
  alternative = "two.sided" # 両側検定
)


# ------------------------------------------------------------------
# 2. 効果量 中 (d = 0.5) の場合
# ------------------------------------------------------------------
power_calc_medium <- pwr.t.test(
  d = 0.5,
  power = target_power,
  sig.level = sig_level,
  type = "two.sample",
  alternative = "two.sided"
)


# ------------------------------------------------------------------
# 3. 効果量 小 (d = 0.2) の場合
# ------------------------------------------------------------------
power_calc_small <- pwr.t.test(
  d = 0.2,
  power = target_power,
  sig.level = sig_level,
  type = "two.sample",
  alternative = "two.sided"
)


# ------------------------------------------------------------------
# 結果の表示
# ------------------------------------------------------------------
cat("--- 有意水準5%、目標検定力80%を達成するために必要な各群のサンプルサイズ ---\n\n")

# 効果量 大 (d = 0.8)
cat("■ 効果量 大 (d = 0.8) の場合:\n")
print(power_calc_large)
# サンプルサイズは整数であるため、計算結果を切り上げる
cat("=> 必要なサンプルサイズ (各群):", ceiling(power_calc_large$n), "人\n\n")

# 効果量 中 (d = 0.5)
cat("■ 効果量 中 (d = 0.5) の場合:\n")
print(power_calc_medium)
cat("=> 必要なサンプルサイズ (各群):", ceiling(power_calc_medium$n), "人\n\n")

# 効果量 小 (d = 0.2)
cat("■ 効果量 小 (d = 0.2) の場合:\n")
print(power_calc_small)
cat("=> 必要なサンプルサイズ (各群):", ceiling(power_calc_small$n), "人\n")