目次

・帰無仮説の棄却は100%正しいわけではない

・第１種の過誤をt分布から眺める

・第２種の過誤と検出力

・適切なサンプルサイズと効果量

・結局はデータに真摯に向き合うこと

帰無仮説の棄却は100%正しいわけではない

　以前の記事でt検定を例とした統計的仮説検定の基本概念について紹介をした。そこでは有意水準を決め、手元の標本で計算して有意水準を下回るようなp値が検出されたときに、帰無仮説を棄却、すなわち2標本の平均値は異なると、結論するのであった。これの根拠は、帰無仮説を正しいとしたときに、有意水準を下回るようなp値が出現する確率は小さい（めったに起こらない）から、そもそも帰無仮説を正しいと考えることに無理がある、と考えているからだ。しかし、めったに起こらないとはいえ、有意水準を下回るp値が出現する可能性は0ではない。だとすれば、有意水準を下回るp値が出現したとしても、真の状態は帰無仮説が正しい可能性も出てくる。このような、真の状態に対して、観測からそれとは異なる判断を下してしまうことを、過誤 errorと呼ぶ。特に、真の状態は帰無仮説が正しいのに、観測から帰無仮説は間違っていると判断してしまうことを、第１種の過誤 Type I errorと呼んでいる。有意水準を5%と定めれば、5%の確率でこの第１種の過誤を犯してしまうことになる。そういう意味で、帰無仮説の棄却は、通常の数学における証明のように、100%正しいと判断を下すことはできないのだ。また、逆に真の状態は帰無仮説は正しくないのに、観測から帰無仮説を棄却しないと判断してしまうことを第２種の過誤 Type II errorと呼ぶ。

第１種の過誤をt分布から眺める

　今回もt分布をシミュレーションで作成して、第１種の過誤と第２種の過誤について考えてみよう。Rでは以下のようにすれば、シミュレーションでt分布を作ることができる。なお、標本の母集団は平均50、標準偏差20の正規分布であり、２標本のサンプルサイズは20と30だ。

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

library(plotn)

#もし、plotnライブラリが欲しい場合は、

#以下のようにコマンドを打つ（#は消す）。devtoolsライブラリが必要。

#library(devtools)

#install_github("bugplant/plotn")

t <- NULL

for (i in 1:10000){ #定義に倣い統計量tを10000個算出

  m1 <- rnorm(20, mean = 50, sd = 20) #上記の正規分布から20個サンプリング

  m2 <- rnorm(30, mean = 50, sd = 20) #上記の正規分布から30個サンプリング

  ue <- ((20 - 1)*(sd(m1)^2) + (30 - 1)*(sd(m2)^2))/(20 + 30 - 2) #2個の分散を統合

  t <- c(t, (mean(m1) - mean(m2))/sqrt(ue*(1/20 + 1/30)))

  #統計量tの算出、分子を見るとわかるがこれは2つの標本の平均値の差に依存した量

}

head(t) #10000個のシミュレートした統計量tを一部表示

## [1] -0.3944782 -2.1391604  0.6972603 -0.1572087 -1.7567524  0.1391997

histn(t, xlab = "統計量t", ylab = "頻度")#統計量tをヒストグラムにして図示 = t分布、hist(t)でも可。

overdraw(abline(v = qt(df = 48, p = 0.975)),

         abline(v = qt(df = 48, p = 0.025)), #95%区間の表示

         axis(side = 1, at = c(qt(df = 48, p = 0.025), qt(df = 48, p = 0.975)))

         )

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

棄却域に含まれるtの個数も計算しよう。以下のコマンドで行う。今回の私のシミュレーションでは、たまたま棄却域にぴったり500個のtが含まれていた！　皆さんが手元で行うときは必ずしも同じに値にならないので注意だ。

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

#95%の境界線にあたるtの計算は以下のコマンドで行う。

#dfは自由度と呼ばれる値で、2標本のt検定なら、（2標本のサンプルサイズの合計　-　2）と計算する。

#pが境界を決める値で、両端から合計5%と設定しているので、半分の2.5%だけ差し引いた97.5%の境界としている。

qt(df = 48, p = 0.975)

## [1] 2.010635

#tの平均から95%の境界線にあたるtを足し引きした範囲に含まれるtの数を計算

sum(t < mean(t) + qt(df = 48, p = 0.975) &

      t > mean(t) - qt(df = 48, p = 0.975)) #95%区間に入っている統計量tの推定値の数

## [1] 9500

#5%区間に入っている統計量tの推定値の数

#95%の範囲に含まれるtの個数を10000個から引けば、両端の5%の範囲に含まれるtの個数となる。

10000 - sum(t < mean(t) + qt(df = 48, p = 0.975) & t > mean(t) - qt(df = 48, p = 0.975))

## [1] 500

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

t分布全体のうち、0からどれくらいの離れたtが出現したら、めったに起こらないこととして処理するか、の基準が有意水準や棄却域であった。定義から当たり前であるが、有意水準を5%とすれば、同じ母集団から2標本を生成したとしても、5%の確率で有意水準以下のp値となる。言い換えれば、5%の確率で、帰無仮説が正しいのに帰無仮説を棄却する危険がある、ということだ。そういう意味で、有意水準は危険率とも呼ばれる。また、そのまま第１種の過誤を犯す確率ともいわれる。有意水準を変えたとしても、状況は変わらない。1%とすれば、1%の確率で、10%とすれば10%の確率で、第１種の過誤を犯すことになる。

第２種の過誤と検出力

　では、第１種の過誤を犯す確率を低くするために、有意水準を厳しくする、5%ではなく、1%や0.1%にすることは最適な判断だろうか？　残念ながらそうではなく、有意水準を厳しくすればするほど第２種の過誤を犯す可能性が高まっていくのである。第２種の過誤を実感するためには、帰無仮説が正しくない、つまり２標本が平均の異なる母集団から生成される状況を考えなくてはいけない。そこで1つの標本は今までと同じく平均50の正規分布から、もう1つの標本は平均70の正規分布から生成された場合を考えてみよう。まずは2つの標本がどんな正規分布から生成されたか見てみる。

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

n1 <- rnorm(10000, mean = 50, sd = 20) #正規分布に従う乱数を10000個生成

n2 <- rnorm(10000, mean = 70, sd = 20) #正規分布に従う乱数を10000個生成

#上記2つは平均値が異なることに注意！

n <- as.data.frame(rep(c("n1","n2"), each = 10000))

n$value <- c(n1,n2)

colnames(n) <- c("pop","value")

histn(n$value ~ n$pop, xlab = "value", ylab = "頻度")#図2の図示

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

さて、平均が20程度異なる正規分布から標本が生成されるのが真の状態だとして、この時のt分布はどうなるだろうか。2標本の平均が異なっていようと、定義にならってtを計算してみよう。今回は、サンプルサイズ20のほうが平均値50の正規分布から生成されたと考えている。

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

t <- NULL

for (i in 1:10000){ #定義に倣い統計量tを10000個算出

  m1 <- rnorm(20, mean = 50, sd = 20) #平均50の正規分布から20個サンプリング

  m2 <- rnorm(30, mean = 70, sd = 20) #平均70の正規分布から30個サンプリング

  ue <- ((20 - 1)*(sd(m1)^2) + (30 - 1)*(sd(m2)^2))/(20 + 30 - 2)

  t <- c(t, (mean(m1) - mean(m2))/sqrt(ue*(1/20 + 1/30)))

}

histn(t, xlab = "統計量t", ylab = "頻度") #図3の図示

overdraw(abline(v = qt(df = 48, p = 0.975)),

         abline(v = qt(df = 48, p = 0.025)), #帰無仮説に基づいた95%区間の表示

         abline(v = qt(df = 48, p = 0.995), col = "red"),

         abline(v = qt(df = 48, p = 0.005), col = "red"), #帰無仮説に基づいた99%区間の表示

         axis(side = 1, at = c(qt(df = 48, p = 0.005), qt(df = 48, p = 0.025), 

                               qt(df = 48, p = 0.975), qt(df = 48, p = 0.995)))

         )

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

すると平均が約3.5 = (母集団の平均値の差)/√(母集団の分散×(1/サンプル数1 + 1/サンプル数2))程度ずれたt分布となる（厳密には、対称性も崩れている）。このとき、有意水準5%で帰無仮説を棄却できないと判断される、つまり-2.01 < t < 2.01の範囲に含まれるtの個数を調べてみよう。

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

sum(t > 0 - qt(df = 48, p = 0.975))#図3の黒線部分より右側に存在するtの個数

## [1] 778

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

今回の場合は2標本は異なる母集団から生成されたわけだが、778/10000 × 100 = 7.78%の確率で、帰無仮説を棄却できない、第２種の過誤を犯すことになる。ではここで、有意水準1%としたらどうだろう。第２種の過誤を犯すtの個数は以下のように計算できる。

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

sum(t > 0 - qt(df = 48, p = 0.995))#図3の黒線部分より右側に存在するtの個数

## [1] 2271

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

第２種の過誤を犯す確率は約23%まで上昇する。つまり、4～5回に1回は帰無仮説は正しくないのに棄却できないのだ。第１種の過誤と同じく、これでは真の状態を正しく推定することはできない。ゆえに、有意水準をただ厳しくすれば、間違いを犯す確率が減少する、というわけではないのだ。そこで、第１種の過誤を犯す確率を減らしつつも、第２種の過誤を犯す確率も減らしたい、というときに登場する概念が、検出力 powerだ。検出力は、帰無仮説を正しく棄却できる確率であり、1 - （第２種の過誤を犯す確率）として計算できる。上記の例なら、有意水準5%の時の検出力は(10000 - 778)/10000 × 100 = 92.2%、1%の時は77.3%だ。有意水準は低くしつつ、検出力は高くする、というのが検定の鉄則となる。有意水準は伝統的に5%とされるが、検出力は伝統的に80%とされることが多いようだ。とりあえず、標本の平均が20程度違う標準偏差20の正規分布から生成されたと考えられるとき、有意水準は5%としても問題ないようだ。

　もう少し平均値が異なる母集団から標本が生成される例を考えてみよう。今度は平均値の差が5程度しかない場合だ。

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

n1 <- rnorm(10000, mean = 50, sd = 20) #正規分布に従う乱数を10000個生成

n2 <- rnorm(10000, mean = 55, sd = 20) #正規分布に従う乱数を10000個生成

#上記2つは平均値が異なることに注意！

n <- as.data.frame(rep(c("n1","n2"), each = 10000))

n$value <- c(n1,n2)

colnames(n) <- c("pop","value")

histn(n$value ~ n$pop, xlab = "value", ylab = "頻度") #図4の図示

t <- NULL

for (i in 1:10000){ #定義に倣い統計量tを10000個算出

  m1 <- rnorm(20, mean = 50, sd = 20) #平均50の正規分布から20個サンプリング

  m2 <- rnorm(30, mean = 55, sd = 20) #平均55の正規分布から30個サンプリング

  ue <- ((20 - 1)*(sd(m1)^2) + (30 - 1)*(sd(m2)^2))/(20 + 30 - 2)

  t <- c(t, (mean(m1) - mean(m2))/sqrt(ue*(1/20 + 1/30)))

}

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

overdraw(abline(v = qt(df = 48, p = 0.975)),

         abline(v = qt(df = 48, p = 0.025)), #帰無仮説に基づいた95%区間の表示

         axis(side = 1, at = c(qt(df = 48, p = 0.025), qt(df = 48, p = 0.975)))

         )

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

すると平均が約0.87 = (母集団の平均値の差)/√(母集団の分散×(1/サンプル数1 + 1/サンプル数2))程度ずれたt分布となる。このとき、有意水準5%で帰無仮説を棄却できないと判断される、つまり-2.01 < t < 2.01の範囲に含まれるtの個数を調べてみよう。

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

sum(t > 0 - qt(df = 48, p = 0.975))#図5の黒線部分より右側に存在するtの個数

## [1] 8594

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

今度は打って変わって、有意水準5%であっても、第２種の過誤を犯す確率は約86%、検出力は約14%となり、極めて心もとない。しかし、これは当然のことで、平均値の差がわずかな標本において差を検出するのは難しいという、直感的な結果だ。このような時は、サンプルサイズを大きくすれば検出力を大きくすることができる。サンプルサイズをともに10倍して、200と300にしてみよう。

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

t <- NULL

for (i in 1:10000){ #定義に倣い統計量tを10000個算出

  m1 <- rnorm(200, mean = 50, sd = 20) #平均50の正規分布から200個サンプリング

  m2 <- rnorm(300, mean = 55, sd = 20) #平均55の正規分布から300個サンプリング

  ue <- ((200 - 1)*(sd(m1)^2) + (300 - 1)*(sd(m2)^2))/(200 + 300 - 2)

  t <- c(t, (mean(m1) - mean(m2))/sqrt(ue*(1/200 + 1/300)))

}

histn(t, xlab = "統計量t", ylab = "頻度") #図6の図示

overdraw(abline(v = qt(df = 498, p = 0.975)),

         abline(v = qt(df = 498, p = 0.025)), #帰無仮説に基づいた95%区間の表示

         axis(side = 1, at = c(qt(df = 498, p = 0.025), qt(df = 498, p = 0.975)))

         )

sum(t > 0 - qt(df = 498, p = 0.975)) #図6の黒線部分より右側に存在するtの個数

## [1] 2174

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

サンプルサイズを10倍すると、有意水準5%で、第２種の過誤を犯す確率は約22%、検出力は約78%となる。検出力80%の基準に到達させるためには、もっとサンプルを集めなければいけないようだ。このように平均値が異なる母集団から標本が生成されたと考える場合、理論上はサンプルサイズを大きくしていけば、どんなに小さな平均値の差でも検出力が上昇する、つまり帰無仮説を棄却しやすくなっていく。

適切なサンプルサイズと効果量

　では、サンプルサイズを極めて大きくすることが検定するうえでの適切な応対か、というとそうではない。第1に、調査に割ける労力や資源は有限だ。生物学の範疇で研究しているとすぐに直面するが、フィールドワークや実験に割ける時間は有限だ。さらには野生生物を捕獲して調べるなら、すべて調査してしまったら自然界に多大な影響を与えてしまいかねない。このあたりは、工場の部品の精度を測る、という状況とはわけが違う。研究者は、コストを念頭に置いた調査が求められている。

　第2に、帰無仮説が棄却されたことと、差に意味があることは別問題だということだ。例えば、図4のように平均値の差が5である標準偏差20の正規分布が真の状態であると信じ、300以上のサンプルサイズをそろえて検出力の問題をクリアして、差を見出すと考えてもよい。しかし、そこまでしないと出ない差は、生物学的に何か影響のある差なのだろうか、ということを常に念頭に置いておく必要がある。一方で、平均値の差は5であり、標準偏差が5の正規分布が真の状態であるときを考えてみよう。

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

n1 <- rnorm(10000, mean = 50, sd = 5) #正規分布に従う乱数を10000個生成

n2 <- rnorm(10000, mean = 55, sd = 5) #正規分布に従う乱数を10000個生成

#上記2つは平均値が異なることに注意！

n <- as.data.frame(rep(c("n1","n2"), each = 10000))

n$value <- c(n1,n2)

colnames(n) <- c("pop","value")

histn(n$value ~ n$pop, xlab = "value", ylab = "頻度") #図7の図示

t <- NULL

for (i in 1:10000){ #定義に倣い統計量tを10000個算出

  m1 <- rnorm(20, mean = 50, sd = 5) #平均50の正規分布から20個サンプリング

  m2 <- rnorm(30, mean = 55, sd = 5) #平均55の正規分布から30個サンプリング

  ue <- ((20 - 1)*(sd(m1)^2) + (30 - 1)*(sd(m2)^2))/(20 + 30 - 2)

  t <- c(t, (mean(m1) - mean(m2))/sqrt(ue*(1/20 + 1/30)))

}

histn(t, xlab = "統計量t", ylab = "頻度") #図8の図示

overdraw(abline(v = qt(df = 48, p = 0.975)),

         abline(v = qt(df = 48, p = 0.025)), #帰無仮説に基づいた95%区間の表示

         axis(side = 1, at = c(qt(df = 48, p = 0.025), qt(df = 48, p = 0.975)))

         )

sum(t > 0 - qt(df = 48, p = 0.975)) #図8の黒線部分より右側に存在するtの個数

## [1] 698

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

この時、サンプルサイズが20と30であっても、有意水準5%での検出力は(10000 - 698)/10000 × 100 = 93%と申し分ない状態だ。こんなときはもちろんサンプルサイズを大きくする意味はあまりない。こんな時でもこの平均値の差には意味があるのか、と考えることは大切だ。とは言えども、サンプルサイズが小さくても出てしまう差には、何か重要な要因が隠れているかもしれない、と一考に値する情報ではある。

　このように有意水準、検出力、帰無仮説を棄却できたかそうでないか、だけでは解釈の落とし穴にはまってしまう。そういう意味で、サンプルサイズの情報は重要であり、適切なサンプルサイズを考えることは必須なのである。とりあえず、ここでは統計の感覚をつかむことに徹してほしいので、具体的にどうサンプルサイズを決めるかは、本稿では申し訳ないが解説しない。

　これまでの議論で、サンプルサイズが大きくなれば小さな差でも帰無仮説を棄却できてしまう、と述べた。では、サンプルサイズに依存せず、2標本にどれくらいの検出しやすい差があるか、を見れる指標を導入しよう。それが効果量 effect sizeだ。いくつか計算方法があり、詳細については各自調べてほしいが、1つの例としては以下のように計算する。

Hedgesのgと呼ばれる指標は、よく見ると検定統計量tとほとんど同じ形をしている。tに√（1/nA + 1/nB）を掛ければ計算できる。つまり、2標本の平均値が大きく、分散が小さいほど効果量は大きくなる。効果量が大きいほど、検出しやすい差、ということになる。今までの、サンプルサイズが20と30の2つの標本が平均値50と55で標準偏差20の正規分布から生成される場合（図5）と、標準偏差5の正規分布から生成される場合（図8）で考えると、t分布に√（1/20 + 1/30）= 0.289くらいの値を掛ければ効果量の分布となる。効果量の平均は大体、平均値50と55で標準偏差20の正規分布から生成される場合でg = 0.79、標準偏差5の正規分布から生成される場合でg = 1.02くらいとなる。母集団の標準偏差が小さいほど検出しやすくなっていることを反映している。

結局はデータに真摯に向き合うこと

　長々と議論してきたが、結局のところ、真の状態を見極めたいと願うのであれば、データに真摯に向き合うしかない。帰無仮説を棄却できたかどうかだけではなく、サンプルサイズや効果量も念頭において解析を進めていこう。p値至上主義、という言葉がある。帰無仮説を棄却できたかできなかったかばかりを重要視する立場を憂慮した表現だ。帰無仮説を棄却できたからと言って、差に意味があるとは限らないし、棄却できなかったからと言って意味のない差とは限らない。検定はあくまで数学的な手続きであり、バックにあるメカニズムに関する正しさを保証してくれるわけではない。差に意味を見出すかどうかは、研究者の手腕の見せどころなのである。