qiitaにも同じ記事があるのですが、試しにzennにも移植してみました。
はじめに
A/Bテストを行う際、施策を打ったことによる効果を推定するだけでなく、推定値がどれだけ信頼出来るものなのかを計算することはとても重要です。
そこで、効果の推定誤差と信頼区間を導出します。
例) ATEの推定値と上下95%信頼区間
問題設定
あるキャンペーンを打つことで、ECサイトの商品購入率がどれだけ変化するかを推定します。
キャンペーンの対象となる介入群をT, キャンペーン対象外となる統制群をC, サンプルサイズをそれぞれn, m、母平均をp_{T}, p_{C}とすると、介入群と統制群におけるユーザーの行動X_{T}, X_{C}はベルヌーイ分布に従います。
ユーザーが商品を購入したなら1, 購入しなかったら0となります。
なお、ユーザーはA/Bを独立して割り付けられ、ユーザーの行動は独立とします。
X_{T1}, X_{T2}, \ldots X_{Tn} \sim Bern(p_{T}) \\
X_{C1}, X_{C2}, \ldots X_{Cm} \sim Bern(p_{C})
効果の分布の推定
キャンペーンを打つことが購入率に与える施策の効果をATEと表記します。
ATEの分布、平均、分散を導出します。
中心極限定理
平均、分散パラメータに\mu, \sigma^2を持つ確率分布Fから独立にサンプリングされた標本があるとすると、
X_1, X_2, \ldots, iid \sim F(\mu, \sigma^2)
中心極限定理より、標本平均\bar{X}について次の分布収束が成り立ちます。
\lim_{n\rightarrow\infty}P \left( \frac{\sqrt n (\bar{X}-\mu)}{\sigma} \leq x \right) = \Phi(x)
つまり、平均と分散パラメータを持つ任意の確率分布から独立に発生した標本の標本平均が従う分布は、漸近的に正規分布に近づいていきます。
ベルヌーイ分布の母平均をpとすれば母分散はp(1-p)で表せるので、サンプルサイズが大きい場合、介入群と統制群の標本平均\bar{X_T}, \bar{X_C}は以下の正規分布に従います。
\bar{X_T} = \frac{1}{n}\sum^n X_{Ti} \\
\bar{X_C} = \frac{1}{m}\sum^m X_{Ci} \\
\bar{X_T} \sim N(p_T, \frac{1}{n}p_{T}(1-p_{T})) \\
\bar{X_C} \sim N(p_C, \frac{1}{m}p_{C}(1-p_{C})) \\
効果の定義
効果ATEは、2つの潜在的な結果変数の期待値の差と定義します。
ATE = E[X_T - X_C] = E[X_T] - E[X_C] = p_T - p_C
「ユーザーの母集団全てが介入群に割り当てられた際の結果」と「ユーザーの母集団全てが統制群に割り当てられた際の結果」の差の平均です。
実際は、ユーザーは介入群か統制群どちらかにしか割り当てられません。
例えば、あるユーザーがキャンペーン対象者になったとき、その後商品を購入したかは観測できます。
しかし、そのユーザーは既にキャンペーン対象に割り当ててしまっているので、そのユーザーがキャンペーン対象にならなかったとき商品を購入したかどうかはわかりません。
なので、ATEは直接データから値を計算することはできません。
そこで、通常A/Bテストは割り当てを完全にランダムに行います。
割り当てがランダムに行われていれば、観測された各群の平均値の差を用いることで効果をバイアスなく(不偏)推定できます。
ATEの推定量を\hat{ATE}とおきます。
\hat{ATE} = \frac{1}{n}\sum^nX_{T,i} + \frac{1}{m}\sum^mX_{C,i} = \bar{X_T} - \bar{X_C} \\
E[\hat{ATE}] = E[\bar{X_T}] - E[\bar{X_C}] = p_T - p_C = ATE \\
効果の推定量が従う確率分布の導出
効果の推定量は介入群の購入率 - 統制群の購入率で計算します。
介入群の購入率と統制群の購入率はそれぞれ正規分布に従うので、それらを元に計算した統計量も正規分布に従います。
\hat{ATE} = \bar{X_T} - \bar{X_C} \sim N(p_T - p_C, \frac{1}{n}p_{T}(1-p_{T})+\frac{1}{m}p_{C}(1-p_{C}))
施策を打ったことによる効果の統計量\hat{ATE}が従う確率分布の母数は
- 平均p_T - p_C
- 分散\frac{1}{n}p_{T}(1-p_{T})+\frac{1}{m}p_{C}(1-p_{C})
- 標準偏差\sqrt{\frac{1}{n}p_{T}(1-p_{T})+\frac{1}{m}p_{C}(1-p_{C})}
となります。
推測統計における単語の整理
ここまで様々な記号を前振りなく使ってきました。
一度様々な用語を整理します。
統計調査と推測統計
統計調査は、母集団の特性を調べるためにデータを抽出し、そのデータを加工して母集団に関する有益な情報を導いて判断や決定の材料にするために行います。
母集団の全員に調査するのが難しい場合、母集団から標本を抽出して母集団に対する推測を行います。
推測統計では、母集団に確率モデルを想定し、その確率分布に従う確率変数の実現値としてデータを捉えます。
ある確率分布に従う標本X_1, X_2, \ldots, X_nは確率変数で、その実現値x_1, x_2, \ldots, x_nがデータです。nを標本の大きさ(sample size)と言います。
母数と統計量・推定量
母集団の平均\muや分散\sigma^2など、母集団の特性を表すものを母数と言います。
標本X_1, X_2, \ldots, X_nに基づいた関数で母数を含んでいないものを統計量といい、その確率分布を標本分布と言います。
例えば、標本X_1や標本平均\bar{X_T}などが統計量です。また、統計量を元に計算する\hat{ATE}も統計量となります。
推測統計においては、統計量を計算することで母数を推測します。
標本から得た統計量をもとに、母集団のパラメータを推定するとき、標本の統計量のことを、推定量といいます。
推定量と推定値について
推定量は確率変数です。確定的に何かしらの値をとるわけではありません。
実際にデータを観測して得られた数字をつかって計算したものが、推定値です。
推定量と推定値の違い
効果の確率分布の母数の推定
実際に観測されたデータを用いて母数を推定します。
母平均の推定
母集団の平均\muは母数ですが、母集団の全てのデータを元に直接計算するのは容易ではありません。
そこで、n個のデータを観測したのち、標本平均\bar{X}=\frac{1}{n}\sum^n X_iを計算します。
標本平均はその平均が母平均と一致します。E[\bar{X}]=\muです。
なので、標本平均を通じて母平均を推測することが可能になります。
データを(x_1, x_2, \ldots, x_n)とすると、推定値は\bar{x}=\sum^n x_iとなります。
母平均の不確実性
ただし、統計量(=推定量)は確率変数なので、不確実性がつきまといます。
10個のデータを元に計算した標本平均と10000個のデータを元にした標本平均では、その値に対する不確実性が大きく変わってきます。
そこで、標準誤差を計算して不確実性を明らかにします。
標準偏差と標準誤差の違い
確率変数Xがあるとき、平均\mu = E[X]、分散\sigma^2 = Var(X) = E[(X-\mu)^2]が定義出来ます。
分散のルートを取ったものを標準偏差\sigma = \sqrt{Var(X)}といいます。
統計量も確率変数なので、同様に標準偏差を定義出来ます。
統計量はその平均が母平均になることを通じて、母数の値を推定することによく使われます。
なので、統計量の標準偏差、すなわち真の値に対する推定値の不確実性は「標準誤差」と呼ばれます。
例えば、標本平均\bar{X}=\frac{1}{n}\sum_{i=1}^n X_iの標準偏差は\frac{\sigma}{\sqrt n}となりますが、\bar{X}は母平均\muを推定するための統計量であるという観点で、その不確実性を表す\frac{\sigma}{\sqrt n}は標準誤差です。
標準偏差の推定
推定量の期待値が母数と等しい場合、不偏推定量と言います。
例えば、標本平均は母平均の不偏推定量です。
母分散\sigma^2の不偏推定量V^2は以下のようになります。
V^2 = \frac{1}{n-1} \sum^n (X_i - \bar{X})^2 \\
E[V^2] = \sigma^2 \\
V^2は不偏分散と呼ばれます。
標本平均は個々のデータの合計をnで割りましたが、不偏分散はn-1で割るのがポイントです。
Pandasの関数df.var()やdf.std()ではデフォルトで不偏分散が計算されるようになっています。
上の定義に従って不偏分散を計算することも出来ますが、ベルヌーイ分布の分散がp(1-p)で計算できることを利用して、標本平均\bar{X}や標本分散\bar{X}(1-\bar{X})を使って不偏分散V^2を計算することも出来ます。
\begin{align}
E[\bar{X}(1-\bar{X})] &= E[\bar{X} - \bar{X}^2] \\
&= E[\bar{X}] - E[\bar{X}^2] \\
&= E[\bar{X}] - (E[\bar{X}]^2 + Var(\bar{X})) \\
&= p - p^2 - \frac{1}{n^2} \sum^n Var(X_i) \\
&= p(1 - p) - \frac{1}{n} p(1 - p) \\
&= \frac{n-1}{n} p (1 - p) \\
E \left[ \frac{n}{n-1}\bar{X}(1-\bar{X}) \right] &= \sigma^2 \\
\end{align} \\
従って、Xの不偏分散はV^2 = \frac{n}{n-1}\bar{X}(1-\bar{X})となります。
介入群と統制群の不偏分散は以下のようになります。
V_T^2 = \frac{n}{n-1} \bar{X_T} (1 - \bar{X_T}) \\
V_C^2 = \frac{m}{m-1} \bar{X_C} (1 - \bar{X_C}) \\
Xの母平均\muの推定値として標本平均\bar{x}を計算します。
Xの母分散\sigma^2=p(1-p)の推定値v^2は以下のようになります。
\bar{x} = \frac{1}{n} \sum^n x_i \\
v^2 = \frac{n}{n-1} \bar{x} (1 - \bar{x}) \\
標準偏差の推定値はvです。
標準誤差の推定
標本平均\bar{X}の分散(標準誤差の2乗)は以下のようになります。
介入群と統制群の標本平均の分散は以下のようになります。
SE_T^2 = \frac{V_T^2}{n} \\
SE_C^2 = \frac{V_C^2}{m} \\
従って、効果の分散は以下のようになります。
V_{ATE}^2 = SE_T^2 + SE_C^2 = \frac{V_T^2}{n} + \frac{V_C^2}{m}
プールした分散
介入群と統制群が従う確率分布の分散がV_{TC}^2で等しいと仮定すると、\hat{ATE}の分散を小さく(=推定精度を高める)ことができます。
等分散性が成り立つか予めチェックが必要です。
V_T^2 = \frac{n}{n-1}\bar{X_T}(1-\bar{X_T}) \\
V_C^2 = \frac{m}{m-1}\bar{X_C}(1-\bar{X_C}) \\
V_{TC}^2 = \frac{(n-1)V_T^2 + (m-1)V_C^2}{(n-1)+(m-1)} \\
= \frac{(n-1)V_T^2 + (m-1)V_C^2}{n+m-2} \\
V_{ATE,pool}^2 = \frac{V_{TC}^2}{n} + \frac{V_{TC}^2}{m} \\
効果の確率分布の母数の推定まとめ
今までに出てきた式をまとめます。
介入群と処置群の確率分布について
X_{T1}, X_{T2}, \ldots X_{Tn} \sim Bern(p_{T}) \\
X_{C1}, X_{C2}, \ldots X_{Cm} \sim Bern(p_{C}) \\
\bar{X_T} \sim N(p_T, \frac{1}{n}p_{T}(1-p_{T})) \\
\bar{X_C} \sim N(p_C, \frac{1}{m}p_{C}(1-p_{C})) \\
\hat{ATE} \sim N(p_T - p_C, \frac{1}{n}p_{T}(1-p_{T})+\frac{1}{m}p_{C}(1-p_{C}))
効果の母数の推定量
平均パラメータの推定量
\bar{X_T} = \frac{1}{n}\sum^n X_{Ti} \\
\bar{X_C} = \frac{1}{m}\sum^n X_{Ci} \\
\hat{ATE} = \bar{X_T} - \bar{X_C}
分散パラメータの推定量
V_T^2 = \frac{n}{n-1}\bar{X_T}(1-\bar{X_T}) \\
V_C^2 = \frac{m}{m-1}\bar{X_C}(1-\bar{X_C}) \\
V_{ATE}^2 = \frac{V_T^2}{n} + \frac{V_C^2}{m} \\
V_{TC}^2 = \frac{(n-1)V_T^2 + (m-1)V_C^2}{n+m-2} \\
V_{ATE,pool}^2 = \frac{V_{TC}^2}{n} + \frac{V_{TC}^2}{m} \\
効果の母数の推定値について
観測されたデータを用いて、統計量の実現値である推定値を計算します。
ほぼ推定量と一緒なので一部省略します。
\hat{ate} = \bar{x_T} - \bar{x_C} \\
v_{ATE}^2 = \frac{v_T^2}{n} + \frac{v_C^2}{m} \\
v_{ATE,pool}^2 = \frac{v_{TC}^2}{n} + \frac{v_{TC}^2}{m} \\
効果の推定値の信頼区間の計算
最後に、効果の推定値の不確実性を測るため、\hat{ATE}の信頼区間を計算します。
仮説検定と信頼区間
仮説検定では、帰無仮説と対立仮説という相対する2つの仮説を立てます。
通常なら母平均\muはある値\mu_0だ、というのが帰無仮説です。
それに対し、母平均の値は\mu_0ではないのではないか、という仮説を対立仮説と言います。
確率変数Xの実現値xを観測することで、帰無仮説が正しい場合にその実現値が観測される確率がどれだけ珍しいのかを計算し、計算結果に基づいて仮説を受容したり棄却します。
ある確率変数が正規分布に従うとします。
分散\sigma_0^2は既知とします。
X \sim N(\mu, \sigma_0^2)
帰無仮説と対立仮説を以下のように設定します。
H_0: \mu = \mu_0 \\
H_1: \mu \neq \mu_0 \\
帰無仮説の棄却域、受容域は以下のようになります。
R = \left\{ x \in \mathbb{R} \middle| \frac{|x - \mu_0|}{\sigma_0} > z_{\alpha/2}\right\} \\
A = \left\{ x \in \mathbb{R} \middle| \frac{|x - \mu_0|}{\sigma_0} \leq z_{\alpha/2}\right\}
z_{\alpha/2}とは標準正規分布のz値です。
標準正規分布\Phi = N(0,1)からz_{\alpha}という値以上が取り出される確率が\alphaとなります。
標準正規分布表の使い方
ここで、検定方式をXから\mu_0に反転させると
A = \left\{ \mu_0 \in \mathbb{R} \middle| \frac{|x - \mu_0|}{\sigma_0} \leq z_{\alpha/2}\right\} \\
= \left\{ \mu_0 \in \mathbb{R} \middle| |x - \mu_0| \leq \sigma_0 z_{\alpha/2}\right\} \\
= \left\{ \mu_0 \in \mathbb{R} \middle| x - \sigma_0 z_{\alpha/2} \leq \mu_0 \leq x + \sigma_0 z_{\alpha/2} \right\}
となります。
帰無仮説が正しい場合、標本Xの実現値がxとして観測される確率が1 - \alpha%以内に収まる区間は[x - \sigma_0 z_{\alpha/2}, x + \sigma_0 z_{\alpha/2}]となります。
これを、信頼係数1 - \alphaの信頼区間と言います。
効果の推定量の信頼区間の計算
信頼区間の計算方法がわかったので、効果の信頼区間を実際に計算します。
効果が従う確率分布の母数
施策を打ったことによる効果の推定量\hat{ATE}が従う確率分布の母数は
\mu = ATE = p_T - p_C \\
\sigma = \sqrt{\frac{1}{n}p_{T}(1-p_{T})+\frac{1}{m}p_{C}(1-p_{C})} \\
でした。
データから母数の推定値を計算
平均の推定値を計算します。
標準偏差の真の値はわからないので、代わりに不偏分散のルートで代用します。
\hat{\mu} = \hat{ate} = \bar{x_T} - \bar{x_C} \\
\sigma_0 \approx v_{ATE} = \sqrt{\frac{1}{n-1}\bar{x_T}(1-\bar{x_T})+\frac{1}{m-1}\bar{x_C}(1-\bar{x_C})} \\
ATEが従う分布が正規分布の場合(※後述)、施策を打つことによる効果ATEの1 - \alpha信頼区間は
[\hat{\mu} - v_{ATE} \, z_{\alpha/2}, \hat{\mu} + v_{ATE} \, z_{\alpha/2}]
となります。
\alphaとz_{\alpha/2}の関係は
\alpha = 0.05 \leftrightarrow z_{\alpha/2} = 1.96 \\
\alpha = 0.1 \leftrightarrow z_{\alpha/2} = 1.65 \\
なので、例えば\alpha=0.05の場合の両側95%信頼区間は
[\hat{\mu} - 1.96v_{ATE}, \hat{\mu} + 1.96v_{ATE}]
です。
「標本平均 ± 2 × 標準誤差」が両側95%信頼区間の近似になります。
実用上はまずこれを使うところから始めるといいと思います。
コード
import numpy as np
import pandas as pd
import plotly.express as px
t_imp = 10000
c_imp = 50000
t_click = 100
c_click = 350
t_ctr = t_click / t_imp
c_ctr = c_click / c_imp
ate_mean = t_ctr - c_ctr
ate_var = (t_ctr * (1 - t_ctr) / (t_imp - 1)) + (c_ctr * (1 - c_ctr) / (c_imp - 1))
ate_se = np.sqrt(ate_var)
t_var = (t_ctr * (1 - t_ctr)) * (t_imp / (t_imp - 1))
c_var = (c_ctr * (1 - c_ctr)) * (c_imp / (c_imp - 1))
tc_var_pool = ((t_imp - 1) * t_var + (c_imp - 1) * c_var) / (t_imp + c_imp - 2)
ate_var_pool = tc_var_pool * (1 / t_imp + 1 / c_imp)
ate_se_pool = np.sqrt(ate_var_pool)
df_plot = pd.DataFrame({
"mean": [ate_mean, ate_mean],
"2se": [2 * ate_se, 2 * ate_se_pool],
"type": ["equal_var=False", "equal_var=True"]
})
px.bar(df_plot, x="type", y="mean", error_y="2se", height=400, width=400)
等分散性を仮定したほうが、推定誤差がわずかに小さいことがわかります。
感想など
色々調べたお陰で推論統計の概念を整理できました。
- ある事象の発生過程に確率分布を仮定し、その形状を決める母数を知りたい。
- 母数は観測不可なので、標本をもとに(不偏)推定量を導出する。
- 推定量は確率変数なので、確からしさを標準偏差や標準誤差として導出する。
- 実際に観測されたデータをもとに平均や標準誤差などの推定値を計算する。
用語がこんがらがってややこしい。。
参考文献
(久保川, 2017) 「現代数理統計学の基礎」
(星野, 2009) 「調査観察データの統計科学」
(安井, 2020) 「効果検証入門〜正しい比較のための因果推論/計量経済学の基礎」
(余談) 正規分布とt分布
信頼区間に正規分布を仮定しましたが、分散が未知の場合\hat{ATE}は正規分布ではなくt分布に従います。
今回の問題設定ではECサイトの商品購入率なので、十分なデータ数が集まり、t分布と正規分布の形状はほぼ一致します。そのため説明を省きました。
t分布を用いてより厳密な検定や信頼区間を推定したい場合は、z_{\alpha/2}をt_{\alpha/2}に置き換えてください。
実際に計算する際は、t値を手計算するのではなく、データから直接stats.ttest_indなどのパッケージを用いた方が楽です。
Discussion