有意水準と検出力について考えてみよう

はじめに

前回の記事では、ランダム化比較試験について架空の事例をもとに考えていきました。

その中で、両側5%の有意水準(Significance level)で検定を行った場合、「本当はグループ間に差がないのに間違って差があると判定する確率は5%になる」という話をしました。
このあたりの説明は本題から逸れるのでかなりすっ飛ばしてしまいましたが...改めて自分なりにまとめるのも良いかも、と思いました。

一方で、サンプルサイズが非常に大きい場合には、グループ間の差がごくわずかであっても差があると判定できてしまう、いわゆる「統計的に有意」と結論できてしまうということもちらっと述べました。
例えば、大変な忍耐を必要とするような勉強法が開発されたとして、それを頑張ってやり遂げると偏差値が1上がることが期待できます、ということが統計的に有意に示されたとしましょう。でも、そんな勉強法を喜んでやる人はいるでしょうか？

f:id:mstour:20210425202346j:plain

本来、検定によってグループ間に差があるかどうかを評価するためには、少なすぎも多すぎもしない適切なサンプルサイズを設定した上で行う必要があります。
この時にカギとなる概念が「検出力(Power)」です。

今回は、有意水準と検出力について、シミュレーション計算も交えてちょっとだけ考えていくことにします。

なお、今回の話はいわゆる「ネイマン・ピアソン(Neyman-Pearson)流の検定」と呼ばれる、我々が学校などで教わることの多い検定方式を前提としています。つまり、否定したい仮説である「帰無仮説(Null hypothesis)」と、帰無仮説と相反する仮説である「対立仮説(Alternative hypothesis)」の2つの仮説を設定し、実験などによって得られたデータからどちらの仮説を採用するかを決定する、という選択問題として検定を考える立場です。
一方で、統計的検定の方法論を最初に確立したフィッシャー(Fisher)の立場はそれとは異なり、対立仮説というものは出てきません。
有意水準や検出力が決定的な意味を持つのは前者のネイマン・ピアソンの考え方で検定を行う場合です*1。

ちょっとややこしい話ですが、統計的有意性やp値に関する混乱の原因の一つとして、検定に関する異なる2つの立場が（多くの場合、明確に意識されることなく）存在していることがあるかと思います。
本題から外れてしまうので、この問題はまたいつか考えることにします。

有意水準と検出力についておさらい

有意水準

検定の手順として、まずは「帰無仮説が本当は正しいのに、間違ってそれを否定する（=対立仮説が正しいと間違って結論づける）」ことを可能な限り避けるように、その確率を適切にコントロールすることを考えます。間違って帰無仮説を否定することを「第一種の過誤(Type1 error)」と呼びます。
「適切にコントロールする」ために、第一種の過誤確率の上限値を設定しておきます。この上限値を有意水準といいます。

なお一般的な検定では、理論的には第一種の過誤確率が有意水準に等しくなるように定義されています*2。
例えば、帰無仮説を $H_0$ 、検定に使用する検定統計量(Test statistic)を $T$ とすると*3、有意水準 $\alpha$ （右辺）と「帰無仮説が本当は正しいのに、間違ってそれを否定する確率」（左辺）は以下のように表現できます：

$\displaystyle P(|T| \geq t_{\alpha} | H_0) = \alpha$

ここで $t_{\alpha}$ は検定の棄却点(Critical value)といい、通常は検定統計量 $T$ が従う確率分布*4において、この点よりも大きな値をとる確率が $\alpha \times 100$ %になるような点、という意味で用います。
この点を検定統計量が超えた場合に帰無仮説を否定することになります。確率を表す記号 $P$ の中の「 $| H_0$ 」というのは、「帰無仮説 $H_0$ が正しかった場合に」ということを表す数学的な書き方です。

つまりこの式が意味しているのは、「帰無仮説が正しい時に、（間違って）帰無仮説を否定する確率は $\alpha$ である」と言うことです。上記の式が成り立つためには、 $\alpha$ の値が変わる時には $t_{\alpha}$ の値も同時に変更する必要があります。

よく言われる「有意水準5%」とは、上の式に出てくる $\alpha$ を0.05と設定したものです*5。
しつこいようですが、有意水準5%の検定とは、「帰無仮説が本当は正しいのに、間違ってそれを否定する確率が5%であるような検定」ということになります。

検出力

検出力は、一言でいうと「正しく帰無仮説を否定できる確率」のことを指します。
逆に言うと、1から検出力を引いたものは「正しく帰無仮説を否定できない確率」、つまり「間違って対立仮説を採用する確率」を意味します。この「正しく帰無仮説を否定できない（=間違って対立仮説を採用する）」という事象は「第二種の過誤(Type2 error)」と呼ばれます。
第二種の過誤を起こす確率は $\beta$ で表すことが多く、その場合検出力は $1 - \beta$ で表されます。

f:id:mstour:20210425202414j:plain

おおざっぱに言うと上記の通りなのですが、検出力は有意水準とは異なる点があります。
有意水準は、あらかじめ5%などの値に固定することが検定の大前提です。一方、検出力とは、正確に言うと「有意水準が $\alpha$ で、サンプルサイズが $N$ で、データのばらつきが $\sigma$ で、（2グループ間の比較の場合）グループ間の真の差が $\Delta \neq 0$ （つまり、グループ間に差がない、という帰無仮説は間違い）である場合に、何%の確率で正しく帰無仮説を否定することができるか」を表しています。
したがって、帰無仮説を否定できる（=いわゆる「有意差がある」と判断する）確率は以下のファクターに影響を受けます。

有意水準 $\alpha$
サンプルサイズ $N$
データのばらつき $\sigma$
（グループ間の比較を行う場合）グループ間の真の差 $\Delta$

有意水準が大きくなる（第一種の過誤確率が高くなる、基準の甘い判定）と、帰無仮説を否定しやすくなるため検出力は上がります。
またサンプルサイズが大きくなると、推定のばらつきが小さくなることを通して、検出力が高くなります。データのばらつきが小さいことも、推定のばらつきが小さくなることに直結するので同様です。
さらにグループ間の差が大きければ大きいほど、帰無仮説に反する結果を得やすくなるため、やはり検出力が高くなります。

つまり、真の差が同じでも、サンプルサイズが大きくなれば検出力が上がり、帰無仮説を否定しやすくなるなるわけです。
極端な話、冒頭で述べたような偏差値を1だけ上げるような勉強法があったとすると、実用上の意味に乏しくてもサンプルサイズを増やしさえすれば、効果がないという帰無仮説を否定し、いわゆる「統計的に有意」という主張ができてしまうことになります。

このような問題があるため、検定で仮説を検証しようとする時には、本当に意味のある差がある場合だけ「統計的に有意」と判定できるよう、事前に「サンプルサイズ計算」を行うことが非常に重要になります。
むしろ適切なサンプルサイズ計算が行われていないような場合には、不用意に検定を行うべきではないと考えた方が良いかもしれません*6。

シミュレーションで理解する

検出力について、より具体的にシミュレーション実験で確認してみましょう。
前回の記事と同じように、ある設定のもとでデータを生成し、検定を行うというプロセスを反復します。
検出力は「正しく帰無仮説を否定できる確率」なので、反復回数を1000回とすると、1000回中何回帰無仮説を否定するという結論になったかを見れば理論的な値に近い数値が得られることが期待できます*7。

冒頭の話を踏まえ、今回は勉強法Aを実施したグループ（グループA）と勉強法Bを実施したグループ（グループB）の「偏差値」の平均*8の差を検定するという場面を想定してみます*9。

まずは前準備を行っておきます。

library(tidyverse)
library(gtsummary)

set.seed(1234)
nsim = 1000

mean_a = 50
mean_b = 45
mean_c = 49
sd = 10

1番目のブロック：使用するパッケージを読み込みます。データフレーム加工のためにtidyverse、結果の要約にgtsummaryを使用しました。
2番目のブロック：乱数のシード値と、シミュレーション回数(1000回に設定)です。
3番目のブロック：mean_aはグループAの平均、mean_bとmean_cはグループBの平均です（いくつかのパターンでシミュレーションするため、2通りの値を用意しました）。またsdは両グループの標準偏差 *10です。

シミュレーション・パターン1

最初に、「狙った検出力を出せるようにサンプルサイズの計算を事前に行った」という想定でやってみます。
例えば、グループAとBの差が5、標準偏差が10であると見込まれる時に両側5%の有意水準で検定を行うと、検出力80%を達成するにはサンプルサイズはいくつ必要か？を事前に計算しておき、算出されたサンプルサイズのデータを得るような実験を計画した、という場面を想定します。

diff = mean_a - mean_b
power.t.test(delta = diff, sd = sd, power = 0.80, sig.level = 0.05,
             alternative = "two.sided", type = "two.sample")

上記のようにすると、この設定で必要なサンプルサイズが計算できます。詳しいことはまた別の機会に記事にしようと思いますので、ここでは計算結果のみ以下に示します。

     Two-sample t test power calculation 

              n = 63.76576
          delta = 5
             sd = 10
      sig.level = 0.05
          power = 0.8
    alternative = two.sided

NOTE: n is number in *each* group

計算の結果、各グループのサンプルサイズとして63.77必要なことがわかりました。
検出力が期待する値を下回ることのないように、通常は小数点以下を切り上げた値を採用します。
したがって、シミュレーション・パターン1では、各グループの人数が64人になるように実験を計画した、という設定にします。

シミュレーションは以下のコードにより行いました。

# シミュレーション1
npergroup1 = 64
pvalue1 = NULL
for (i in 1:nsim){
  group = c(rep("A", npergroup1), c(rep("B", npergroup1)))
  val_a = rnorm(n = npergroup1, mean = mean_a, sd = sd)
  val_b = rnorm(n = npergroup1, mean = mean_b, sd = sd)
  val = c(val_a, val_b)
  simdata = data.frame(group = factor(group), val = val)
  
  ttest = t.test(val ~ group, data = simdata, var.equal = TRUE)
  pvalue1 = c(pvalue1, ttest$p.value)
}

result_df1 = data.frame(pvalue = pvalue1) %>%
  mutate(reject = case_when(
    pvalue < 0.05 ~ 1,
    TRUE ~ 0)
  )

for内の1番目のブロック：正規分布に従う乱数を用いて模擬データを生成しています。先のサンプルサイズ計算で想定した通りにグループAの平均が50、グループBの平均が45、標準偏差が両グループともに10となるようにしました。
for内の2番目のブロック：2グループ間の平均差についてt検定を行い、p値を計算しています。
for後のブロック：1000回の結果のうち、p値が0.05未満だったものに1をつけます。

結果はgtsummaryを用いて以下のようにまとめました。計算の過程で使ったp値についても、おまけとして中央値と最小値・最大値を載せました。

tbl_summary(result_df1,
            type = list(reject ~ "continuous"),
            statistic = list(pvalue ~ "{median}[{min}, {max}]",
                             reject ~ "{mean}"),
            digits = list(pvalue ~ 4,
                          reject ~ 2),
            label = list(pvalue ~ "p-Value",
                         reject ~ "Power")
            ) %>%
  modify_header(
    update = list(
      label ~ "",
      stat_0 ~ "**Simulation1 : Delta=5, N = 64**"
    )
  )

シミュレーション結果は以下のようになりました。

f:id:mstour:20210425202434j:plain

「Power」の欄に、p値が0.05未満になった割合(=両側有意水準5%で検定した時に帰無仮説を否定できた割合)を示しました。これがシミュレーションから推定される検出力を表しています。
検出力が80%になるように計算されたサンプルサイズで実験を行った場合、期待通りの検出力が確保できそうだということが分かります。
「p-Value」の欄を見ると、p値の中央値は0.0069と小さな値となっていますが、最大値は0.9358であり非常に大きな値が出ることもありました。実際、2グループに差があるような場合、p値は右側に裾の長い分布を示すことが知られています。

シミュレーション・パターン2

次に、各グループの平均・標準偏差は同じだが、サンプルサイズを各グループ2000に増やした場合をやってみます。サンプルサイズが増えれば検出力は高くなるはずです。
コードはほぼ同様です。

# シミュレーション2
npergroup2 = 2000
pvalue2 = NULL
for (i in 1:nsim){
  group = c(rep("A", npergroup2), c(rep("B", npergroup2)))
  val_a = rnorm(n = npergroup2, mean = mean_a, sd = sd)
  val_b = rnorm(n = npergroup2, mean = mean_b, sd = sd)
  val = c(val_a, val_b)
  simdata = data.frame(group = factor(group), val = val)

  ttest = t.test(val ~ group, data = simdata, var.equal = TRUE)
  pvalue2 = c(pvalue2, ttest$p.value)
}

result_df2 = data.frame(pvalue = pvalue2) %>%
  mutate(reject = case_when(
    pvalue < 0.05 ~ 1,
    TRUE ~ 0)
  )

tbl_summary(result_df2,
            type = list(reject ~ "continuous"),
            statistic = list(pvalue ~ "{median}[{min}, {max}]",
                             reject ~ "{mean}"),
            digits = list(pvalue ~ 4,
                          reject ~ 2),
            label = list(pvalue ~ "p-Value",
                         reject ~ "Power")
            ) %>%
  modify_header(
    update = list(
      label ~ "",
      stat_0 ~ "**Simulation2 : Delta=5, N = 2000**"
    )
  )

f:id:mstour:20210425202449j:plain

予想通り検出力は非常に高くなり、1000回のシミュレーションの全てで帰無仮説を否定できる結果になりました。
ではサンプルサイズを増やしさえすればよいかと言うと、そうとも限らない問題があります。

シミュレーション・パターン3

今度は、サンプルサイズをシミュレーション・パターン2と同じく各グループ2000とし、グループAの平均を50、グループBの平均を49に設定してみます（標準偏差はこれまでと同様に各グループで10）。
2グループの平均の差は1であり、当初よりかなり小さい設定です。

# シミュレーション3
npergroup3 = 2000
pvalue3 = NULL
for (i in 1:nsim){
  group = c(rep("A", npergroup3), c(rep("B", npergroup3)))
  val_a = rnorm(n = npergroup3, mean = mean_a, sd = sd)
  val_b = rnorm(n = npergroup3, mean = mean_c, sd = sd)
  val = c(val_a, val_b)
  simdata = data.frame(group = factor(group), val = val)
  
  ttest = t.test(val ~ group, data = simdata, var.equal = TRUE)
  pvalue3 = c(pvalue3, ttest$p.value)
}

result_df3 = data.frame(pvalue = pvalue3) %>%
  mutate(reject = case_when(
    pvalue < 0.05 ~ 1,
    TRUE ~ 0)
  )

tbl_summary(result_df3,
            type = list(reject ~ "continuous"),
            statistic = list(pvalue ~ "{median}[{min}, {max}]",
                             reject ~ "{mean}"),
            digits = list(pvalue ~ 4,
                          reject ~ 2),
            label = list(pvalue ~ "p-Value",
                         reject ~ "Power")
) %>%
  modify_header(
    update = list(
      label ~ "",
      stat_0 ~ "**Simulation3 : Delta=1, N = 2000**"
    )
  )

f:id:mstour:20210425202503j:plain

検出力の推定値は0.87であり、かなりの確率で帰無仮説が否定されることになります。

繰り返しになりますが、このように必要以上にサンプルサイズが増えすぎると、実用上あまり意味がないような差であっても「統計的に有意」と判定する可能性が高くなります。

おわりに

検定にまつわる話は、個人的な感覚ですが、統計モデリングの手法と比べて正しい理解がなかなか難しいように思います（私自身そうです）。
また誤用や解釈に関する議論の最も多いトピックの1つであることは間違いありません。
少なくとも、検定やp値「だけ」で何かを判断しようとしないようにするのが無難な態度ではないか、と考えています。

今回の記事は以下の資料を参考にしました。
[1] Gibson E.W.(2021) "The Role of p-Values in Judging the Strength of Evidence and Realistic Replication Expectations", Statistics in Biopharmaceutical Research, 13, 6-18.
https://www.tandfonline.com/doi/full/10.1080/19466315.2020.1724560
[2] 柳川尭(2017) "p値は臨床研究データ解析結果報告に有用な優れたモノサシである", 計量生物学, 38, 153-161.
https://www.jstage.jst.go.jp/article/jjb/38/2/38_153/_pdf
[3] 汪金芳, 桜井裕仁(2011), "ブートストラップ入門", 共立出版.
[4] John M. Lachin(2020), "医薬データのための統計解析", 共立出版.

*1:後で書くように、有意水準が「間違って対立仮説を選択する確率」、検出力が「正しく対立仮説を選択する確率」に対応します。フィッシャーの立場でも有意水準は考えますが、検出力という概念はありません。

*2:理論的にはそうなのですが、必要な仮定を満たしていない等の状況によって第一種の過誤確率は設定した有意水準より大きくなったり小さくなったりします。

*3:例えば2グループの平均に関するt検定であれば、2グループの平均の差をその標準誤差で割ったものが $T$ に相当します。

*4:t検定の場合は、t分布です。

*5:ここでは一般的に使用されることの多い両側検定を考えています。片側検定の場合は $|T|$ の絶対値を取らず、しかるべき方向の不等号を設定すれば同様に有意水準が定義されます。