📊

比率の差の検定におけるサンプルサイズの計算

2021/05/29に公開

概要

A/B テストなど、比率の差を検定する場合に、検定を行うのに十分なサンプルサイズを算出する以下式については、ネットを検索すると日本語のページもいくつか見つかります[1] (各変数の意味については後述します)。

n = \frac{2p (1 - p) (z_{\alpha / 2} + z_{\beta})^2}{(p_1 - p_2)^2}

ですが、この式の導出方法についてはネットを検索してもなかなか見つけることができなかったため、導出方法についてまとめ、また実際に Python で計算してみました。

動機

A/B テストでは、二つの施策 A と B のいずれかをユーザに提示し、A と B のどちらでよりユーザが反応するか検証します[2]。両施策でどちらがより効果があったかは、両施策における反応率の差の検定で確認することができます。

A/B テストを実施するときに、検定において十分な検出力 (施策 A と B の反応率に差があったときに、それを検知できる確率) を保証するためには、全部で何人のユーザに施策を提示すればよいかを知る必要があります。

これは比率の差の検定を実施する場合の、所望の検出力を達成するために必要なサンプルサイズを求めることと同義となります。

検定の設定

今二つの二項母集団があり、各々は母比率 p_1p_2 を持つとします。又、これらの標本比率を \hat{p_1} = X_1 / n_1\hat{p}_2 = X_2 / n_2 とします。両者の比率の差 d = \hat{p_1} - \hat{p_2} について、その期待値と分散は、

\begin{aligned} E[d] &= E[\hat{p_1} - \hat{p_2}] \\ &= E[\hat{p_1}] - E[\hat{p_2}] \\ &= p_1 - p_2 \\ \\ V[d] &= V[\hat{p_1} - \hat{p_2}] \\ &= V[\hat{p_1}] + V[\hat{p_2}] \\ &= \frac{p_1 (1 - p_1)}{n_1} + \frac{p_2 (1 - p_2)}{n_2} \end{aligned}

であり、n_1n_2 が十分に大きければ d は正規分布に従います。

比率の差の検定では帰無仮説を H_0: p_1 = p_2 (両群の比率に差はなし)、対立仮説を H_1: p_1 \neq p_2 (両群の比率に差はある) として両側検定を行います。この際の帰無仮説の元での検定統計量は、

z_0 = \frac{\hat{p_1} - \hat{p_2}}{\sqrt{\frac{\hat{p_1}(1 - \hat{p_1})}{n_1} + \frac{\hat{p_2}(1 - \hat{p_2})}{n_1}}}

で与えられ、これは標準正規分布します。ただし今回は帰無仮説を H_0: p_1 = p_2 としているため、代わりに標本比率をプールした p を使った以下の式を用います [3]

\begin{aligned} p &= \frac{x_1 + x_2}{n_1 + n_2} \\ \\ z_0 &= \frac{\hat{p_1} - \hat{p_2}}{\sqrt{p (1 - p)(\frac{1}{n_1} + \frac{1}{n_2})}} \end{aligned}

ここで A/B テストの場合、両群のサンプルサイズが等しい (n = n_1 = n_2) と仮定できるので、前出の式は以下となります、

\begin{aligned} z_0 &= \frac{\hat{p_1} - \hat{p_2}}{\sqrt{p (1 - p)\frac{2}{n}}} \\ &= \frac{\sqrt{n}(\hat{p_1} - \hat{p_2})}{\sqrt{2p (1 - p)}} \\ \\ p &= \frac{x_1 + x_2}{2n} \\ &= \frac{n\hat{p_1} + n\hat{p_2}}{2n} \\ &= \frac{\hat{p_1} + \hat{p_2}}{2} \end{aligned}

サンプルサイズの算出

この検定では有意水準 \alpha のとき |z_0| > z_{\alpha / 2} で帰無仮説 H_0 は棄却されます。

また、対立仮説 H_1 が正しいときこの検定の検出力は以下で求めることができます (\Phi は累積分布関数)[4]

\Phi(\frac{\sqrt{n}(p_1 - p_2)}{\sqrt{2p(1 - p)}} - z_{\alpha / 2})
p = \frac{p_1 + p_2}{2}

従って検出力 1 - \beta を達成するのに必要となるサンプルサイズは以下式を n について解くことで近似的に求めることができます。

\frac{\sqrt{n}(p_1 - p_2)}{\sqrt{2p(1 - p)}} - z_{\alpha / 2} = z_{\beta}
n = \frac{2p (1 - p) (z_{\alpha / 2} + z_{\beta})^2}{(p_1 - p_2)^2}

サンプルサイズを求める関数を Python で実装してみます[5]

import math
from scipy.stats import norm
    
def sample_power_probtest(
       p1: float,
       p2: float,
       power: float = 0.8,
       sig: float = 0.05
) -> int:
   z_half_alpha = norm.isf([sig / 2])
   z_beta = -1 * norm.isf([power])
   p = (p1 + p2) / 2
   return math.ceil(
       (2 * p * (1 - p) * ((z_half_alpha + z_beta) ** 2))
       / ((p1 - p2) ** 2)
   )

ここで scipy.stats.norm.isf は生存関数の逆関数 (Inverse survival function) です。生存関数は P(X > x) = 1 - F(x) を表わすので、norm.isf([sig / 2]) により P(X > z_{\alpha / 2}) = \alpha / 2 を満たす点 z_{\alpha / 2} が得られます。

同様に nprm.isf([power]) は、P(X > z_{1 - \beta}) となる点 z_{1 - \beta} が得られるので、この符号を逆転して z_{\beta} を得ています。

A/B テストを実施する場合に、A 群の比率が 0.1、B 群の比率が 0.11 と想定されるとき、有意水準 \alpha = 0.05、検出力を 1 - \beta = 0.8 としたい場合に必要なサンプルサイズを、この関数を使って求めてみます。

sample_power_probtest(0.1, 0.11)  # >>>> 14752

よって A/B テストを仮定した場合、14752 \times 2 = 29504 人に無作為に施策 A と B のいずれかを提示して、その結果に対し検定を行えばよいことになります。

仮に B 群の比率が 0.12 と想定される場合は、

sample_power_probtest(0.1, 0.12)  # >>>> 3843 

3843 \times 2 = 7686 人となり、想定される両比率の差が大きければ大きい程必要となるサンプルサイズは小さくなります。

statsmodels.stats.power での計算

statsmodels.stats.power には同等の計算をする関数があります。

import math
from statsmodels.stats.power import normal_sample_size_one_tail
    
p1 = 0.1
p2 = 0.11
p = (p1 + p2) / 2
std = math.sqrt(2 * p * (1 - p))
    
normal_sample_size_one_tail(
   diff=0.1 - 0.11,
   power=0.8,
   alpha=0.05 / 2,
   std_null=std,
   std_alternative=std,
)  # >>>> 14751.969460709131
脚注
  1. 統計 WEB - 幾つデータが必要か?―比率の差の検定 ↩︎

  2. Wikipedia - A/B testing ↩︎

  3. 生物科学研究所 井口研究室 - 比率の差 Z 検定の注意点:統合比率を使う理由 ↩︎

  4. Wang, Hansheng, and Shein‐Chung Chow. “Sample size calculation for comparing proportions.” Encyclopedia of Statistical Sciences (2004): 1-14. ↩︎

  5. stack overflow - Is there a python (scipy) function to determine parameters needed to obtain a target power? ↩︎

Discussion