目次

・Studentのt検定が使える条件

・独立でもないし、正規分布でもない

・Wilcoxonの符号順位検定

・検出力

・危険率

・Rで行う符号順位検定

・まとめ

・追記1：検定統計量W+の平均と分散―同順がないとき―

・追記2：検定統計量W+の平均と分散―0と同順があるとき―

　　・Wilcoxon法

　　・Pratt法

Studentのt検定が使える条件

　2標本のStudentのt検定は、母集団に次のような過程が成り立つときに有効に使うことができる。

1. 等分散性……2標本の分散が等しい

2. 独立性……2標本は独立にサンプリングされている

3. 正規性……2標本は正規分布に従う

では上記の3条件が崩れたときには、どれくらい信頼をおけない解析となってしまうのだろうか？　本稿では、シミュレーションを通して、t検定の頑健性と、仮定が満たされないときの対策を紹介する。

独立でもないし、正規分布でもない

　これまで、Studentのt検定を行うために必要な3要素、等分散性、独立性および正規性について検討し、それぞれに対する対応策を紹介した。本稿では、もっと極端な例として、標本が独立でもないし、正規分布でもない状況を考える。

Wilcoxonの符号順位検定

　標本が独立ではなく、正規分布ではないので、普通の2標本のt検定の要請は満たされていない。また、Mann-WhitneyのU検定は標本の独立性を仮定しているため、使えない。そこで、対応のあるt検定のノンパラメトリック版としてWilcoxon（ウィルコクソン）の符号順位検定 Wilcoxon signed-rank testを紹介しよう。Man-WhitneyのU検定の別名、Wilcoxonの順位和検定とよく似た名前だが、別物なので注意しよう。しかしながら、先に言うとこのWilcoxonの符号順位検定は使いどころが限られる。多くのデータの場合、対応のあるt検定で事足りる可能性が高い。さて、Wilcoxonの符号順位検定の検定統計量W+は以下のように計算する。0に順位をつけない方法はWilcoxon法と呼ばれる（0に順位をつける方法はPratt法、追記、こちらも参照）。

　このW+は、以下のような平均と分散をもつ（証明は追記）。

検出力

では、検討していきたいのだが、独立性の検討の時に、独立ではない標本m22は、正規分布に従う標本に、さらに正規分布を足す、という方法で作成した。

------------------------------------------------------

m11 <- rnorm(30, mean = 50, sd = 20)

m22 <- m11 + rnorm(length(m11), mean = 10, sd = 2)

------------------------------------------------------

この方法で作成した、m11とm22はともにほとんど正規分布だったため、対応のあるt検定で解析したのだった。もし、独立ではない標本において、正規分布でなくするとしたらどこだろう？　まずは、足す前の標本m11を正規分布ではない形にすることを考える。例えば、正規性で紹介したコーシー分布にするとしたら下記のようになる。

------------------------------------------------------

m11 <- rcauthy(30, location = 50)

m22 <- m11 + rnorm(length(m11), mean = 10, sd = 2)

------------------------------------------------------

確かにこの２標本m11とm22はコーシー分布に従い、正規性はない（自分で図を描いてみよう）。一見すると、対応のあるt検定に出番はないように見える。ところで、対応のあるt検定では対応する観測値で差をとる。すると、上記の2標本の差は正規分布になり、普通の対応のあるt検定と同じになってしまう。つまり、足す前の分布が正規分布でなくても、差が正規分布に従いそうなら、対応のあるt検定で問題ない。ゆえに問題となるのは、差が正規分布に従わないときだ。なので、片方の標本は正規分布とし、独立ではない標本はその標本にコーシー分布を足すことにする。つまり、差の分布がコーシー分布になる、ということだ。

------------------------------------------------------

library(plotn)

library(exactRankTests)

#exactRankTestsパッケージがない場合は、下記でインストール

#install.packages("exactRankTests", repos="http://cran.ism.ac.jp/")

n1 <- 30

t2 <- NULL

w3 <- NULL

p2 <- NULL

p3 <- NULL

for (i in 1:10000){

  m1 <- rnorm(n1, mean = 50, sd = 20)

  m2 <- m1 + rcauchy(n1, location = 10)

  m <- m1 - m2

  t2 <- c(t2, (mean(m) - 0)/(sd(m)/sqrt(30)))

  p2 <- c(p2, t.test(m1 - m2)$p.value)

  d <- data.frame(m = m)

  d <- d[order(abs(d$m)), , drop = F]

  d$rank <- 1:nrow(d)

  w1 <- sum(d[d$m > 0,"rank"])

  w2 <- sum(d[d$m < 0,"rank"])

  w3 <- c(w3, min(w1,w2))

  p3 <- c(p3, wilcox.exact(m1,m2, paired = T)$p.value)

}

histn(t2, xlab = "統計量t", ylab = "頻度")#図1の描画

histn(p2, xlab = "P value", ylab = "頻度")#図2の描画

histn(w3, xlab = "統計量W", ylab = "頻度")#図3の描画

histn(p3, xlab = "P value", ylab = "頻度")#図4の描画

sum(p2 < 0.05)

## [1] 8716

sum(p3 < 0.05)

## [1] 10000

------------------------------------------------------

t検定の検出力は87%、符号検定の検出力は100%となった。正規性を議論した時と同様に、図1を見ると統計量tの分布が明らかにt分布ではない。したがって、このような状況の場合、t検定は全くあてにならない。差の分布に外れ値が生じやすい状況では、解析方法に注意が必要となる（その一方で、統計量Wの分布もかなり変な分布ではあるが……）。

危険率

　危険率についても同様に検討してみよう。

------------------------------------------------------

n1 <- 30

t2 <- NULL

w3 <- NULL

p2 <- NULL

p3 <- NULL

for (i in 1:10000){

  m1 <- rnorm(n1, mean = 50, sd = 20)

  m2 <- m1 + rcauchy(n1, location = 0)

  m <- m1 - m2

  t2 <- c(t2, (mean(m) - 0)/(sd(m)/sqrt(30)))

  p2 <- c(p2, t.test(m1 - m2)$p.value)

  d <- data.frame(m = m)

  d <- d[order(abs(d$m)), , drop = F]

  d$rank <- 1:nrow(d)

  w1 <- sum(d[d$m > 0,"rank"])

  w2 <- sum(d[d$m < 0,"rank"])

  w3 <- c(w3, min(w1,w2))

  p3 <- c(p3, wilcox.exact(m1,m2, paired = T)$p.value)

}

histn(t2, xlab = "統計量t", ylab = "頻度")#図5の描画

histn(p2, xlab = "P value", ylab = "頻度")#図6の描画

histn(w3, xlab = "統計量W", ylab = "頻度")#図7の描画

histn(p3, xlab = "P value", ylab = "頻度")#図8の描画

sum(p2 < 0.05)

## [1] 197

sum(p3 < 0.05)

## [1] 481

------------------------------------------------------

t検定の危険率は2%、U検定は5%程度であった。検出力の場合と同様、統計量tの分布はt分布となっていない。ゆえに、このような状況ではt検定は符号検定と比較してあてにならないということになる。

　特に、危険率を確認する場合、p値の分布に関する累積分布関数(CDF)を確認するのが良い。標本に差がないとき、期待されるp値のCDFはy=xの直線となる。それに対するt検定とU検定のCDFも確認しよう。

------------------------------------------------------

cdf <- function(p, v){

  sapply(v, function(x) sum(p < x))/length(p)

}

v <- seq(0, 1, length = 200)

p2c <- cdf(p2,v)

p3c <- cdf(p3,v)

cs <- data.frame(v,v,p2c,p3c)

plotn(cs[,1],cs[,2:4], mode = "m", xlab = "P-value", ylab = "CDF", type = "l",

      col.line = c("black","red","blue"), 

      pt.col.leg = c("black","red","blue"), 

      legend = T, line = T, pt.cex.leg = 0, 

      leg.lab = c("expected", "Student","Wilcoxon"), leg.sp = 8) #追加図1の描画

------------------------------------------------------

t検定は、期待されるp値の分布から大きく外れている。全体的に5%周りの頻度が低く、より高いp値が出やすい状況だ。つまり、過剰に保守的であり、解析に向かない。これは、平均値の検定であるt検定は、外れ値の影響を強く受け、平均値も分散もかなりぶれてしまう。それ反映して検出力が落ちていると期待される。一方で、符号順位検定は期待されるp値の分布に近いため、より適切であることがわかる。符号順位検定はデータを順位に変換して解析するので、外れ値の影響を受けにくい。

Rで行う符号順位検定

　上記のような背景のもと、Wilcoxonの符号順位検定を行ってみよう。以下のようなサンプルサイズ20の独立ではないデータがあるとする。有意水準は5%としよう。まずは、データを描画してみて様子を確認する。

------------------------------------------------------

m1 <- c(46, 82, 28, 52, 50, 12, 29, -1, 56, 50, 

        47, 45, 41, 37, 63, 64, 52, 84, 65, 68)

m2 <- c(-31, 236, 27, 89, 34, 6, 14, 44, 66, 

        -19, 90, 41, 48, 23, 65, 55, 82, 81, 57, 266)

m <- m1 - m2

d <- data.frame(d = c(m1, m2), g = c(rep("m1", length(m1)), rep("m2", length(m2))))

boxplotn(d$d ~ d$g, ylab = "data", jitter.method = "swarm")#図9の描画

------------------------------------------------------

今回は独立ではないデータなので2次元平面にプロットもしてみよう。独立かどうかの判断は、自分の実験設計と相談して考えよう。

------------------------------------------------------

plotn(m2 ~ m1) #図10の描画

------------------------------------------------------

独立ではなく、外れ値も存在しそうなので、Wilcoxonの符号順位検定を行う。今回は、wilcox.test関数を用いるが、0や同順が存在するときは正確なp値の算出のために、exactRankTestsパッケージのwilcox.exact関数を用いる。下記のようにpaired = Tとすることで、対応のある2標本の検定を行ってくれるが、あらかじめm <- m1 - m2と差をとり、それを1標本の検定としてもよい。

------------------------------------------------------

fit <- wilcox.test(m1, m2, paired = T, exact = T)

#fit <- wilcox.exact(m1, m2, paired = T, exact = T)

#fit <- wilcox.exact(m, mu = 0, exact = T)でも可

fit

## 

##  Wilcoxon signed rank exact test

## 

## data:  m1 and m2

## V = 96, p-value = 0.7562

## alternative hypothesis: true location shift is not equal to 0

------------------------------------------------------

Vがここでいう統計量W+のことである。下記のように計算する。

------------------------------------------------------

n <- length(m1)

m_0 <- m[!m == 0]

si <- sign(m_0)

abm <- abs(m_0)

r <- rank(abm)

w1 <- sum(r[si > 0])

w2 <- sum(r[si < 0])

w3 <- w1 #定義上ここで最小値をとるのであるが、上記関数ではVは最小値側を表示しない

w3

## [1] 96

------------------------------------------------------

統計量W+は、符号順位和の小さい方、と定義されているのであるが、なぜか上記関数のVは最小値側ではなく、正の順位和が表示される。別に検定上問題はないのであるが、こういう微妙な差異が混乱を生んでいる気がする。p値であるが正確な値を計算するのは難しい。なので下記のようにシミュレーションでW+分布を作って計算してみる。

------------------------------------------------------

n <- 20

w <- NULL

for(i in 1:300000){

  s <- rnorm(n, mean = 0, sd = 20)

  rs <- rank(abs(s[!s == 0]))

  ss <- sign(s[!s == 0])

  w <- c(w, sum(rs[ss > 0]))

}

histn(w, freq = F, xlab = "統計量W") #図11の描画

xx <- 0:400

yy <- dnorm(0:400, mean = n*(n+1)/4, sd = sqrt(n*(n+1)*(2*n+1)/24))

overdraw(points(xx, yy, type = "l", col = "red"))

sum(w < w3)

## [1] 108907

------------------------------------------------------

上記のように標本から計算された統計量W+よりも小さくなるW+は108907となり、この値の2倍を30万で割った値が近似的なp値となる。つまり、108907*2/300000 = 0.726である。やや、wilcox.test関数よりも小さな値だ。また、統計量W+はサンプルサイズが大きければ平均n×(n+1)/4、分散n×(n+1)×(2n+1)/24の正規分布に近似できることが知られており、これを使ってp値は次のように計算できる。

------------------------------------------------------

wilcox.test(m1, m2, paired = T, exact = F, correct = F)

## 

##  Wilcoxon signed rank test

## 

## data:  m1 and m2

## V = 96, p-value = 0.7369

## alternative hypothesis: true location shift is not equal to 0

pnorm(w3, mean = n*(n+1)/4, sd = sqrt(n*(n+1)*(2*n+1)/24),lower.tail = T) * 2

## [1] 0.7368754

------------------------------------------------------

p値の近似値計算は正規分布を使ったほうとはよく一致している。一般的に言われていることとしては、上記の値はexact = Fとしたときの漸近検定のp値と一致するはずだが、今回はexact = TでもFでも値は大きく変わらない。私もまだまだ勉強不足だ……。この値は、wilcox.test関数で、漸近検定（exact = F）とし、連続値の補正もなし（correct = F）としたときのp値と一致する。

まとめ

　このように対応のある標本の場合は、差の分布に注目して解析する必要がある。これまでの標本自体の分布ではないので注意が必要だ。一般的な量的データであれば、基本的には差の分布が極端に正規分布から外れるというのはない気がするが、カテゴリカルデータなら十分にあり得る。例えば、５段階評価を下してもらう実験などがそうだ。自分のデータに合わせた解析方法を見つけるのは、始めのうちはなかなか難しい。ちょっとずつ時間をかけて、「こういう解析ができたはず」というような解決辞書を豊かにしていくことで、徐々に状況に応じた対応策ができてくる。辛抱強く、自分と付き合っていこう。

追記1：検定統計量W+の平均と分散―同順がないとき―

　同順がないとき、統計量W+の平均と分散は、以下のように求める。

追記2：検定統計量W+の平均と分散―0と同順があるとき―

　実は、これらの取り扱いには2通りある。取り扱いとしては、Wilcoxon法のほうがシンプルだが、Pratt法は0になったサンプルも考慮する、という点でよりデータを活用している。

・Wilcoxon法

　まず、0が存在するときは、その分だけサンプルサイズから除く。そして、同順があるとき、平均は変わることがない。統計量W+の分散は、以下のように求める。

・Pratt法

　Pratt (1959, Journal of the American Statistical Association)は0を除くことによる検定結果の変化について指摘をし、0にも順位をつけて検定を行うことを提案した。

この時のW+の平均と分散は以下の通り。