Studentのt検定が使える条件

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

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

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

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

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

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

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

検出力

　標本が独立ではなく、正規分布ではないので、普通の2標本のt検定の要請は満たされていない。また、Mann-WhitneyのU検定は標本の独立性を仮定しているため、使えない。そこで、対応のあるt検定のノンパラメトリック版としてWilcoxonの符号順位検定を紹介しよう。Man-WhitneyのU検定の別名、Wilcoxonの順位和検定とよく似た名前だが、別物なので注意しよう。しかしながら、先に言うとこのWilcoxonの符号順位検定は使いどころが限られる。多くのデータの場合、対応のあるt検定で事足りる可能性が高い。さて、Wilcoxonの符号順位検定の検定統計量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検定は符号検定と比較してあてにならないということになる。

Rで行う符号順位検定

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

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

m1 <- c(42, 34, 63, 37, 86, 80, 78, 38, 49, 55,

49, 40, 56, 59, 62, 93, 72, 7, 76, 42)

m2 <- c(54, 47, 63, 44, 88, 137, 83, 51, 32, 50,

51, -36, 60, 74, 66, 99, 69, 5, 77, 40)

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の符号順位検定を行う。正確なp値の算出のために、exactRankTestsパッケージのwilcoxon.exact()関数を用いるのがおすすめだ。下記のようにpaired = Tとすることで、対応のある2標本の検定を行ってくれるが、あらかじめm <- m1 - m2と差をとり、それを1標本の検定としてもよい。

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

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

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

fit

##

## Exact Wilcoxon signed rank test

##

## data: m1 and m2

## V = 58.5, p-value = 0.1467

## alternative hypothesis: true mu 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] 58.5

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

統計量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] 12285

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

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

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

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

## [1] 0.08256929

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

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

まとめ

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