Rでウェルチのt検定(Welch’s t test)

Rで ウェルチのt検定(Welch’s t test) を試みます。

始めに平均値も分散も異なるAとBの2つのサンプルを作成します。

seed <- 20250410
set.seed(seed)

n_a <- 100
n_b <- 200

# サンプルA: 平均50, 標準偏差10 の正規分布に従うデータ
sample_a <- rnorm(n_a, mean = 50, sd = 10)

# サンプルB: 平均55, 標準偏差15 の正規分布に従うデータ
sample_b <- rnorm(n_b, mean = 55, sd = 15)

# 生成したデータの基本統計量を表示
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_a <- mean(sample_a)
mean_b <- mean(sample_b)
var_a <- var(sample_a) # サンプルAの不偏分散(n-1)
var_b <- var(sample_b) # サンプルBの不偏分散(n-1)

# Welchのt統計量の計算
# t = (mean_a - mean_b) / sqrt(var_a/n_a + var_b/n_b)
t_numerator <- mean_a - mean_b
t_denominator <- sqrt(var_a / n_a + var_b / n_b)
welch_t_stat <- t_numerator / t_denominator

# 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) }
df_numerator <- (var_a / n_a + var_b / n_b)^2
df_denominator_term1 <- (var_a / n_a)^2 / (n_a - 1)
df_denominator_term2 <- (var_b / n_b)^2 / (n_b - 1)
welch_df <- df_numerator / (df_denominator_term1 + df_denominator_term2)

# p値の計算 (両側検定)
# t分布の累積分布関数pt()を使用
# P(T > |t|) を計算し、2倍する
p_value <- 2 * pt(abs(welch_t_stat), df = welch_df, lower.tail = FALSE)

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 

以上です。