時系列データを持っているとき、時間経過で値が上昇/減少しているものがほしい、とする。
そういうとき、ヨンクヒール・タプストラ検定というものがあるらしい。
install.package("SAGx") library(SAGx) JT.test
でやってみよう。
マイクロアレイデータ統計解析プロトコール―Excelを中心としたデータの標準化から有意差解析、クラスタリング、ネットワーク解析法のすべて (実験医学別冊 23)
- 作者: 藤渕航,堀本勝久
- 出版社/メーカー: 羊土社
- 発売日: 2008/05/01
- メディア: 単行本
- 購入: 4人 クリック: 92回
- この商品を含むブログ (3件) を見る
#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])