利用基于回归的因果推断方法进行产品实验：使用Python及statsmodels库来估算大语言模型各项特征对模型效果的影响

以下的平衡性检查可以验证随机化处理是否在各个观测变量上得到了有效应用。在“失败模式”部分中，会指出每个实际问题违反了上述四个假设中的哪一项。

当随机化处理得当时，回归分析能够高效地得出因果效应的估计值；而一旦某个假设被违反，回归分析结果就会反映出现实的错误情况，而无法准确体现处理措施的实际效果。如果平衡性分析结果显示某些协变量存在系统性差异，在继续进行估计之前，必须立即停止当前的操作并深入调查数据分配的过程。

先决条件

本教程中的每一个代码块都可以在09_regression/regression_demo.ipynb这个配套笔记本中端到端运行。

您需要使用Python 3.11或更高版本，并且需要对pandas和统计学有基本的了解。在这里，statsmodels这个库可能对您来说比较陌生：它能够通过一次调用就计算出HC3统计量以及基于聚类分析得出的稳健标准误差，而这些是scipy.stats单独使用无法实现的。

请安装以下所需的软件包：

pip install numpy pandas statsmodels scipy

请克隆配套代码仓库以获取合成数据集：

git clone https://github.com/RudrenduPaul/product-experimentation-causal-inference-genai-llm.git
cd product-experimentation-causal-inference-genai-llm
python data/generate_data.py --seed 42 --n-users 50000 --out data/synthetic_llm_logs.csv

设置工作示例环境

这个数据集模拟了50,000名用户分布在50个不同的工作环境中。prompt_variant这一列记录了每个用户被分配到了哪个实验组：1代表新模板组，0代表对照组。

用户的分配是通过对其用户ID进行哈希运算来完成的，因此这种分配方式是随机的，并且与数据中的其他所有信息都无关。

task_completed这一列表示实验结果是否完成。在生成数据的过程中，预设的因果效应为效果增加4个百分点。

在建立任何模型之前，请先确认随机化处理是否使得各组在所有观测变量上达到了平衡状态。一个设计合理的随机化实验应该确保各个实验组在所有被测量的特征上具有相近的平均值。

import pandas as pd
import numpy as np

df = pd.read_csv("data/synthetic_llm_logs.csv")

print("数据集的结构：", df.shape)
print("\n不同实验组的提示选项分布情况：")
print(df.prompt_variant.value_counts().to_dict())

# 随机化处理平衡性检查：各组中各协变量的平均值
check_cols = ["query_confidence", "session_minutes", "cost_usd"]
balance_table = (
    df.groupby("promptVariant")[checkcols]
    .mean()
    .round(4)
    .T
)
balance_table.columns = ["对照组 (variant=0)", "实验组 (variant=1)"]
balance_table["差异"] = (
    balance_table["实验组 (variant=1)"]
    - balance_table["对照组 (variant=0)"]
)
print("\n各协变量的平衡性检查结果：")
print(balance_table)

# 用户参与程度的分布情况
print("\n不同实验组中用户参与程度的比例：")
print(
    df.groupby("prompt_variant")
    .engagement_tier.value_counts(normalize=True)
    .unstack()
    .round(3)
)

预期输出：

[提示：在包含50,000条数据的数据集上运行regression_demo.py程序，以便获取实际数值]

具体操作过程如下：首先加载这50,000条数据，然后计算其中各组的数据数量（大致每组有25,000条记录）。接下来，分别计算对照组和实验组中三个连续变量（query_confidence、session_minutes和cost_usd）的平均值。

这些列所反映的数据是在为样本分配不同处理方案之前收集到的，因此从本质上讲它们属于“处理前的数据”。因此，在每一行中，“Difference”这一列的值都应该非常小。

你还需要检查分类变量“参与程度”（高、中、低）在两组中的分布比例是否相似。轻微的差异属于正常的抽样误差，但如果某个协变量的分布存在系统性偏差，那就说明基于哈希值的分配方式出现了问题，或者数据处理流程在随机化之后引入了选择性偏倚。如果发现这种明显的偏差，请立即停止当前操作，先仔细检查数据分配流程，然后再继续进行后续分析。

在这个数据集上，所有协变量之间的绝对差异都小于0.01，且各组中“参与程度”变量的比例相差也在两个百分点以内，这说明随机化分配过程是有效的。

图2： 在25,000名对照组用户和25,000名实验组用户中，query_confidence这一变量的分布情况。两条曲线几乎完全重合（平均差异为-0.0013），这证明了基于哈希值的随机分配确实能够确保各组之间的协变量分布平衡。这就是对实际数据集进行的诊断分析结果。将其与图1中的示意图进行对比即可。

步骤1：简单均值差异分析法

首先使用最简单的统计方法来进行分析：用实验组的平均结果减去对照组的平均结果。

from scipy import stats

mean_control = df[df.prompt_variant == 0].task_completed.mean()
mean_treatment = df[df.prompt VARIANT == 1].task_completed.mean()

naive_effect = mean_treatment - mean_control

print(f"对照组的平均值：{mean_control:.4f}")
print(f>实验组的平均值：{mean_treatment:.4f}")
print(f>简单效应：{naive_effect:+.4f})

# 手动进行双样本t检验
n0 = (df.prompt_variant == 0).sum()
n1 = (df.prompt VARIANT == 1).sum()
var0 = df[df.prompt.variant == 0].task_completed.var()
var1 = df[df.prompt_variant == 1].task_completed.var()
se = np.sqrt(var0 / n0 + var1 / n1)
t_stat = naive_effect / se

p_val = 2 * stats.t.sf(abs(t_stat), df=n0 + n1 - 2)

print(f"双样本t检验的标准误差：{se:.4f}")
print(f>t统计量的值：{t_stat:.3f}")
print(f>p值：{p_val:.4f}")

预期输出：

[占位符——在包含5万条数据的数据集上运行regression_demo.py程序，以便获取实际数值]

具体操作过程如下：首先计算每个实验组中的平均任务完成率，然后计算两者之间的差异，并使用双样本t检验的合并方差公式来计算标准误差。由于该实验是随机设计的，因此这种简单的差异值其实是一个有效的因果估计结果。

所得到的估计值可能与预先设定的+4个百分点的真实值相差一两个百分点。在这种数据量下，这种差异属于正常的抽样波动，并非估计偏差。在下一步进行OLS回归分析时，如果不加入任何协变量，得到的结果会与这个数值完全一致；而一旦加入协变量，标准误差就会变小。

简单的t检验将所有观测值视为相互独立的。在这个研究中，这种假设是合理的，但在第三步中就不再适用了——因为处于同一工作环境中的用户之间存在相关性，因此使用简单的标准误差会低估实际的不确定性。

步骤2：具有异方差性鲁棒性的OLS分析（HC3）

当没有协变量时，使用二元处理变量对二元结果进行普通最小二乘回归分析，得到的点估计值与计算平均差异所得的结果是相同的。加入协变量后，可以吸收残差方差，从而降低标准误差。

与简单的t检验相比，HC3标准误差具有更大的优势：即使误差项的方差在不同观测值之间存在变化，HC3标准误差仍然有效。

对于有限样本来说，HC3比HC0到HC2更优，因为它会对那些具有较高影响力的观测值施加更严格的惩罚机制，从而在样本量适中时使置信区间的覆盖范围更加准确。

import statsmodels.formula.api as smf

# 不加入协变量的OLS分析：结果应与简单差异值一致
m1 = smf.ols(
    "task_completed ~ prompt_variant",
    data=df
).fit(cov_type="HC3")

print("=== 不加入协变量的OLS分析（HC3） ===")
print(m1.summary().tables[1])
print(f"\n系数：{m1.params['promptVariant']:+.4f}")
print(f"HC3标准误差：{m1.bse['prompt_variant']:.4f}")
print(f"p值：{m1.pvalues['prompt_variant']:.4f}")

预期输出：

[请在50k数据集上运行regression_demo.py程序以获取实际结果]

具体操作过程是：使用HC3鲁棒标准误差且不加入任何协变量来进行OLS分析。结果显示，prompt_variant这个变量的系数与计算平均差异所得的结果在四位小数范围内是一致的，这证明了OLS分析实际上就是一种用于估计平均差异的回归方法。

由于HC3标准误差能够纠正异方差性问题，因此其数值通常会比传统的OLS标准误差略大一些。不过在实际应用中，这种差异往往很小；但无论如何，都应该优先选择HC3方法。如果不需要使用HC3，就不会有任何损失；而如果确实需要它，那么使用它所带来的好处将是显而易见的。

现在，让我们加入协变量吧：

# 定义包含协变量的回归公式
formula = (
    "task_completed ~ promptVariant + query_confidence + "
    "session_minutes + C(engagement_tier)"
)

# 包含协变量的OLS分析：估计值相同，但标准误差更小
m2 = smf.ols(formula, data=df).fit(cov_type="HC3")

print("=== 含有协变量的OLS分析（HC3） ===")
print(m2.summary().tables[1])
print(f"\n系数：{m2.params['prompt_variant']:+.4f}")
print(f"HC3标准误差：{m2.bse['promptVariant']:.4f}")
print(f"p值：{m2.pvalues['prompt_variant']:.4f}")

# 比较两种情况下的标准误差
print("\n--- 标准误差比较 ---")
print(f"不包含协变量时：{m1.bse['prompt_variant']:.4f}")
print(f"包含协变量时：{m2.bse['promptVariant']:.4f}")
print(f"R平方值（包含协变量时）：{m2.rsquared:.4f}")

预期输出：

[提示：在包含50,000条数据的数据集上运行regression_demo.py程序，以便获得真实的分析结果]

具体操作过程如下：你将query_confidence、session_minutes和engagement_tier这些变量作为控制变量纳入分析。这三者都属于预处理阶段测量的变量，在应用不同处理方案之前就已经被记录下来，因此将它们包含在分析中不会导致“碰撞偏差”。

prompt_variant这一变量的系数值与朴素估计值相差不大，因为随机化设计确保了这些协变量与处理方案的分配结果之间不存在相关性。因此，该系数的点估计值保持不变，只是其不确定性有所降低而已。

当纳入这些协变量后，R平方值会从接近零上升到几个百分点，这说明这些协变量确实能够解释部分任务完成情况的变化。随着标准误差的减小，prompt_variant这一变量的HC3 p值也会变得更小。

这就是随机实验中利用协变量进行调整所带来的“免费午餐”效应：只要是任何能够预测研究结果的预处理变量，都可以被纳入分析范围——无论是基线参与度、历史任务完成率，还是注册用户群体等等。不过，必须只选择在处理方案实施之前就已经确定的变量，因为那些可能因处理方案而发生变化的变量并不适合用于此类分析。

步骤3：群组稳健标准误差

在第二步中采用的HC3方法虽然能够处理异方差性问题，但仍然将每个观测值视为独立的。实际上，同一工作空间内的用户会共享同一个支持团队、相同的产品等级、相同的IT政策，而且他们的使用场景往往也相似，因此他们的结果之间存在相关性。

如果新的提示模板在工作空间12中的效果很好，但在工作空间37中效果很差，那么无论是否进行了某种处理，这两个工作空间内的结果仍然会存在相关性。如果忽略这种相关性，计算出的标准误差就会偏小，从而导致t统计量的值被夸大，使分析结果看起来比实际情况更显著。

群组稳健标准误差通过将每个工作空间视为一个独立的信息单元来纠正这一问题。这样一来，处理效应系数的方差就反映了50个工作空间的综合数据情况，而不是50,000次独立的随机实验结果。

# 朴素标准误差（假设各工作空间内的数据是独立的）
m3_naive = smf.ols(formula, data=df).fit(cov_type="HC3")

# 群组稳健标准误差（考虑了各工作空间内的相关性）
m3_cluster = smf.ols(formula, data=df).fit(
    cov_type="cluster",
    cov_kwds={"groups": df["workspace_id"]}
)

print("=== 标准误差比较：HC3与群组稳健标准误差 ===")
print(f"系数值（两者均相同）：{m3_cluster.params['prompt_variant']:+.4f}")
print(f"HC3标准误差：                  {m3_naive.bse['prompt_variant']:.4f}")
print(f"群组稳健标准误差：       {m3_cluster.bse['promptVariant']:.4f}")
print(f"HC3 p值：             {m3_naive.pvalues['promptARIANT']:.4f}")
print(f"群组稳健p值：         {m3_cluster.pvalues['prompt_variant']:.4f}")

# 查明共有多少个工作空间
print(f"\n工作空间的数量：{df.workspace_id.nunique()}")
print(f"每个工作空间的用户数（平均值）：{len(df) / df_workspace_id.nunique():.0f}")

预期输出：

[提示：在包含50,000条数据的数据集上运行regression_demo.py程序，以获得准确的数值结果]

具体操作过程如下：你两次使用了相同的、经过协变量调整的OLS模型进行分析——一次使用HC3方法，另一次则按照workspace_id对数据进行分组处理，然后使用群集稳健误差估计方法。两次分析得出的点估计值是相同的，因为标准误差的选择并不会影响系数本身的数值，只会影响其不确定性。在这个包含50个工作区、每个工作区有1,000名用户的数据集上，群集稳健标准误差会略大于HC3方法计算得到的标准误差，这是因为你的有效样本规模实际上是50个工作区的整体数据，而不是50,000条单独的记录。

需要记住的一条规则是：如果你的实验是在个体层面分配处理组的，但你的数据具有聚类结构（用户分布在不同的工作区内，用户的会话记录存在关联，不同产品之间的使用时间也存在规律），那么就应该在具有自然相关性的单元层面进行聚类分析。如果聚类程度过低，会导致分析结果出现过度自信的情况；而如果聚类的粒度超过实际的数据关联结构，虽然不会影响点估计值的准确性，但会使得标准误差被夸大，从而降低分析的精确度。

当不确定该使用何种聚类方法时，最好选择聚类程度较高的方式。如果聚类的数量少于30个，那么群集稳健标准误差就会变得不可靠，此时应该改用排列检验来进行分析。

步骤4：通过交互作用分析处理效应的异质性

在步骤2和步骤3中使用的OLS模型计算得出的系数，反映的是所有用户平均而言所受到的处理效果。然而，平均值往往会掩盖一些重要的数据特征。新的提示模板可能对使用频率较高的用户有效，但对使用频率较低的用户则毫无作用；或者它可能会为所有类型的用户产生相同的效果。要检测这种异质性，就需要在处理变量与调节变量之间添加交互作用项。

# 交互作用模型：prompt_variant与engagement_tier的交互作用
interaction_formula = (
    "task_completed ~ promptVariant * C(engagement_tier) + "
    "query_confidence + session_minutes"
)

m4 = smf.ols(interaction_formula, data=df).fit(
    cov_type="cluster",
    cov_kwds={"groups": df["workspace_id"]}
)

print("=== 交互作用模型（群集稳健方法）===")
print(m4.summary().tables[1])

# 分析不同层次的用户所受到的处理效应
print("\n=== 不同参与程度的用户所受到的处理效应 ===")
baseline_effect = m4.params["prompt_variant"]
tiers = ["medium", "heavy"]  # 'light'为参考类别

effects = {"light": baseline_effect}
for tier in tiers:
    interaction_key = f"promptVariant:C(engagement_tier)[T.{tier}]"
    if interaction_key in m4.params:
        effects[tier] = baseline_effect + m4.params[interaction_key]
    else:
        effects[tier] = baseline_effect

for tier, eff in effects.items():
    print(f"  {tier:8s}: {eff:+.4f}")

# 进行联合F检验：这些交互作用项是否都具有显著性？
interaction_terms = [k for k in m4.params.index if "prompt_variant:C" in k]
if interaction_terms:
    f_test = m4.f_test([f"({t} = 0)" for t in interaction_terms])
    print(f"\n交互作用项的联合F检验结果：p值 = {f_test.pvalue:.4f}")

预期输出：

[占位符——在包含5万条数据的数据集上运行regression_demo.py程序，以获得实际结果]

具体操作过程如下：你在prompt_variant与C(engagement_tier)之间添加了交互项。其中light层级被设为参考类别，因此promptVariant的系数实际上反映的是针对“轻度用户”这一群体的具体效果；而当加入medium或heavy层级的交互系数后，就能得到这些相应层级中的处理效应。

对所有交互项进行联合F检验，目的是确定这些效应在不同层级之间是否存在差异，而这些差异并非由抽样误差造成的。如果检验结果不显著，说明提示模板的效果在各个参与程度不同的群体中都大体一致；而如果结果显著，那就意味着应该分别报告针对不同层级的具体效果，并优先在那些效果提升最为明显的层级中推广该提示模板。

要正确运用交互效应分析模型，就必须严格遵守研究规范。在开始分析数据之前，必须事先确定好自己打算测试哪些变量。如果进行10次这样的分析，却只报告其中一次在p值小于0.05的情况下显示出显著性的结果，这种做法实际上属于多重比较错误，或者说是一种伪装成“亚组分析”的统计欺诈行为。

如果你是在没有事先进行注册登记的情况下探索新的数据集，那么应该使用Bonferroni校正方法，或者采用虚假发现率控制程序，并且明确说明自己的分析属于探索性研究范畴。

步骤5：使用自助法计算置信区间

普通最小二乘法得出的点估计值虽然效率较高，但使用自助法计算的置信区间则不会受到分布假设的影响。具体操作方法是：重复进行500次抽样分析——每次都从原始数据集中有放回地抽取样本，然后重新拟合集群稳健模型，并记录每次计算得到的处理效应系数。最终，该分布中第2.5百分位和第97.5百分位的值就是所求的95%置信区间。

rng = np.random.default_rng(seed=7)
n_boot = 500
boot_coefs = []

for _ in range(n_boot):
    idx = rng.integers(0, len(df), size=len(df))
    boot_df = df.iloc[idx].reset_index(drop=True)
    boot_model = smf.ols(
        formula,
        data=boot_df
    ).fit(
        cov_type="cluster",
        cov_kwds={"groups": boot_df["workspace_id"]}
    )
    boot_coefs.append(boot_model.params["prompt_variant"])

boot_coefs = np.array(boot_coefs)
ci_low, ci_high = nppercentile(boot_coefs, [2.5, 97.5])

print(f"自助法计算的95%置信区间：[{ci_low:+.4f}, {ci_high:+.4f}]")
print(f"自助法计算的平均值：{boot_coefs.mean():+.4f}")
print(f"基于分析模型的标准误差：{m3_cluster.bse['prompt_variant']:.4f}")

预期输出：

[占位符——在包含5万条数据的数据集上运行regression_demo.py程序，以获得实际结果]

具体操作过程是：对整个数据集进行500次有放回的抽样，每次抽样后都重新拟合经过协变量调整后的集群稳健模型，并记录每次计算得到的处理效应系数。通过这种方式得到的分布能够同时反映抽样误差以及数据中的群体结构特征。一个有效的自助法置信区间应该既能包含真实的效应范围，又不会包含零值；而自助法计算得出的平均值也应该与基于分析模型的点估计值非常接近。如果两者之间存在较大差距，那就说明分析模型对某些特定观测值过于敏感。

当单独使用回归分析不够时

在随机化环境下进行的回归分析能够清晰地揭示因果关系，因为随机化措施切断了处理组与混杂因素之间的关联。然而，实际开发的大型语言模型系统很少会进行纯粹的实验设计。下面列出的每种故障模式都对应于前面提到的四个假设中的某一个。

观察性数据中未测量的混杂因素

假设你的团队从未对提示模板进行随机化处理，而是让高置信度的查询默认使用新的模板。此时，prompt_variant与query_confidence之间存在很强的相关性，而query_confidence本身又可以预测task_completed的结果。

这种情况违反了随机分配的假设（E[ε|D] = 0）：误差项不再与处理组无关了。

使用普通最小二乘法进行分析时，会将部分置信效应归因于模板本身，从而导致对处理效应的高估。只有当你已经测量并正确识别了这些混杂因素时，将query_confidence作为控制变量才能消除这种偏差。

任何未被测量的、既影响分配方式又影响结果的因素都会直接反映在普通最小二乘法的计算结果中。因此，必须对这些混杂因素进行测量并将其纳入控制范围，或者采用工具变量分析或不连续设计来恢复随机分配的效果。

SUTVA假设的违反与溢出效应

普通最小二乘法假定每个用户的最终结果仅取决于他们所被分配到的处理组（即SUTVA假设，也就是上面列出的第三个识别假设）。

在多用户协作环境中，这一假设很容易被打破。如果某个工作空间中的核心用户开始使用新的提示模板，并影响其他成员的查询方式，那么同一工作空间中的普通用户也会间接受到这种影响，从而产生处理效应。此时，你的结果就不再仅仅取决于你自己所被分配到的处理组了。

虽然聚类稳健标准误差能够考虑这种相关性，但系数仍然无法区分直接效应与溢出效应。要检测溢出效应，就需要采用两级随机化设计：首先将所有工作空间分为处理组和对照组，然后再分别测量每个组内所有用户的最终结果。

随时间变化的混杂因素

如果提示模板是在某个特定时间点分配给用户的，但由于产品更新、技术支持问题或季节性使用习惯的变化，用户在分析期间的使用行为发生了改变，那么处理组与结果之间的关联关系可能会发生波动，而普通最小二乘法无法将这些波动与真正的因果效应区分开来。

这种情况违反了随机分配的假设——当协变量分布在分配后发生变化时，处理组的分配方式就不再与最终结果无关了。

在这种情况下，你需要采用面板数据分析方法，并设置特定时期的控制变量，或者使用工具变量来解释这些随时间变化的因素对结果的影响。

二元结果与线性概率模型

任务完成的结果只能是0或1。对于二元结局变量，普通最小二乘法实际上就是线性概率模型，这种模型适用于估算平均处理效应，在A/B测试的背景下，其结果也比逻辑回归更易于理解。

这种方法的局限性在于它依赖于线性假设：当某些协变量的取值极端时，线性条件期望计算得出的预测概率可能会超出[0, 1]的范围。虽然这并不会影响平均处理效应的估算结果，但会导致针对单个个体的预测结果不可靠。因此，当你需要得到经过校准的概率分数时，应该使用逻辑回归；而当你需要一个易于理解的平均处理效应时，则应选择普通最小二乘法。

下一步该怎么做

当实验设计符合要求且四个基本假设都成立时，通过以下四个步骤就能得到完整的研究结果：首先计算简单均值差异，然后使用HC3方法进行分析，接着进行群组鲁棒性检验，最后研究预先设定的交互作用效应。务必确保随机化分配过程准确无误，运行平衡性检查，并按照变量之间的相关关系自然地将样本分组。每完成一步，置信区间都会变得更窄，这样在做出推广决策时，你就能清楚地知道自己的数据能够提供多高的精确度。

但如果实验设计存在缺陷，所需的分析方法也会相应改变。对于那些在参与过程中存在选择性的观察性数据，就需要使用倾向得分方法或对丰富的协变量集进行回归调整；而当分配依据是连续阈值时，则需要运用回归不连续性分析；如果实验是在不同时间段内分阶段进行的，那么差分分析法就是必要的工具。

每种方法都是针对普通最小二乘法无法解决的特定混杂因素设计的，同时也能帮助我们判断实验设计违反了哪一条识别假设。

本教程配套的笔记本文件位于github.com/RudrenduPaul/product-experimentation-causal-inference-genai-llm/tree/main/09_regression。你可以克隆这个仓库，生成合成数据集，然后运行regression_demo.py文件，从而完整地体验本教程中每一步的代码实现过程。