2014-04-16

順序とかのための正確検定法パッケージ

R exactRankTests

量的データだけれど、分布を仮定していいのかわからない
同じ値がある(Tie(タイ)がある)
グレード別のデータ(タイばっかり)
というような場合
しかもサンプル数が、え、これだけ？
いくつかの方法を正確検定中心にまとめたらしいパッケージがexactRankTestsパッケージ
wilcox.exact() ウィルコクソン・ランク和テストのexact版
- サンプル数が少ないとき(50個)は正確検定、多いときはそうでないのがデフォルトであったり、タイのサンプル(同じ値のサンプル)があるときはどうする、２群に対応があるときは…、など、デフォルトの設定を覚えるのが面倒くさければ、引数指定するのが無難

# サンプル数が少ない、タイあり
A <- sample(1:10,18,replace=TRUE) # 必ずタイがある
B <- sample(5:18,15,replace=TRUE)
wilcox.exact(A,B,paired=FALSE)
wilcox.exact(A,B,paired=FALSE,exact=TRUE)
wilcox.exact(A,B,paired=FALSE,exact=FALSE)
# サンプル数が少ない、タイなし
A <- sample(1:20,10)*2 # タイなし　偶数
B <- sample(5:18,12)*2 + 1 # タイなし　奇数
wilcox.exact(A,B,paired=FALSE)
wilcox.exact(A,B,paired=FALSE,exact=TRUE)
wilcox.exact(A,B,paired=FALSE,exact=FALSE)

# サンプル数が多い、タイあり
A2 <- sample(1:20,30,replace=TRUE) # タイあり
B2 <- sample(5:25,25,replace=TRUE) # タイあり
wilcox.exact(A2,B2,paired=FALSE)
wilcox.exact(A2,B2,paired=FALSE,exact=TRUE)
wilcox.exact(A2,B2,paired=FALSE,exact=FALSE)
# サンプル数が多い、タイなし
A2 <- sample(1:50,30)*2 # タイなし　偶数
B2 <- sample(5:55,20)*2 + 1 # タイなし　奇数
wilcox.exact(A2,B2,paired=FALSE)
wilcox.exact(A2,B2,paired=FALSE,exact=TRUE)
wilcox.exact(A2,B2,paired=FALSE,exact=FALSE)

perm.test()　パーミュテーションテスト、exactオプションで正確テストに
- 値は整数値
- サンプル数が少ないとき(50個)は正確検定、多いときはそうでないのがデフォルトであったり、タイのサンプル(同じ値のサンプル)があるときはどうする、２群に対応があるときは…、など、デフォルトの設定を覚えるのが面倒くさければ、引数指定するのが無難
- 指定するべき引数は、２群の値ベクトル(もちろんだが)の他、paired ={TRUE,FALSE}, exact={TRUE,FALSE},alternative={"less","greater","two.sided"}、は指定しておけば安心

A <- sample(1:10,18,replace=TRUE) # 必ずタイがある
B <- sample(5:18,15,replace=TRUE)
perm.test(A,B,paired=FALSE,exact=TRUE,alternative="less")
perm.test(A,B,paired=FALSE,exact=TRUE,alternative="greater")
perm.test(A,B,paired=FALSE,exact=TRUE,alternative="two.sided")

- A群の方が小さいから"less"で小さいp値、"greater"で大きいp値

> perm.test(A,B,paired=FALSE,exact=TRUE,alternative="less")

        2-sample Permutation Test

data:  A and B
T = 99, p-value = 0.0001732
alternative hypothesis: true mu is less than 0

> perm.test(A,B,paired=FALSE,exact=TRUE,alternative="greater")

        2-sample Permutation Test

data:  A and B
T = 99, p-value = 0.9999
alternative hypothesis: true mu is greater than 0

> perm.test(A,B,paired=FALSE,exact=TRUE,alternative="two.sided")

        2-sample Permutation Test

data:  A and B
T = 99, p-value = 0.0002749
alternative hypothesis: true mu is not equal to 0

ランク和統計量の分布関数に関する関数　dperm(),pperm(),qperm(),rperm()
- サンプルをシャッフルし、そのたびに「はじめのm個のランクの和」＝ランク和統計量を算出する。
- 「はじめの」というとき、シャッフルしていれば、ただ、m個のどれが選ばれるか、ということになる。
- x,yと２群の値を与えて、その「はじめのx群の個数のランクの和」をとれば、それは「ランク和統計量」なので、その値(標本のランク和)とシャッフルしたときのランク和の値とを比較してp値を出す(pperm())。ランク和の分布を取るのがdperm()。p値に相当するランク和値を出すのがqperm()、ランク和統計量の乱数を出すのがrperm()

library(exactRankTests)
# ２群の値
x <- c(0.5, 0.5, 0.6, 0.6, 0.7, 0.8, 0.9)
y <- c(0.5, 1.0, 1.2, 1.2, 1.4, 1.5, 1.9, 2.0)
# それをwilcoxonランク和するように値変換
r <- cscores(c(x,y), type="Wilcoxon")
# 標本のランク和がsum(r([seq(along=x)])、シャッフルするベクトルがr、xの要素数が7なので、rから7要素を選んでそのランク和は？という処理
pperm(sum(r[seq(along=x)]), r, 7)
hist(rperm(1000,r,7))

ansari.exact()
- Ansari-Bradleyテスト(こちらのMATLABの記事がよい)の正確版
- ２サンプルが同じ分布由来かどうかを確かめたい。ただし、中央値は両者で同じで分布の形だけ(ばらつきを含めて)が違うかどうかを調べる

ramsay <- c(111, 107, 100, 99, 102, 106, 109, 108, 104, 99,
101, 96, 97, 102, 107, 113, 116, 113, 110, 98)
jung.parekh <- c(107, 108, 106, 98, 105, 103, 110, 105, 104,
100, 96, 108, 103, 104, 114, 114, 113, 108, 106, 99)
h <- hist(c(ramsay,jung.parekh))
hist(ramsay,density=20,breaks=h$breaks,ylim=c(0,max(h$counts)))
hist(jung.parekh,density=17,col=2,breaks=h$breaks,add=TRUE)

ansari.test(ramsay, jung.parekh)
ansari.exact(ramsay, jung.parekh)