OKIYUKI99 Blog

データ分析や日常に関するブログ

逐次的に検定を行っていると第一種過誤が増幅してしまう例(KDD 2017 : Peeking at A/B Tests より)

最近KDD 2017のPeeking at A/B Testsの論文を読んでいました。

www.kdd.org

その論文の問題設定に関わるこの図が少し気になっていました。

f:id:gingi99:20190703211931p:plain

これは何かというと、

このシミュレーションでは、サンプルサイズ10,000のA/B Test(仮にAが5,000、Bが5,000の対象が集まるまで続ける)を行うという設定で、データはAとBともに平均0、標準偏差1の正規分布からサンプリングしたデータを使います。 これは全く同じ分布から生成されるデータですので、平均が同じである帰無仮説を棄却してしまうことは、False Positive = 第一種過誤が生じていることを意味します。 そのような設定のもと、対象が集まる度に逐次的に検定を行い、X個のサンプルサイズ(横軸)が集まるまでに設定した有意水準(0.05とか)で有意差が 一度でも 検出される(= P値が有意水準を下回った = 第一種過誤が生じてしまった)確率をこの図の縦軸Yが表現しています。

一例を述べると、有意水準0.05の場合(緑の線)、サンプルが集まる度に検定をしていると、10,000サンプルサイズ集まるころには、約60%の確率で一度はFalse Positiveが発生してしまうことを意味しています。

有意水準0.05なので、そもそも5%の確率でFalse Positiveが発生するのですが、逐次的に検定を行っていると、それが何倍にも増幅してしまうことを示すとても良い一例だと思いました。

論文ではこのシミュレーションの細かい設定までは書いていなかったのですが、自分の理解の確認のため、Rで再現してみました。

パッケージの読み込み等は省略し、以下のシミュレーションを回してみます(結構時間がかかります…)。

sim_func <- function(simulations, sample_size, times, alpha){
  obs <- sapply(1:simulations, function(j){
    data_a <- c()
    data_b <- c()
    n <- sample_size / times
    for(i in 1:(times+1)){
      a <- rnorm(n / 2, mean = 0, sd = 1)
      b <- rnorm(n / 2, mean = 0, sd = 1)
      data_a <- c(data_a, a)
      data_b <- c(data_b, b)
      pval <- t.test(data_a, data_b, alternative = "two.sided", paired = F, var.equal = F)$p.value
      if(pval < alpha){
        break
      }
    }
    return(length(data_a) * 2)
  })
  
  data.frame(alpha = as.character(alpha), obs = obs, stringsAsFactors = F) %>%
    dplyr::group_by(alpha, obs) %>%
    dplyr::count() %>%
    ungroup() %>%
    mutate(false_positive_prob = cumsum(n) / simulations) -> df
  df
}

simulations <- 10000
sample_size <- 10000
times <- 1000
df_010 <- sim_func(simulations, sample_size, times, alpha = 0.1)
df_005 <- sim_func(simulations, sample_size, times, alpha = 0.05)
df_001 <- sim_func(simulations, sample_size, times, alpha = 0.01)

df <- bind_rows(df_010, df_005, df_001) 

1回のシミュレーションの中で、10個のサンプルが集まる(すでに集まってるサンプルに新しいサンプルが追加されるイメージ)度にt検定を行い、設定した有意水準のもと有意差が検出されるかをチェックします。 検出されたときのサンプルサイズを返す(スクリプトの中の変数obs)ことで、10,000個のサンプルが集まるまでに初めて検出されてしまうのはどれくらいサンプルが集まったときなのか?についてのデータを集めます。 これを3つの有意水準でそれぞれ10,000回ずつシミュレーションします。論文と同じ図を描くために累積数を計算して、可視化します。可視化のコードは以下です。

df %>%
  filter(obs <= sample_size) %>%
  mutate(alpha = format(as.numeric(alpha), nsmall = 2)) %>%
  ggplot(aes(x = obs, y = false_positive_prob, group = alpha, color = alpha, linetype = alpha)) + 
    geom_line(size = 1) + 
    scale_y_continuous(breaks = seq(0.0, 0.8, by=0.2), limits = c(0,0.8)) +
    theme_bw() + 
    labs(x = "# of observations", y = "False positive probability")

結果は以下のようになり、だいたい同じような図が描けました。

f:id:gingi99:20190703215943p:plain

論文よりも、False positiveの確率が低くなっているのは、私のシミュレーションでは計算時間を減らす&簡略化するために10個のサンプルが集まる度に検定を行っているためだと思います。 これを1個ずつにすれば、検出される確率は大きくなるため、論文ではそのようにしているのかと思います。

手を動かすことで、論文の問題設定が少しクリアに理解できた気がします(先は長い)