2017.11.14

数理統計学演習

ノンパラメトリック検定

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

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

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

(1)

(2)

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

2-1．概要

(5)

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

(7)

(8)

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

(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．連検定（run test）

　なお，R で連検定を行うには，"tseries" (Time series analysis and computational finance) ライブラリーをインストールしておけば利用できる．

library(tseries)
# y， x， x， y， x， y， y， y， x
a <- c(0,1,1,0,1,0,0,0,1)
a <- factor(a)
runs.test(a)

3-3．Mann-Whitney の U 検定

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

3-4．解析例

データ分布

boxplot(mse$before, mse$after, names= c("対照区","処理区"), ylab="MSE 得点", cex.axis=0.8)

連検定

nrow(mse)
score <- c(mse$before, mse$after)
treat <- factor(rep(c(0,1), each=nrow(mse)))
o <- order(score)
score[o]
treat[o]
library(tseries)
runs.test(treat[o])

Mann-Whitney の U 検定

wilcox.test(mse$after, mse$before)
wilcox.test(score ~ treat)
library(coin)
wilcox_test(score ~ treat)

t 検定

t.test(mse$after, mse$before, var.equal=T)
t.test(score ~ treat, var.equal=T)

Kolmogorov-Smirnov 検定

plot(ecdf(mse$before), do.points=F, verticals=T,xlab="MSE 得点",
  ylab="累積確率", xlim=range(mse$after,mse$before), main="")
plot(ecdf(mse$after), do.points=F, verticals=T, add=T, col.h='red', col.v='red')
title(main="MSE 得点の経験累積分布")
legend(locator(1), legend=c("処理","対照"), lty=c(1,1),col=c("red","black"))
ks.test(mse$after, mse$before)

「若年層」のデータ分布

boxplot(mse$before[young], mse$after[young], names= c("対照区","処理区"), ylab="MSE 得点", cex.axis=0.8)
title(main="高齢でない層の MSE の得点分布", cex.main=1.0)

連検定

young <- which(mse$age<86)
yscore <- c(mse$before[young], mse$after[young])
ytreat <- factor(rep(c(0,1), each=length(young)))
oy <- order(yscore)
yscore[oy]
treat[oy]
library(tseries)
runs.test(treat[oy])

Mann-Whitney の U 検定

wilcox.test(mse$after[young], mse$before[young])
library(coin)
wilcox_test(yscore ~ ytreat)

t 検定

t.test(mse$after[young], mse$before[young], var.equal=T)
t.test(yscore ~ ytreat, var.equal=T)

Kolmogorov-Smirnov 検定

plot(ecdf(mse$before[young]), do.points=F, verticals=T,xlab="MSE 得点",
  ylab="累積確率", xlim=range(mse$after[young],mse$before[young]), main="")
plot(ecdf(mse$after[young]), do.points=F, verticals=T, add=T, col.h='red', col.v='red')
title(main="高齢でない層の MSE 得点の経験累積分布", cex.main=0.8)
legend(locator(1), legend=c("処理","対照"), lty=c(1,1),col=c("red","black"))
ks.test(mse$after[young], mse$before[young])

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

• 品種によるコメの収量の違い
  水稲の９品種をそれぞれ，６区画の水田で栽培したときのアール当たりの玄米重量は以下のようであった． このうち，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. 石井進(1971)．"生物統計学入門"，培風館．
4. 鈴木修（編），大森宏，児玉良治，渡辺俊之，矢野広，山根健治（2004）． 専修学校における園芸療法士教育育成システムの研究開発 （文部科学省委託平成15年度「専修学校先進的教育研究開発事業」）
5. 鈴木修，渡辺俊之，矢野広，山根健治，大森宏，伊東正信，最上正秀，山下容子，小泉力， 児玉良治，頭士智美，細井薫，水口聡子，遠藤久子，樋田奈穂子，小島ユリ，郡司敏幸（2005）． 福祉サービス提供者に対する園芸療法教育システムの研究開発 （平成16年度文部科学省「専修学校社会人キャリアアップ教育推進事業」）
6. 竹村彰通(1991)．"現代数理統計学"，創文社．