読者です 読者をやめる 読者になる 読者になる

多群間で差のある変数を特定する

r stat

 前回の二群間で差のある変数を特定するという記事に続き、今回は、三群以上(多群)のデータで差のある変数を特定する方法を考えます。*1
 群間の中央値を比較する場合は、クラスカル・ウォリス検定を用います。以下は、データ中の全ての変数にクラスカル・ウォリス検定を行い、その結果として得られる統計量が大きい順に出力する関数です。*2

# 全ての変数にクラスカル・ウォリス検定を行う関数
kruskalRep <- function(dat){
	x <- ncol(dat)
	y <- x - 1
	n <- 1 : y 
	mat <- matrix(rep(0, y), ncol = 1)
	colnames(mat) <- "chi-squared"
	rownames(mat) <- colnames(dat[, -x])
	for (i in n){
		a <- kruskal.test(dat[, i] ~ dat[, x])
		options(warn = -1)
		mat[i, 1] <- a$statistic
	}
	kw.chi <- mat[order(mat[, 1], decreasing = TRUE), ]
	kw.chi.mat <- matrix(kw.chi, ncol = 1, dimnames = list(names(kw.chi), "chi-squared"))
	print(kw.chi.mat)
}

 例えば、irisデータに対して上記の関数を用いるには、以下のような処理を行います。*3

# 関数の実行(データの指定)
kruskalRep(iris)

 上記のスクリプトを実行すると、以下のような結果が表示されます。

             chi-squared
Petal.Width    131.18538
Petal.Length   130.41105
Sepal.Length    96.93744
Sepal.Width     63.57115

 また、群間の平均値を比較する場合は、以下のように、一元配置の分散分析を用います。

# 全ての変数に分散分析を行う関数
anovaRep <- function(dat){
	x <- ncol(dat)
	y <- x - 1
	n <- 1 : y 
	mat <- matrix(rep(0, y), ncol = 1)
	colnames(mat) <- "F"
	rownames(mat) <- colnames(dat[, -x])
	for (i in n){
		a <- oneway.test(dat[, i] ~ dat[, x])
		options(warn = -1)
		mat[i, 1] <- a$statistic
	}
	fv <- mat[order(mat[, 1], decreasing = TRUE), ]
	fv.mat <- matrix(fv, ncol = 1, dimnames = list(names(fv), "F"))
	print(fv.mat)
}

 上記の関数を用いるには、以下のような処理を行います。

# 関数の実行(データの指定)
anovaRep(iris)

 上記のスクリプトを実行すると、以下のような結果が表示されます。

                      F
Petal.Length 1828.09195
Petal.Width  1276.88456
Sepal.Length  138.90829
Sepal.Width    45.01204

(追記)裏RjpWikiの方が、こちらの記事で、1つの関数でwilcox.test、t.test、kruskal.test、oneway.testを実行する方法を教えてくださいました。どうもありがとうございます。以下に、その関数を引用させて頂きます。

# 全ての変数に wilcox.test / t.test / kruskal.test / oneway.test を行う関数
Rep = function(data, group, method = c("wilcox", "t", "kruskal", "oneway"), sort = TRUE) {
  name = method = match.arg(method)
  if (method == "wilcox") {
    func = wilcox.test
  } else if (method == "t") {
    func = t.test
  } else if (method == "kruskal") {
    func = kruskal.test; name = "chi-squared"
  } else { # if (method == "oneway")
    func = oneway.test; name = "F"
  }
  old = options(warn = -1)
  ans = sapply(data, function(d) func(d ~ group)$statistic)
  options(old)
  ans = cbind(ans)
  rownames(ans) = colnames(data)
  if (sort) ans = cbind(ans[order(ans[,1], decreasing = TRUE), ])
  colnames(ans) = name
  ans
}

# 関数の実行
# library(kernlab)
# data(spam)
# Rep(spam[c(1:10, 16, 19)], spam$type)
# Rep(spam[c(1:10, 16, 19)], spam$type, sort=FALSE)
# Rep(spam[c(1:10, 16, 19)], spam$type, method = "t")
# Rep(spam[c(1:10, 16, 19)], spam$type, method = "t", sort=FALSE)

Rep(iris[1:4], iris$Species, method="k") # method は省略形でも指定できる
Rep(iris[1:4], iris$Species, method="k", sort=FALSE)
Rep(iris[1:4], iris[,5], method="o") # group はベクトルでなければならない
Rep(iris[1:4], iris[,5], method="o", sort=FALSE)

*1:前回の記事に関して、裏RjpWikiの方がこちらで2つの関数を1つにまとめる方法を紹介してくださいました。その方法は、今回扱った多群の検定にも応用することが可能です。

*2:統計量が大きい順に並び替えない場合は、kw.chi.matではなく、matをprintしてください。

*3:この関数では、irisデータのように、データ中の最後の列に群情報があることを想定しています。