本教程探讨了协变量如何影响随机实验中的A/B测试精度。正确随机化的A/B测试通过比较治疗组和对照组的平均结果来计算提升。然而，除了治疗因素以外的其他特征对结果的影响决定了A/B测试的统计特性。例如，在测试提升计算中省略重要特征可能导致提升估计非常不精确，即使随着样本量增加，其收敛到真实值。
通过生成模拟数据并进行蒙特卡洛实验，您将了解RMSE、偏差和测试规模是什么，并理解A/B测试的性能。这种工作有助于了解数据生成过程（DGP）的特性如何影响A/B测试的性能，并帮助您在真实数据上运行A/B测试。首先，我们讨论一些估计器的基本统计特性。
估计器的统计特性
均方根误差（RMSE）
均方根误差（Root Mean Square Error，RMSE）：RMSE是一种经常用于衡量模型或估计器预测值与观察值之间差异的指标。它是预测值与实际观测值之间差异的平方的平均值的平方根。RMSE的公式如下：
RMSE = sqrt[(1/n) * Σ(实际值 – 预测值)²]
RMSE对大误差给予相对较高的权重，因为在求平均值之前它们被平方，这意味着当大误差不可取时，RMSE应该更有用。
偏差
在统计学中，估计器的偏差是该估计器的期望值与估计参数的真实值之间的差异。期望值为零的估计器或决策规则称为无偏估计器；否则，估计器被称为有偏。换句话说，当算法始终学习相同的错误事物而未能看到准确的基础关系时，就会出现偏差。
例如，如果你试图根据房屋特征来预测房价，并且你的预测始终比实际价格低10万美元，那么你的模型就有偏差。
规模
在统计学中的假设检验中，“检验规模”指的是检验的显著性水平，通常用希腊字母α（alpha）表示。显著性水平或检验规模是一个测试统计量必须超过才能拒绝一个假设的阈值。
它表示当原假设为真时，拒绝原假设的概率，这是一种称为第一类错误或假阳性的错误类型。
例如，如果一个检验设置为5%的显著性水平（α = 0.05），这意味着当原假设为真时，有5%的风险拒绝原假设。这个水平0.05是α的常见选择，尽管根据上下文和研究领域的不同，也可以使用其他水平，如0.01或0.10。
检验规模越小，拒绝原假设所需的证据越强，降低了第一类错误的可能性，但可能增加了第二类错误（当原假设为假时未能拒绝原假设）的机会。在任何统计检验的设计中，第一类和第二类错误之间的平衡是一个重要的考虑因素。
经验规模
在通过蒙特卡洛模拟进行假设检验的背景下，经验规模指的是在原假设为真时，在模拟中错误地拒绝原假设的比例。这本质上是一种模拟版本的第一类错误率。
以下是您将执行此操作的一般过程：
1. 建立您的零假设并为您的测试选择一个显著性水平（例如，α = 0.05）。
2. 假设零假设为真，生成大量样本。样本数量通常很大，例如10,000或100,000，以确保结果的稳定性。
3. 对于每个样本，执行假设检验并记录是否拒绝零假设（您可以将其记录为1表示拒绝，0表示未拒绝）。
4. 计算经验规模，即拒绝零假设的模拟比例。这估计了在给定的检验过程下，当零假设为真时拒绝零假设的概率。
下面的代码展示了如何实现这个。
import numpy as npfrom scipy.stats import ttest_1sampimport randomrandom.seed(10)def calculate_empirical_size(num_simulations: int, sample_size: int, true_mean: float, significance_level: float) -> float: """ 模拟一组样本并对每个样本进行假设检验，然后计算经验大小。 参数： num_simulations（int）：要运行的模拟次数。 sample_size（int）：每个模拟样本的大小。 true_mean（float）：零假设下的真实均值。 significance_level（float）：假设检验的显著性水平。 返回： float：经验大小，即拒绝零假设的测试比例。 """ import numpy as np from scipy.stats import ttest_1samp # 初始化零假设拒绝计数器 rejections = 0 # 运行模拟 np.random.seed(0) # 为了重现性 for _ in range(num_simulations): sample = np.random.normal(loc=true_mean, scale=1, size=sample_size) t_stat, p_value = ttest_1samp(sample, popmean=true_mean) if p_value < significance_level: rejections += 1 # 计算经验大小 empirical_size = rejections / num_simulations return empirical_sizecalculate_empirical_size(1000, 1000, 0, 0.05)对于每个1000次模拟，从均值为0，标准差为1的正态分布中随机抽取大小为1000的样本。进行单样本t检验，检验样本均值是否显著不同于真实均值（在本例中为0）。如果检验的p值小于显著性水平（0.05），则拒绝零假设。
经验大小计算为拒绝零假设的次数（假阳性的次数）除以总模拟次数的比例。在良好校准的测试中，该值应接近名义显著性水平。在这种情况下，函数返回经验大小，这给出了在实际应用中，当零假设实际上为真且条件与模拟相同的情况下，错误拒绝零假设的频率。
由于随机波动，经验大小可能不完全匹配名义显著性水平，但如果样本量足够大且测试假设得到满足，则它们应该接近。经验大小与名义大小之间的差异是为什么进行此类模拟研究的原因，以查看名义大小与实际情况的匹配程度。
A/B测试估计的统
在DGP中不考虑协变量的实验
下面我们模拟数据，这些数据遵循一个只有处理和随机误差影响的DGP。
y_i = tau*T_i+e_i
如果结果只受处理影响，即使样本量相对较小，处理效应参数的估计也是精确的，并且随着样本量的增加而迅速收敛到真实参数值。在下面的代码中，参数tau的值设为2。
tau = 2corr = .5p = 10p0 = 0 # 在DGP中使用的协变量数目Nrange = range(10,1000,2) # 循环N值(nvalues,tauhats,sehats,lb,ub) = fn_run_experiments(tau,Nrange,p,p0,corr)caption = """在没有协变量的随机实验中，处理效应参数的估计"""fn_plot_with_ci(nvalues,tauhats,tau,lb,ub,caption)
对于给定的样本量，验证这与带有截距的回归的结果是相同的。
通过在截距和处理之间运行OLS回归，您可以验证您获得相同的结果。
N = 100Yexp,T = fn_generate_data(tau,N,10,0,corr)Yt = Yexp[np.where(T==1)[0],:]Yc = Yexp[np.where(T==0)[0],:]tauhat,se_tauhat = fn_tauhat_means(Yt,Yc)# n_values = n_values + [N]# tauhats = tauhats + [tauhat]lb = lb + [tauhat-1.96*se_tauhat]ub = ub + [tauhat+1.96*se_tauhat]print(f"通过计算均值差计算的参数估计和标准误差:{tauhat:.5f},{se_tauhat:.5f}")const = np.ones([N,1])model = sm.OLS(Yexp,np.concatenate([T,const],axis = 1))res = model.fit()print(f"通过带有截距的OLS回归计算的参数估计和标准误差:{res.params[0]:.5f},{ res.HC1_se[0]:.5f}")
通过计算均值差计算的参数估计和标准误差:1.91756,0.21187通过带有截距的OLS回归计算的参数估计和标准误差:1.91756,0.21187
运行R次蒙特卡洛迭代并计算偏差、RMSE和大小
现在，您将通过循环一个N参数列表来运行蒙特卡洛模拟，从而增加样本量。您将计算每次迭代的测试RMSE、偏差和实际大小。
这个Python脚本进行了一项实验模拟，研究了样本量（N）对在不考虑协变量的情况下A/B测试性能的偏差、RMSE和大小的影响。让我们一步一步来看：
1. estDict = {}初始化一个空字典，用于存储实验结果。
2. R=2000设置实验的重复次数为2000。
3. for N in [10,50,100,500,1000]循环不同的样本量。
4. 在这个循环内，tauhats=[], sehats=[]初始化为空列表，用于存储每个实验的估计处理效应tauhat及其相应的标准误差se_tauhat。
5. for r in tqdm(range(R)):循环R次实验，进度条由tqdm提供。
6. Yexp,T = fn_generate_data(tau,N,10,0,corr)为每个实验生成合成数据，其中预定义的处理效应tau，观测数量N，10个协变量，没有非零系数的协变量和预定义的相关性。
7. Yt=Yexp[np.where(T==1)[0],:]和 Yc=Yexp[np.where(T==0)[0],;] 将合成数据分成了处理组和对照组。
8. tauhat,se_tauhat=fn_tauhat_means(Yt,Yc) 计算了处理效应估计值和其标准误。
9. tauhats=tauhats+[tauhat]和 sehats=sehats+[se_tauhat] 将处理效应估计值和其标准误添加到相应的列表中。
10. estDict[N]={‘tauhat':np.array(tauhats).reshape([len(tauahts),1]),’sehat':np.array(sehats).reshape([len(sehats),1])} 将估计值存储在字典中，以样本大小为键。
11. tau0 = tau*np.ones([R,1]) 创建了一个大小为 R 的数组，其中所有元素均等于真实处理效应。
12. 对于 estDict 中的每个样本大小，脚本使用 fn_bias_rmse_size() 函数计算并打印处理效应估计值的偏差、RMSE 和大小。
正如预期的那样，随着样本大小的增加，偏差和 RMSE 会降低，而大小会接近真实大小 0.05。
estDict = {}R = 2000for N in [10,50,100,500,1000]: tauhats = [] sehats = [] for r in tqdm(range(R)): Yexp,T = fn_generate_data(tau,N,10,0,corr) Yt = Yexp[np.where(T==1)[0],:] Yc = Yexp[np.where(T==0)[0],:] tauhat,se_tauhat = fn_tauhat_means(Yt,Yc) tauhats = tauhats + [tauhat] sehats = sehats + [se_tauhat] estDict[N] = { 'tauhat':np.array(tauhats).reshape([len(tauhats),1]), 'sehat':np.array(sehats).reshape([len(sehats),1]) } tau0 = tau*np.ones([R,1])for N, results in estDict.items(): (bias,rmse,size) = fn_bias_rmse_size(tau0,results['tauhat'], results['sehat']) print(f'N={N}: bias={bias}, RMSE={rmse}, size={size}')
100%|██████████| 2000/2000 [00:00<00:00, 3182.81it/s]100%|██████████| 2000/2000 [00:00<00:00, 2729.99it/s]100%|██████████| 2000/2000 [00:00<00:00, 2238.62it/s]100%|██████████| 2000/2000 [00:04<00:00, 479.67it/s]100%|██████████| 2000/2000 [02:16<00:00, 14.67it/s]N=10: bias=0.038139125088721144, RMSE=0.6593256331782233, size=0.084N=50: bias=0.002694446014687934, RMSE=0.29664599979723183, size=0.0635N=100: bias=-0.0006785229668018156, RMSE=0.20246779253127453, size=0.0615N=500: bias=-0.0009696751953095926, RMSE=0.08985542730497854, size=0.062N=1000: bias=-0.0011137216061364087, RMSE=0.06156258265280801, size=0.047
带协变量的实验
接下来，您将在 DGP 中添加协变量。现在感兴趣的结果不仅依赖于处理，还依赖于一些其他变量 X。下面的代码模拟了包含在 DGP 中的50个协变量的数据。使用与之前没有协变量的模拟中相同的样本大小和处理效应参数，您可以看到这次估计结果更加嘈杂，但仍然收敛到正确的解决方案。
y_i = tau*T_i + beta*x_i + e_i
下面的图比较了两个DGP的估计值，可以看出当引入协变量时，估计值的噪声要多得多。
tau = 2corr = .5p = 100p0 = 50 # 在DGP中使用的协变量数Nrange = range(10,1000,2) # 循环N值(nvalues_x,tauhats_x,sehats_x,lb_x,ub_x) = fn_run_experiments(tau,Nrange,p,p0,corr)caption = """在DGP中使用X的随机实验的处理效果参数估计，但估计器中不包括X"""fn_plot_with_ci(nvalues_x,tauhats_x,tau,lb_x,ub_x,caption)# 重新运行没有协变量的实验p0 = 0 # 在DGP中使用的协变量数Nrange = range(10,1000,2) # 循环N值(nvalues_x0,tauhats_x0,sehats_x0,lb_x0,ub_x0) = fn_run_experiments(tau,Nrange,p,p0,corr)fig = plt.figure(figsize = (10,6))plt.title("""在DGP中使用X的随机实验的处理效果参数估计，但估计器中不包括X，放大查看大N""")plt.plot(nvalues_x[400:],tauhats_x[400:],label = '$\hat{\\tau}^{(x)}$')plt.plot(nvalues_x[400:],tauhats_x0[400:],label = '$\hat{\\tau}$',color = 'green')plt.legend()fig = plt.figure(figsize = (10,6))plt.title("""带有和不带有协变量的DGP的处理效果估计，放大查看大N""")plt.plot(nvalues_x[400:],tauhats_x[400:],label = '$\hat{\\tau}^{(x)}$')plt.plot(nvalues_x[400:],tauhats_x0[400:],label = '$\hat{\\tau}$',color = 'green')plt.legend()
100%|██████████| 495/495 [00:41<00:00, 12.06it/s]100%|██████████| 495/495 [00:42<00:00, 11.70it/s]
通过使用更大的样本量来重复实验是否能解决这个问题？并非一定。尽管样本量增加了，但估计值仍然显得非常嘈杂。
tau = 2corr = .5p = 100p0 = 50 # 在DGP中使用的协变量数Nrange = range(1000,50000,10000) # 循环N值(nvalues_x2,tauhats_x2,sehats_x2,lb_x2,ub_x2) = fn_run_experiments(tau,Nrange,p,p0,corr)fn_plot_with_ci(nvalues_x2,tauhats_x2,tau,lb_x2,ub_x2,caption)
带有X的DGP — 在回归中加入协变量
在这部分中，您将使用与之前相同的DGP：
y_i = tau*T_i + beta*x_i + e_i
现在，您将在回归模型中包括这些协变量X。您将发现这显著提高了估计值的精确性。但请记住，这有点“作弊” — 在这种情况下，您从一开始就包含了正确的协变量。
在真实世界的情景中，您可能不知道应该包含哪些协变量，并且可能需要尝试不同的模型和协变量。
tau = 2corr = .5p = 100p0 = 50 # 在DGP中使用的协变量数量range = range(100,1000,2) # 循环N的值# 我们需要比p更多的观测标志X = 1(nvalues2,tauhats2,sehats2,lb2,ub2) = fn_run_experiments(tau,Nrange,p,p0,corr,flagX)caption = """使用正确的X进行回归获得的随机试验治疗效应参数的估计"""fn_plot_with_ci(nvalues2,tauhats2,tau,lb2,ub2,caption)
如果同时使用影响结果和不影响结果的一些X会发生什么？
本节将检查在回归模型中包含一些相关和一些不相关的协变量。这模拟了一个真实世界的情况，其中可能不清楚哪些协变量影响结果。
即使包含了一些不相关的变量，可以观察到整体估计与不包含任何协变量的情况相比有所改善。然而，包含无关变量可能会引入一些噪音和不确定性到估计中，它们可能不像只包含相关协变量时那样精确。
总之，了解数据中协变量的影响对于改进A/B测试结果的精确性和可靠性至关重要。本教程通过蒙特卡洛模拟探讨了估计器的均方根误差、偏差和大小等统计特性，并突出了包括DGP和回归模型中的协变量的影响，强调了在实践中进行谨慎的模型选择和假设检验的重要性。
# 使用与之前相同的DGPtau = 2corr = .5p = 100p0 = 50 # 在DGP中使用的协变量数量a = 0.9b = 0.1Nrange = range(100,1000,2) # 循环N的值# 我们需要比p更多的观测标志X = 2(nvalues3,tauhats3,sehats3,lb3,ub3) = fn_run_experiments(tau,Nrange,p,p0,corr,flagX,a,b)caption = f"""使用{100*a:.1f}%正确的X和{100*b:.1f}%不相关的X进行回归获得的随机试验治疗效应参数的估计"""fn_plot_with_ci(nvalues3,tauhats3,tau,lb3,ub3,caption)
除非另有说明，所有图片由作者提供。
y_i = tau*T_i+e_itau的值设为2。tau = 2corr = .5p = 10p0 = 0 # 在DGP中使用的协变量数目Nrange = range(10,1000,2) # 循环N值(nvalues,tauhats,sehats,lb,ub) = fn_run_experiments(tau,Nrange,p,p0,corr)caption = """在没有协变量的随机实验中，处理效应参数的估计"""fn_plot_with_ci(nvalues,tauhats,tau,lb,ub,caption)N = 100Yexp,T = fn_generate_data(tau,N,10,0,corr)Yt = Yexp[np.where(T==1)[0],:]Yc = Yexp[np.where(T==0)[0],:]tauhat,se_tauhat = fn_tauhat_means(Yt,Yc)# n_values = n_values + [N]# tauhats = tauhats + [tauhat]lb = lb + [tauhat-1.96*se_tauhat]ub = ub + [tauhat+1.96*se_tauhat]print(f"通过计算均值差计算的参数估计和标准误差:{tauhat:.5f},{se_tauhat:.5f}")const = np.ones([N,1])model = sm.OLS(Yexp,np.concatenate([T,const],axis = 1))res = model.fit()print(f"通过带有截距的OLS回归计算的参数估计和标准误差:{res.params[0]:.5f},{ res.HC1_se[0]:.5f}")
通过计算均值差计算的参数估计和标准误差:1.91756,0.21187通过带有截距的OLS回归计算的参数估计和标准误差:1.91756,0.21187N参数列表来运行蒙特卡洛模拟，从而增加样本量。您将计算每次迭代的测试RMSE、偏差和实际大小。N）对在不考虑协变量的情况下A/B测试性能的偏差、RMSE和大小的影响。让我们一步一步来看：estDict = {}初始化一个空字典，用于存储实验结果。R=2000设置实验的重复次数为2000。for N in [10,50,100,500,1000]循环不同的样本量。tauhats=[], sehats=[]初始化为空列表，用于存储每个实验的估计处理效应tauhat及其相应的标准误差se_tauhat。for r in tqdm(range(R)):循环R次实验，进度条由tqdm提供。Yexp,T = fn_generate_data(tau,N,10,0,corr)为每个实验生成合成数据，其中预定义的处理效应tau，观测数量N，10个协变量，没有非零系数的协变量和预定义的相关性。Yt=Yexp[np.where(T==1)[0],:]和 Yc=Yexp[np.where(T==0)[0],;] 将合成数据分成了处理组和对照组。tauhat,se_tauhat=fn_tauhat_means(Yt,Yc) 计算了处理效应估计值和其标准误。tauhats=tauhats+[tauhat]和 sehats=sehats+[se_tauhat] 将处理效应估计值和其标准误添加到相应的列表中。estDict[N]={‘tauhat':np.array(tauhats).reshape([len(tauahts),1]),’sehat':np.array(sehats).reshape([len(sehats),1])} 将估计值存储在字典中，以样本大小为键。tau0 = tau*np.ones([R,1]) 创建了一个大小为 R 的数组，其中所有元素均等于真实处理效应。estDict 中的每个样本大小，脚本使用 fn_bias_rmse_size() 函数计算并打印处理效应估计值的偏差、RMSE 和大小。estDict = {}R = 2000for N in [10,50,100,500,1000]: tauhats = [] sehats = [] for r in tqdm(range(R)): Yexp,T = fn_generate_data(tau,N,10,0,corr) Yt = Yexp[np.where(T==1)[0],:] Yc = Yexp[np.where(T==0)[0],:] tauhat,se_tauhat = fn_tauhat_means(Yt,Yc) tauhats = tauhats + [tauhat] sehats = sehats + [se_tauhat] estDict[N] = { 'tauhat':np.array(tauhats).reshape([len(tauhats),1]), 'sehat':np.array(sehats).reshape([len(sehats),1]) } tau0 = tau*np.ones([R,1])for N, results in estDict.items(): (bias,rmse,size) = fn_bias_rmse_size(tau0,results['tauhat'], results['sehat']) print(f'N={N}: bias={bias}, RMSE={rmse}, size={size}')
100%|██████████| 2000/2000 [00:00<00:00, 3182.81it/s]100%|██████████| 2000/2000 [00:00<00:00, 2729.99it/s]100%|██████████| 2000/2000 [00:00<00:00, 2238.62it/s]100%|██████████| 2000/2000 [00:04<00:00, 479.67it/s]100%|██████████| 2000/2000 [02:16<00:00, 14.67it/s]N=10: bias=0.038139125088721144, RMSE=0.6593256331782233, size=0.084N=50: bias=0.002694446014687934, RMSE=0.29664599979723183, size=0.0635N=100: bias=-0.0006785229668018156, RMSE=0.20246779253127453, size=0.0615N=500: bias=-0.0009696751953095926, RMSE=0.08985542730497854, size=0.062N=1000: bias=-0.0011137216061364087, RMSE=0.06156258265280801, size=0.047X。下面的代码模拟了包含在 DGP 中的50个协变量的数据。使用与之前没有协变量的模拟中相同的样本大小和处理效应参数，您可以看到这次估计结果更加嘈杂，但仍然收敛到正确的解决方案。y_i = tau*T_i + beta*x_i + e_itau = 2corr = .5p = 100p0 = 50 # 在DGP中使用的协变量数Nrange = range(10,1000,2) # 循环N值(nvalues_x,tauhats_x,sehats_x,lb_x,ub_x) = fn_run_experiments(tau,Nrange,p,p0,corr)caption = """在DGP中使用X的随机实验的处理效果参数估计，但估计器中不包括X"""fn_plot_with_ci(nvalues_x,tauhats_x,tau,lb_x,ub_x,caption)# 重新运行没有协变量的实验p0 = 0 # 在DGP中使用的协变量数Nrange = range(10,1000,2) # 循环N值(nvalues_x0,tauhats_x0,sehats_x0,lb_x0,ub_x0) = fn_run_experiments(tau,Nrange,p,p0,corr)fig = plt.figure(figsize = (10,6))plt.title("""在DGP中使用X的随机实验的处理效果参数估计，但估计器中不包括X，放大查看大N""")plt.plot(nvalues_x[400:],tauhats_x[400:],label = '$\hat{\\tau}^{(x)}$')plt.plot(nvalues_x[400:],tauhats_x0[400:],label = '$\hat{\\tau}$',color = 'green')plt.legend()fig = plt.figure(figsize = (10,6))plt.title("""带有和不带有协变量的DGP的处理效果估计，放大查看大N""")plt.plot(nvalues_x[400:],tauhats_x[400:],label = '$\hat{\\tau}^{(x)}$')plt.plot(nvalues_x[400:],tauhats_x0[400:],label = '$\hat{\\tau}$',color = 'green')plt.legend()
100%|██████████| 495/495 [00:41<00:00, 12.06it/s]100%|██████████| 495/495 [00:42<00:00, 11.70it/s]tau = 2corr = .5p = 100p0 = 50 # 在DGP中使用的协变量数Nrange = range(1000,50000,10000) # 循环N值(nvalues_x2,tauhats_x2,sehats_x2,lb_x2,ub_x2) = fn_run_experiments(tau,Nrange,p,p0,corr)fn_plot_with_ci(nvalues_x2,tauhats_x2,tau,lb_x2,ub_x2,caption)y_i = tau*T_i + beta*x_i + e_iX。您将发现这显著提高了估计值的精确性。但请记住，这有点“作弊” — 在这种情况下，您从一开始就包含了正确的协变量。tau = 2corr = .5p = 100p0 = 50 # 在DGP中使用的协变量数量range = range(100,1000,2) # 循环N的值# 我们需要比p更多的观测标志X = 1(nvalues2,tauhats2,sehats2,lb2,ub2) = fn_run_experiments(tau,Nrange,p,p0,corr,flagX)caption = """使用正确的X进行回归获得的随机试验治疗效应参数的估计"""fn_plot_with_ci(nvalues2,tauhats2,tau,lb2,ub2,caption)# 使用与之前相同的DGPtau = 2corr = .5p = 100p0 = 50 # 在DGP中使用的协变量数量a = 0.9b = 0.1Nrange = range(100,1000,2) # 循环N的值# 我们需要比p更多的观测标志X = 2(nvalues3,tauhats3,sehats3,lb3,ub3) = fn_run_experiments(tau,Nrange,p,p0,corr,flagX,a,b)caption = f"""使用{100*a:.1f}%正确的X和{100*b:.1f}%不相关的X进行回归获得的随机试验治疗效应参数的估计"""fn_plot_with_ci(nvalues3,tauhats3,tau,lb3,ub3,caption)