同じ母集団からサンプリングしても100回に1回くらいは統計学的有意差が出てしまう話

同じ母集団から2群をランダムにサンプリングしても、100回に1回くらいは同じ母集団なので本当は統計学的有意な差がでないはずなのに統計学的有意差（p<0.05）が出てしまうというという話を聞いたことがあるでしょうか。

例をあげるなら、平均的な身長の成人男子100人（つまり平均値で正規分布しているはず）を集めてその人達に適当にナンバリングして1から100まで並んでもらって前50人をA群、後ろ50人をB群としたときに2群の身長に統計学的有意差があるか検証したときに、A群、B群ともに同じ集団からサンプリングされた2群なので統計学的有意差はでないのが2群の性質上正しいのだけれど、ランダムサンプリングされたときにたまたま両群の身長に統計学的有意差が出てしまうことが100回に1回くらいはありえるよ、という話でした。

私はこの話（上記の例）を統計の授業で、確か過誤の話のおまけで聞いたと記憶していますが、配布されたレジュメなどを見てもそんな話は載っていないのでランダムサンプリングして統計学的有意差がでてもそのことの合理性が説明できないときは要注意だよ、というおまけ話以上でもそれ以下でもないのかもしれません。

たしかに、たまたまランダムサンプリングの結果で身長が小さい順から並べてA群50人、51番目から100番目をB群とすると統計学的有意差が出そうです。

しかし「100回に1回くらい」というところの根拠がいまいちわからなくて「ふーん」程度に思っていました。

数学が得意だったら手計算でもこの根拠を求めることができるかもしれませんが、数学力の乏しいので、Rのちからを借りてこの「100回に1回くらい」という条件を検証してみました。

ちなみに、結論としては「100回に1回くらい」は本当は統計学的有意差が出ないはずの同一の母集団から2群をランダムサンプリングしたときに統計学的有意差はでるだろう、という結果でした。

さらにいうと多分ランダムサンプリング20～60回くらいに1回くらいで結構な確率ででてきそうという結果でした（もちろんランダムサンプリングなので1～数回ででることもある）。

検証方法

検証法としては、

• ランダム発生させたサンプルデータを用いてシミュレーションする
• 力技で全サンプリング方法についてp値を計算して確認する

の2つの方法が考えついたので、両方で検証してみたいと思います。

ランダム発生させたサンプルデータを用いてシミュレーションする

方法その1のランダム発生させたサンプルデータを用いてシミュレーションしてみる方法をやってみます。

シミュレーションのイメージとしては1～10000の数字が書かれたボールが入っている袋からボールを順番に取り出して、最初の50番目までに取り出されたボールをA群、51番目から100番目までに取り出されたボールをB群としてその数字を得点としたときに2群間に統計学的有意差が出るまでそれを繰り返しその回数をカウントする、というものです。

2群とも母集団の平均値は同じ（同じ母集団からサンプリングしている）なので、ランダムに50個取り出したボールの得点も平均値に近づく……はずですが、ランダムサンプリングなので前述したようにランダム故に得点に無作為な偏りが出てしまう可能性があります。

そこまで至るまでの回数をカウントします。

検定には母集団について前提条件を厳しく求めないBrunner-Munzel検定（lawstatパッケージのbrunner.manzel.test()）を用いました。

他の検定方法でも検証しましたほぼ同じような結果でした。

検証用自作関数

library("lawstat")
RS.p.bmt <- function(n = 10, times = 100){
  p.bmt <- NULL
  p0.05times <- NULL
  for (i in 1:times) {
    p.bmt[i] <- 1
    p0.05times[i] <- 0
    while (p.bmt[i] >= 0.05) {
      z <- sample(1:(n^2), n*2, replace = F)
      x <- z[c(1:n)]
      y <- z[-c(1:n)]
      p.bmt[i] <- brunner.munzel.test(x,y)$p.value
      p0.05times[i] <- p0.05times[i]+1 
    }
  }
  return(p0.05times)
}

これで（デフォルトではデータ数10個から2群抽出）pが0.05未満になるまでwhileでp.0.05timesをカウンターにしてループし続けて指定回数までシミュレーションし続ける関数ができました。

これを1000回シミュレーションしてp<0.05に至るまでの回数を求めてヒストグラムを描くと

library("magrittr")
 RS.p.bmt(n=100, times = 1000) %>% hist(breaks = seq(0,140,1))

fig.1

となり、おおよそ少なくとも100回に1回くらいは統計学的有意差がでる様子。

ランダムに発生させているサンプルなので実行するたび結果は変わるけれどおおよそ少なくとも100回に1回くらいは統計学的有意差が出る傾向にある感じ。

ちなみに、母集団を1～100でサンプルを50個ずつにした条件でもほぼ同じ結果になる。

ヒストグラムを見るかぎり、印象として40回くらいランダムサンプリングを繰り返せば1回は統計学的有意差が出そうな感じがする。

ちなみに上の自作関数はたまにエラーを出す（シミュレーション回数を増やすと結構でる）けど、多分brunner-munzel検定でエラーが出ているからだと思う。

t検定（Welch）でも、ウィルコクソンの順位和検定でもほぼ同じような結果だった。

力技で全サンプリング方法についてp値を計算して確認する

力技で、ランダムサンプリングしたp値を求めまくって、p<0.05となる数を数える方法は、条件として順列バージョンと組み合わせバージョンを作ってみたけど、結局どちらも結果として一緒だった。

順列バージョンは尋常じゃない計算回数でサンプル数は容易に増やせない。

組み合わせバージョンの方は計算早いが、後述する理由で大きなサンプルで計算できない。

あと、サンプル数をちょっと増やしただけで劇的に計算回数が増えてしまうので1～10から5個ずつで2群抽出とした。

順列バージョン

1～10を並べ替えて前半と後半で2群に分けて、全例でt検定（Welch）を行いp<0.05の数を数える。

ウィルコクソンの順位和検定の方が良かったかもしれないけど、10の階乗が

factorial(10)
[1] 3628800

で362万回も計算をもう一回やるのが辛いから……

検証用自作関数

ALL.RS.p.tt <- function(n = 10){
  N <- n*2
  dat <- permutations(N)
  times <- nrow(dat)
  p.tt <- NULL
  for(i in 1:times){
    x <- dat[i,c(1:n)]
    y <- dat[i,-c(1:n)]
    p.tt[i] <- t.test(x,y)$p.value
  }
  return(p.tt)
}

デフォルトでは1～20を並び替えて10ずつの2群にしようと考えたが、計算過程で3.6Gbメモリ不足だよと怒られてそんなに割り当てられなかったので、1～10で5個ずつの2群にした。

p.tt <- ALL.RS.p.tt(n=5)
length(p.tt[p.tt<0.05])
[1] 201600
length(p.tt[p.tt<0.05])/length(p.tt)
[1] 0.05555556

362万回もt検定を計算するのでかなり時間がかかる。

結果として約20万パターンで統計学的有意差がでて、それは全サンプリングパターンの5.6％に相当するとのこと。

少なくともi回サンプリングして1回統計学的有意差がでる確率を求めてプロットすると

x <- NULL
for(i in 1:100){
  x[i] <- 1-(1-length(p.tt[p.tt<0.05])/length(p.tt))^i
}
plot(x)

fig.2

となり100回もやれば少なくとも1回は統計学的有意差がでそう。

ちなみに60回で96.8％の確率で1回は統計学的有意差がでる。

後述する問題のため自作関数では組み合わせバージョンでは大きなサンプルを扱えないので、計算資源に余裕があれば順列バージョンでなら大きなサンプルで計算可能……

組み合わせバージョン

順列バージョン回している最中（かなり長い）に「今回の条件（対応なし）なら順番関係ないから組み合わせでも同じじゃん！？（時既に遅し）」と思った。

1～10を並べ替えて前半と後半で2群に分けて、全例でウィルコクソンの順位和検定を行いp<0.05の数を数える。

検証用自作関数

AC.RS.p.wt <- function(n = 10){
  N <- n*2
  dat <- combn(c(1:N),n)
  dat2 <- c(1:N)
  times <- ncol(dat)
  c.p.wt <- NULL
  for(i in 1:times){
    x <- dat[c(1:n),i]
    y <- dat2[-c(x)]
    c.p.wt[i] <- t.test(x,y)$p.value
  }
  return(c.p.wt)
}

計算回数が少ないので順列よりだいぶ早い。

あとは順列バージョンと同じ

c.p.wt <- AC.RS.p.wt(n=5)
length(c.p.wt[c.p.wt<0.05])
[1] 14
length(c.p.wt[c.p.wt<0.05])/length(c.p.wt)
[1] 0.05555556
 x <- NULL
for(i in 1:100){
   x[i] <- 1-(1-length(c.p.wt[c.p.wt<0.05])/length(c.p.wt))^i
 }
plot(x)

fig.3

で順列バージョンと同じ結果が出る。

つまりランダムサンプリングを100回くらいやると少なくとも1回は統計学的有意差がでる。

ついでに60回で98.6％の確率で1回は統計学的有意差が出る。

サンプル数を増やしてみてみたいが、combnに大きな数を入れるとエラーを返すので大きいサンプルバージョンはこの自作関数のままだと求められない。