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 だとどちらの検定もほぼ正確でしたから、一致してます。
教訓
他人の意見を鵜呑みにすることなかれ。たとえシミュレーションがやってあったとしてもだ。