t分布とt検定
統計的仮説検定の考え方
統計的仮説検定の考え方
生物統計学を履修すると、真っ先に習うのが、（Studentの）t検定と呼ばれる仮説検定法だ。t検定の「t」とはこの検定で用いる検定統計量がtと呼ばれていることに由来する。また、ある仮説を設定したときに検定統計量tが従う分布をt分布と呼んでいる。基本的な統計的仮説検定の流れは以下のとおりである。図中の検定の流れ（1-i、2など）と以降の言葉の説明は対応させてある。まずは、以下の流れを押さえよう。本ページの最後で、なぜこのような考え方をしているのか、解説する。
1. 仮説の設定
仮説はいくつか設定する必要がある。統計的仮説検定の中で最も重要であるといって過言ではなく、また、検定を行う前に設定しておく必要がある。
i. 検証したい命題の設定
「ある標本の平均値は0と等しいか」「2つの標本の平均値は等しいか」といった問題の設定。統計的仮説検定において設定される命題は「△△の平均値は○○と等しい」や「○○と△△の平均値は等しい」というような「差がない」仮説である。この差がない仮説を帰無仮説と呼ぶ。この差がない仮説を棄却（否定）することで、「差がある」ことを言うのが統計的仮説検定である。このように、「差がある」ことを言うために棄却してしまう仮説なので、「無に帰する」仮説、帰無仮説と呼ばれる。
ii. 標本が従う分布（母集団）の仮定
標本が従う分布は正規分布なのか、ポアソン分布なのか……といった設定。帰無仮説と標本の従う分布を設定すると、帰無仮説の下での検定統計量の分布を計算できる。この検定統計量の分布を帰無分布と呼ぶ。
iii. 有意水準（棄却域）の設定
標本から得られた検定統計量が、帰無分布の端の何％を占める領域に存在すれば、差があるとみなすか、の設定。何%の部分を有意水準と呼ぶ。伝統的に5%とされることが多い。
2. 標本をもとに検定統計量の計算
検定方法に従って、検定統計量を計算。例えば、t検定なら統計量tを算出するし、F検定なら統計量Fを算出する。
3. 検定統計量と帰無分布との比較・判断
標本から得られた検定統計量が帰無分布のどの場所に位置するか確認。棄却域に検定統計量が存在する場合は帰無仮説を棄却する。そうでなければ、帰無仮説を保留とする。検定統計量より極端な値（帰無分布の左側なら検定統計量以下の値、右側なら以上の値）をとるときの帰無分布に対する割合をp値と呼ぶ。実際の検定ではこのp値と有意水準を比較して、棄却するかを判断する。
t検定の考え方
上記の統計的仮説検定の流れを2標本のt検定を例に考えてみよう。
1. 仮説の設定
i. 検証したい命題の設定
「2つの標本（AおよびB）の平均値は等しい」=帰無仮説
ii. 標本が従う分布（母集団）の仮定
2つの標本が従う分布は正規分布=正規分布から2つの標本は生成したと考える
⇒iとiiを合わせて、平均値が等しい正規分布=単一の正規分布から2標本が生成したと考える
正規分布から2標本が生成され、その2標本を使って検定統計量tが計算される。この過程を無限回繰り返すことで、帰無仮説の下での検定統計量のt分布が決定する（帰無分布）。tの計算は以下のとおりである。定義を見ればわかるが、tは2標本の平均値の差が大きく、分散が小さいほど、大きな値をとる。
iii. 有意水準（棄却域）の設定
両側検定の場合、t分布（帰無分布）の面積のうち、両端から合計5％に相当する部分を棄却域とする。（片側検定の考え方については省略）
2. 標本をもとに検定統計量の計算
帰無分布の手続きと同様に、自分の手元の標本について検定統計量t（t0とする）を計算する。
3. 検定統計量と帰無分布との比較・判断
t0がt分布のどの位置にいるかを確認する。棄却域にt0が存在する場合、帰無仮説を棄却する。t0が5%の棄却域に存在すれば、p値は5%を下回っている。
t分布をシミュレートする
さて、ここからはRを使ってシミュレートし、もっと具体的にt検定について考えていく。今、平均=50、標準偏差=20の正規分布があるとする。正規分布は無限の数の集合であるが、そこから十分な個数、例えば10000個数字を取り出しても、ほとんどその形は変わらない。Rでは以下のようになる。
------------------------------------------------------
library(plotn)
#このライブラリは私オリジナルの作図ライブラリである。
#もし、plotnライブラリが欲しい場合は、
#以下のようにコマンドを打つ（#は消す）。devtoolsライブラリが必要。
#library(devtools)
#install_github("bugplant/plotn")
n <- rnorm(10000, mean = 50, sd = 20) #平均50、標準偏差20の正規分布に従う乱数を10000個生成
histn(n, xlab = "value", ylab = "頻度") #図1の作図。シミュレーションなので乱数を生成するごとに少しずつ図は異なってくる。hist(n)でも可。
------------------------------------------------------
図1 平均50、標準偏差20の正規分布
次に同じ正規分布からサンプルサイズ20と30の標本が生成されたと考える。同じ分布から標本が生成されているので、これら2つの標本の平均値はともに50近くになるはずだ。併せて、検定統計量tも計算してみよう。
------------------------------------------------------
m1 <- rnorm(20, mean = 50, sd = 20) #上記の正規分布から20個サンプリング
m2 <- rnorm(30, mean = 50, sd = 20) #上記の正規分布から30個サンプリング
m1 #m1の中身、シミュレーションするごとに中身は変わる
## [1] 55.984449 -6.190775 59.108486 64.134176 60.107704 36.713796 22.630454
## [8] 42.320423 52.978341 44.691751 38.597839 40.694713 85.325101 44.382624
## [15] 19.621492 60.551150 45.246971 64.784980 37.893546 59.837865
m2 #m2の中身、シミュレーションするごとに中身は変わる
## [1] 50.96129 42.62763 50.38982 74.22224 52.55635 44.02777 53.30684 64.36510
## [9] 52.23433 48.19708 44.53958 36.08335 28.97115 39.16197 70.57005 31.77959
## [17] 44.70954 62.20962 70.65534 43.73504 50.10141 40.53722 61.63951 48.69854
## [25] 12.90490 60.16799 50.10335 39.61425 48.63319 59.75502
mean(m1) #m1の平均、シミュレーションするごとに変わる
## [1] 46.47075
mean(m2) #m2の平均、シミュレーションするごとに変わる
## [1] 49.24863
#多くの場合は、m1とm2の平均は近くなる。しかし、かけ離れた値となることもある。
ue <- ((20 - 1)*(sd(m1)^2) + (30 - 1)*(sd(m2)^2))/(20 + 30 - 2) #2個の分散を統合(合併不偏分散)
t <- (mean(m1) - mean(m2))/sqrt(ue*(1/20 + 1/30)) #検定統計量t
t #tの中身
## [1] -0.5999021
------------------------------------------------------
試してみるとわかるが、m1とm2の平均は元の正規分布の平均50近くになり、近い値を示すことが多い。しかし、全くかけ離れた値となることもある。今回の場合は、t = -0.60だ。では、この正規分布から標本を生成し、tを計算することを、10000回繰り返したらどうなるだろう。2つの標本ペアが10000個でき、それに応じたtが10000個できる。これによって、t分布をシミュレーションで作ることができる。
------------------------------------------------------
t <- NULL
for (i in 1:10000){ #定義に倣い統計量tを10000個算出
m1 <- rnorm(20, mean = 50, sd = 20) #上記の正規分布から20個サンプリング
m2 <- rnorm(30, mean = 50, sd = 20) #上記の正規分布から30個サンプリング
ue <- ((20 - 1)*(sd(m1)^2) + (30 - 1)*(sd(m2)^2))/(20 + 30 - 2) #2個の分散を統合
t <- c(t, (mean(m1) - mean(m2))/sqrt(ue*(1/20 + 1/30)))
#統計量tの算出、分子を見るとわかるがこれは2つの標本の平均値の差に依存した量
}
head(t) #10000個のシミュレートした統計量tを一部表示
## [1] -0.3944782 -2.1391604 0.6972603 -0.1572087 -1.7567524 0.1391997
histn(t, xlab = "統計量t", ylab = "頻度") #統計量tをヒストグラムにして図示 = t分布、hist(t)でも可。（図2）
overdraw(abline(v = qt(df = 48, p = 0.975)),
abline(v = qt(df = 48, p = 0.025)), #95%区間の表示
axis(side = 1, at = c(qt(df = 48, p = 0.025), qt(df = 48, p = 0.975)))
)
------------------------------------------------------
図2 平均50、標準偏差20の正規分布からサンプルサイズ20および30の標本が生成されたときのt分布。-2.01～2.01の範囲に、約95%のtが含まれる。
t分布は0を中心とした釣り鐘型となる。すなわち、tは0付近の値を示すことが多いということだ。これは、2標本は同じ正規分布から生じているので、平均値が近くなりやすい、ということを示している。tの定義には分子に平均値の差があるため、2標本の平均値が近いと、tは0に近い値をとる。
帰無仮説を棄却するとはどういう状況か？
帰無分布にあたるt分布をシミュレーションで作成できるところまできた。次に考えるのは、「帰無仮説が棄却される状況」だ。図2には10000個のtのうち、両端から約5%=約500個のtが含まれるであろう境界線も示している。実際に、含まれているか計算してみよう。
------------------------------------------------------
#95%の境界線にあたるtの計算は以下のコマンドで行う。
#dfは自由度と呼ばれる値で、2標本のt検定なら、（2標本のサンプルサイズの合計 - 2）と計算する。
#pが境界を決める値で、両端から合計5%と設定しているので、半分の2.5%だけ差し引いた97.5%の境界としている。
qt(df = 48, p = 0.975)
## [1] 2.010635
#tの平均から95%の境界線にあたるtを足し引きした範囲に含まれるtの数を計算
sum(t < mean(t) + qt(df = 48, p = 0.975) &
t > mean(t) - qt(df = 48, p = 0.975)) #95%区間に入っている統計量tの推定値の数
## [1] 9500
#5%区間に入っている統計量tの推定値の数
#95%の範囲に含まれるtの個数を10000個から引けば、両端の5%の範囲に含まれるtの個数となる。
10000 - sum(t < mean(t) + qt(df = 48, p = 0.975) & t > mean(t) - qt(df = 48, p = 0.975))
## [1] 500
------------------------------------------------------
シミュレートするごとにこの個数は微妙に変わってくるが、境界線の外側には約500個のtが含まれている（今回はたまたま500個ぴったりだった）。もっとシミュレートするtの個数を増やせば上記の境界線の外側に含まれるtの個数は全体の5%に近づいていく。なお、今回はシミュレーションによってt分布を求めているが、シミュレーションせずとも標本が生成される正規分布の平均と標準偏差、2標本のサンプルサイズがわかれば、自動的に収束先のt分布は決定できる。そして、収束したt分布では上記の境界線の外側には正確に全体の5%のtが含まれている。このページでは、具体的にtの個数がわかる方がイメージしやすいと考えて、シミュレーションによるt分布の再現を採用している。
帰無仮説が正しいとき、同じ正規分布から生成した2標本の平均値は近い値になりやすく、それゆえにtは0付近に多く分布するのであった。逆に、tが0からかけ離れた値になることはめったにない。ここで、有意水準を5%と定めたとしよう。手元の標本から得られた統計量t、これをt0とする、が帰無分布の両端から5%の領域に含まれる、つまり、t0が0からかけ離れた値となれば、帰無仮説を棄却するのであった（図2も参照）。tが0からかけ離れた値になることはめったにないはずなのに、手元のt0は0から逸脱した値を示しているなら、「そもそも帰無仮説が間違っていた」、つまり「2標本は同じ正規分布から生成されたと考えるのは間違い」ではないか、と考えることが帰無仮説の棄却なのである。ちょうど、高校数学における背理法に近い考え方だと言える。t0が0から離れた値をとるというのは、2標本の平均値が大きく異なるとき、に相当している。
なぜ帰無仮説の棄却と「保留」なのか？ 「採択」ではないのか？
統計的仮説検定では、帰無仮説を棄却できないときは、帰無仮説を保留とする。これは差がないことを積極的に肯定しているのではなく、また否定もしない。という立場だ。つまり、帰無仮説が棄却できなかったときは、差があるとも、ないとも主張できないということになる。なぜ、帰無仮説が棄却できなかったときは、帰無仮説を採択、差がないことを積極的に主張することができないのだろうか。これは論理学的な問題である。まずは本当の背理法を使って説明しよう。背理法は下記の図のような論理構造をとっている。例として、京都大学の有名な入試問題「tan1°は有理数か？」を使っている。
tan1°は無理数であることが知られているが、直接的に示すのは難しい。そこで、「tan1°は有理数である」と仮定して矛盾を導く、つまり「tan1°でありながら、有理数である数は存在しない」ことを示す。そうすれば、上の図で言うtan1°の円が、無理数の円の中に取り込まれて、無事tan1°は無理数であると示すことができる。簡単に証明の流れを紹介すると、tan1°が有理数と仮定して、tanの倍角の公式や角度の和の公式を用いてtan30° = 1/√3（無理数）まで計算する。倍角の公式や角度の和の公式は無理数を用いない有理式なのでtan1°を有理数と仮定するとtan30°も有理数となり、本当はtan30°が無理数であることと矛盾するので、tan1°は無理数なのだ。さて、この証明は落としどころを間違えると大失敗する。上の証明では最後の落としどころをtan30° = 1/√3にしたのでうまくいったが、tan45° = 1（有理数）だったらどうだろう。これだとtan1°は有理数ということになってしまうのだろうか？ しかし、tan1°は無理数であるのは事実であり、そうはならないことは明白だ。上の図で言えば、tan1°の円の灰色部分が存在する可能性が示唆されただけで、さらに言えば無理数に取り込まれている部分に関しては何も証明できていない。ゆえに、tan1°は有理数か、無理数か、判断がつかないことになる。これと同じことは統計的仮説検定にも言えてしまう。
帰無仮説の棄却は背理法とほとんど同じ論理的流れだ。上記の図で言う灰色部分を否定することで、異なる母集団に由来した標本の集合に、手元のある標本が取り込まれる、ということだ。一方で、灰色部分を否定できなかった場合、ある標本が同じ母集団に由来しているとは言えないし、異なる母集団に由来しているとも言えないのだ。これが、帰無仮説を採択することができず、保留するしかできいない理由なのだ。
なぜ帰無仮説は「差のない」仮説とするのか？
統計的仮説検定では、差のない仮説=帰無仮説を否定することで、差があることを主張する手続きだ。なぜ、「差のある」仮説を採択する、というような直接的な方法ではないのだろうか？ 理由の一つは、前述の通り、帰無仮説を採択するという結論は取れない点だ。そして、もう一つは差のない仮説はたった一つの仮説であるが、差のある帰無仮説は無限に作ることができる点だ。
差のない仮説は、2標本の平均が同じという仮説なので、言い換えれば、平均値の差は0、という仮説のことだ。一方、差のある仮説は、平均値の差が1だろうが、100だろうが、0以外のあらゆる数で作れてしまう。そして、この差がある仮説を棄却することで言えるのは、平均値の差が1の可能性が低いことや、100の可能性が低いことだけで、もしかしたら平均値の差は10000かもしれないわけで、差がないとは言えないのである。
挑戦的なことを言うのであれば、具体的な差に関する仮説、例えば「平均値の差が10以上でなければ利用するうえで無意味なので、それ以上の差があるか知りたい」とか、「理論上、平均値の差は100以上になるはずだから、本当にそれだけの差が生まれているか知りたい」というような場合は、帰無仮説を差がある仮説とできる場合もあるだろう。ただし、差があるかだけ知れればいい場合が多いし、具体的な○○以上の差、という理論的な仮説を立てることができない場合が多いだろう。
Rで行うStudentのt検定
上記のような背景のもと、t検定を行ってみよう。以下のようなサンプルサイズ20と30のデータがあるとする。ともに同じ標準偏差の正規分布からデータを生成していると考えてもらってよい。有意水準は5%としよう。まずは、データを描画してみて様子を確認する。解析をあれこれする前に、生データを眺めてみるというのは、偏ったデータ解釈をしないためにも重要な行為だ。
------------------------------------------------------
m1 <- c(41, 58, 42, 47, 51, 70, 36, 49, 74, 32,
39, 38, 81, 17, 45, 32, 72, 50, 61, 61)
m2 <- c(101, 69, 50, 69, 88, 55, 100, 60, 68, 46,
80, 74, 51, 60, 75, 70, 75, 63, 83, 46, 75,
53, 53, 56, 103, 47, 85, 50, 91, 82)
d <- data.frame(d = c(m1, m2), g = c(rep("1", 30), rep("2", 20)))
boxplotn(d$d ~ d$g, ylab = "data") #図3の描画
------------------------------------------------------
図3 2標本の描画。平均値に違いがありそうな様子がわかる。
確認してみると、標本2のほうが平均値が大きそうであることがわかる。では本当に違いがありそうかをt検定してみよう。
------------------------------------------------------
fit <- t.test(m1, m2, var.equal = T)#Rでのt検定のコマンド
#2標本の標準偏差（分散）が等しいと仮定できるときはvar.equal = Tとする。（Studentのt検定と呼ばれる）
#2標本の標準偏差（分散）が等しくないと仮定できるときはvar.equal = Fとする。（Welchのt検定と呼ばれる）
fit #t検定の解析結果の確認
##
## Two Sample t-test
##
## data: m1 and m2
## t = -4.0019, df = 48, p-value = 0.0002167
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -29.247030 -9.686303
## sample estimates:
## mean of x mean of y
## 49.80000 69.26667
------------------------------------------------------
fitの中身を確認してみると、いくつかのデータが並んでいる。tは検定統計量tのことで以下のコマンドで計算できるものと同じだ。
------------------------------------------------------
ue <- ((20 - 1)*(sd(m1)^2) + (30 - 1)*(sd(m2)^2))/(20 + 30 - 2)
t0 <- (mean(m1) - mean(m2))/sqrt(ue*(1/20 + 1/30))
t0
## [1] -4.001933
------------------------------------------------------
次にdfは自由度（degree of freedom）のことで、今回は(サンプルサイズの合計 - 2)と計算する。次にp-valueで、あらかじめ定めた有意水準を下回っていれば、帰無仮説を棄却する。p-valueはt分布を使って次のように計算できる。
------------------------------------------------------
df <- 20 + 30 - 2 #自由度の計算
pt(df = df, q = t0, lower.tail = T) * 2
## [1] 0.0002166834
------------------------------------------------------
今回は0.0002167なので、0.05を下回っており、帰無仮説を棄却する。このように、帰無仮説を棄却することで、差があることを示せた場合は有意差がある(significant)と表現する。つまり、この2標本の平均値は違うと、結論を下すわけだ。次の行の”alternative hypothesis: true difference in means is not equal to 0”は「対立仮説(alternative hypothesis)は平均値の差が0ではないこと」を示している。対立仮説とは、帰無仮説ではないほうの仮説、つまり差がある方の仮説だ。95 percent confidence intervalは平均値の差に関する95%信頼区間である。これらの2標本をもとにしたt分布は、平均0のt分布からt0 = -4.001933ずれた形となっている。このときの、95％のtが含まれる範囲は、下記のように計算できる。
------------------------------------------------------
t0 - qt(df = 48, p = 0.975)
## [1] -6.012567
t0 + qt(df = 48, p = 0.975)
## [1] -1.991298
------------------------------------------------------
これはt分布における範囲なのでもとの平均値の差の分布に戻すときは下記のように計算できる。
------------------------------------------------------
(t0 - qt(df = 48, p = 0.975)) * sqrt(ue*(1/20 + 1/30))
## [1] -29.24703
t0 * sqrt(ue*(1/20 + 1/30))
## [1] -19.46667
(t0 + qt(df = 48, p = 0.975)) * sqrt(ue*(1/20 + 1/30))
## [1] -9.686303
------------------------------------------------------
確かにt検定を行った時の95%信頼区間と一致している。また、t0を平均値の差の分布に戻したときの値はちょうど、2標本の平均値の差となっている。
------------------------------------------------------
mean(m1) - mean(m2)
## [1] -19.46667
------------------------------------------------------
最後のmean of xやmean of yはそれぞれ標本の平均値なので以下の値と一致する。
------------------------------------------------------
mean(m1)
## [1] 49.8
mean(m2)
## [1] 69.26667
------------------------------------------------------