廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

permutation テスト入門

p 値ってなに?

p 値ってなに? ってことをてっとり早く理解するためには permutation テストをやってみるといいと思う。日本語に訳すと「置換検定」とかになるだろうか。と思ったけどやっぱ、permutation順列ですね。「順列組み合わせ検定」とかにしたほうが実態にあってるような気もするけど。(2014年2月2日)

で、p 値ってなに?

帰無仮説を棄却できる最小の優位水準のこと。

せっかくなのでもうちょっと説明する。

2 つのグループ A と B を考える。A と B は男と女とか、右利ききと左利きとか、喫煙者と非喫煙者とか、赤い薬を飲んだ人と青い薬を飲んだ人とかだ。いま、この 2 つのグループ A と B には違いがある、ということをいいたいとしよう。でも「違いがある」っていうことを直接証明するのは難しい。
0.01mm 違ってるだけだって違いは違いだって考える人もいるけど、そんなのは違いのうちに入らないと思う人もいるからだ。

そこでちょっとひねりを加える。もし A と B に違いがないとしたら、いま観測された 0.01mm の差は偶然生じたものだ。その偶然はどのくらいの確率で起こるだろうか。その確率を計算して推定しよう。
推定された確率が低かったら、そんな低い確率のことがたまたま起こるって考えるのは不自然だよね、じゃあやっぱもともと A と B は違うものだったって考えたほうがいいんじゃないの、ということになる。

その「推定された確率」のことを p 値と呼ぶ。

R による permutation テスト入門

例はこのページから拝借。
Permutation test

2 つのグループ、 A チームと B チームの体重がこんなふうに得られた。

A <- c(60,54,75,48,55)
B <- c(49,76,65,58,62)

B の平均は 58.4、A の平均は 62.0、A と B の差は 62.0-58.4=3.6。

Dobs <- mean(B)-mean(A)

さてこの差 (Dobs) は違いのうちに入るか入らないか。

こういうときは「もし A と B に違いがないとしたら……」ということを考えるんだった。A と B に違いがないとしたら、グループ A とグループ B のメンバーをバラバラに入れ替えても同じくらいの差が出ておかしくない。

「グループ A とグループ B をバラバラに入れ替える」のは 10 個のものから 5 個選ぶ場合の数
 _{10} \rm{C}_5 = 252
通りある。

(len <-choose(10,5))

252くらいだったら、全パターンの組み合わせを試せる。

パッケージ gtools に組み合わせを列挙してくれる関数、combinations がある。

library(gtools)
perm1 <- combinations(10,5)
> head(perm1)
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    2    3    4    5
[2,]    1    2    3    4    6
[3,]    1    2    3    4    7
[4,]    1    2    3    4    8
[5,]    1    2    3    4    9
[6,]    1    2    3    4   10

これを使おう。A と B をくっつけて、一個のベクトルにすれば、

> c(A,B)
 [1] 60 54 75 48 55 49 76 65 58 62

perm1 を使って入れ替えられた新しいグループ A ができる。

> c(A,B)[perm1[2,]]
[1] 60 54 75 48 49
> c(A,B)[perm1[3,]]
[1] 60 54 75 48 76

グループ B は A で選ばれなかった方をとればいいから、差集合 setdiff

> setdiff(1:10, perm1[2,])
[1]  5  7  8  9 10
> setdiff(1:10, perm1[3,])
[1]  5  6  8  9 10

これを使って

> c(A,B)[setdiff(1:10, perm1[2,])]
[1] 55 76 65 58 62
> c(A,B)[setdiff(1:10, perm1[3,])]
[1] 55 49 65 58 62

と書ける。

for 文で全パターン回す。

Ds <- numeric()
for(i in 1:len){
  pA <- c(A,B)[perm1[i,]]
  pB <- c(A,B)[setdiff(1:10, perm1[i,])]
  Ds[i] <- mean(pB)-mean(pA)
} 

結果として 3.6 くらいの差は何%あるだろうか。3.6 以上の差が何回出たかは、

 sum(Ds >= Dobs)
[1] 75

75回。まあまま多そうだ。これを組み合わせ全パターンの数、252でわると、

> sum(Ds >= Dobs)/len
[1] 0.297619

29.75%。これが目標としていた p 値であった。上記のページだと 3.57%になっているが、ちゃんと計算すると29.75%だ。

permutation テストはこういうふうに「全パターンやってみる」という素朴な方法なので、平均の差じゃなくてもなんでもいい。メジアンの差でもいい 。
関数化しておく。

library(gtools)
permtest <- function(A,B, average=mean){
Dobs <- average(B)-average(A)
m <- length(A)
n <- length(B)
len <-choose(m+n,m)
perm1 <- combinations(m+n,n)
head(perm1)
Ds <- numeric()
for(i in 1:len){
  pA <- c(A,B)[perm1[i,]]
  pB <- c(A,B)[setdiff(1:10, perm1[i,])]
  Ds[i] <-average(pB)-average(pA)
}
sum(Ds >= Dobs)/len
}

> permtest(A,B)
[1] 0.297619
> permtest(A,B, median)
[1] 0.1904762

coin パッケージを使うと permutation テストができるらしい。

参考文献はこちら。
Ferry Butar Butar, Ph.D. & Jae-Wan Park. Permutation Tests for Comparing Two Populations