小標本問題と t検定 - ほくそ笑む

統計を学び始めると「t検定」というのが最初のほうで出てくると思います。
t検定は、20世紀前半に活躍した統計学者、ウィリアム・ゴセットによって「小標本問題」というのを解決するために考案されました。
小標本問題とは、正規分布の平均値の検定に正規分布を用いると、サンプルサイズが小さい場合にαエラーを過小評価してしまうという問題です。
今日はこの小標本問題とそれを解決する t検定について R によるシミュレーションを使って説明してみたいと思います。

正規分布の平均値の検定

確率変数 $X$ が正規分布に従うとき、その平均値もまた、正規分布に従います。
数式で書くと、
　 $X \sim {\rm Norm}(\mu, \sigma^2)$
　 $\overline{X} \sim {\rm Norm}(\mu, \sigma^2/n)$
となります。（分散が $1/n$ されていることに注意）
なので、正規分布の平均値の検定には正規分布を使用すれば良いように思われます。
これを R でシミュレートしてみましょう。

# 正規分布を使用して平均値が 0 と等しいかの p値を求める
norm.dist.p <- function(x) {
  n <- length(x)
  mean <- mean(x)
  sd <- sd(x) / sqrt(n)
  p <- pnorm(-abs(mean), mean=0, sd=sd) * 2
  p
}
x <- rnorm(10, mean=0)
p <- norm.dist.p(x)
cat("p =", p, "\n")

p = 0.974622

このようにして、p値を求めることができます。
上では、有意差の無い場合、つまり平均値 0 の正規分布から得られる x に対してこの検定を行っているので、p は十分大きな値になります。
ただし、x はランダムに取られるので、何度かやってみると、偶然による偏りが生じてしまい、p値が小さくなることがあります。
これが、αエラーです。
αエラーは、検定を何度も繰り返したときに、本当は有意差の無いものが有意差ありと判断されてしまうことを指しています。
統計的検定では、このαエラーが生じる確率を有意水準αによって制御します。
有意水準αを 0.05 と決め、有意差の無い場合に対して検定を何度も繰り返すと、αエラーが生じる確率は 0.05 程度に制御されるはずです。
これをシミュレーションして確かめてみましょう。

# 有意差の無いランダムデータに対して検定を繰り返す
alpha = 0.05
ps <- sapply(1:10000, function(i) {
  x <- rnorm(10)
  p <- norm.dist.p(x)
  p
})
fp <- sum(ps < alpha) / length(ps)
cat("alpha error rate =", fp, "\n")

alpha error rate = 0.0812

なんと、αエラー率が 0.08 程度になってしまいました。
α=0.05 としたのに、実際は 0.08 ということは、αエラーが低く見積もられていたということになります。
これは、有意でないものを有意とする危険性が想定より高くなることを示しています。これはいけません。
この傾向は、サンプルサイズ(標本サイズ)が小さい場合に顕著にみられるため、小標本問題と呼ばれます。

小標本問題はなぜ起こるのか？

小標本問題はなぜ起こるのでしょうか？
これは、正規分布の平均値の検定において、分散が推定値にすぎないということに問題があります。
上で示した式
　 $\overline{X} \sim {\rm Norm}(\mu, \sigma^2/n)$
において、分散 $\sigma^2/n$ は、真の値が分かっていません。
したがって、 $\sigma$ の推定値である s = sd(x) が使用されています。
しかし、s 自体にもぶれがあり、x の偏りによって大きくなったり小さくなったりします。

sds <- sapply(1:10000, function(i) {
  x <- rnorm(3, mean=0, sd=1)
  sd(x)
})
hist(sds)
quantile(sds)

       0%       25%       50%       75%      100% 
0.0048260 0.5302899 0.8344423 1.1935087 3.1106698

このように、ヒストグラムが左側に偏っているのと、中央値(50%点)が 0.83 であり、真値である 1 よりも小さい値が出やすいということがわかります。
分散が小さく見積もられやすいということは、それだけ検定が有意に出やすくなるということであり、これが小標本問題の原因となるのです。

小標本問題の解決

この小標本問題を解決するために、ウィリアム・ゴセットは t統計量というのを考えました。*1
　 $t = \frac{\overline{X} - \mu}{s/ \sqrt{n}}$
この t統計量が従う、t分布という確率分布を発明したのです。*2
それでは、この t分布を使って、上と同じく正規分布の平均値の差の検定をシミュレーションしてみましょう。

# t分布を使用して平均値が 0 と等しいかの p値を求める
t.dist.p <- function(x) {
  n <- length(x)
  mean <- mean(x)
  sd <- sd(x) / sqrt(n)
  t <- mean / sd
  p <- pt(-abs(t), df=n-1) * 2
  p
}
# 有意差の無いランダムデータに対して検定を繰り返す
alpha = 0.05
ps <- sapply(1:10000, function(i) {
  x <- rnorm(10)
  p <- t.dist.p(x)
  p
})
fp <- sum(ps < alpha) / length(ps)
cat("alpha error rate =", fp, "\n")

alpha error rate = 0.049

おおー、αエラー率がちゃんと 0.05 近くになりましたね。
ウィリアム・ゴセットのこの仕事によって、正規分布の平均値の差の検定における小標本問題は見事に解決したのでした。

比較

これだけでは物足りないので、小標本に対する正規分布での検定と t分布での検定で、αエラー率がどのように異なるのかシミュレーションしてみましょう。
コードと結果を下記に示します。

calc <- function(sample.nums, calc.p.func, rep.num=100, num.of.ps=rep.num, num.of.fps=rep.num, alpha=0.05) {
  data.frame(t(sapply(sample.nums, function(sample.num) {
    fps <- sapply(1:num.of.fps, function(i) {
      ps <- sapply(1:num.of.ps, function(i) {
        x <- rnorm(sample.num)
        p <- calc.p.func(x)
        p
      })
      fp <- sum(ps < alpha) / length(ps)
      fp
    })
    c(sample.num = sample.num, alpha.error = mean(fps), sd = sd(fps))
  })))
}

sample.nums <- seq(2, 15, by=1)
# sample.nums <- seq(10, 100, by=10)
rep.num <- 250
res.norm <- calc(sample.nums, norm.dist.p, rep.num)
res.t <- calc(sample.nums, t.dist.p, rep.num)
res <- rbind(cbind(res.norm, Dist.="Normal dist."), cbind(res.t, Dist.="Student's t"))

library(ggplot2)
limits <- aes(ymax = res$alpha.error + res$sd, ymin = res$alpha.error - res$sd)
ggplot(res, aes(y=alpha.error, x=sample.num, colour=Dist.)) + ylim(0, 0.35) +
  geom_line(aes(group=Dist.), size=2) + geom_errorbar(limits, alpha=0.3, size=1) +
  ylab("Alpha Error Rate") + xlab("Sample Num")

小さいサンプルサイズ(2〜15)に対してαエラー率をプロットしてみました。

n = 2 の場合、正規分布を使った検定はα=0.05 と設定しているにもかかわらず実際は 0.3 程度と非常に大きくなります。
つまり、設定されたαは実際よりも小さく見積もられているということです。
この差はサンプルサイズが増えるごとに縮まっていきます。

対して t検定のほうはα=0.05 程度で安定しています。素晴らしい！

サンプルサイズ 10〜100 をプロットしてみたものがこちらです。
n = 100 の場合、正規分布を使った検定ではα=0.05 の設定に対して実際は 0.051 程度となり、あまり差はありません。
つまり、正規分布を使った検定でも、サンプルサイズが十分大きければ使えるということです。*3
逆に言うと、t検定が発明されなければ、サンプルサイズ100程度は無いと検定できなかったのかなーとか考えると、ゴセットの偉大さがわかります。
n = 3 でひいひい言ってる生物学者はイギリスに足向けて寝れませんね。