Rで ウェルチのt検定(Welch’s t test) を試みます。
始めに平均値も分散も異なるAとBの2つのサンプルを作成します。
<- 20250410
seed set.seed(seed)
<- 100
n_a <- 200
n_b
# サンプルA: 平均50, 標準偏差10 の正規分布に従うデータ
<- rnorm(n_a, mean = 50, sd = 10)
sample_a
# サンプルB: 平均55, 標準偏差15 の正規分布に従うデータ
<- rnorm(n_b, mean = 55, sd = 15)
sample_b
# 生成したデータの基本統計量を表示
summary(sample_a)
summary(sample_b)
sd(sample_a)
sd(sample_b)
Min. 1st Qu. Median Mean 3rd Qu. Max.
23.46 44.41 50.41 49.77 57.01 70.20
Min. 1st Qu. Median Mean 3rd Qu. Max.
15.39 44.40 54.28 55.64 67.21 96.28
[1] 9.906609
[1] 15.45127
両サンプルの平均値が有意に異なっているか否かを検定(両側検定)します。
Welchのt検定 の検定統計量と自由度は以下の通りです。
検定統計量
\[ t = \dfrac{\bar{A} - \bar{B}}{\sqrt{\dfrac{S^2_A}{N_A} +\dfrac{S^2_B}{N_B}}} \]
通常のt検定と異なり Pooled Variance(\(S^2_p\)) は利用していません。
\[\begin{eqnarray}S_p^2 &=& \dfrac{ (N_A - 1)S_A^2 + (N_B - 1)S_B^2 }{ N_A + N_B - 2}\\t&=&\dfrac{\bar{A} - \bar{B}}{\sqrt{ S_p^2 \left(\dfrac{1}{N_A} + \dfrac{1}{N_B}\right) }} \end{eqnarray}\]
Welch-Satterthwaiteの自由度
\[df = \dfrac{\left(\dfrac{S^2_A}{N_A}+\dfrac{S^2_B}{N_B}\right)^2}{ {\dfrac{\left(\dfrac{S^2_A}{N_A}\right)^2}{N_A-1} + \dfrac{\left(\dfrac{S^2_B}{N_B}\right)^2}{N_B-1} }} \]
通常のt検定の単純な合計とは異なります。
\[df=N_A + N_B - 2\]
# 必要な統計量を計算
<- mean(sample_a)
mean_a <- mean(sample_b)
mean_b <- var(sample_a) # サンプルAの不偏分散(n-1)
var_a <- var(sample_b) # サンプルBの不偏分散(n-1)
var_b
# Welchのt統計量の計算
# t = (mean_a - mean_b) / sqrt(var_a/n_a + var_b/n_b)
<- mean_a - mean_b
t_numerator <- sqrt(var_a / n_a + var_b / n_b)
t_denominator <- t_numerator / t_denominator
welch_t_stat
# Welch-Satterthwaiteの自由度の計算
# df = (var_a/n_a + var_b/n_b)^2 / { (var_a/n_a)^2/(n_a-1) + (var_b/n_b)^2/(n_b-1) }
<- (var_a / n_a + var_b / n_b)^2
df_numerator <- (var_a / n_a)^2 / (n_a - 1)
df_denominator_term1 <- (var_b / n_b)^2 / (n_b - 1)
df_denominator_term2 <- df_numerator / (df_denominator_term1 + df_denominator_term2)
welch_df
# p値の計算 (両側検定)
# t分布の累積分布関数pt()を使用
# P(T > |t|) を計算し、2倍する
<- 2 * pt(abs(welch_t_stat), df = welch_df, lower.tail = FALSE)
p_value
list(welch_t_stat = welch_t_stat, welch_df = welch_df, p_value = p_value)
$welch_t_stat
[1] -3.974897
$welch_df
[1] 280.1242
$p_value
[1] 8.965141e-05
よって、有意水準5%でサンプルAとサンプルBの平均値には統計的に有意な差があると言えます。
Rの関数t.test {stats}による結果を確認します。
t.test(sample_a, sample_b, var.equal = FALSE, paired = F, alternative = "two.sided")
Welch Two Sample t-test
data: sample_a and sample_b
t = -3.9749, df = 280.12, p-value = 8.965e-05
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-8.765444 -2.959138
sample estimates:
mean of x mean of y
49.77284 55.63513
以上です。