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

 データ中のどの変数に群間の差が見られるかが知りたいとします。群間の中央値を比較する場合は、ウィルコクソンの順位和検定を用います。以下は、データ中の全ての変数にウィルコクソンの順位和検定を行い、その結果として得られる統計量が大きい順に出力する関数です。*1

# 全ての変数にウィルコクソンの順位和検定を行う関数
wilcoxRep <- function(dat){
	x <- ncol(dat)
	y <- x - 1
	n <- 1 : y 
	mat <- matrix(rep(0, y), ncol = 1)
	colnames(mat) <- "W"
	rownames(mat) <- colnames(dat[, -x])
	for (i in n){
		a <- wilcox.test(dat[, i] ~ dat[, x], data = dat, paired = F)
		options(warn = -1)
		mat[i, 1] <- a$statistic
	}
	wc <- mat[order(mat[, 1], decreasing = TRUE), ]
	wc.mat <- matrix(wc, ncol = 1, dimnames = list(names(wc), "W"))
	print(wc.mat)
}

 例えば、kernlabパッケージのspamデータに対して上記の関数を用いるには、以下のような処理を行います。*2

# データの読み込み
library(kernlab)
data(spam)
# 関数の実行(データの指定)
wilcoxRep(spam)

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

                          W
hp                3417266.5
george            3218129.5
hpl               3205327.0
num1999           3048972.0
labs              2912719.5
num650            2879927.5
num85             2868122.5
edu               2849913.5
lab               2838851.0
technology        2818578.5
meeting           2792788.5
telnet            2786002.0
data              2758544.5
pm                2738175.5
project           2719578.0
charSquarebracket 2711287.0
num857            2708625.5
num415            2699893.5
re                2698012.5
original          2680076.5
conference        2675045.0
cs                2659248.0
charSemicolon     2638835.5
charRoundbracket  2621359.5
table             2540946.5
parts             2529060.5
direct            2483895.0
num3d             2480110.5
font              2414968.0
report            2319855.0
addresses         2170161.0
will              2111655.0
people            2096653.5
credit            2041585.0
charHash          2015146.5
make              1998442.5
order             1942506.0
address           1912670.0
email             1873687.5
receive           1861383.5
internet          1841103.0
over              1841056.5
mail              1823039.0
business          1777956.5
num000            1742893.5
all               1641171.0
money             1628113.0
remove            1500570.5
you               1463480.5
our               1461133.5
free              1350962.5
capitalTotal      1201589.5
charDollar        1151754.0
your              1111183.0
capitalAve        1070571.5
capitalLong        989888.5
charExclamation    864111.0

 また、群間の平均値を比較する場合は、以下のように、t検定を用います。

# 全ての変数にt検定を行う関数
tRep <- function(dat){
	x <- ncol(dat)
	y <- x - 1
	n <- 1 : y 
	mat <- matrix(rep(0, y), ncol = 1)
	colnames(mat) <- "t"
	rownames(mat) <- colnames(dat[, -x])
	for (i in n){
		a <- t.test(dat[, i] ~ dat[, x], data = dat, paired = F, var.equal = F)
		options(warn = -1)
		mat[i, 1] <- a$statistic
	}
	tv <- mat[order(mat[, 1], decreasing = TRUE), ]
	tv.mat <- matrix(tv, ncol = 1, dimnames = list(names(tv), "t"))
	print(tv.mat)
}

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

# 関数の実行(データの指定)
tRep(spam)

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

                            t
hp                 22.2798617
hpl                20.0753758
george             15.6899835
labs               14.4017342
num1999            13.7587247
num85              12.6387644
num650             12.4706106
edu                12.3504927
re                 11.6068291
meeting            11.5940095
original           11.3600678
lab                11.3293427
technology         11.1527540
telnet             10.7373572
pm                 10.2903496
data               10.0486150
num857              9.6627262
num415              9.5126392
cs                  8.2282261
project             7.9702950
conference          7.0742747
charRoundbracket    6.0040846
direct              5.2024130
charSquarebracket   5.1959098
charSemicolon       4.8421622
table               3.7029262
parts               2.5755406
address             2.4943663
will               -0.5714892
num3d              -3.1425915
charHash           -3.7867131
report             -4.1453434
font               -5.3923937
capitalAve         -6.0799143
people             -8.5466943
make               -8.5512703
mail               -9.5536316
credit            -10.6385427
addresses         -11.1982116
capitalLong       -12.1932457
internet          -12.4558833
email             -12.7857414
money             -13.0374454
all               -13.7570298
order             -14.4406553
receive           -14.7885890
over              -15.0511733
capitalTotal      -15.0691498
business          -15.6505238
our               -16.4153699
free              -16.7803447
charExclamation   -17.2598988
charDollar        -19.0039911
remove            -19.5629979
num000            -19.5732386
you               -19.8518352
your              -27.0946879

 次回は、多群間で差のある変数を特定する例を扱う予定です。*3

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

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

*3:裏RjpWikiの方がこちらで「中央値の検定」と「平均値の検定」を1つの関数にまとめる方法を紹介してくださいました。いつもありがとうございます。この記事をアップした段階で、forをapply系に直されるかな、と予想はしていました。。。