上昇トレンドの検定 ヨンクヒール・タプストラ検定

時系列データを持っているとき、時間経過で値が上昇/減少しているものがほしい、とする。
そういうとき、ヨンクヒール・タプストラ検定というものがあるらしい。

install.package("SAGx")
library(SAGx)
JT.test

でやってみよう。
 

この本の例をやってみる。

#3サンプルのマイクロアレイデータ的な
y <- rbind(c(1.06, 1.03, 1.10, 1.66, 1.72, 1.82, 1.91),
				c(0.98, 1.02, 1.08, 1.58, 1.68, 1.78, 1.89),
				c(1.09, 1.10, 1.11, NA, 1.71, 1.83, 1.92))
JT.test(y[1, ], rep(1, 7), alternative="increasing")
matplot(t(y), type="l")

class was not an ordered factor.  Redefined to be one.

	Jonckheere-Terpstra

data:  y[1, ] by rep(1, 7) 
J = 26.5, p-value = 1
alternative hypothesis: increasing: 1 

 警告メッセージ: 
In cor(rank(x), rank(class)) :  標準偏差が0です 

本ではものすごい小さいp値になるらしいけど1になる。
保留

 
JT.testのexampleをやってみる。

# Enter the data as a vector
A <- as.matrix(c(99,114,116,127,146,111, 125,143,148,157,133,139, 149, 160, 184))
# create the class labels
g <- c(rep(1,5),rep(2,5),rep(3,5))
# The groups have the medians
tapply(A, g, median)
# JT.test indicates that this trend is significant at the 5
JT.test(data = A, class = g, labs = c("GRP 1", "GRP 2", "GRP 3"), alternative = "two-sided")
matplot(matrix(A, nr=5), type="l")

class was not an ordered factor.  Redefined to be one.

	Jonckheere-Terpstra

data:  A by g 
J = 16, p-value = 0.02311
alternative hypothesis: two-sided: GRP 1 two-sided: GRP 2 two-sided: GRP 3 


 
マイクロアレイの解析っぽくやってみる。

m <- matrix(rnorm(60000), 2000, 30)
b <- ordered(rep(1:3, each=10))
p <- JT.test(m, b)$p.value
#ランダムサンプリングなp値が一様分布になっている
plot(sort(p))
#マイクロアレイっぽい色付けにしたいだけ
library(gplots)
heatmap(m[p <0.05,], col=greenred(7), scale="none", ColSide=c("red", "blue", "green")[b])