「続・高校生の統計学」のプログラムコードを掲載しています。

  • RstudioやPosit Cloud.
  • GoogleColabのR環境.

においてコピー＆ペーストで利用できます。
初めは指導者のもとで行うのが良いでしょう。
（電子書籍はデジタル著作権管理（DRM）仕様のためコピーできません）

インストールせずにブラウザでプログラムを実行するには、次の方法があります。

以下、電子書籍との章番号を対応させるため、項目名は残してあります。

  • 必要なパッケージを用いるのに必要なコードで書籍に表示されていないものもありますが、以下では表示してあります。
    例えば、library(mvtnorm)、library(scatterplot3d).

1 はじめに.

2 パラメトリック検定.

2.1 ２標本 t 検定.

＜例＞
２つの学校A,Bの学力に差があるかを調べたいとします。
各校から大きさ40の標本を抽出し比較します。
シミュレーションでは、共通の平均・標準偏差を持つ正規乱数を用いてデータを出力し、boxplot（箱ヒゲ図）に表します。

  • 群G1（平均50, 標準偏差10）の正規分布からの大きさ40の標本 \(x_1\)
  • 群G2（平均50, 標準偏差10）の正規分布からの大きさ40の標本 \(x_2\)
x1 <- rnorm(40,50,10);x2 <- rnorm(40,50,10)
list=list(x1=x1, x2=x2)
boxplot(list,ylim=c(20,80))

母平均と母分散が同じであっても、誤差変動のため boxplot（箱ヒゲ図）にはズレが生じています。

今、平均差は.

mean(x1)-mean(x2)
## [1] -2.56542

同じ母集団からとった２群ですが、平均に差があります。

2.1.1 マハラノビスの距離.

2.1.1.1 標準正規分布.

標準正規分布は、平均０, 分散１の正規分布です。
標準正規分布に従うデータを1000個出力しヒストグラムを描きます。

x <- rnorm(1000,0,1)
hist(x,prob=T,breaks = "Scott",col="lightblue")
curve(dnorm(x,0,1),add=T,col="blue")

2.1.1.2 正規分布.

2.1.1.3 検定との関係.

2.1.2 平均差の確率分布.

2.1.2.1 期待値・分散の性質.

分散の性質についてシミュレーションで確かめてみましょう。

　　　x1：平均 3 、標準偏差 1 の正規乱数からの1000個のデータ.
　　　x2：平均 2 、標準偏差 1 の正規乱数からの1000個のデータ.
　　　x3=x1-x2.

x1 = rnorm(1000,mean=3,sd=1)
x2 = rnorm(1000,mean=2,sd=1)
x3 = x1-x2
par(family = "HiraKakuProN-W3") 
hist(x1,prob=T,col="pink",xlim=c(-2,12),main="差の確率分布",xlab="")
hist(x2,prob=T,col="lightblue",add=T)
hist(x3,prob=T,col="lightgreen",add=T)

\(x3=x1-x2\) のヒストグラムの幅は広がり、分散は大きくなっているようです。 分散の性質を確認します。

v <- c(var(x3),var(x3/2));v
## [1] 2.0714815 0.5178704

2.1.2.2 標本平均の分布.

2.1.2.3 標本平均の差の分布.

2.1.2.4 ２群の平均差の検定.

2.1.3 シミュレーション①.

2.1.3.1 母分散が既知.

2.1.3.1.1 標本数40.

＜帰無仮説＞.

  • 群G1（平均5, 標準偏差1）からの大きさ40の標本 \(x_1\)
  • 群G2（平均5, 標準偏差1）からの大きさ40の標本 \(x_2\)

のもとで、
　　標本平均の差 　 \(\overline{X_1}-\overline{X_2}~~\sim N\left(0,~\frac{1^2}{40}+\frac{1^2}{40}\right)\)

\(\overline{X_1}-\overline{X_2}\) の母標準偏差を尺度としたマハラノビスの座標は、
標準正規分布 \(N(0,1)\) に従うはずです。
標準正規分布の理論曲線を赤で重ね書きして確かめましょう。

XZ=c()
for (i in 1:1000){
  x1 <- rnorm(40,5,1);x2 <- rnorm(40,5,1)
  x <- mean(x1)-mean(x2)
  xz <- x/sqrt(1^2/40+1^2/40)
  XZ <- c(XZ,xz)
}
hist(XZ,prob=T,breaks = "Scott",main="",xlab="",col="lightblue")
curve(dnorm(x, 0, 1),add=T,col="red")

マハラノビスの座標は、標準正規分布に従っているようです。

＜対立仮説＞.

  • 群G1（平均5.5, 標準偏差1）からの大きさ40の標本x1.
  • 群G2（平均5, 標準偏差1）からの大きさ40の標本x2.

のもとで、
　　標本平均の差 　 \(\overline{X_1}-\overline{X_2}~~\sim N\left(0.5,~\frac{1^2}{40}+\frac{1^2}{40}\right)\)

このように、２群に平均差があるとしたら、\(\overline{X_1}-\overline{X_2}~~\) の 0 からのマハラノビスの距離はどうなるでしょうか。

帰無仮説と対立仮説の中心間の距離は 0.5 ですが、

　　　マハラノビス距離は \(\frac{0.5}{\sqrt{\frac{1}{40}+\frac{1}{40}}}\fallingdotseq 2.24~\)

Z検定では、有意差が検知される可能性が高いようです。
ヒストグラムで視覚化してみましょう。

  • 帰無仮説のもとでの理論曲線を青で示します。
  • 赤の縦線は帰無仮説のもとでの棄却境界であり、
    ヒストグラムで赤の縦線の右側の値が得られたとき帰無仮説は棄却されます。 　　
YZ=c()
for (i in 1:1000){
  y1 <- rnorm(40,5.5,1);y2 <- rnorm(40,5,1)
  y <- mean(y1)-mean(y2)
  yz <- y/sqrt(1^2/40+1^2/40)
  YZ <- c(YZ,yz)
}
hist(YZ,prob=T,breaks = "Scott",main="",xlab="")
curve(dnorm(x, 0, 1),col="blue",add=T)
abline(v = qnorm(0.95,0,1),col="red")

ヒストグラムの赤い縦線より右の範囲では帰無仮説が棄却されることになります。

2.1.3.1.2 標本数10.

＜帰無仮説＞.

  • 群G1（平均5, 標準偏差1）からの大きさ10の標本 \(x_1\)
  • 群G2（平均5, 標準偏差1）からの大きさ10の標本 \(x_2\)

のもとで、
　　標本平均の差 　 \(\overline{X_1}-\overline{X_2}~~\sim N\left(0,~\frac{1^2}{10}+\frac{1^2}{10}\right)\)

\(\overline{X_1}-\overline{X_2}\) の母標準偏差を尺度としたマハラノビスの座標は、
標準正規分布 \(N(0,1)\) に従うはずです。
標準正規分布の理論曲線を赤で重ね書きして確かめましょう。

XZ=c()
for (i in 1:1000){
  x1 <- rnorm(10,5,1);x2 <- rnorm(10,5,1)
  x <- mean(x1)-mean(x2)
  xz <- x/sqrt(1^2/10+1^2/10)
  XZ <- c(XZ,xz)
}
hist(XZ,prob=T,breaks = "Scott",main="",xlab="",col="lightblue")
curve(dnorm(x, 0, 1),add=T,col="red")

マハラノビスの座標は、標準正規分布に従っているようです。

＜対立仮説＞.

  • 群G1（平均5.5, 標準偏差1）からの大きさ10の標本 \(x_1\)
  • 群G2（平均5, 標準偏差1）からの大きさ10の標本 \(x_2\)

のもとで、
　　標本平均の差 　 \(\overline{X_1}-\overline{X_2}~~\sim N\left(0.5,~\frac{1^2}{10}+\frac{1^2}{10}\right)\)

このように、２群に平均差があるとしたら、\(\overline{X_1}-\overline{X_2}~~\) の 0 からのマハラノビスの距離はどうなるでしょうか。

帰無仮説と対立仮説の中心間の距離は 0.5 ですが、

　　　マハラノビス距離は \(\frac{0.5}{\sqrt{\frac{1}{10}+\frac{1}{10}}}\fallingdotseq 1.12~\)

標本数が小さいと、尺度である標準偏差が大きくなり、Z検定では棄却されない可能性が高いようです。
ヒストグラムで視覚化してみましょう。

  • 帰無仮説のもとでの理論曲線を青で示します。
  • 赤の縦線は帰無仮説のもとでの棄却境界であり、
    ヒストグラムで赤の縦線の右側の値が得られたとき帰無仮説は棄却されます。 　　
YZ=c()
for (i in 1:1000){
  y1 <- rnorm(10,5.5,1);y2 <- rnorm(10,5,1)
  y <- mean(y1)-mean(y2)
  yz <- y/sqrt(1^2/10+1^2/10)
  YZ <- c(YZ,yz)
}
hist(YZ,prob=T,breaks = "Scott",main="",xlab="")
curve(dnorm(x, 0, 1),col="blue",add=T)
abline(v = qnorm(0.95,0,1),col="red")

ヒストグラムの赤い縦線より右の範囲は少なく、有意差が検出される可能性は低いと言えます。

2.1.3.2 母分散が未知.

2.1.3.2.1 標本数20.

＜帰無仮説＞.

  • 群G1（平均5, 標準偏差1）からの大きさ20の標本 \(x_1\)
  • 群G2（平均5, 標準偏差1）からの大きさ20の標本 \(x_2\)

母標準偏差は未知であるとし、母標準偏差の推定値 s を用いて標準化します。
t分布の確率密度関数を青で、標準正規分布を赤で描き比較します。

XZ=c()
n1 <- 20
n2 <- 20
for (i in 1:1000){
  x1 <- rnorm(n1,5,1);x2 <- rnorm(n2,5,1)
  x <- mean(x1)-mean(x2)
  s <- (var(x1)*(n1-1)+var(x2)*(n2-1))/(n1+n1-2)
  xz <- x/sqrt(s/n1+s/n2)
  XZ <- c(XZ,xz)
}
hist(XZ,prob=T,breaks = "Scott",main="",xlab="",col="lightblue")
curve(dnorm(x, 0, 1),add=T,col="red")
abline(v = qnorm(0.95,0,1),col="red")
curve(dt(x,(n1+n2-2)),add=T,col="blue")
abline(v = qt(0.95,n1+n2-2),col="blue")

t検定の棄却境界を青の縦線で、z検定の棄却境界を赤の縦線で表示していますが、ほとんど差がないようです。（標本数40にすると完全に重なって見えます）.

＜対立仮説＞.

  • 群G1（平均5.5, 標準偏差1）からの大きさ20の標本 \(x_1\)
  • 群G2（平均5, 標準偏差1）からの大きさ20の標本 \(x_2\)

20個の標本で、平均差 0.5 を検知できるでしょうか。
対立仮説のもとでの出力をヒストグラムで、帰無仮説のもとでの理論曲線を青で表示します。

YZ=c()
n1 <- 20
n2 <- 20
for (i in 1:1000){
  y1 <- rnorm(n1,5.5,1);y2 <- rnorm(n2,5,1)
  y <- mean(y1)-mean(y2)
  s <- (var(y1)*(n1-1)+var(y2)*(n2-1))/(n1+n1-2)
  yz <- y/sqrt(s/n1+s/n2)
  YZ <- c(YZ,yz)
}
hist(YZ,prob=T,breaks = "Scott",main="",xlab="")
curve(dt(x,(n1+n2-2)),add=T,col="blue")
abline(v = qt(0.95,n1+n2-2),col="red")

検出力（ヒストグラムで赤縦線の右の割合）は 50% ほどで低いですが、これは標本数が20 で小さいため有意差が検出しにくいからです。
標本数を40とすれば、Z検定とほぼ同じ検出力になります。

2.1.3.2.2 標本数5.

＜帰無仮説＞.

  • 群G1（平均5, 標準偏差1）からの大きさ5の標本 \(x_1\)
  • 群G2（平均5, 標準偏差1）からの大きさ5の標本 \(x_2\)
XZ=c()
n1 <- 5
n2 <- 5
for (i in 1:1000){
  x1 <- rnorm(n1,5,1);x2 <- rnorm(n2,5,1)
  x <- mean(x1)-mean(x2)
  s <- (var(x1)*(n1-1)+var(x2)*(n2-1))/(n1+n1-2)
  xz <- x/sqrt(s/n1+s/n2)
  XZ <- c(XZ,xz)
}
hist(XZ,prob=T,breaks = "Scott",main="",xlab="",col="lightblue")
curve(dnorm(x, 0, 1),add=T,col="red")
abline(v = qnorm(0.95,0,1),col="red")
curve(dt(x,(n1+n2-2)),add=T,col="blue")
abline(v = qt(0.95,n1+n2-2),col="blue")

標本数が小さいと、t分布の裾は標準正規分布よりも広がり、棄却境界もZ検定から離れてきます。
この状態でZ検定を行えば検定は甘くなります。推定や検定を精密にするには、t検定がより適切と考えられます。

＜対立仮説＞.

  • 群G1（平均5.5, 標準偏差1）からの大きさ5の標本 \(x_1\)
  • 群G2（平均5, 標準偏差1）からの大きさ5の標本 \(x_2\)

小標本で、平均差 0.5 が検出できるでしょうか？

YZ=c()
n1 <- 10
n2 <- 10
for (i in 1:1000){
  y1 <- rnorm(n1,5.5,1);y2 <- rnorm(n2,5,1)
  y <- mean(y1)-mean(y2)
  s <- (var(y1)*(n1-1)+var(y2)*(n2-1))/(n1+n2-2)
  yz <- y/sqrt(s/n1+s/n2)
  YZ <- c(YZ,yz)
}
hist(YZ,prob=T,breaks = "Scott",main="",xlab="")
curve(dt(x,(n1+n2-2)),col="blue",add=T)
abline(v = qt(0.95,n1+n2-2),col="red")

t分布への当てはまりは良いですが、小標本のため有意差の検出が難しくなっています。

2.2 ウェルチのt検定.

例えば、「母平均5、母分散100の母集団から大きさ10の標本を抽出し不偏分散を算出する」という操作を1000回繰り返し、不偏分散の変動を見ると、

XV=c()
n1 <- 10
for (i in 1:1000){
  x1 <- rnorm(n1,5,10)
  XV <- c(XV,var(x1))
}
mv <- c(mean(XV),sd(XV));mv
## [1] 99.27459 49.58168

不偏分散の平均は100に近いですが、変動幅である標準偏差は 50 に近く、かなり大きな変動となっています。

2.2.1 ウェルチの標準化.

2.2.1.1 不等分散・母分散既知

＜\(n_1=5,\sigma_1=1、 n_2=15,\sigma_2=10\)＞.

XW=c()
n1 <- 5;s1 <- 1
n2 <- 15;s2 <- 10
for (i in 1:1000){
  x1 <- rnorm(n1,5,s1);x2 <- rnorm(n2,5,s2)
  x <- mean(x1)-mean(x2)
  xw <- x/sqrt(s1^2/n1+s2^2/n2)
  XW <- c(XW,xw)
}
hist(XW,prob=T,breaks = "Scott",main="",xlab="",col=0,xlim=c(-3,3))
curve(dnorm(x, 0, 1),add=T,col="red")
abline(v = qnorm(0.95,0,1),col="red")
curve(dt(x,(n1+n2-2)),add=T,col="blue")
abline(v = qt(0.95,n1+n2-2),col="blue")

母分散が既知であれば、

  • 分散が大きく異なっていても、
  • 標本数が大きく異なっていても、

標本差の標準化は標準正規分布に従うと考えて良さそうです。

2.2.1.2 不等分散・母分散未知

2.2.1.3 ウェルチの自由度.

赤の縦線はプール分散の自由度、青の縦線はウェルチの自由度の平均です。

＜2群の分散が１と４＞.

FF=c()
n1 <- 4
n2 <- 7
for (i in 1:1000){
  x1 <- rnorm(n1,5,1);x2 <- rnorm(n2,5,2)
  f <- ((var(x1)/n1+var(x2)/n2))^2/((var(x1)^2/(n1^2*(n1-1))+var(x2)^2/(n2^2*(n2-1))))
  FF <- c(FF,f)
}
par(family = "HiraKakuProN-W3") 
hist(FF,prob=T,main="ウェルチの自由度",xlab="")
abline(v = n1+n2-2,col="red")
abline(v = mean(FF),col="blue")

次に、２群の分散の差が大きい場合を考えます。

＜2群の分散が１と９＞.

FF=c()
n1 <- 4
n2 <- 7
for (i in 1:1000){
  x1 <- rnorm(n1,5,1);x2 <- rnorm(n2,5,3)
  f <- ((var(x1)/n1+var(x2)/n2))^2/((var(x1)^2/(n1^2*(n1-1))+var(x2)^2/(n2^2*(n2-1))))
  FF <- c(FF,f)
}
par(family = "HiraKakuProN-W3") 
hist(FF,prob=T,main="ウェルチの自由度",xlab="")
abline(v = n1+n2-2,col="red")
abline(v = mean(FF),col="blue")

補正自由度はプール分散の自由度より小さく、２群の分散の差が大きくなければ、プール分散の自由度を補正して小さくしたものと考えることができそうです。シミュレーションでは、変動する自由度の平均を自由度としたt分布を考えます。

2.2.2 シミュレーション②.

2.2.2.1 不等分散かつ等平均.

2.2.2.1.1 \(n_1=4,\sigma_1=1、 n_2=7,\sigma_2=5\).

＜帰無仮説＞.

  • 群G1（平均5, 標準偏差1）からの大きさ4の標本 \(x_1\)
  • 群G2（平均5, 標準偏差5）からの大きさ7の標本 \(x_2\)

ウェルチの標準化の自由度は f とし、その平均を自由度とした t分布曲線を表示します。

XW=c()
XP=c()
FF=c()
n1 <- 4
n2 <- 7
for (i in 1:1000){
  x1 <- rnorm(n1,5,1);x2 <- rnorm(n2,5,5)
  x <- mean(x1)-mean(x2)
  sp <- (var(x1)*(n1-1)+var(x2)*(n2-1))/(n1+n1-2)#プール分散
  xp <- x/sqrt(sp/n1+sp/n2)
  XP <- c(XP,xp)
  xw <- x/sqrt(var(x1)/n1+var(x2)/n2)
  XW <- c(XW,xw)
  f <- ((var(x1)/n1+var(x2)/n2))^2/((var(x1)^2/(n1^2*(n1-1))+var(x2)^2/(n2^2*(n2-1))))
  FF <- c(FF,f)
}
par(family = "HiraKakuProN-W3") 
par(mfrow=c(1,2))
hist(XP,prob=T,breaks ="Scott",col=0,xlim=c(-4,4),main="プール分散による標準化")
curve(dt(x,(n1+n2-2)),add=T,col="red")
hist(XW,prob=T,breaks ="Scott",col=0,xlim=c(-4,4),main="ウェルチの標準化")
curve(dt(x,mean(FF)),add=T,col="blue")

#補正自由度、プール自由度
free <- c(mean(FF),n1+n2-2);free
## [1] 6.947691 9.000000

プール分散による標準化は、理論分布から乖離しています。

2.2.2.1.2 \(n_1=10,\sigma_1=1、 n_2=4,\sigma_2=5\).

今度は、分散の小さい群の標本数を大きくしてみます。

＜帰無仮説＞.

  • 群G1（平均5, 標準偏差1）からの大きさ10の標本 \(x_1\)
  • 群G2（平均5, 標準偏差5）からの大きさ4の標本 \(x_2\)
XW=c()
XP=c()
FF=c()
n1 <- 10
n2 <- 4
for (i in 1:1000){
  x1 <- rnorm(n1,5,1);x2 <- rnorm(n2,5,5)
  x <- mean(x1)-mean(x2)
  sp <- (var(x1)*(n1-1)+var(x2)*(n2-1))/(n1+n1-2)#プール分散
  xp <- x/sqrt(sp/n1+sp/n2)
  XP <- c(XP,xp)
  xw <- x/sqrt(var(x1)/n1+var(x2)/n2)
  XW <- c(XW,xw)
  f <- ((var(x1)/n1+var(x2)/n2))^2/((var(x1)^2/(n1^2*(n1-1))+var(x2)^2/(n2^2*(n2-1))))
  FF <- c(FF,f)
}
par(family = "HiraKakuProN-W3") 
par(mfrow=c(1,2))
hist(XP,prob=T,breaks = "Scott",xlab="",col=0,xlim=c(-8,8),main="プール分散による標準化")
curve(dt(x,(n1+n2-2)),add=T,col="red")
hist(XW,prob=T,breaks = "Scott",xlab="",col=0,xlim=c(-8,8),main="ウェルチの標準化")
curve(dt(x,mean(FF)),add=T,col="blue")

#補正自由度、プール自由度
free <- c(mean(FF),n1+n2-2);free
## [1]  3.275582 12.000000

分散の小さい群の標本数が大きいと、標準化は大きな値になるようです。

2.2.3 検定の有効性.

＜帰無仮説＞.

  • 群G1（平均5, 標準偏差1）からの大きさ15の標本 \(x_1\)
  • 群G2（平均5, 標準偏差3）からの大きさ5の標本 \(x_2\)
XW=c()
FF=c()
n1 <- 15
n2 <- 5
for (i in 1:1000){
  x1 <- rnorm(n1,5,1);x2 <- rnorm(n2,5,3)
  x <- mean(x1)-mean(x2)
  xw <- x/sqrt(var(x1)/n1+var(x2)/n2)
  XW <- c(XW,xw)
  f <- ((var(x1)/n1+var(x2)/n2))^2/((var(x1)^2/(n1^2*(n1-1))+var(x2)^2/(n2^2*(n2-1))))
  FF <- c(FF,f)
}
hist(XW,prob=T,breaks = "Scott",main="",xlab="",col=0,xlim=c(-3,3))
curve(dnorm(x, 0, 1),add=T,col="red")
abline(v = qnorm(0.95,0,1),col="red")
curve(dt(x,(n1+n2-2)),add=T,col="blue")
abline(v = qt(0.95,n1+n2-2),col="blue")
curve(dt(x,mean(FF)),add=T,col="green")
abline(v = qt(0.95,f),col="green")

free <- c(mean(FF),n1+n2-2);free
## [1]  4.586026 18.000000

ヒストグラムへの当てはまりは、ウェルチの自由度による曲線のほうが良く、棄却境界もかなり異なるので、ウェルチの検定の有効性が高いと考えられます。

以下、ウェルチの自由度による t分布とその右側棄却境界のみ表示します。

＜対立仮説＞.

  • 群G1（平均8, 標準偏差1）からの大きさ15の標本 \(x_1\)
  • 群G2（平均5, 標準偏差3）からの大きさ5の標本 \(x_2\)
XW=c()
FF=c()
n1 <- 15
n2 <- 5
for (i in 1:1000){
  x1 <- rnorm(n1,8,1);x2 <- rnorm(n2,5,3)
  x <- mean(x1)-mean(x2)
  xw <- x/sqrt(var(x1)/n1+var(x2)/n2)
  XW <- c(XW,xw)
  f <- ((var(x1)/n1+var(x2)/n2))^2/((var(x1)^2/(n1^2*(n1-1))+var(x2)^2/(n2^2*(n2-1))))
  FF <- c(FF,f)
}
free <- c(mean(FF),n1+n2-2);free
## [1]  4.686231 18.000000
hist(XW,prob=T,breaks = "Scott",main="",xlab="",col=0,xlim=c(-3,3))
curve(dt(x,mean(FF)),add=T,col="green")
abline(v = qt(0.95,f),col="green")

検出力は高いとは言えないようです。 　

分散が大きい群の標本を大きくした対立仮説でシミュレーションします。

＜対立仮説＞.

  • 群G1（平均8, 標準偏差1）からの大きさ5の標本 \(x_1\)
  • 群G2（平均5, 標準偏差3）からの大きさ15の標本 \(x_2\)
XW=c()
FF=c()
n1 <- 5
n2 <- 15
for (i in 1:1000){
  x1 <- rnorm(n1,8,1);x2 <- rnorm(n2,5,3)
  x <- mean(x1)-mean(x2)
  xw <- x/sqrt(var(x1)/n1+var(x2)/n2)
  XW <- c(XW,xw)
  f <- ((var(x1)/n1+var(x2)/n2))^2/((var(x1)^2/(n1^2*(n1-1))+var(x2)^2/(n2^2*(n2-1))))
  FF <- c(FF,f)
}
#自由度
f <- ((var(x1)/n1+var(x2)/n2))^2/((var(x1)^2/(n1^2*(n1-1))+var(x2)^2/(n2^2*(n2-1))))
free <- c(mean(FF),n1+n2-2);free
## [1] 16.45355 18.00000
hist(XW,prob=T,breaks = "Scott",main="",xlab="",col=0,xlim=c(-3,3))
curve(dt(x,mean(FF)),add=T,col="green")
abline(v = qt(0.95,f),col="green")

分散が大きい群の標本数が大きいほうが、検出力は高くなるようです。

2.2.4 P値の分布の比較

＜\(n_1=5,\sigma_1=1、 n_2=15,\sigma_2=3\)＞.

＜対立仮説＞.

  • 群G1（平均8, 標準偏差1）からの大きさ15の標本 \(x_1\)
  • 群G2（平均5, 標準偏差3）からの大きさ5の標本 \(x_2\)
PP=c()
PW=c()
n1 <- 5
n2 <- 15
for (i in 1:1000){
  x1 <- rnorm(n1,8,1);x2 <- rnorm(n2,5,3)
  x <- mean(x1)-mean(x2)
  sp <- (var(x1)*(n1-1)+var(x2)*(n2-1))/(n1+n1-2)#プール分散
  xp <- x/sqrt(sp/n1+sp/n2)
  pp <- 1- pt(xp,n1+n2-2)
  PP <- c(PP,pp)  
  xw <- x/sqrt(var(x1)/n1+var(x2)/n2)
  f <- ((var(x1)/n1+var(x2)/n2))^2/((var(x1)^2/(n1^2*(n1-1))+var(x2)^2/(n2^2*(n2-1))))
  pw <- 1- pt(xw,f)
  PW <- c(PW,pw)
}
par(mfrow=c(1,2))
par(family = "HiraKakuProN-W3") 
hist(PP,prob=T,breaks = "Scott",main="等分散仮定でのP値",xlab="",ylim=c(0,12))
hist(PW,prob=T,breaks = "Scott",main="ウェルチt検定のP値",xlab="",ylim=c(0,12))

mean <- c(mean(PP),mean(PW));mean
## [1] 0.09896764 0.01355102

ウェルチt検定の方が、P値が小さい傾向があり、検出力が高いと言えそうです。

次に、分散が大きい群の標本数を大きくしてみます。

＜対立仮説＞.

  • 群G1（平均7, 標準偏差1）からの大きさ5の標本 \(x_1\)
  • 群G2（平均5, 標準偏差4）からの大きさ15の標本 \(x_2\)
PP=c()
PW=c()
n1 <- 5
n2 <- 15
for (i in 1:1000){
  x1 <- rnorm(n1,7,1);x2 <- rnorm(n2,5,4)
  x <- mean(x1)-mean(x2)
  sp <- (var(x1)*(n1-1)+var(x2)*(n2-1))/(n1+n1-2)#プール分散
  xp <- x/sqrt(sp/n1+sp/n2)
  pp <- 1- pt(xp,n1+n2-2)
  PP <- c(PP,pp)  
  xw <- x/sqrt(var(x1)/n1+var(x2)/n2)
  f <- ((var(x1)/n1+var(x2)/n2))^2/((var(x1)^2/(n1^2*(n1-1))+var(x2)^2/(n2^2*(n2-1))))
  pw <- 1- pt(xw,f)
  PW <- c(PW,pw)
}
par(mfrow=c(1,2))
par(family = "HiraKakuProN-W3") 
hist(PP,prob=T,breaks = "Scott",main="等分散仮定でのP値",xlab="",ylim=c(0,12))
hist(PW,prob=T,breaks = "Scott",main="ウェルチt検定のP値",xlab="",ylim=c(0,12))

mean <- c(mean(PP),mean(PW));mean
## [1] 0.2417339 0.1007390

やはり、ウェルチのt検定の方が、有意差を検出しやすいようです。

次に、分散の大きい群の標本数を小さくしてみます。

＜\(n_1=15,\sigma_1=1、 n_2=5,\sigma_2=2\)＞.

PP=c()
PW=c()
n1 <- 10
n2 <- 15
for (i in 1:1000){
  x1 <- rnorm(n1,7,1);x2 <- rnorm(n2,5,2)
  x <- mean(x1)-mean(x2)
  sp <- (var(x1)*(n1-1)+var(x2)*(n2-1))/(n1+n1-2)#プール分散
  xp <- x/sqrt(sp/n1+sp/n2)
  pp <- 1- pt(xp,n1+n2-2)
  PP <- c(PP,pp)  
  xw <- x/sqrt(var(x1)/n1+var(x2)/n2)
  f <- ((var(x1)/n1+var(x2)/n2))^2/((var(x1)^2/(n1^2*(n1-1))+var(x2)^2/(n2^2*(n2-1))))
  pw <- 1- pt(xw,f)
  PW <- c(PW,pw)
}
par(mfrow=c(1,2))
par(family = "HiraKakuProN-W3") 
hist(PP,prob=T,breaks = "Scott",main="等分散仮定でのP値",xlab="",ylim=c(0,12))
hist(PW,prob=T,breaks = "Scott",main="ウェルチt検定のP値",xlab="",ylim=c(0,12))

mean <- c(mean(PP),mean(PW));mean
## [1] 0.02567482 0.01272526

ウェルチのt検定の方が安定しているようです。

2.3 正規性の問題.

2.3.1 片方が非正規分布.

カイ二乗分布と正規分布の場合、平均差の分布はどのようになるでしょうか。

＜帰無仮説＞.

  • 群G1（平均8, 標準偏差4）のカイ二乗分布からの大きさ5の標本 \(x_1\)
  • 群G2（平均8, 標準偏差4）の正規分布からの大きさ10の標本 \(x_2\)
XW=c()
FF=c()
n1 <- 5
n2 <- 10
for (i in 1:1000){
  x1 <- rchisq(n1,8);x2 <- rnorm(n2,8,4)
  x <- mean(x1)-mean(x2)
  xw <- x/sqrt(var(x1)/n1+var(x2)/n2)
  XW <- c(XW,xw)
  f <- ((var(x1)/n1+var(x2)/n2))^2/((var(x1)^2/(n1^2*(n1-1))+var(x2)^2/(n2^2*(n2-1))))
  FF <- c(FF,f)
}
hist(XW,prob=T,breaks = "Scott",main="",xlab="",col="lightblue")
curve(dt(x,mean(FF)),add=T,col="blue")#自由度

ヴェルチの標準化は、t分布にほぼ一致しているようです。
検定の有効性はどうでしょうか。

＜対立仮説＞.

  • 群G1（平均8, 標準偏差4）のカイ二乗分布からの大きさ5の標本 \(x_1\)
  • 群G2（平均6, 標準偏差4）からの大きさ10の標本 \(x_2\)
XW=c()
FF=c()
n1 <- 5
n2 <- 10
for (i in 1:1000){
  x1 <- rchisq(n1,8);x2 <- rnorm(n2,6,4)
  x <- mean(x1)-mean(x2)
  xw <- x/sqrt(var(x1)/n1+var(x2)/n2)
  XW <- c(XW,xw)
  f <- ((var(x1)/n1+var(x2)/n2))^2/((var(x1)^2/(n1^2*(n1-1))+var(x2)^2/(n2^2*(n2-1))))
  FF <- c(FF,f)
}
hist(XW,prob=T,breaks = "Scott",main="",xlab="",col="lightblue")
curve(dt(x,mean(FF)),add=T,col="blue")#自由度
abline(v = qt(0.95,mean(FF)),col="red")

検出力は高くありませんが、これは小標本のためとも考えられます。

P値の分布と平均を調べます。.

PW=c()
n1 <- 5
n2 <- 10
for (i in 1:1000){
  x1 <- rchisq(n1,8);x2 <- rnorm(n2,6,4)
  x <- mean(x1)-mean(x2)
  xw <- x/sqrt(var(x1)/n1+var(x2)/n2)
  f <- ((var(x1)/n1+var(x2)/n2))^2/((var(x1)^2/(n1^2*(n1-1))+var(x2)^2/(n2^2*(n2-1))))
  pw <- 1- pt(xw,f)
  PW <- c(PW,pw)
}
par(family = "HiraKakuProN-W3") 
hist(PW,prob=T,breaks = "Scott",main="P値の分布",xlab="")

mean(PW)
## [1] 0.2912395

P値は大きい傾向があり、平均でも0.05を超えています。
「平均差2」は検出されにくいようです。

2.3.2 外れ値の影響.

一方の群からの標本に外れ値を付け加えます。
平均は、ほぼ一致するよう操作します。

＜帰無仮説＞.

  • 群G1（平均5, 標準偏差1）の正規分布からの大きさ4の標本、及び20 \(x_1\)
  • 群G2（平均8, 標準偏差2）の正規分布からの大きさ10の標本 \(x_2\)

各群の平均は 8 となり、２群の平均差は 0 の近くに分布します。
ヴェルチの標準化の分布と併置し、ヴェルチのt分布への当てはまりを見ます。

X=c();XW=c()
n1 <- 5
n2 <- 10
for (i in 1:1000){
  x1 <- c(rnorm(n1-1,5,1),20);x2 <- rnorm(n2,8,2)
  x <- mean(x1)-mean(x2)
  X <- c(X,x)
  xw <- x/sqrt(var(x1)/n1+var(x2)/n2)
  XW <- c(XW,xw)
  f <- ((var(x1)/n1+var(x2)/n2))^2/((var(x1)^2/(n1^2*(n1-1))+var(x2)^2/(n2^2*(n2-1))))
  FF <- c(FF,f)
}
par(family = "HiraKakuProN-W3") 
par(mfrow=c(1,2))
hist(X,prob=T,breaks = "Scott",main="平均差の分布",xlab="")
hist(XW,prob=T,breaks = "Scott",main="ヴェルチの標準化",xlab="")
curve(dt(x,mean(FF)),add=T,col="blue")#自由度
abline(v = qt(0.95,mean(FF)),col="red")

外れ値のため標準偏差が過大となり、ヴェルチ標準化は過小となるので、0の近くに分布が集中しています。
また、fitする理論曲線がないので検定は困難になります。

＜対立仮説＞.

  • 群G1（平均5, 標準偏差1）の正規分布からの大きさ4の標本、及び20 \(x_1\)
  • 群G2（平均6, 標準偏差2）の正規分布からの大きさ10の標本 \(x_2\)

群G1の平均は 8 、群G2の平均は 6 となりますが、ヴェルチの標準化の分布を青のヒストグラムで重ねます。

Y=c();YW=c()
n1 <- 5
n2 <- 10
for (i in 1:1000){
  y1 <- c(rnorm(n1-1,5,1),20);y2 <- rnorm(n2,6,3)
  y <- mean(y1)-mean(y2)
  Y <- c(Y,y)
  yw <- y/sqrt(var(y1)/n1+var(y2)/n2)
  YW <- c(YW,yw)
  f <- ((var(y1)/n1+var(y2)/n2))^2/((var(y1)^2/(n1^2*(n1-1))+var(y2)^2/(n2^2*(n2-1))))
  FF <- c(FF,f)
}
par(family = "HiraKakuProN-W3") 
hist(XW,prob=T,breaks = "Scott",main="",xlab="",col=0)
hist(YW,prob=T,breaks = "Scott",main="",xlab="",col="lightblue",add=T)

平均差があることは、ヒストグラムのズレから分かりますが、理論曲線が不明なのでP値を求めることができません。

このように外れ値がある場合は、平均差の統計的検定は不可能になります。

3 ノンパラメトリック検定.

3.1 Wilcoxon順位和検定.

3.2 順位分布

3.2.1 順位と一様分布.

大きさ10の乱数を発生させ、rank( )関数で順位をつけます。

  x <- rnorm(10,5,1)
  round(x,2)
##  [1] 6.56 5.48 3.91 4.51 5.72 5.13 3.50 4.53 5.99 5.28
  r <- rank(x, ties.method = "min");r
##  [1] 10  7  2  3  8  5  1  4  9  6

1番目のデータの順位 r1は確率変数と考えられます。
rank( )関数が昇順の順位を返すことが分かります。
以下、順位は断らない限り昇順として話を進めます。

r1 を1000回出力したデータを S1 として分布を調べます。

S1=c()
n <- 10
for (i in 1:1000){
  x <- rnorm(n,5,1)
  r <- rank(x, ties.method = "min")
  r1 <- r[1]
  S1<- c(S1,r1)
}
hist(S1,prob=T,breaks=seq(0.5, 10.5, by = 1),main="",xlab="")

一様分布になっており、どの順位も同様に確からしいということになります。

3.2.2 平均と分散

n 個から１個選んだものの順位をRとすると、Rは一様分布に従う確率変数で、

　　\(P(R=i)=\frac{1}{n}\)　　　\(\displaystyle E(R)=\sum_{i=1}^n i\cdot \frac{1}{n}=\frac{1}{n}\cdot\frac{n(n+1)}{2}=\frac{n+1}{2}\)

　　\(\displaystyle E(R^2)=\sum_{i=1}^n i^2\cdot \frac{1}{n}=\frac{1}{n}\cdot\frac{n(n+1)(2n+1)}{6}=\frac{(n+1)(2n+1)}{6}\)

　\(\therefore V(R)=\frac{(n+1)(2n+1)}{6}-\left(\frac{n+1}{2}\right)^2=\frac{(n+1)(n-1)}{12}\)

シミュレーションデータで確かめてみましょう。

m <- c(mean(S1),(n+1)/2)
m
## [1] 5.509 5.500
v <- c(var(S1),(n+1)*(n-1)/12)
v
## [1] 8.444363 8.250000

確かに一致しているようです。

3.2.3 順位和統計量

3.3 小標本の場合

\(n_1=4,n_2=6\) の場合について、２群の順位分布を比較します。

3.3.1 同一分布の２群.

  • 群G1（平均3, 標準偏差1）からの大きさ4の標本 \(x_1\)
  • 群G2（平均3, 標準偏差1）からの大きさ6の標本 \(x_2\)

２つの標本を合わせて順位をつけ、順位データを r1,r2 とします。
同じ母集団からの標本なので、２群の順位分布に差はないはずです。

n1 <- 4
n2 <- 6
  x1 <- rnorm(n1,3,1);x2 <- rnorm(n2,3,1)
  x <- c(x1,x2)
  r <- rank(x, ties.method = "average")
  r1 <- head(r,n1)
  r2 <- tail(r,n2)

ties.method = “average” は、同じ順位のとき平均順位にするというオプションです。
例えば、５位と６位が同じ点数なら 5.5位とします。

print(r1)
## [1]  4  9 10  3
print(r2)
## [1] 5 1 2 6 8 7

２群の順位分布を箱ひげ図で表示し比較します。

library(beeswarm)
list=list(A=r1, B=r2)
boxplot(list)
beeswarm(list,  pch = 16, add = TRUE)

同じ順位分布のはずですが、小標本だと変動が見られます。

繰り返しプログラムで、順位和 \(W_1,~W_2\) を1000回出力します。
２群の順位和の分布を確かめましょう。

S1=c()
S2=c()
n1 <- 4
n2 <- 6
for (i in 1:1000){
  x1 <- rnorm(n1,3,1);x2 <- rnorm(n2,3,1)
  x <- c(x1,x2)
  r <- rank(x, ties.method = "average")
  r1 <- head(r,n1)
  r2 <- tail(r,n2)
  S1<- c(S1,sum(r1))
  S2<- c(S2,sum(r2))
}
list=list(A=S1, B=S2)
par(family = "HiraKakuProN-W3") 
boxplot(list,main="順位和")

２群の標本数が異なるので、順位和の期待値は異なりますが、分散は同じになっています。

次は、順位平均の分布を比較します。

list=list(A=S1/n1, B=S2/n2)
par(family = "HiraKakuProN-W3") 
boxplot(list,main="順位平均")

順位平均の中央値は 5.5 で一致しています。
また、範囲や四分位範囲は相似形となっていることが観察されます。

　　　\(V(\bar{r_1})=\frac{V(r_1)}{4^2}=\frac{22}{16}\)　　　\(V(\bar{r_2})=\frac{V(r_1)}{6^2}=\frac{22}{36}\)

分散は \(\frac{36}{16}\)倍、標準偏差は \(\frac{6}{4}=1.5\)倍 となり、範囲も 1.5 倍です。

max(S1/n1)-min(S1/n1)
## [1] 6
max(S2/n2)-min(S2/n2)
## [1] 4

3.3.2 ２群の順位和分布.

２群の順位和のヒストグラムを並置して比較します。

par(mfrow=c(1,2))
par(family = "HiraKakuProN-W3") 
hist(S1,breaks = "Scott",main="W1",xlab="",col=0)
hist(S2,breaks = "Scott",main="W2",xlab="",col=0)

ヒストグラムはシンメトリー（対称）であることが分かります。

hist(S1,main="",xlab="",col=0,xlim=c(8,50),breaks=seq(8.5, 47.5, by = 2))
hist(S2,main="",xlab="",col="lightblue",breaks=seq(8.5, 47.5, by = 2), add=T)
abline(v = 55/2,col="red")

赤の縦線は \(x=27.5\)であり、２つのヒストグラムは赤の縦線について左右対称になっています。

3.3.3 理論分布.

小標本なので、\(W_1\)の理論分布は組み合わせによる計算が可能です。

Rの関数 combn(x=c(1:10),m=4) は、1から10までの順位から4個の順位を取り出す組合せのすべてを出力します。
はじめの6個を表示します。

combn(x=c(1:10),m=4)[,1:6]
##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,]    1    1    1    1    1    1
## [2,]    2    2    2    2    2    2
## [3,]    3    3    3    3    3    3
## [4,]    4    5    6    7    8    9

列の和が順位和ということになります。

列数は \({}_{10} C_4\) で、Rの関数choose(10,4) で計算されます。

N <- choose(10,4)
N
## [1] 210

次のコードで、210個の順位和が出力されます。
head(S) で 210個のはじめの6個を表示します。

data <- combn(x=c(1:10),m=4)
N <- choose(10,4)
S=c()
for (i in 1:N){
  s <- sum(data[,i])
  S <- c(S,s)
}
head(S)
## [1] 10 11 12 13 14 15

210個の順位和を table( ) 関数で整理し、相対度数を求めれば確率分布が分かります。
表は、はじめの６列を表示し、全体はグラフで表示します。

head(table(S))
## S
## 10 11 12 13 14 15 
##  1  1  2  3  5  6
#sum(table(S)/N)
head(round(table(S)/N,3))
## S
##    10    11    12    13    14    15 
## 0.005 0.005 0.010 0.014 0.024 0.029
plot(table(S)/N,type="l",col="lightblue")

＜累積分布＞

分布関数は、cumsum( ) 関数で相対度数の累積和を求めます。

head(round(cumsum(table(S)/N),3))
##    10    11    12    13    14    15 
## 0.005 0.010 0.019 0.033 0.057 0.086

barplot( ) 関数を用いて、棒グラフで図示できます。
赤い横線は両側 5% 検定の棄却境界です。

par(family = "HiraKakuProN-W3") 
barplot(cumsum(table(S)/N),main="P(W1≦x)")
abline(h = c(0.025,0.975),col="red")

barplot(1-cumsum(table(S)/N),main="1-P(W1≦x)")
abline(h = c(0.025,0.975),col="red")

3.3.4 帰無仮説.

＜帰無仮説＞.

  • 群G1（平均3, 標準偏差1）の正規分布からの大きさ4の標本 \(x_1\)
  • 群G2（平均3, 標準偏差1）の正規分布からの大きさ6の標本 \(x_2\)
S1=c()
S2=c()
n1 <- 4
n2 <- 6
for (i in 1:1000){
  x1 <- rnorm(n1,3,1);x2 <- rnorm(n2,3,1)
  x <- c(x1,x2)
  r <- rank(x, ties.method = "average")
  r1 <- head(r,n1)
  r2 <- tail(r,n2)
  S1<- c(S1,sum(r1))
  S2<- c(S2,sum(r2))
}
par(family = "HiraKakuProN-W3") 
hist(S1,prob=T,main="",ylab="",col=0,xlab="")
par(new=T)
plot(table(S)/N,type="l",col="blue", yaxt="n",ylab="", xaxt="n",xlab="")
curve(dnorm(x,22,sqrt(22)),add=T,col="green")#自由度

シミュレーションによる確率分布と理論分布（青）はほぼ一致しています。

3.3.5 対立仮説.

＜平均差がある場合＞.

２群の平均に差があるとき、２群の順位分布にも差が出てきます。

＜対立仮説＞.

  • 群G1（平均3, 標準偏差1）の正規分布からの大きさ4の標本 \(x_1\)
  • 群G2（平均4, 標準偏差1）の正規分布からの大きさ6の標本 \(x_2\)
n1 <- 4
n2 <- 6
  x1 <- rnorm(n1,3,1);x2 <- rnorm(n2,4,1)
  x <- c(x1,x2)
  r <- rank(x, ties.method = "average")
  r1 <- head(r,n1)
  r2 <- tail(r,n2)
  list=list(A=r1, B=r2)
boxplot(list)
beeswarm(list,  pch = 16, add = TRUE)

「２群の平均に差がある」（平均差1）という対立仮説のもとで \(W_1\) を1000回出力し、有意差を検出できるかどうかシミュレーションします。

＜対立仮説＞.

  • 群G1（平均3, 標準偏差1）の正規分布からの大きさ4の標本 \(x_1\)
  • 群G2（平均2, 標準偏差1）の正規分布からの大きさ6の標本 \(x_2\)
S1=c()
S2=c()
n1 <- 4
n2 <- 6
for (i in 1:1000){
  x1 <- rnorm(n1,3,1);x2 <- rnorm(n2,2,1)
  x <- c(x1,x2)
  r <- rank(x, ties.method = "average")
  r1 <- head(r,n1)
  r2 <- tail(r,n2)
  S1<- c(S1,sum(r1))
  S2<- c(S2,sum(r2))
}
par(family = "HiraKakuProN-W3") 
hist(S1,prob=T,main="",ylab="",col=0,xlab="",xaxt="n")
par(new=T)
plot(table(S)/N,type="l",col="blue", yaxt="n",ylab="",cex.axis=0.5,xlab="")
abline(v = 31,col="blue")
curve(dnorm(x,22,sqrt(22)),add=T,col="green")
abline(v = qnorm(0.95,22,sqrt(22)),col="red")

赤の縦線は正規近似による右側棄却境界、青の縦線は理論分布による右側棄却境界です。赤と青の棄却境界が離れていますが、整数値補正を行わなかったためとも考えられます。

3.3.6 外れ値がある場合.

\(n_1=7,n_2=7\) の場合について、順位和の理論分布を求めます。

組合せ計算による理論分布を青で、正規近似による分布を赤で表示します。

data <- combn(x=c(1:14),m=7)
N <- choose(14,7)
S=c()
for (i in 1:N){
  s <- sum(data[,i])
  S <- c(S,s)
}
par(family = "HiraKakuProN-W3") 
plot(table(S)/N,type="l",col="lightblue",main="順位和の理論分布（青）と正規近似（赤）",xlab="")
curve(dnorm(x,52.5,sqrt(61.25)),add=T,col="red")

標本数が増えると曲線は滑らかになり、正規分布との乖離も小さくなっています。

3.3.6.1 帰無仮説.

＜帰無仮説＞.

  • 平均5, 標準偏差2 の正規乱数6個に20を付加した \(x_1\)
  • 平均5, 標準偏差2 の正規乱数7個 \(x_2\)

ともに大きさ7の標本ですが、\(x_1\) に外れ値があります。

n1 <- 7
n2 <- 7
x1 <- c(rnorm(n1-1,5,2),20);x2 <- rnorm(n2,5,2)
round(x1,2);round(x2,2)
## [1]  7.89  7.37  2.38  4.62  3.14  4.23 20.00
## [1] 5.23 4.62 7.36 5.27 4.76 4.13 4.46
par(family = "HiraKakuProN-W3") 
list=list(x1=x1, x2=x2)
boxplot(list,main="生データの分布")
beeswarm(list,  pch = 16, add = TRUE)

x <- c(x1,x2)
r <- rank(x, ties.method = "average")
r1 <- head(r,n1);r1
## [1] 13 12  1  7  2  4 14
r2 <- tail(r,n2);r2
## [1]  9  6 11 10  8  3  5
list=list(r1=r1, r2=r2)
boxplot(list,main="順位データの分布")
beeswarm(list,  pch = 16, add = TRUE)

順位分布の差はあまりないようです。

順位和（昇順）の分布をみてみましょう。
\(x_1\) に外れ値がある状況では、昇順の順位和は \(W_1\) が少し大きくなるはずです。
外れ値の影響がどれくらいあるか見てみましょう。
同じ条件で1000回順位和を出力し、 \(W_1\) のヒストグラムと理論分布を比べます。

S1=c()
S2=c()
Y=c()
n1 <- 7
n2 <- 7
for (i in 1:1000){
  x1 <- c(rnorm(n1-1,5,2),20);x2 <- rnorm(n2,5,2)
  x <- c(x1,x2)
  r <- rank(x, ties.method = "average")
  r1 <- head(r,n1)
  r2 <- tail(r,n2)
  S1<- c(S1,sum(r1))
  S2<- c(S2,sum(r2))
}
par(family = "HiraKakuProN-W3") 
hist(S1,prob=T,col=0, xaxt="n", yaxt="n",xlab="",main="順位和の分布")
par(new=T)
plot(table(S)/N,type="l",col="lightblue",ylab="",xlab="")

外れ値の影響は少なく、順位和の分布は理論分布にほぼ当てはまっています。

3.3.6.2 対立仮説.

＜外れ値を除くと平均差がある＞.

  • 平均3, 標準偏差2 の正規乱数6個に20を付加した \(x_1\)
  • 平均5, 標準偏差2 の正規乱数7個 \(x_2\)
n1 <- 7
n2 <- 7
x1 <- c(rnorm(n1-1,3,2),20);x2 <- rnorm(n2,5,2)
round(x1,2);round(x2,2)
## [1]  4.75  2.98  5.91  4.54 -1.19  3.56 20.00
## [1] 2.68 4.66 4.90 4.34 3.94 6.30 7.32
par(family = "HiraKakuProN-W3") 
list=list(x1=x1, x2=x2)
boxplot(list,main="生データの分布")
beeswarm(list,  pch = 16, add = TRUE)

x <- c(x1,x2)
r <- rank(x, ties.method = "average")
r1 <- head(r,n1);r1
## [1]  9  3 11  7  1  4 14
r2 <- tail(r,n2);r2
## [1]  2  8 10  6  5 12 13
list=list(r1=r1, r2=r2)
boxplot(list,main="順位データの分布")
beeswarm(list,  pch = 16, add = TRUE)

あまり差がないように見えるので、検出力は高くないかもしれません。
同じ条件で1000回順位和を出力し、ヒストグラムと理論曲線を重ねます。

S1=c()
S2=c()
n1 <- 7
n2 <- 7
for (i in 1:1000){
  x1 <- c(rnorm(n1-1,3,2),20);x2 <- rnorm(n2,5,2)
  x <- c(x1,x2)
  r <- rank(x, ties.method = "average")
  r1 <- head(r,n1)
  r2 <- tail(r,n2)
  S1<- c(S1,sum(r1))
  S2<- c(S2,sum(r2))
}
par(family = "HiraKakuProN-W3") 
hist(S1,prob=T,col=0, xaxt="n", yaxt="n",xlab="",main="順位和の分布")
par(new=T)
plot(table(S)/N,type="l",col="lightblue",ylab="",xlab="")
abline(v = qnorm(0.05,52.5,sqrt(61.25)),col="red")

赤の縦線は左側5%の棄却境界であり、その左側が棄却されます。
ある程度は平均差が検出されますが、検出力は高くありません。

3.3.6.3 t検定との比較

＜外れ値を含めると平均差がない場合＞.

G1に1番の生徒がいるとして、
外れ値を含めてクラス平均が同じであるときの２群の学力差を考えます。
外れ値を除くと平均差があることになるのですが、学力差は検出されるでしょうか。

  • 平均4, 標準偏差2 の正規乱数6個に11を付加した \(x_1\)
  • 平均5, 標準偏差2 の正規乱数7個 \(x_2\)
n1 <- 20
n2 <- 20
  x1 <- c(rnorm(n1-1,4,2),24);x2 <- rnorm(n2,5,2)
  mean <- c(mean(x1),mean(x2));mean
## [1] 4.432449 4.996990

標本数が小さいので、出力した平均には若干の差があります。

簡単のため Rの wilcox.exact()、t.test() を用いてP値を算出し比較します。

library(exactRankTests)
wilcox.exact(x=x1,y=x2,paired=F)$p.value
## [1] 0.06343245
t.test(x1,x2)$ p.value
## [1] 0.6431074

外れ値がある場合は、順位和統計量のほうが P値は小さいようです。
同じ条件で1000回シミュレーションを行い、P値の出力分布を調べます。

＜1000回繰り返し＞

PWは順位和統計量によるP値、PVはヴェルチのt検定によるP値です。

S1=c()
S2=c()
PV=c()
PW=c()
n1 <- 20
n2 <- 20
for (i in 1:1000){
  x1 <- c(rnorm(n1-1,4,2),24);x2 <- rnorm(n2,5,2)

  PW=c(PW,wilcox.exact(x=x1,y=x2,paired=F)$p.value)
  PV=c(PV,t.test(x1,x2)$ p.value)
}
par(mfrow=c(1,2))
hist(PW,prob=T)
hist(PV,prob=T)

外れ値がある場合は、順位和検定の方が有利なようです。　　

次に、外れ値がない場合を調べます。

  • 平均4, 標準偏差2 の正規乱数10個 \(x_1\)
  • 平均5, 標準偏差2 の正規乱数10個 \(x_2\)
PV=c()
PW=c()
n1 <- 10
n2 <- 10
for (i in 1:1000){
  x1 <- c(rnorm(n1,4,2));x2 <- rnorm(n2,5,2)
  PW=c(PW,wilcox.exact(x=x1,y=x2,paired=F)$p.value)
  PV=c(PV,t.test(x1,x2)$ p.value)
}
par(mfrow=c(1,2))
hist(PW,prob=T)
hist(PV,prob=T)

外れ値のない場合は、若干、t検定の方が有利なことが分かります。

3.4 大標本の場合.

3.4.1 帰無仮説

\(n_1=20,n_2=20\)、2群に差がない場合を考えます。.

  • 群G1（平均4, 標準偏差1）の正規分布からの大きさ20の標本 \(x_1\)
  • 群G2（平均4, 標準偏差1）の正規分布からの大きさ20の標本 \(x_2\)

２群の分布に差がないので、\(z = \frac{W_1-E(W_1)}{\sqrt{Var(W_1)}}\)　は標準正規分布に従うはずです。

S1=c()
S2=c()
n1 <- 20
n2 <- 20
for (i in 1:1000){
  x1 <- rnorm(n1,4,1);x2 <- rnorm(n2,4,1)
  x <- c(x1,x2)
  r <- rank(x, ties.method = "average")
  r1 <- head(r,n1)
  r2 <- tail(r,n2)
  S1<- c(S1,sum(r1))
  S2<- c(S2,sum(r2))
}
Z1  <- (S1-n1*(n1+n2+1)/2)/sqrt(n1*n2*(n1+n2+1)/12)
par(family = "HiraKakuProN-W3") 
hist(Z1,prob=T,main="",ylab="",col=0,xlab="")
curve(dnorm(x,0,1),add=T,col="blue")
abline(v = qnorm(0.95),col="red")

青は標準正規分布のグラフ、赤の縦線は右側95%の棄却境界です。
標準正規分布で近似できると考えてよさそうです。

次に、対立仮説を設定し、平均差を検出できるかシミュレーションします。

3.4.2 対立仮説.

  • 群G1（平均4, 標準偏差1）の正規分布からの大きさ20の標本 \(x_1\)
  • 群G2（平均3, 標準偏差1）の正規分布からの大きさ20の標本 \(x_2\)

順位和の出力値がどの程度ずれるか、1000回シミュレーションして確かめます。

S1=c()
S2=c()
n1 <- 20
n2 <- 20
for (i in 1:1000){
  x1 <- rnorm(n1,4,1);x2 <- rnorm(n2,3,1)
  x <- c(x1,x2)
  r <- rank(x, ties.method = "average")
  r1 <- head(r,n1)
  r2 <- tail(r,n2)
  S1<- c(S1,sum(r1))
  S2<- c(S2,sum(r2))
}
Z1  <- (S1-n1*(n1+n2+1)/2)/sqrt(n1*n2*(n1+n2+1)/12)
par(family = "HiraKakuProN-W3") 
hist(Z1,prob=T,main="",ylab="",col=0,xlab="")
curve(dnorm(x,0,1),add=T,col="blue")
abline(v = qnorm(0.95),col="red")

平均差が１のとき、順位和の正規化を統計量としたときの検出力はかなり大きくなっています。

実際に、検出力を計算してみましょう。
sum(Z > a) は Zの中でaより大きい要素の個数をカウントします。
赤の縦線より右にくる Z1 の要素の割合を計算します。

qnorm(0.95)
## [1] 1.644854
sum(Z1>qnorm(0.95))/1000
## [1] 0.912

検出力は90%を超えています。

3.4.3 標本設計

標本設計の重要性

標本数を \(n_1=10,n_2=10\) に減ずると検出力がどう変わるか見ておきましょう。.

S1=c()
S2=c()
n1 <- 10
n2 <- 10
for (i in 1:1000){
  x1 <- rnorm(n1,4,1);x2 <- rnorm(n2,3,1)
  x <- c(x1,x2)
  r <- rank(x, ties.method = "average")
  r1 <- head(r,n1)
  r2 <- tail(r,n2)
  S1<- c(S1,sum(r1))
  S2<- c(S2,sum(r2))
}
Z1  <- (S1-n1*(n1+n2+1)/2)/sqrt(n1*n2*(n1+n2+1)/12)
par(family = "HiraKakuProN-W3") 
hist(Z1,prob=T,main="",ylab="",col=0,xlab="")
curve(dnorm(x,0,1),add=T,col="blue")
abline(v = qnorm(0.95),col="red")

sum(Z1>qnorm(0.95))/1000
## [1] 0.68

検出力は、標本数が20ずつの場合に比べて小さいようです。 検出力が60%台では検定の信頼性に問題が残るでしょう。

3.5 マン・ホイットニーU検定.

3.5.1 U統計量

3.5.2 計算例.

＜Wilcox.test( )関数＞

Rには wilcox.test( ) 関数があり、２群を比較するノンパラメトリック統計量が出力されます。
どのような計算をしているのか見ておきましょう。

n1 <- 4
n2 <- 6
x1 <- rnorm(n1,3,1);x2 <- rnorm(n2,3,1)
x <- c(x1,x2)
r <- rank(x, ties.method = "average")
r1 <- head(r,n1)
r2 <- tail(r,n2)
wilcox.test(x1, x2)$ statistic
## W 
## 8
wilcox.exact(x=x1,y=x2,paired=F)$ statistic
## W 
## 8

＜順位和統計量＞

W1=sum(r1)
W2=sum(r2)
c(W1,W2)
## [1] 18 37

Wilcoxonの提案した順位和統計量は、２群のどちらも wilcox.test( )の統計量と一致しません。
wilcox.test( ) で出力されるのは、進化型のU統計量です。

＜U統計量＞

  • \(U_1=W_1-\frac{n_1(n_1+1)}{2}\).
  • \(U_2=W_2-\frac{n_2(n_2+1)}{2}\).
U1=sum(r1)-n1*(n1+1)/2
U2=sum(r2)-n2*(n2+1)/2
c(U1,U2)
## [1]  8 16

wilcox.test( ) では、U1が出力されています。

＜降順の場合＞.

デフォルトの順位は昇順です。降順にするとどうなるか見てみましょう。
rank(-x) とすれば、xの降順の順位が出力されます。

r <- rank(-x, ties.method = "average")
r1 <- head(r,n1)
r2 <- tail(r,n2)
U1=sum(r1)-n1*(n1+1)/2
U2=sum(r2)-n2*(n2+1)/2
c(U1,U2)
## [1] 16  8

降順では U1,U2の値が入れ替わるだけです。

3.5.3 U1,U2 の関係

組合せ計算を用いた理論分布の図において、\(U_1\)を青の縦線、\(U_2\) を赤の縦線で表示します。

data <- combn(x=c(1:10),m=4)
N <- choose(10,4)
S=c()
for (i in 1:N){
  s <- sum(data[,i])
  S <- c(S,s)
}
plot(table(S-4*5/2)/N,type="l",col="lightblue")
abline(v = 4*6/2,col="lightblue")
abline(v = U1,col="blue")
abline(v = U2,col="red")

\(U_1\)が右側有意なら、\(U_2\) は左側有意で、P値は等く、どちらも統計量として使用可能です。

3.5.4 もう一つの公式.

3.5.5 順位差による計算.

\(X_1\) の中での順位を \(r_{1i}\) 、２群をあわせた中での順位を \(R_{1i}\) とすると、

　　　\(U_1=\displaystyle\sum_{i=1}^{n_1}(R_{1i}-r_{1i})\)　　が成立します。

\(\sum\) の中の \((R_{1i}-r_{1i})\) については、

  • 順位が降順のとき、\(x_{1i}\)より大きい他方の群のデータ数
  • 順位が昇順のとき、\(x_{1i}\)より小さい他方の群のデータ数

を意味します。

このことを用いると、順位和計算は簡単になります。

3.5.5.1 順位が降順のとき

　　　

n1 <- 4
n2 <- 6
  x1 <- rnorm(n1,3,1);x2 <- rnorm(n2,3,1)
  list=list(A=x1, B=x2)
boxplot(list)
beeswarm(list,  pch = 16, add = TRUE)

箱ひげ図のドットを見て、順位差を根拠とした

  • \(U_1= \displaystyle \sum_i(x_{1i}\) より大きい他方の群のデータ数）.

による計算をしてください。
ほとんど黙視だけで、順位和が分かります。

同じ計算をRで行うと、

  s = 0
  for (j in 1:n1){
    s <- s+sum(x2 > x1[j])
  }
s
## [1] 4

集合Xと定数aに対し、sum(X > a) は、aより大きい X の要素の数をカウントします。

次に、U統計量の定義により、

  • \(U_1=W_1-\frac{n_1(n_1+1)}{2}\) を計算します。.
    順位はr(-x)により、降順に直します。
  x <- c(x1,x2)
  r <- rank(-x, ties.method = "average")
  r1 <- head(r,n1)
  r2 <- tail(r,n2)
  sum(r1)-n1*(n1+1)/2
## [1] 4

どちらの方法でも同じ値になることが確認されました。 　

3.5.5.2 順位が昇順のとき.

昇順の順位でもU統計量は計算できます。

箱ひげ図のドットを見て、順位差を根拠とした

  • \(U_1= \displaystyle \sum_i(x_{1i}\) より小さい他方の群のデータ数）.

を計算をしてください。

boxplot(list)
beeswarm(list,  pch = 16, add = TRUE)

同じ計算をRで行うと、

  s = 0
  for (j in 1:n1){
    s <- s+sum(x2 < x1[j])
  }
s
## [1] 20

集合Xと定数aに対し、sum(X < a) は、aより小さい X の要素の数をカウントします。

次に、U統計量の定義により、

  • \(U_1=W_1-\frac{n_1(n_1+1)}{2}\) を計算します。.
    順位はr(x)で昇順のままです。
r <- rank(x, ties.method = "average")#昇順
r1 <- head(r,n1)
r2 <- tail(r,n2)
U1=sum(r1)-n1*(n1+1)/2
U2=sum(r2)-n2*(n2+1)/2
U1
## [1] 20

やはり、どちらの方法でも同じ値になることが確認されました。 　

［注意］
昇順と降順では、\(U_1,U_2\)の値が入れ替割っていることに注意してください。

U1=sum(r1)-n1*(n1+1)/2
U2=sum(r2)-n2*(n2+1)/2
c(U1,U2)
## [1] 20  4

どちらを統計量としても問題ありませんが、 このような仕組みを手計算せずブラックボックスのままにしておいたのでは心配になります。

3.5.6 参考.

\(U_1,U_2\) の対称性は箱ひげ図からも分かります。

S1=c()
S2=c()
n1 <- 4
n2 <- 6
for (i in 1:10){
  x1 <- rnorm(n1,3,1);x2 <- rnorm(n2,3,1)
  x <- c(x1,x2)
  r <- rank(x, ties.method = "average")
  r1 <- head(r,n1)
  r2 <- tail(r,n2)
  S1<- c(S1,sum(r1-rank(x1)))
  S2<- c(S2,sum(r2-rank(x2)))
}
list=list(A=S1, B=S2)
par(family = "HiraKakuProN-W3") 
boxplot(list,main="順位差の和")
beeswarm(list,  pch = 16, add = TRUE)

\(U_1,U_2\) の分布は、12を中心に上下反転した形状です。.

１つ１つの公式を自分で確かめる。
順位和検定は、高校生の探究課題として格好の教材と考えます。