Student t検定と Welch t検定の前提条件

今更ですが... Student t検定は母分散が等しいことが前提ですから、等しいと言い切れなければWelch t検定を使いましょう

告白します。
そんな基本的なこと知ってましたけど、リアルワールドではそもそも等分散なんてことないのにみんな使ってるし、だいたい正規性を満たすという条件もあやふやなんだから、少しくらい分散が違ったって Student t検定でいいじゃん！ どうせたいして違わないだろうし！
くらいに思ってました。...確かめもしないで。

シミュレーションやってみた。まずは結果のグラフから

検定が正しければ、棒グラフの高さが赤の点線に揃うはずです。

Student t検定がダメ使いづらいというのは本当だった

正直、ここまで違うとは驚きでした。
これ見たら、全部 Welch t検定でいいんじゃね？！ってなりました。

「前提条件」を満たしていれば正しいのだから、ダメは言いすぎですね。

こんなに違うんだ...

Student t検定は、n が等しくない時に分散が違うとむちゃくちゃ有意差が出やすくなってた！
要するに、誤って有意差があるとしてしまう可能性 (第一種の過誤) が想定以上に高いことがあるということ。
p=0.05 だと思ってたのに実は p=0.25 だったとか。目を疑いました。
~~発表するときには有意差がありますって言いやすくて助かるけどね。~~

Welch t検定がダメという理由は少なくともこのグラフにはない

保守的なので有意差が出にくいかと思い込んでいたが、そんなことはなかったです。

R でのシミュレーション

以下がシミュレーションのコードです。
https://rdrr.io/snippets/
で試せます。
負荷がかかるので遠慮して n = 2000 にしてます。
ローカル環境なら適当に増やしてください。

# 必要なライブラリの読み込み
library(ggplot2)
library(tidyr)

# シミュレーション関数の定義
simulate_t_tests <- function(n1, n2, mean1, mean2, sd1, sd2, num_simulations = 2000, alpha = 0.05) {
  student_t_rejections <- 0
  welch_t_rejections <- 0
  
  for (i in 1:num_simulations) {
    group1 <- rnorm(n1, mean = mean1, sd = sd1)
    group2 <- rnorm(n2, mean = mean2, sd = sd2)
    
    # Student's t-test
    student_t_result <- t.test(group1, group2, var.equal = TRUE)
    if (student_t_result$p.value < alpha) {
      student_t_rejections <- student_t_rejections + 1
    }
    
    # Welch's t-test
    welch_t_result <- t.test(group1, group2, var.equal = FALSE)
    if (welch_t_result$p.value < alpha) {
      welch_t_rejections <- welch_t_rejections + 1
    }
  }
  
  # Type I エラー率の計算
  student_t_error_rate <- student_t_rejections / num_simulations
  welch_t_error_rate <- welch_t_rejections / num_simulations
  
  return(list(student_t = student_t_error_rate, welch_t = welch_t_error_rate))
}

# サンプルサイズ、平均、標準偏差の組み合わせ
conditions <- list(
  # 等分散、等サンプルサイズ
  list(n1 = 50, n2 = 50, mean1 = 0, mean2 = 0, sd1 = 1, sd2 = 1),
  list(n1 = 100, n2 = 100, mean1 = 0, mean2 = 0, sd1 = 1, sd2 = 1),
  list(n1 = 500, n2 = 500, mean1 = 0, mean2 = 0, sd1 = 1, sd2 = 1),
  
  # 等分散、不等サンプルサイズ
  list(n1 = 50, n2 = 100, mean1 = 0, mean2 = 0, sd1 = 1, sd2 = 1),
  list(n1 = 50, n2 = 500, mean1 = 0, mean2 = 0, sd1 = 1, sd2 = 1),
  list(n1 = 100, n2 = 500, mean1 = 0, mean2 = 0, sd1 = 1, sd2 = 1),
  
  # 不等分散、等サンプルサイズ
  list(n1 = 100, n2 = 100, mean1 = 0, mean2 = 0, sd1 = 1, sd2 = 2),
  list(n1 = 100, n2 = 100, mean1 = 0, mean2 = 0, sd1 = 1, sd2 = 3),
  list(n1 = 100, n2 = 100, mean1 = 0, mean2 = 0, sd1 = 1, sd2 = 4),
  
  # 不等分散、不等サンプルサイズ（小サンプル、大分散）
  list(n1 = 50, n2 = 100, mean1 = 0, mean2 = 0, sd1 = 2, sd2 = 1),
  list(n1 = 50, n2 = 500, mean1 = 0, mean2 = 0, sd1 = 2, sd2 = 1),
  list(n1 = 100, n2 = 500, mean1 = 0, mean2 = 0, sd1 = 2, sd2 = 1),
  
  # 不等分散、不等サンプルサイズ（大サンプル、大分散）
  list(n1 = 100, n2 = 50, mean1 = 0, mean2 = 0, sd1 = 1, sd2 = 2),
  list(n1 = 500, n2 = 50, mean1 = 0, mean2 = 0, sd1 = 1, sd2 = 2),
  list(n1 = 500, n2 = 100, mean1 = 0, mean2 = 0, sd1 = 1, sd2 = 2)
)

# シミュレーションの実行
results <- lapply(conditions, function(cond) {
  result <- simulate_t_tests(cond$n1, cond$n2, cond$mean1, cond$mean2, cond$sd1, cond$sd2)
  c(n1 = cond$n1, n2 = cond$n2, mean1 = cond$mean1, mean2 = cond$mean2, 
    sd1 = cond$sd1, sd2 = cond$sd2, student_t = result$student_t, welch_t = result$welch_t)
})

# 結果の表示
results_df <- do.call(rbind, results)
print(results_df)

# 結果のプロット
plot_data <- as.data.frame(results_df)
plot_data$condition <- paste("n1 =", plot_data$n1, ", n2 =", plot_data$n2, 
                             "\nsd1 =", plot_data$sd1, ", sd2 =", plot_data$sd2)

# 条件に番号を付加
plot_data$condition_num <- 1:nrow(plot_data)
plot_data$condition <- paste(plot_data$condition_num, plot_data$condition, sep=". ")

plot_data_long <- pivot_longer(plot_data, 
                               cols = c("student_t", "welch_t"),
                               names_to = "Test", 
                               values_to = "Type_I_Error_Rate")

# factor型に変換し、レベルを設定することで順序を制御
plot_data_long$condition <- factor(plot_data_long$condition, 
                                   levels = unique(plot_data_long$condition))

ggplot(plot_data_long, aes(x = condition, y = Type_I_Error_Rate, fill = Test)) +
  geom_bar(stat = "identity", position = position_dodge(width = 0.9), width = 0.8) +
  geom_hline(yintercept = 0.05, linetype = "dashed", color = "red") +
  labs(title = "Type I Error Rates: Student's t-test vs Welch's t-test",
       subtitle = "With varying sample sizes and standard deviations",
       x = "Conditions",
       y = "Type I Error Rate",
       fill = "Test Type") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  coord_cartesian(ylim = c(0, 0.3)) +  # y軸の範囲を0から0.3に設定
  scale_x_discrete(labels = function(x) gsub("^\\d+\\. ", "", x))  # 番号を除去してラベルを簡潔に

結果

グラフは冒頭に示してあります。

       n1  n2 mean1 mean2 sd1 sd2 student_t welch_t
[1,]  50  50     0     0   1   1    0.0430  0.0430
[2,] 100 100     0     0   1   1    0.0490  0.0490
[3,] 500 500     0     0   1   1    0.0555  0.0555
[4,]  50 100     0     0   1   1    0.0575  0.0555
[5,]  50 500     0     0   1   1    0.0480  0.0485
[6,] 100 500     0     0   1   1    0.0505  0.0510
[7,] 100 100     0     0   1   2    0.0465  0.0460
[8,] 100 100     0     0   1   3    0.0475  0.0465
[9,] 100 100     0     0   1   4    0.0515  0.0500
[10,]  50 100     0     0   2   1    0.1095  0.0470
[11,]  50 500     0     0   2   1    0.2655  0.0485
[12,] 100 500     0     0   2   1    0.1905  0.0465
[13,] 100  50     0     0   1   2    0.1030  0.0470
[14,] 500  50     0     0   1   2    0.2360  0.0515
[15,] 500 100     0     0   1   2    0.1935  0.0395

等分散でなく、n1とn2が異なる条件では、Student t検定は第一種の過誤を大幅に増加させていることがわかります。

有意差がないのにあると結論してしまうことは、特に医療統計においては重大な結果を引き起こしかねないわけで、ある程度保守的であるほうが医療者の態度としては正しいと思う (個人的見解)。

シミュレーションの結果をもとに、Student t検定と Welch t検定でほとんど違わないという主張を見たことがあったが...

ネットで医療者向けの教育資料として公開されていて、信じていましたが、どうやら間違っていたようです。

なぜそうなった？

見直してみると、すべて n1 = n2 でやってありました。
上記の今回の結果も、n1 = n2 だとどちらの検定もほぼ正確でしたから、一致してます。

教訓

他人の意見を鵜呑みにすることなかれ。たとえシミュレーションがやってあったとしてもだ。

等分散の仮定が成り立たないことが多い現実では、Welch t検定の方がより安全で信頼できる選択肢です！