【Python】層化A/Bテストを実行し、その信頼性をチェックする
この記事は、CA Tech Lounge Advent Calendar 2024 7日目の記事となっています。
はじめに
共変量の値に基づいて層を作り、層ごとに割り当てを行う層化A/Bテスト(層化RCT)について以前書いた記事を、多くの方に読んでいただきました。
この記事では共変量である連続値の年齢や会員歴を層化するにあたって、母集団における因果関係と同じ基準で層を作ったため、かなり理想的な推定結果が算出されました。
そこで今回は、層化基準と異なる基準で母集団の説明変数-被説明変数間の因果関係が構成されていると仮定し、層化A/Bテストを試行しました。加えて、A/Bテストの信頼性をチェックするために、各層でA/Aテストのシミュレーションを実施するPythonコードを作成しました。
なお、本記事は安井、伊藤、金子（2024）「Pythonで学ぶ効果検証入門」を参考とさせていただきました。
層化A/Bテストとは
通常のA/Bテストでは、Treatment GroupとControl Groupがランダムに割り振られることで処置を行ったことによる平均効果(平均処置効果：ATE)の推定量を算出することができます。一方で、問題となるのはTreatment GroupとControl Groupに属するサンプルに偏りが生じうるという点です。これを解決し推定量の分散を小さくするために、層化A/Bテストでは共変量に基づいた層化を行い、層ごとに独立したRCTを実行して処置効果を推定し、それらを適切に重み付けして全体の平均処置効果（ATE）を算出します。
層化RCTにおいて、ATEを算出する方法を二通り紹介します。
その1 ~層ごとの平均の差をとる方法~
個体の共変量に基づき、層を1,2,…,Sに分けます。各層で同一の処置割り当て確率を置き、層sにおける平均処置効果を
-
\overline{Y}_{t,s} -
\overline{Y}_{c,s}
このとき、全体の平均処置効果
-
n_s -
N
この
その2 ~回帰~
Yを層ダミーと、層ダミーと処置ダミーの交差項とに回帰することで得られる交差項の係数が、各層の平均処置効果の推定量、すなわち「その1」の
回帰式は以下の通り
-
S_{is} -
Z_i
上記の回帰式で得られた
A/Aテストのシミュレーションとは
今回は、A/Bテストの設計の妥当性を判断するために、A/Aテストのシミュレーションを以下のような手順で実行します。
- Treatment Group/Control Groupの割り当てを行う
- 推定量を算出し、p値を算出する
- 1,2を300回繰り返す。このとき、毎回シード値を変更する
- 3を全ての層で実施する
このとき、割り当てごとのGroup間で差がないと想定されるため、P値は偏った値をとらないと考えられます。一方で、P値が偏る場合は、設計に問題があると考えられます。
今回は共変量に基づいて層化を行ったうえで各層で割り当てを実施しているため、従来のA/Bテストと比較するとA/Aテストが失敗する、即ちP値が0付近に偏ることは起こりにくいと考えられます。
詳細に関しては安井、伊藤、金子（2024）「Pythonで学ぶ効果検証入門」P85-110をご参照ください。
データ設定
実行場面の設定
ECサイト会員に向けて購入プロセスの簡素化施策を実施したことが、実際の購買につながったかどうかの効果検証を行うケースを想定します。このECサイトは、カゴ落ち（購入放棄）に悩んでいると仮定します。そのためサイト運営会社は、購入までのプロセスの一部を簡素化することで、かご落ちを防ぎ購買の完結を促したいと考えています。今回は、この効果を検証するためにA/Bテストを実行する場面を想定します。
結果変数は、ユーザーごとの、購入かごに入れた全アイテムに対する実際に購入されたアイテムの割合とします。
共変量は、ユーザーの年齢と会員歴(年)とします。また、上記の結果変数は、共変量の一部とTreatmentに依存すると考えます。
人工データ作成
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import stats
import statsmodels.formula.api as smf
from scipy.stats import mannwhitneyu
def make_covariate(n_samples, seed):
np.random.seed(seed)
#共変量1:年齢を正規分布（平均35歳、標準偏差10歳）から生成
ages_continuous = np.random.normal(35, 10, size=n_samples).astype(int)
#共変量2:会員歴(簡素化のため、1~9の値の出現回数はある程度一致する仕様とする)
membership_years = np.random.randint(1, 10, size=n_samples)
df = pd.DataFrame({
'age': ages_continuous,
'membership_years': membership_years
})
#ユーザーを21歳以上50歳未満に絞る
df = df[(df['age'] >= 21) & (df['age'] < 51)]
#年層を3年区切りで作成
df['year_group'] = pd.cut(df['membership_years'], bins=np.arange(1, 11, 3), labels=[f'{i}-{i+2}' for i in range(1, 8, 3)], right=False)
#年齢層を10年区切りで作成
df['age_group'] = pd.cut(df['age'], bins=[20, 30, 40, 50], labels=[f'{i}0s' for i in range(2, 5)], right=False)
#層を作成
df['Group1'] = df.groupby(['age_group', 'year_group']).ngroup()
#層ごとに割り当てを実施
ratio = 0.3 #処置割当比率
rng = np.random.default_rng(seed)
df['Treatment'] = 0
for group, group_df in df.groupby('Group1'):
num_samples = round(len(group_df) * ratio)
treat_indices = rng.choice(group_df.index, num_samples, replace=False)
df.loc[treat_indices, 'Treatment'] = 1
#ノイズとYスコア(目的変数、購入かごに入れた全アイテムに対する実際に購入されたアイテムの割合)の生成
df['Y_noise'] = np.random.normal(0, 0.3, len(df))
def calculate_age_score(age):
if 20 <= age <= 25:
return 0.2
elif 25 < age <= 43:
return (1 - 0.02 * (age - 25)) * 0.2
else:
return 0 # 条件外の年齢はスコアに影響しない
#年齢スコア列の作成
df['age_score'] = df['age'].apply(calculate_age_score)
# membership_years に基づくスコア増加
df['membership_score'] = df['membership_years'] * 0.02
# Y_score(目的変数)の算出
df['Y_score'] = (
df['age_score'] + # 年齢スコア
df['membership_score'] + # 会員年数スコア
0.3 * df['Treatment'] + # 処置の影響
df['Y_noise'] # ノイズ
)
return df
df = make_covariate(n_samples=1000, seed=4)
このとき、各層のサイズは以下のようになります。
また、被説明変数は、以下のように決定されると仮定します。
- ageが20歳以上25歳以下：0.2上昇
- ageが26歳以上43歳以下：(0.2-0.02*(age-25))上昇
- age = 30のとき、0.2 - 0.02*(30 - 25) = 0.1分上昇
- membership_years(会員歴)が1年増加すると0.02上昇
- Treatmentに割り当てられると0.3上昇
- ノイズ(0から0.3までの間の値がランダムに与えられる)
A/Aテストのシミュレーション実行
def aa_test_simulation(df, n_simulations=300, treatment_prob=0.3, random_seed=4):
np.random.seed(random_seed)
results = []
#各層ごとに処理
for group in df['Group1'].unique():
group_data = df[df['Group1'] == group].copy()
group_size = len(group_data)
#シミュレーションを実行
for sim in range(n_simulations):
#各シミュレーションのシードを設定
sim_seed = random_seed + sim
rng = np.random.default_rng(sim_seed)
#処置群と対照群をランダムに割り当て
shuffled_treatment = rng.choice(
[0, 1], size=group_size, replace=True, p=[1 - treatment_prob, treatment_prob]
)
group_data['aa_Treatment'] = shuffled_treatment
#Mann-WhitneyのU検定を実行
treatment_scores = group_data[group_data['aa_Treatment'] == 1]['Y_score']
control_scores = group_data[group_data['aa_Treatment'] == 0]['Y_score']
#検定実施
if len(treatment_scores) > 0 and len(control_scores) > 0:
result = mannwhitneyu(treatment_scores, control_scores, alternative="two-sided")
results.append({
'Group1': group,
'Simulation': sim,
'p_value': result.pvalue
})
results_df = pd.DataFrame(results)
return results_df
n_simulations = 300
treatment_prob = 0.3
random_seed = 42
simulation_results = aa_test_simulation(df, n_simulations=n_simulations, treatment_prob=treatment_prob, random_seed=random_seed)
#各層ごとのp値の分布を確認
plt.figure(figsize=(12, 6))
for group in simulation_results['Group1'].unique():
p_values = simulation_results[simulation_results['Group1'] == group]['p_value']
plt.hist(p_values, bins=20, alpha=0.5, label=f'Group {group}')
plt.title('P-value Distribution for Each Group')
plt.xlabel('P-value')
plt.ylabel('Frequency')
plt.legend()
plt.show()
300回のシミュレーションを各層で回した結果、各試行のP値の出現回数は以下のようになりました。
層ごとのシミュレーション結果は以下の通りです。
unique_groups = simulation_results['Group1'].unique()
num_groups = len(unique_groups)
fig, axes = plt.subplots(3, 3, figsize=(20, 10)) #3行3列
axes = axes.flatten()
#層ごとにヒストグラムを作成
for i, group in enumerate(unique_groups):
p_values = simulation_results[simulation_results['Group1'] == group]['p_value']
axes[i].hist(p_values, bins=20, alpha=0.75, color='blue')
axes[i].set_title(f'Group {group}')
axes[i].set_xlabel('P-value')
axes[i].set_ylabel('Frequency')
plt.tight_layout()
plt.show()
この結果から、P値に偏りがなく、一様分布に近い分布となっていることから、A/Bテストが成功していることが分かります。
処置の実行結果
上で示したATEの算出方法のうち、今回は回帰によって推定量を算出する方法を試行しました。
コードは以下の通りです。
#回帰分析
model = smf.ols(formula='Y_score ~ C(Group1) + C(Group1):C(Treatment)', data=df).fit()
print(model.summary())
各層の推定量を重みづけし、全体のATEを算出します。
#1.交差項の係数を取得
coefficients = [model.params[f'C(Group1)[{i}.0]:C(Treatment)[T.1]'] for i in range(9)]
#2.Groupごとのサンプル数をカウントして全体比を計算
group_counts = df['Group1'].value_counts().sort_index()[:9] # グループ0〜8の人数
group_ratios = group_counts / group_counts.sum() # 全体に対する割合
#3.交差項の係数に人数の割合を重みづけして総和を計算
weighted_sum = np.dot(coefficients, group_ratios)
print(f'重み付けした総和: {weighted_sum}')
結果として、「重み付けした総和: 0.31799145189195127」という、0.3に近い値を算出することができました。
おわりに
層化A/Bテストを紹介しました。推定量の分散を小さくすることができることに加え、共変量を踏まえた割り当てが可能になることからA/Bテストが失敗する可能性を小さくすることができるのではと考えられます。また、今回のケースでは各層のA/BテストのATEの推定量が1%有意水準で有意であったため、層に基づいて条件付き平均処置効果を得ることができるとも考えられます。事前の設計コストはかかるものの、A/Bテストの成功確率を上げることができるという点で実施の検討の余地は大いにあるのではないかと思います。
Discussion