t 検定

1標本の平均の検定

確率変数 X の分布が正規分布 N(μ, σ2) であり，そこから n 個取り出した標本の平均値を

\[ \bar{X} = \frac{X_1 + X_2 + \cdots + X_n}{n} \]

とし，標本分散を

\[ s^2 = \frac{1}{n-1} \sum_{i=1}^n (X_i - \bar{X})^2 \]

とすると，

\[ t = \frac{\bar{X} - \mu}{\sqrt{s^2 \! / n}} \]

は，σ2 の値にかかわらず，自由度 n - 1 の t 分布になります。このことを使って，母集団が平均 μ の正規分布かどうかを検定できます。

実際には X が正規分布でなくても，n が大きければ中心極限定理により X は正規分布に近づくので，この検定は母集団が正規分布かどうかには鈍感です。母平均が μ かどうかを調べるだけの検定と考えて大きな間違いはありません。

1標本の平均の検定の例

10人の患者にある睡眠薬を飲ませたところ，睡眠時間がそれぞれ次の時間だけ増えました (Arthur R. Cushny and A. Roy Peebles, The Journal of Physiology 32, 501-510 (1905)):

1.9, 0.8, 1.1, 0.1, -0.1, 4.4, 5.5, 1.6, 4.6, 3.4

> x = c(1.9, 0.8, 1.1, 0.1, -0.1, 4.4, 5.5, 1.6, 4.6, 3.4)

検定に走る前に，まずはグラフを描いてみるべきです。どんなグラフがいいでしょうか。一例です：

> stripchart(x, pch=16, method="stack")  # または "jitter"

母集団が平均 0（の正規分布）であるという帰無仮説をRで検定してみましょう。

> t.test(x)

	One Sample t-test

data:  x 
t = 3.6799, df = 9, p-value = 0.005076
alternative hypothesis: true mean is not equal to 0 
95 percent confidence interval:
 0.8976775 3.7623225 
sample estimates:
mean of x 
     2.33

つまり，t 値の絶対値がこの 3.6799 以上になる確率（p 値）は 0.005076 と小さく，データは帰無仮説を支持していないようです。つまり，睡眠薬の効果はあったと言えます。

この検定は正規分布かどうかには鈍感ですが，まったく正規分布を仮定しない検定法として符号検定とWilcoxonの符号付き順位検定がありますので，正規分布でない場合はそちらを積極的に使いましょう。

2標本の平均の差の検定

確率変数 X から引き出した m 個の値 X1, X2, ..., Xm と，確率変数 Y から引き出した n 個の値 Y1, Y2, ..., Yn があったとします。どちらの分布も正規分布と仮定したとき，それらの母平均が等しいという仮説を検定してみましょう。

それぞれの標本平均を X, Y, 標本分散を sX2, sY2 とします。

分散が等しいと仮定できる場合

二つの標本をプールしたものから次の式で分散の推定値 s2 を求めます：

\begin{align*} s^2 &= \frac{(X_1 - \bar{X})^2 + \cdots + (X_m - \bar{X})^2 + (Y_1 - \bar{Y})^2 + \cdots + (Y_n - \bar{Y})^2}{m + n - 2} \\ &= \frac{(m - 1) s_X^2 + (n - 1) s_Y^2}{m + n - 2} \end{align*}

このとき，

\[ t = \frac{\bar{X} - \bar{Y}}{\sqrt{s^2 ( \frac{1}{m} + \frac{1}{n} ) }} \]

は自由度 m + n - 2 の t 分布に従います。

分散が等しいと仮定できない場合

一般には，X の母分散も Y の母分散も不明であるなら，それらが等しいかどうかも不明のはずです。この場合によく用いられる方法は，

\[ t = \frac{\bar{X} - \bar{Y}}{\sqrt{s_X^2 \big/ m + s_Y^2 \big/ n}} \]

が近似的に自由度

\[ \nu = \frac{(s_X^2 \big/ m + s_Y^2 \big/ n)^2}{\dfrac{(s_X^2 \big/ m)^2}{m-1} + \dfrac{(s_Y^2 \big/ n)^2}{n-1}} \]

の t 分布に従うことを使うものです。これはWelch（ウェルチ）が次の論文で提案した方法の一つで，広く用いられています。

この自由度 ν は整数になりませんが，t 分布は ν が整数でない場合にも拡張できます。

2標本の平均の差の検定の例

2組の患者たちに痛みのレベルを報告してもらったところ，Rの書き方で次のような結果を得ました (T. Lumley, Biometrics 52, 354-361 (1996))。

x = c(1,2,1,1,1,1,1,1,1,1,2,4,1,1)
y = c(3,3,4,3,1,2,3,1,1,5,4)

この平均値の違いをRで t 検定してみます。

> t.test(x,y)

	Welch Two Sample t-test

data:  x and y 
t = -2.9486, df = 15.916, p-value = 0.009479
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:
 -2.3556088 -0.3846509 
sample estimates:
mean of x mean of y 
 1.357143  2.727273

これがRの t.test() のデフォルトのWelchの t 検定で，両群の分散が等しいことを仮定しないものです。分散が等しいことを仮定するものは，Rでは次のようにして計算します。

> t.test(x,y,var.equal=TRUE)

	Two Sample t-test

data:  x and y 
t = -3.1158, df = 23, p-value = 0.004862
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:
 -2.2797891 -0.4604706 
sample estimates:
mean of x mean of y 
 1.357143  2.727273

これも正規分布かどうかには鈍感ですが，まったく正規分布を仮定しない検定法としてWilcoxon-Mann-Whitney検定（分布が同じ場合）やBrunner-Munzel検定（分布が異なる場合）がありますので，それらを積極的に使いましょう。

平均と標準偏差を与えられた場合の検定

Rの t.test() はデータそのものを与えますが，個数1，平均1，標準偏差1，個数2，平均2，標準偏差2を与えて等分散・非等分散の t 検定をする関数は次のように書くことができます。標準偏差は n - 1 で割る方式であることに注意してください。

ttest = function(n1, x1, s1, n2, x2, s2) {
  n = n1 + n2 - 2
  u = ((n1 - 1) * s1^2 + (n2 - 1) * s2^2) / n
  t = (x1 - x2) / sqrt(u / n1 + u / n2)
  r = cat("Equal variance:\n\t", sep="")
  r = cat(r, "t = ", t, ", df = ", n, ", p = ", 2 * pt(-abs(t), n), "\n", sep="")
  t = (x1 - x2) / sqrt(s1^2 / n1 + s2^2 / n2)
  n = (s1^2 / n1 + s2^2 / n2)^2 / ((s1^2 / n1)^2 / (n1-1) + (s2^2 / n2)^2 / (n2-1))
  r = cat(r, "Unequal variance:\n\t", sep="")
  cat(r, "t = ", t, ", df = ", n, ", p = ", 2 * pt(-abs(t), n), "\n", sep="")
}

分散が等しいかどうか不明の場合

分散が等しいかどうかわからないときは，まず等分散かどうかの検定をして，それによって振り分ける，という方法が過去において薦められたことがありました。これは論理的にまずい方法です。まずいことはシミュレーションを使っても確かめられます。

ni = 1000000                     # 100万回のシミュレーション
p = numeric(ni)                  # p[1]..p[ni] というni個の変数を作る
for (i in 1:ni) {                # iが1からniまでループ
  x = rnorm(10, mean=0, sd=1.5)  # 平均0，標準偏差1.5の正規乱数を10個作る
  y = rnorm(30, mean=0, sd=1.0)  # 平均0，標準偏差1.0の正規乱数を30個作る
  p[i] = t.test(x, y)$p.value    # Welchのt検定によるp値をp[i]に代入
}
mean(p <= 0.05)                  # p≦0.05となる割合を出力
mean(p <= 0.01)                  # p≦0.01となる割合を出力

上を実行すると，論理的に辻褄が合っていれば，p ≦ 0.05 となる割合は 0.05 と出力されるはずです。同様に，t.test(x, y) を t.test(x, y, var.equal=TRUE) とすれば，等分散を仮定する方法についての p 値が求められます。

等分散かどうか検定してから振り分ける方法は，t.test() の行を次の2行にします。

  vp = var.test(x, y)$p.value
  p[i] = t.test(x, y, var.equal=(vp > 0.05))$p.value

論理的に辻褄が合っていれば出力される値と実際に出力された値は次のようになりました。2段階の方法については，等分散かどうかの判断は水準0.05で行う場合と0.2で行う場合について調べました（それ以外については下に書きます）。

 0.050.01
等分散を仮定0.1074690.033762
2段階，0.050.0801980.024796
2段階，0.20.0642140.018550
Welchだけ0.0515150.011337
順位をWelch0.0549440.014522
WMW0.0785910.019815
Brunner-Munzel0.0569380.017402
並べ替えBrunner-Munzel0.0510.012

このように，実際に等分散でない場合は，等分散かどうかの検定を使わず，Welchの方法だけで行うのが正しい方法です。

ついでに順位についてのWelchの t 検定や，WMW検定，Brunner-Munzel検定も試してみました。 WMW検定も等分散でない場合に甘くなります。Brunner-Munzel検定は小さい p 値については少し甘くなるようです。

最後の並べ替えBrunner-Munzel検定については，あまりたくさん試せませんでしたが，下のようなプログラムでやってみました。

library(lawstat)
ni = 1000
nj = 10000
p = numeric(ni)
q = numeric(nj)
for (i in 1:ni) {
  x = rnorm(10, mean=0, sd=1.5)
  y = rnorm(30, mean=0, sd=1.0)
  z = c(x, y)
  d = abs(brunner.munzel.test(x,y)$statistic)
  for (j in 1:nj) {
    s = sample(1:40, 10)
    e = abs(brunner.munzel.test(z[s],z[-s])$statistic)
    q[j] = (e > d);
  }
  p[i] = mean(q)
}
mean(p <= 0.05)
mean(p <= 0.01)

対応がある場合の検定

同じ n 人について X1, X2, ..., Xn と Y1, Y2, ..., Yn を測定した場合，2標本の平均の差を検定するのではなく，差 di = Xi - Yi について1標本の平均の検定をします。どうしてそうするか考えてみましょう。

Rでは t.test(X-Y) としてもいいのですが t.test(X, Y, paired=TRUE) のようにすることもできます。

奥村晴彦

Last modified: 2012-01-26 08:29:59