統計的仮説検定（t-検定と多重検定補正）

まずはホームディレクトリを設定します。

home <- [配布したフォルダをコピーしたパス]
setwd(home)

library(data.table)

Warning message:
In scan(file = file, what = what, sep = sep, quote = quote, dec = dec,  :
  EOF within quoted string

t検定について説明します。以下のように、10個の遺伝子に対して二つの条件下（例えば、Tumor / Normal）で5サンプルずつの遺伝子発現データ考えます。目標は、条件間で発現量の異なる遺伝子群を同定することです。

set.seed(100)
n <- 10
p <- 5
cond1 <- matrix(rnbinom(n*p,2,mu = 300),ncol=p,dimnames=list(paste0("Gene",1:n),paste0("Condition",rep(1,p))))
cond2 <- matrix(rnbinom(n*p,2,mu = 100),ncol=p,dimnames=list(paste0("Gene",1:n),paste0("Condition",rep(2,p))))
data_example <- cbind(cond1,cond2)

	Condition1	Condition1	Condition1	Condition1	Condition1	Condition2	Condition2	Condition2	Condition2	Condition2
Gene1	123	91	178	86	327	92	86	75	46	329
Gene2	419	250	648	210	264	210	93	43	89	62
Gene3	277	119	397	148	538	87	58	87	46	33
Gene4	306	198	530	139	428	163	246	54	96	135
Gene5	243	211	172	152	47	102	3	233	58	49
Gene6	205	299	620	202	206	127	64	128	52	284
Gene7	125	56	361	549	80	345	62	189	39	45
Gene8	513	14	217	451	179	89	37	294	125	169
Gene9	71	525	36	276	83	24	141	70	29	166
Gene10	288	151	1005	108	275	37	96	228	7	195

二つの条件間の違いを見る方法はいくつか考えられます。一つは、log fold change (LFC)を計算することです。遺伝子\(i\)に対するLFCとは

\[\mbox{logFC}_i=\mbox{log2}(M_{i,1} / M_{i,2})\]

(\(M_{i,j}\)は遺伝子\(i\)条件\(j\)の発現量の平均です。

M1 <- apply(cond1,1,function(x)mean(x))
M2 <- apply(cond2,1,function(x)mean(x))
(logFC <- log2((M1 + 1) / (M2 + 1)))

    Gene1     Gene2     Gene3     Gene4     Gene5     Gene6     Gene7     Gene8     Gene9    Gene10 
0.3557164 1.8390281 2.2314946 1.2001075 0.8831863 1.2195792 0.7797122 0.9395588 1.1951303 1.6894567 

LFCはどの程度、違うかを表す指標ですが、どの閾値以上であれば差があるかの明確な基準はありません。また遺伝子ごとにデータのばらつき方は異なり、得られたLFCが本当に偶然かどうかも見定める必要がありますがそれはできません。そのために、統計的仮説検定という枠組みを用いて、ばらつきを確率的に考えて、本当にその差があるのかを評価します。t-検定は、二群間の平均の差を検定する代表的な検定手法です。

早速、t-検定を実行してみましょう。遺伝子1に対してt-検定を実行してみます。

x <- cond1[1,]
y <- cond2[1,]
(result_t <- t.test(x,y,var.equal = TRUE))

    Two Sample t-test

data:  x and y
t = 0.51975, df = 8, p-value = 0.6173
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -121.6617  192.4617
sample estimates:
mean of x mean of y 
    161.0     125.6 

p値はいくらかというと、

result_t$p.value

[1] 0.6173102

となり、有意水準5%で考えると有意ではありません。Rでt-検定を行うのは非常に簡単です。ただし、その仕組みを知っていないと、あまり意味も分からずにp値を振り回したり、振り回されてしまいます。さきほどの数行で何をやったのかの意味を理解してみようと思います。

　統計的仮説検定では、まず帰無仮説と対立仮説を立てます。帰無仮説とは最終的に否定したい事柄を、対立仮説にはそれと対比をなす事柄を考えます。今回の場合だと、条件1と条件2で、発現量が異なる遺伝子群を同定したいので、仮説を以下のようにします。 \[\mbox{帰無仮説}:M_{i,1}=M_{i,2}\] \[\mbox{対立仮説}:M_{i,1}\neq M_{i,2}\]

次に、帰無仮説が本当ではないことを確かめるために、帰無仮説を定量化した検定統計量を求めます。

t-統計量

\[\mbox{t-統計量} = \frac{\mbox{群間の平均の違い}}{\mbox{各群のばらつき}}\]

先ほどのt統計量は、

result_t$statistic

        t 
0.5197481 

です。

このt-統計量は単なる数値ではありません。自分の確率的なばらつきかたを表す確率分布を背負っています。特に帰無仮説の確率分布なので、帰無分布と呼びます。仮説検定のポイントはこの帰無分布を正しく理解することです。t統計量の従う分布はt分布という名前がついており、その形はサンプル数のみによって決まります。Rでは簡単にt分布からデータを生成させることができるので、実際に見てみましょう。

data_t <- rt(10000,df = n - 2)
hist(data_t,freq = F,xlab="t",ylab="Probability")

上図はt分布からt統計量を10000個、仮想的に生成したときの頻度分布です。全ての頻度の和を取ると、1になるように正規化しています。これが先ほど実行したt検定の帰無分布です。帰無仮説は二群間の平均の違い=0、すなわち帰無仮説に対応するt統計量が0なので、理論上はゼロになるはずですが、分布上は、確率的なばらつきを持っています。別の言い方をすると、もしも、10000回（実際にはそんなことはできませんが）、繰り返しt統計量を計算した時のばらつき方を表しています。

繰り返しになりますが、この分布は帰無仮説に対応したt統計量のばらつきを表しています。では、対立仮説に対応するt統計量の値(t=-2.026503)はどのあたりかというと、

hist(data_t,freq = F,xlab="t",ylab="Probability")
abline(v=result_t$statistic,col="red")

確かに帰無仮説に対応する0ではないですが、誤差の範囲内と言えそうです。これが確率的に誤差の範囲かどうかを判定するには、有意水準という閾値を設けます。具体的に有意水準5%といったときには、

図中青色の線が有意水準に対応する閾値を表しており、その両側に対応する確率の和が有意水準に対応する5%となっています。有意水準5%とは、偶然0ではない（帰無仮説が正しい）とする閾値以上となる確率が高々5%であると言っています。仮に間違って帰無仮説が間違っているといったとしても、100回中5回程度だということです。裏を返せば、もしもその閾値を越えていれば、それは帰無仮説の下では、滅多に起こりえないことが起こったということになります。したがって、5%という誤差はあるけれども、今回の帰無仮説は成り立たないすなわち有意ではないと考えるわけです。

今回は、5%の有意水準を越えていないので、帰無仮説は5%の有意水準で棄却できません。

ではp値とは何かというと、図中の赤色の線の右側の確率の和です（もし負であれば左側）。先ほど、p値は0.613でしたが、右側だけの有意水準に対応する領域は0.025(2.5%)なので、この値を用いても先ほどと同様に赤色の線が青色の線を越えていないということができます。つまり、p値は有意水準に対応する閾値を越えているかを判断するための値です。

全ての遺伝子に対して、検定してみましょう。

pvals <- NULL
for(i in 1:n){
x <- cond1[i,]
y <- cond2[i,]
result_t <- t.test(x,y,var.equal = TRUE)
pvals[i] <- result_t$p.value
}
pvals

 [1] 0.61731016 0.01664043 0.01819819 0.05069772 0.17883478 0.08861926 0.40741858 0.22987814 0.27743056 0.17355748

次に有意水準5%の説明が本当かどうかを確かめてみましょう。確かめ方は簡単です。同じ確率分布から二群の遺伝子発現量を生成させて、5%有意水準でのt検定を10000回繰り返したときに、間違って帰無仮説を棄却する回数がおよそ500回であれば5%の意味は正しいと言えます。

set.seed(100)
n <- 10000
p <- 5
cond1_many <- matrix(rnbinom(n*p,2,mu = 100),ncol=p,dimnames=list(paste0("Gene",1:n),paste0("Condition",rep(1,p))))
cond2_many <- matrix(rnbinom(n*p,2,mu = 100),ncol=p,dimnames=list(paste0("Gene",1:n),paste0("Condition",rep(2,p))))
pvals_many <- NULL
for(i in 1:n){
  x <- cond1_many[i,]
  y <- cond2_many[i,]
  result_t <- t.test(x,y,alternative = "two.sided",var.equal = TRUE)
  pvals_many[i] <- result_t$p.value
}
sum(pvals_many <= 0.05) / n

[1] 0.044

確かに、5%くらいの確率で帰無仮説を誤って棄却しているようです。

今回は、両側検定と呼ばれる方法で検定を行いました。それがalternative="two.sidedの部分です。両側の意味は、対立仮説が\(M_{i,1} < M_{i,2}\)と\(M_{i,1}> M_{i,2}\)(つまり\(M_{i,1}\neq M_{i,2}\))の両方を含めているということです。片側検定の場合は、どちらかを検定することになります。つまり、

\[\mbox{帰無仮説}:M_{i,1}>M_{i,2}\] \[\mbox{対立仮説}:M_{i,1}\neq M_{i,2}\]

を検定します。このときの注意は、有意水準の考え方です。両側検定で有意水準5%といった場合には、二つのしきい値は、外側の和が5%になるようなしきい値となりますが、片側検定ですと、片方のしきい値しか用いなので、片方のしきい値の外側の和が5%となるようなしきい値となります。

alternative="greater"とすれば、対立仮説は\(M_{i,1} > M_{i,2}\)となり、alternative="less"とすれば対立仮説は\(M_{i,1} < M_{i,2}\)となります。t.test関数の基本設定ではalternative="two.sided"となっています。

x <- cond1[1,]
y <- cond2[1,]
(result_t <- t.test(x,y,greater="greater",var.equal = TRUE))

    Two Sample t-test

data:  x and y
t = 0.51975, df = 8, p-value = 0.6173
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -121.6617  192.4617
sample estimates:
mean of x mean of y 
    161.0     125.6 

var.equalは2群のばらつきを同じと仮定するか否かを設定します。今回は同じと仮定したのでvar.equal=TRUEとしています。もしvar.equal=FALSEとすると、ウェルチの検定という方法になります。状況に応じて使い分けてください。

ところで、有意水準5%で検定を100回繰り返すと、有意水準の定義から、確率的には5回、間違って帰無仮説を棄却することになります。遺伝子数は２万近くありますから、遺伝子ごとに検定を行うと20000回の検定を行うことになり、理論上は$20000\times 0.05=1000$回は間違って、発現の差異がないのに、あると判断してしまうことになります。つまり、擬陽性の遺伝子を1000個程度含んでしまうことになります。統計的仮説検定を繰り返し行うことを多重検定と言い、擬陽性をいかにして減らすかという問題は統計学では多重検定補正といって、重要な問題です。擬陽性率をコントロールして擬陽性を減らすために、いくつかの方法が提案されています。研究論文で広く用いられている方法の一つはBenjamini and HochbergによるFalse discovery rate (FDR)です。FDRとは間違って棄却した帰無仮説のうち、本当は帰無仮説が正しい割合です。FDR=5%といった場合には、本当は発現差異がないのに、あると判断してしまった遺伝子1000個のうち、\(1000\times 0.05=50\)個の遺伝子は本当は発現差異がないという仮定で、p値を補正します。詳しい原理はここでは省略します。FDRによる多重検定補正を行うには、p.adjust関数を用います。まず、10個の遺伝子に対するp値が以下のように与えられたとします。

pvals

 [1] 0.61731016 0.01664043 0.01819819 0.05069772 0.17883478 0.08861926 0.40741858 0.22987814 0.27743056 0.17355748

5%有意な遺伝子は

pvals<=0.05

 [1] FALSE  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE

です。

これに対して、FDRによる多重検定補正を行うには以下のように実行します。

(qvals <- p.adjust(pvals,method="BH"))

 [1] 0.61731016 0.09099094 0.09099094 0.16899241 0.29805796 0.22154816 0.45268731 0.32839735 0.34678820 0.29805796

FDR5%の結果は、以下のようになります。

qvals <= 0.05

 [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE

p値による結果では有意だったものが、FDRによる多重検定補正で有意ではなくなりました。多重検定補正により有意となる基準が厳しくなるので、実際の解析では擬陽性の遺伝子の数を減らすことができます。