累積確率分布の「開き」の大きさ - ryamadaの遺伝学・遺伝統計学メモ

Gene-set enrichment analysisの一環としてこちらでコルモゴロフ=スミルノフ検定(KS検定)の話がちらっと出た
KS検定は１標本と２標本の区別がある
- 観察された標本の元の分布が、ある特定の分布ではない、と考えた方がよいかを判断する(１標本のKS検定)
- もしくは観察された２つの標本の元の分布が同一ではない、と考えた方がよいかを判断する(２標本のKS検定)
KS検定について書くためには、経験分布・経験分布関数というのを知っておくのがよい
経験分布・経験分布関数
- 標本の１番小さい値より小さい範囲の「累積確率」はゼロとする。標本の１番大きい値より大きい範囲の「累積確率」は１とする。途中は、値が来るたびに、標本数nについて $\frac{1}{n}$ だけジャンプして階段状にすることにする。
- 標本の値を $y_1,...,y_n$ とすると $F(x)=\frac{1}{n}\sum_{i}^n \delta(i,x)$ ( $\delta(i,x)$ は $y_i \le x$ のときに1で $y_i > x$ のときは0とする)とも書くが(こちら参照)、式を見ると恐ろしげだけれどやっていることは単純
- ちなみにRなら

y <- rnorm(10)
ecdf.y <- ecdf(y)
plot(ecdf.y)

２標本のKSテストをなぞってみよう
- もちろんRならks.test()関数を使えばよいだけではあるが
- 黒いプロットが赤いプロットより小さめ(左)に位置する

# 2標本を作る
n1 <- 20
n2 <- 30
m1 <- 0
m2 <- 1
sd1 <- 1
sd2 <- 1.5
x <- rnorm(n1,m1,sd1)
y <- rnorm(n2,m2,sd2)
# xの方がyより大きいというのが「正しい」ところ
ks.test(x,y,alternative="less")
ks.test(x,y,alternative="two.sided")
ks.test(x,y,alternative="greater")
ecdf1 <- ecdf(x)
ecdf2 <- ecdf(y)
plot(ecdf(x), xlim=range(c(x, y)))
plot(ecdf(y), add=TRUE, col=2)

- KSテストでは、２つの経験分布関数について、標本x,yの両方の観測値において、２つの経験分布関数の値の差を調べる
- 赤プロットが黒プロットより上にあるかどうかを知りたいとき("less")には、赤の経験分布値から黒の経験分布値を引いた値の最大値を問題にする
- 赤プロットが黒プロットより下にあるかどうかを知りたいとき("greater")には、黒の経験分布値から赤の経験分布値を引いた値の最大値を問題にする
- 赤・黒のどちらが上か下かは気にせず、とにかく違うのかどうかを知りたいときは、赤・黒の値の絶対値の最大値を問題にする
- この差は次のように、各観測ポイントでの値が算出できる
  - 階段関数なので、その段差について考慮するとそうなる

n.x <- length(x)
n.y <- length(y)
w <- c(x, y)
z <- cumsum(ifelse(order(w) <= n.x, 1/n.x, -1/n.y))
plot(w,z)

- ２つのecdfプロットと「差」を上下に並べるとよくわかる

- 検定するためには、この「最大値」が「帰無仮説(２標本が同じ分布に由来する)」のときにどのような分布を取るかと引き比べる必要がある
- これは元の分布の形によらずブラウン運動(ウィーナー過程)から導き出されるので、それに照らしてｐ値化する。
では、次の記事でこの統計量がどんな具合かを調べてみよう