独立２群の差の検定
1 二標本 t 検定
このページは、市原清志先生の「バイオサイエンスの統計学」の問題を、統計ソフト”R”で解説するものです。
詳細については、市原先生のテキスト（p72〜P78）を参照ください。
数学的説明も補足してありますが、市原先生のテキストにない部分は飛ばしてもよいです。
1.1 ［基本］
２群の等分散性、正規性、独立性を前提とする。（等分散が仮定できなければ
welch検定）
実際には、次の手順で行うことが多い。
\(X_1\sim
N(\mu_1,\sigma_1^2)\:,\:X_2\sim N(\mu_2,\sigma_2^2)\::i.i.d.\)
は独立な２群を表す確率変数とする。
独立な２群の標本数を \(n_1,n_2\)
、母分散を \(\sigma_1^2,\sigma_2^2\)
、不偏分散を \(s_1^2,s_2^2\)
とする。
（1）等分散性の検定を行う。
＜帰無仮説＞ \(\sigma_1^2=\sigma_2^2\) のもとでは、 \(F=\dfrac{s_1^2}{s_2^2}\sim
F(n_1-1,n_2-1)\)
（市原p75によると、不偏分散の大きい方を分子とした片側検定である）
（2）等分散が棄却されなければ、２群の平均に差があるかどうか t
検定を行う。
\(X_1\sim N(\mu_1,\sigma^2)\:,\:X_2\sim
N(\mu_2,\sigma^2)\) 、ただし \(\sigma^2=\sigma_1^2=\sigma_2^2\)
\(\overline{X_1}\sim N\left(\mu_1,\dfrac{\sigma^2}{n_1}\right)\:,\:\overline{X_2}\sim N\left(\mu_2,\dfrac{\sigma^2}{n_2}\right)\) より
\(\overline{X_1}-\overline{X_2}\sim
N\left(\mu_1-\mu_2\:,\:\dfrac{\sigma^2}{n_1}+\dfrac{\sigma^2}{n_2}\right)\)
＜帰無仮説＞ ２群の平均に差はないのもとでは、
\(\overline{X_1}-\overline{X_2}\sim
N\left(0\:,\:\dfrac{\sigma^2}{n_1}+\dfrac{\sigma^2}{n_2}\right)\)
２群をプールした偏差平方和は \(s_1^2(n_1-1)+s_2^2(n_2-1)\)
自由度は \((n_1-1)+(n_2-1)=n_1+n_2-2\)
２群をプールした不偏分散は \(s^2=\dfrac{s_1^2(n_1-1)+s_2^2(n_2-1)}{n_1+n_2-2}\) となるので
検定量 \(z=\dfrac{\overline{X_1}-\overline{X_2}}{\sigma\sqrt{\dfrac{1}{n_1}+\dfrac{1}{n_2}}}\sim
N(0\:,\:1)\) \(t=\dfrac{\overline{X_1}-\overline{X_2}}{s\sqrt{\dfrac{1}{n_1}+\dfrac{1}{n_2}}}\sim
t(n_1+n_2-2)\)
母分散は未知なので、自由度 (\((n_1+n_2-2)\) の
t検定を行う。（両側検定）
1.2 ［例3-1］バセドウ病患者（p78）
健常人5名（N群）とバセドウ病患者7名（B群）について、糖負荷検査の際の血糖上昇値に差があるかどうか検定する。
Normal <- c(54,49,42,40,35)
Basedow <- c(68,65,60,56,52,47,44)
boxplot(Normal,Basedow)
箱ひげ図（boxplot）によって2群を比較すると、2群には差があるように見える。
この差が有意かどうかRで調べる。
1.2.1 等分散性検定
帰無仮説 「分散が等しい」のもとで、
F=(Basedowの不偏分散)/(Normalの不偏分散) は自由度(6,4)のF分布に従う。
Fの値が１から離れていれば、等分散の仮説は棄却される。
（市原p80では、大きい分散を分子にとり、F>1の状態で片側検定を行う）
var.test(Basedow,Normal) ##
## F test to compare two variances
##
## data: Basedow and Normal
## F = 1.4218, num df = 6, denom df = 4, p-value = 0.7646
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.1545918 8.8539578
## sample estimates:
## ratio of variances
## 1.421829出力３行目の numはnumerator（分子）、denomはdenominator（分母）
分子の自由度が6、分母の自由度が4であり、F値=1.4218
は１に近いので分散が異なるとは考えにくい。
F値が1.4218以上になる確率は 0.7646 で
0.05より大きいので、帰無仮説は棄却されない。
従って、等分散という仮説を棄却しないまま進む。
即ち、等分散を前提として進むことになる。
［注意］
もし、var.test(Normal,Basedow)
とすれば、分母と分子が逆に計算され、F値は１より小さくなるが、検定の結果は同じになる。
正規性については、市原p77「t分布の頑強性」にあるように、厳しく検定する必要はない。
（実用上、問題は起こらないということである）
1.2.2 ２群のt検定
等分散の仮説のもとで、不偏分散として2群をプールした分散 \(s^2\)を求め、
\(t=\dfrac{\overline{X_1}-\overline{X_2}}{s\sqrt{\dfrac{1}{n_1}+\dfrac{1}{n_2}}}\sim
t(n_1+n_2-2)\) の値が 0 から離れているかで検定を行う。
t.test(Normal, Basedow, paired=FALSE, alternative ="two.sided",var.equal=T)##
## Two Sample t-test
##
## data: Normal and Basedow
## t = -2.4356, df = 10, p-value = 0.03511
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -22.977804 -1.022196
## sample estimates:
## mean of x mean of y
## 44 56“paired=FALSE” は対応のないt検定であることを表し、”two.sided”
は両側検定を表す。
“var.equal=T” は等分散を仮定した検定であることを意味する。
- 自由度は \((n_1+n_2-2)=5+7-2=10\)
- t値=-2.4356 は0から離れた異常値と見える。
t 値が、この値かそれ以上に0から離れる確率は、p-value = 0.03511 この確率は基準とする0.05より小さく、t = -2.4356は異常値と考えられる。
帰無仮説「２群の平均に差はない」は棄却され、
対立仮説（alternative hypothesis）「 平均の差が0ではない」を真と判断する。
1.3 例3-3 Welchの検定 (p81)
等分散を仮定しない場合は、Welchの検定を行う。
これは、市原p73によれば、自由度を少し小さめに見積もることによる微調整と考えられる。
簡便のため２群の標準偏差の差を大きめにとり、
男 sd=20 女 sd=35 とする。
male<-rnorm(30,120,20)
female<-rnorm(40,135,35)
boxplot(male,female)
まず、等分散仮説を検定する。
var.test(female,male) ##
## F test to compare two variances
##
## data: female and male
## F = 2.3511, num df = 39, denom df = 29, p-value = 0.01889
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 1.156672 4.612605
## sample estimates:
## ratio of variances
## 2.351128等分散仮説は棄却される。
次にWelchの検定を行う。
Rのコードで、”var.equal=F” とすればWelchの検定を行ってくれる。
t.test(male, female, paired=FALSE, alternative ="two.sided",var.equal=F)##
## Welch Two Sample t-test
##
## data: male and female
## t = -0.55345, df = 66.86, p-value = 0.5818
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -16.176323 9.153262
## sample estimates:
## mean of x mean of y
## 121.8392 125.3507帰無仮説 「男女に差はない」 は棄却され、
対立仮説 「男女に差がある」 が採択される。
以上は正規性を前提条件とするパラメトリック検定である。
データの正規性が保証されない場合はノンパラメトリック検定を用いる。
次に紹介するのはノンパラメトリックな方法である。
市原先生のテキストではMann-Whitney検定が紹介されているが、ウイルコクスン順位和検定と同等である。
ここではまずウイルコクスン順位和検定の方法を説明する。
2 ウイルコクスン順位和検定
2.1 Wilcoxon rank sum test
比較したい２標本の観測値に全体順位（Rank）をつけて、標本ごとの合計（順位和）を算出。 観測値に同一順位がある場合は、順位の平均値をつける。
標本 \(A,B\) の順位和を \(R_A,R_B\) とすると、統計量は
\(W = min \left( R_{A},\ R_{B}
\right)\) （実際はどちらをWとしてもよい）
2.1.1 シミュレーション
平均に差がないデータ X1,X2 を1000回出力し、ウイルコクスン順位和 W
がどのような分布になるかシミュレーションする。
平均が10、標準偏差が1.5 の2つの群 X1,X2 を出力し
Wを計算するという試行を1000回繰り返す。
X1の順位和を W として計算する。
\(X1\:,\:X2\) のデータ数を \(n_1,n_2\) とすると第１群の順位和\(W_1\)は、
\(\displaystyle
E(W_1)=\dfrac{n_1(n_1+n_2+1)}{2}\) \(V(W_1)=\dfrac{n_1n_2(n_1+n_2+1)}{12}\)
導出は当サイト Wilcoxon
を参照（数列の計算がやや難）。
X1 = c()
X2 = c()
X=c()
for (i in 1:1000) {
x1<-rnorm(5,10,1.5);x2<-rnorm(4,10,1.5)
x3 <- c(x1,x2)
x <- sum(rank(x3)[1:5])
X <- c(X,x)
}期待値と分散は \(\displaystyle E(W_1)=\dfrac{5\cdot10}{2}=25\) \(V(W_1)=\dfrac{5\cdot4\cdot10}{12}=16.66\) （理論値）
mean(X)## [1] 24.938var(X)## [1] 17.32148hist(X,prob=T)
curve(dnorm(x,mean(X),sqrt(var(X))),add=TRUE,col="red")
1000回のヒストグラムと理論曲線は、ほぼ一致している。
２群の平均が異なる場合を1回出力すると、
x1<-rnorm(5,10,1.5);x2<-rnorm(4,8,1)
x3 <- c(x1,x2)
sum(rank(x3)[1:5])## [1] 30Wの値は棄却域に近くなってくる。
3 Mann-Whitney検定
Wilcoxonの順位和検定は、「マン・ホイットニーU検定」と同等であり、検定量Uは両群で
U=\(\displaystyle \sum_i\) ( \(x_i\) より大きい他方の群のデータ数 )
を考え値の小さい方を選ぶ。(市原p90)
\(U=W_1-\dfrac{n_{1}\left(
n_{1}+1\right) }{2}\)
であることが導かれる。
導出は当サイト Wilcoxon
参照。
例3-4を想定してシミュレーションする。
3.1 例3-4（p92）
X1 = c();X2 = c()
X=c()
for (i in 1:1000) {
x1<-rnorm(5,10,1.5);x2<-rnorm(4,10,1.5)
x3 <- c(x1,x2)
x <- sum(rank(x3)[1:5])-5*6/2
X <- c(X,x)
}理論値は \(E[ U] =\dfrac{n_{1} n_{2}}{2}\:\:,\:\:V(U)=\dfrac{n_{1}n_{2}( n_1+n_2+1) }{12}\)
mean(X)## [1] 9.993var(X)## [1] 16.87182hist(X,prob=T)
curve(dnorm(x,mean(X),sqrt(var(X))),add=TRUE,col="red")
ヒストグラムはテキストp92の図とほぼ一致している。
3.2 例3-5（p93）
例3-5を想定しシミュレーションを行う。
テキストでは個数の少ない２群を用いているが、
\(E[ U] =\dfrac{n_{1}
n_{2}}{2}\:\:,\:\:V(U)=\dfrac{n_{1}n_{2}( n_1+n_2+1)
}{12}\)
から分かるように２群のどちらの順位和を用いてもよい。
理論値は
\(E[ U] =\dfrac{n_{1}
n_{2}}{2}=\dfrac{13\cdot7}{2}=45\:\:,\:\:V(U)=\dfrac{n_{1}n_{2}(
n_1+n_2+1) }{12}=151.6\)
ここでは２群を用いて計算する。
\(U=W_2-\dfrac{n_{2}\left( n_{2}+1\right) }{2}\)
X1 = c();X2 = c()
X=c()
for (i in 1:1000) {
x1<-rnorm(13,1.7,0.5);x2<-rnorm(7,1.7,0.5)
x3 <- c(x1,x2)
x <- sum(rank(x3)[14:20])-7*8/2
X <- c(X,x)
}
mean(X)## [1] 45.774var(X)## [1] 141.6145hist(X,prob=T)
curve(dnorm(x,mean(X),sqrt(var(X))),add=TRUE,col="red")
これは２群に差がない場合のヒストグラムである。
２群に差がある場合を1回出力すると、
x1<-rnorm(13,1.7,0.5);x2<-rnorm(7,3,0.8)
x3 <- c(x1,x2)
sum(rank(x3)[14:20])-7*8/2## [1] 76有意な値が出力される。
１群を用いた場合も出力する。
X1 = c();X2 = c()
X=c()
for (i in 1:1000) {
x1<-rnorm(13,1.7,0.5);x2<-rnorm(7,1.7,0.5)
x3 <- c(x1,x2)
x <- sum(rank(x3)[1:13])-13*14/2
X <- c(X,x)
}
mean(X)## [1] 45.675var(X)## [1] 160.1856hist(X,prob=T)
curve(dnorm(x,mean(X),sqrt(var(X))),add=TRUE,col="red")
x1<-rnorm(13,1.7,0.5);x2<-rnorm(7,3,0.8)
x3 <- c(x1,x2)
sum(rank(x3)[1:13])-13*14/2## [1] 0やはり有意な値が出力される。
Rの関数を使うには wilcox.test を用いる。
wilcox.test(x1, x2)##
## Wilcoxon rank sum exact test
##
## data: x1 and x2
## W = 0, p-value = 2.58e-05
## alternative hypothesis: true location shift is not equal to 0あまりに簡単すぎて何をやっているのか分からない。
やはり学習するときにはシミュレーションでイメージを掴むと良いと思う。