Rでサンプルサイズ設計:2変数の相関分析

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

# パラメータ設定
effect_size_r <- 0.3 # 検出したい相関係数r
sig_level <- 0.05 # 有意水準 alpha
power_level <- 0.80 # 検出力 1-beta

# サンプルサイズの計算
# nをNULLにすることで、必要なサンプルサイズを計算させる
sample_size_calc <- pwr.r.test(
  r = effect_size_r,
  sig.level = sig_level,
  power = power_level,
  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 # シミュレーション回数

# シナリオの母集団パラメータ
# 平均ベクトル（スマホ利用時間6時間、睡眠時間7時間と仮定）
mu_vector <- c(6, 7)
# 標準偏差（スマホ2時間、睡眠1時間と仮定）
sd_vector <- c(2, 1)
# 相関係数
rho <- -0.3 # 負の相関を想定

# 共分散行列を作成
sigma_matrix <- matrix(
  c(
    sd_vector[1]^2, rho * sd_vector[1] * sd_vector[2],
    rho * sd_vector[1] * sd_vector[2], sd_vector[2]^2
  ),
  nrow = 2
)

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

# シミュレーションループ
for (i in 1:n_sim) {
  # 1. 対立仮説が真である状況の相関のあるデータを生成
  correlated_data <- mvrnorm(n = n_required, mu = mu_vector, Sigma = sigma_matrix)

  # 2. 相関係数の検定を実行
  cor_test_result <- cor.test(correlated_data[, 1], correlated_data[, 2])

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

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

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

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

# 帰無仮説が棄却されない場合 (相関係数 rho = 0)
rho_h0 <- 0
sigma_matrix_h0 <- matrix(
  c(
    sd_vector[1]^2, rho_h0 * sd_vector[1] * sd_vector[2],
    rho_h0 * sd_vector[1] * sd_vector[2], sd_vector[2]^2
  ),
  nrow = 2
)

p_values_h0 <- numeric(n_sim)

for (i in 1:n_sim) {
  # 1. 帰無仮説が棄却されない状況の無相関のデータを生成
  uncorrelated_data <- mvrnorm(n = n_required, mu = mu_vector, Sigma = sigma_matrix_h0)

  # 2. 相関係数の検定を実行
  cor_test_result_h0 <- cor.test(uncorrelated_data[, 1], uncorrelated_data[, 2])

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

type1_error_rate <- mean(p_values_h0 < sig_level)

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