R

Rのprop.testを読んで、jsで書いてみる

A/Bテストの評価で用いる二群の比率の差の検定のprop.testが実際にどのようなアルゴリズムで動いているか調べてみました。

結論

高校生で勉強する母比率の信頼区間推定にイェイツの修正を加えて実装されています。

p-Valueの計算

prop.testでは、A/Bテストで用いる2つの母集団の比較に加えて、1も2より大きいものもサポートされています。
ただ今回は簡略化のために、2母集団比較のみ考慮して計算しています。

そのままprop.testを実行すると以下のようになります。

> x=c(20, 44)
> n=c(40, 80)
> prop.test(x, n)

    2-sample test for equality of proportions with continuity correction

data:  x out of n
X-squared = 0.10463, df = 1, p-value = 0.7463
alternative hypothesis: two.sided
95 percent confidence interval:
 -0.2582061  0.1582061
sample estimates:
prop 1 prop 2 
  0.50   0.55 

では、実際にはどのような処理をしているかと言うと、以下のような処理をしています。

p <- sum(x)/sum(n)

x <- cbind(x, n - x)
E <- cbind(n * p, n * (1 - p))

ESTIMATE <- x/n

DELTA <- ESTIMATE[1L] - ESTIMATE[2L]
YATES = 0.5
YATES <- min(YATES, abs(DELTA)/sum(1/n))

STATISTIC <- sum((abs(x - E) - YATES)^2 / E)
PVAL <- pchisq(STATISTIC, 1, lower.tail = FALSE)

ざっくり言うと、実際のデータは、理論的に計算するよりも信頼区間をある程度大きめにとった方が正しく補正できると言うことで処理しています。

実際にどのくらい補正されているのか

YATES-correctionの補正はcorrectと言う引数を使い分けることで違いを確認できます。
以下の例では、p-valueは0.24と0.20と0.04ほど変化がありました。
仮に母集団の数で比較するとだいたい8%ほど補正することで多くのnが必要になるくらいの補正でした。

x = c(80, 50)
n = c(2000, 1000)

prop.test(x, n)
=> p-value = 0.2408

> YATES => 0.5
> prop.test(x, n, correct=FALSE)
=> p-value = 0.2048

xx=2.4;prop.test(x*xx, n*xx, correct=FALSE)
=> p-value = 0.04946
xx=2.6;prop.test(x*xx, n*xx)
=> p-value = 0.04705

BigQuery上で実行する方法

今回の実装をそのままjsで書き換え、BigQuery上で実行するには、以下のようになります。

CREATE TEMPORARY FUNCTION prop_test(x_1 FLOAT64, x_2 FLOAT64, n_1 FLOAT64, n_2 FLOAT64)
RETURNS FLOAT64
LANGUAGE js AS """

  p = (x_1 + x_2)/(n_1 + n_2)
  yates = 0.5
  diff = Math.pow((Math.abs(x_1 - n_1 * p) - yates), 2)
  statistic = diff/(n_1*p) + diff/(n_1*(1-p)) + diff/(n_2*p) + diff/(n_2*(1-p))

  return 1 - jStat.chisquare.cdf(statistic, 1);
"""
OPTIONS (
  library="gs://#{project_name}/jstat.js"
);

SELECT
  prop_test(20, 44, 40, 80)

=> 0.24078852577936793

参考URL

高等学校数学C/統計処理
https://ja.wikibooks.org/wiki/%E9%AB%98%E7%AD%89%E5%AD%A6%E6%A0%A1%E6%95%B0%E5%AD%A6C/%E7%B5%B1%E8%A8%88%E5%87%A6%E7%90%86
イェイツのカイ二乗検定
https://bellcurve.jp/statistics/course/9122.html