2018.11.13

数理統計学演習

ノンパラメトリック検定

東京大学大学院農学生命科学研究科　大森宏

http://lbm.ab.a.u-tokyo.ac.jp/~omori/kensyu/nonpara18.htm

１．ノンパラメトリック検定とは

(1)

(2)

２．１標本問題（one sample problem）

2-1．概要

(5)

2-2．符号検定（sign test）

(7)

(8)

2-3．Wilcoxon の符号順位和検定（signed rank sum test）

表１．順位のあり得る組み合わせ（n＝4）

１	２	３	４	R+	R-	T
＋	＋	＋	＋	10	0	0
−	＋	＋	＋	9	1	1
＋	−	＋	＋	8	2	2
＋	＋	−	＋	7	3	3
＋	＋	＋	−	6	4	4
−	−	＋	＋	7	3	3
−	＋	−	＋	6	4	4
−	＋	＋	−	5	5	5
＋	−	−	＋	5	5	5
＋	−	＋	−	4	6	4
＋	＋	−	−	3	7	3
−	−	−	＋	4	6	4
−	−	＋	−	3	7	3
−	＋	−	−	2	8	2
＋	−	−	−	1	9	1
−	−	−	−	0	10	0

(9)

2-4．解析例

データダウンロード

符号検定

mse <- read.csv("mse.csv")	#データ読み込み
d <- mse$after - mse$before	#処方後得点−処方前得点
np <- length(d[d > 0])
nm <- length(d[d < 0])
binom.test(c(np,nm))

Wilcoxon の順位和検定

wilcox.test(d)	# Wilcoxon 検定を行う関数
wilcox.test(mse$after, mse$before, paired=TRUE)	# このように記述しても良い

 警告メッセージ: 
1: In wilcox.test.default(mse$before, mse$after, paired = TRUE) :
   タイがあるため、正確な p 値を計算することができません  
2: In wilcox.test.default(mse$before, mse$after, paired = TRUE) :
   ゼロ値のため、正確な p 値を計算することができません 

という表示が出て何となくすわりが悪い．タイがあるときでも Wilcoxon 検定で 近似的な p 値を出すには，"coin" (Conditional Inference Procedures in a Permutation Test) を”パッケージのインストール”を用いてインストールすれば良い．

library(coin)			# coin ライブラリー使用宣言
wilcoxsign_test(mse$after ~ mse$before)		# タイがあることに対応
wilcoxsign_test(mse$after ~ mse$before, zero.method="Pratt") # 0 に対応

t 検定

t.test(d)　		# t 検定を行う関数
t.test(mse$after, mse$before, paired=TRUE)	# こう書いてもよい．

それ程高齢でない層の抽出

符号検定

young <- which(mse$age<86)	 # 年齢86歳未満のデータ番号を取得
d2 <- d[young]
np <- length(d2[d2 > 0])
nm <- length(d2[d2 < 0])
binom.test(c(np,nm))

Wilcoxon の順位和検定

wilcox.test(d2)
wilcoxsign_test(mse$after[young] ~ mse$before[young], zero.method="Pratt")

t 検定

t.test(d2)
t.test(mse$after[young], mse$before[young], paired=TRUE)

2-5．検出力（power）の違い

　まず母集団分布として， 標準正規分布 N(0, 1)，自由度 10 の t 分布 t(10)，自由度 2 の t 分布 t(2) を考える．t(10) は正規分布より 少し裾が重いが，t(2) は裾が大変重く，極端に大きな値や小さな値が出やすくなっている．

シミュレーションの概要

まず，データ数が n = 10 のときを考える．帰無仮説が真であるとする．すなわち，母集団平均 m = 0 である． この母集団から n 個の正規乱数 x を抽出し，そのデータに対し符号検定と Wilcoxon の順位和検定， t 検定を行い，それぞれの p 値を出す．

n <- 10	# サンプルサイズ
m <- 0	# 真の平均
x <- rnorm(n, mean=m)	# 標準正規乱数 n 個
np <- length(x[x > 0])	# x > 0 となる個数
nm <- length(x[x < 0])	# x < 0 となる個数
x	# データ
binom.test(c(np,nm))$p.value	# 符号検定の p 値
wilcox.test(x)$p.value	# Wilcoxon の順位和検定 p 値
t.test(x)$p.value	# t 検定の p 値

次に，データが正規分布に従っていない場合を考える．自由度 d = 10 と 2 の t 分布 とした場合の符号検定と Wilcoxon の順位和検定， t 検定を行い，それぞれの p 値を出す．

x0 <- x	# 標準正規分布
x <- rt(n, df=2) + m	# 自由度 2 の t 分布
xx <- data.frame(cbind(x0, x))	#
names(xx) <- c("N(0, 1)", "t(2)")	# 変数名の定義
stripchart(xx, ylim=c(0.5,2.5))	# 乱数の分布範囲の表示
title(main="乱数標本の例")	#

これを N = 10000 回行い，p 値が 5％や 1％以下になる個数を出す．これより， 検定方式による検出力（検定力）がシミュレーションにより求められる．
　以下に，乱数のサンプルサイズ n と平均 m と t 分布のときの自由度 d から N = 10000回のシミュレーションにより， 検定ごとに p 値が 5％や 1％以下になる頻度を算出する関数を定義する．

srst.f <- function(m, n, d=NULL){	# 関数定義：平均 m，サンプル数 n，自由度 d
N <- 10000	# シミュレーション回数
yb <- NULL; yw <- NULL; yt <- NULL	# p 値格納ベクトルの定義
for(i in 1:N){	#
if(is.null(d)) x <- rnorm(n, mean=m)	# d の値がないときは正規乱数
else x <- rt(n, df=d) + m	# そうでないときは自由度 d の t 分布乱数 + m
np <- length(x[x > 0])	# 乱数で x > 0 の個数
nm <- length(x[x < 0])	# 乱数で x < 0 の個数
yb <- c(yb, binom.test(c(np,nm))$p.value)	# 符号検定の p 値を順に格納
yw <- c(yw, wilcox.test(x)$p.value)	# Wilcoxon の順位和検定 p 値を順に格納
yt <- c(yt, t.test(x)$p.value)	# t 検定の p 値を順に格納
}	#
b <- length(yb[yb<0.05])/N	# 符号検定の p 値が 0.05 未満の個数
w <- length(yw[yw<0.05])/N	# Wilcoxon の順位和検定 p 値が 0.05 未満の個数
t <- length(yt[yt<0.05])/N	# t 検定の p 値が 0.05 未満の個数
q5 <- c(t, w, b)	#
b <- length(yb[yb<0.01])/N	# 符号検定の p 値が 0.01 未満の個数
w <- length(yw[yw<0.01])/N	# Wilcoxon の順位和検定 p 値が 0.01 未満の個数
t <- length(yt[yt<0.01])/N	# t 検定の p 値が 0.01 未満の個数
q1 <- c(t, w, b)	#
return(rbind(q5, q1))	#
}	#

　まず，母集団が正規分布のとき，平均を 0 から 0.2 ずつ 2 まで変化させて，検定 方式ごとに有意水準 5 ％検定と 1 ％検定を行ったときの棄却確率をグラフ化すると， 検出力の変化と違いがよくわかる．以下のスクリプトを少し変えて，引き続き母集団が t 分布のときの検出力の グラフを出してみよう．

mv <- seq(0, 2, by=0.2)	# 母集団平均の定義
resnorm5 <- NULL	# 5 ％検定を行ったときの棄却確率行列
resnorm1 <- NULL	# 1 ％検定を行ったときの棄却確率行列
for(i in 1:length(mv)){	#
s <- srst.f(mv[i], 10)	# 平均 mv[i] の正規乱数（t 分布のときは，srst.f(mv[i], 10, d) とする）
resnorm5 <- rbind(resnorm5, s[1,])	#
resnorm1 <- rbind(resnorm1, s[2,])	#
}	# 行列データの同時プロット
matplot(mv, resnorm5, type="b", pch=1:3, ylim=c(0,1), xlab="真の平均", ylab="棄却確率")
abline(h=0.05); abline(h=1)	# 有意水準（0.05）と確率 1 のライン
legend(1.1,0.4,c("t 検定","W 検定","符号検定"), pch=1:3,col=c("black","red","green"))
title(main="5 ％検定検出力の比較，N(0, 1)")	#

真の平均 m = 1 のときの母集団分布と検定法による検出力の比較

	N(0, 1)	t(10)	t(5)	t(2)
t 検定	0.7989	0.7296	0.6395	0.4135
Wilcoxon の順位和検定	0.7830	0.7187	0.6391	0.4616
符号検定	0.5121	0.4761	0.4370	0.3428

　シミュレーションの結果，符号検定の検出力は予想されたように，かなり低い．一方， t 検定と Wilcoxon の順位和検定では検出力がほとんど変わらない．Wilcoxon の順位和検定のロバスト性が 生きてくるのは自由度 2 の t 分布という相当裾の重い分布になってからである．裾がそれほど重くない ときは，t 検定の方が検出力がわずかとはいえ高いので，t 検定を使用してよいように思える． 標本の大きさ n を増やしたときも調べてみるとよいだろう．
　データ分布の左右の対称性の「ずれ」に対する検出力の違いは，ガンマ分布などを使用したシミュレーションが 必要であろう．

３．２標本問題（two sample problem）

3-1．概要

3-2．Mann-Whitney の U 検定

3-3．Kolmogorov-Smirnov の 2 標本検定

表３．Kolmogorov-Smirnov 2 標本両側検定の臨界値

確率：γ	0.99	0.95	0.90	0.85	0.80
臨界値：kγ	1.63	1.36	1.22	1.14	1.07

3-4. Behrens-Fisher 問題

3-5. Brunner-Munzel 検定

3-6．解析例

データダウンロード

unagi <- read.csv("unagi2.csv"); unagi  	#データ入力
boxplot(unagi)
x <- unagi$noNa; y <- unagi$Na

Mann-Whitney の U 検定

wilcox.test(x, y)

Welch 検定

t.test(x, y)

plot(ecdf(x), do.points=F, verticals=T,xlab="生理活性",
  ylab="累積確率", xlim=range(x, y), main="")
plot(ecdf(y), do.points=F, verticals=T, add=T, col.h='red', col.v='red')
title(main="ウナギ生理活性の経験累積分布")
legend("bottomright", legend=c("Na","noNa"), lty=c(1,1),col=c("red","black"))
ks.test(x, y)

Brunner-Munzel 検定

library(lawstat)
brunner.munzel.test(x, y)

3-7. 検定方法の比較

母分散が等しい場合

X 〜 N(μ, σ²), Y 〜 N(μ, σ²)

N <- 10000			# シミュレーション回数
Welch <- U <- KS <- BM <- numeric(0)
for(i in 1:N){
	x <- rnorm(10, sd=1)
	y <- rnorm(30, sd=1)
	Welch <- c(Welch, t.test(x, y)$p.value)
	U <- c(U, wilcox.test(x ,y)$p.value)
	KS <- c(KS, ks.test(x ,y)$p.value)
	BM <- c(BM, brunner.munzel.test(x, y)$p.value)
}
length(Welch[Welch<0.05])
length(U[U<0.05])
length(KS[KS<0.05])
length(BM[BM<0.05])

10000回の検定で有意になった回数

Welch 検定	U 検定	KS 検定	BM 検定
518	494	336	564

母分散が異なる場合

X 〜 N(μ, σ_x²), Y 〜 N(μ, σ_y²)

N <- 10000			# シミュレーション回数
Welch <- U <- KS <- BM <- numeric(0)
for(i in 1:N){
	x <- rnorm(10, sd=2)
	y <- rnorm(30, sd=1)
	Welch <- c(Welch, t.test(x, y)$p.value)
	U <- c(U, wilcox.test(x ,y)$p.value)
	KS <- c(KS, ks.test(x ,y)$p.value)
	BM <- c(BM, brunner.munzel.test(x, y)$p.value)
}
length(Welch[Welch<0.05])
length(U[U<0.05])
length(KS[KS<0.05])
length(BM[BM<0.05])

10000回の検定で有意になった回数

Welch 検定	U 検定	KS 検定	BM 検定
498	980	1357	533

母平均が異なる場合

X 〜 N(μ_x, σ_x²), Y 〜 N(μ_y, σ_y²)

N <- 10000			# シミュレーション回数
Welch <- U <- KS <- BM <- numeric(0)
for(i in 1:N){
	x <- rnorm(10, mean=1, sd=2)
	y <- rnorm(30, sd=1)
	Welch <- c(Welch, t.test(x, y)$p.value)
	U <- c(U, wilcox.test(x ,y)$p.value)
	KS <- c(KS, ks.test(x ,y)$p.value)
	BM <- c(BM, brunner.munzel.test(x, y)$p.value)
}
length(Welch[Welch<0.05])
length(U[U<0.05])
length(KS[KS<0.05])
length(BM[BM<0.05])

10000回の検定で有意になった回数

Welch 検定	U 検定	KS 検定	BM 検定
2818	4012	4248	2773

１元配置のノンパラメトリック検定

• 品種によるコメの収量の違い
  水稲の９品種をそれぞれ，６区画の水田で栽培したときのアール当たりの玄米重量は以下のようであった． このうち，A，B，D，それぞれ，同じ母本から育成された品種であり，C は標準（対照 control）品種である．
  このデータは一元配置分散分析で解析できる．処理が品種で，9 水準からなっている．帰無仮説は，
  H₀：収量はどの品種も同じである
  である．

  データダウンロード

  yy <- read.csv("hinsyu.csv"); yy
  n <- nrow(yy)		# 繰り返し数
  a <- ncol(yy)		# 処理数 
  y <- as.vector(as.matrix(yy))		# データのベクトル化
  x <- factor(rep(1:a, c(rep(n,a))))	# データのラベル
  cbind(y, x)
  summary(aov(y ~ x))	# 一元配置分散分析
  #
  kruskal.test(y, x)	# クラスカル・ウォリス順位和検定                            

参考文献

1. Mood A. M., Graybill F. A. & Bose D. C. (1974). Introduction to the Theory of Statistics, McGraw-Hill.
2. Sokal R. R., & Rohlf F. J. (1981). Biometry, 2nd. ed., W. H. Freeman & Company.
3. Brunner E. and Munzel U. (2000). The nonparametric Behrens-Fisher problem: Asmptotic theory and a small-sample approximation, Biometrical Journal, Vol. 42, 17-25.
4. 石井進(1971)．"生物統計学入門"，培風館．
5. 鈴木修（編），大森宏，児玉良治，渡辺俊之，矢野広，山根健治（2004）． 専修学校における園芸療法士教育育成システムの研究開発 （文部科学省委託平成15年度「専修学校先進的教育研究開発事業」）
6. 鈴木修，渡辺俊之，矢野広，山根健治，大森宏，伊東正信，最上正秀，山下容子，小泉力， 児玉良治，頭士智美，細井薫，水口聡子，遠藤久子，樋田奈穂子，小島ユリ，郡司敏幸（2005）． 福祉サービス提供者に対する園芸療法教育システムの研究開発 （平成16年度文部科学省「専修学校社会人キャリアアップ教育推進事業」）
7. 竹村彰通(1991)．"現代数理統計学"，創文社．