如何用Python统计方法估计两变量间的因果关系及因果效应？

首先要明确：统计方法无法绝对“证明”因果关系，只能通过合理的分析设计结合领域知识，强化x是y原因的证据，并量化潜在的因果效应。以下是针对你的数据集可操作的思路和Python实现：

一、因果推断的核心前提 #

要建立x到y的因果关系，必须满足三个核心条件：

• 相关性：你已经通过test.corr()得到0.96的高相关系数，满足这一点
• 时间先后：x的变化必须发生在y的变化之前（如果你的数据是时序数据，可通过观测顺序确认；如果是横截面数据，需要领域知识判断）
• 排除混淆：不存在其他未观测/未控制的变量同时影响x和y（这是观测数据因果推断的核心难点）

二、Python中的分析方法 #

1. 线性回归（初步量化边际效应） #

线性回归的系数可以给出x对y的边际关联效应，如果能确保无混淆变量，这个系数可近似为因果效应：

import statsmodels.api as sm

# 给自变量添加截距项
X = sm.add_constant(test['x'])
# 拟合OLS模型
model = sm.OLS(test['y'], X).fit()
# 输出模型结果
print(model.summary())

• 结果中x对应的系数就是：x每增加1单位，y平均变化的数值。但注意：如果存在未控制的混淆变量，这个系数只是关联，不是因果。

2. 格兰杰因果检验（时序数据的时间先后证据） #

如果你的数据是按时间顺序收集的（从行号0到62的趋势看，y随时间递增，x也整体递增），可以用格兰杰检验判断x是否能显著预测y的未来值，这能提供“x先于y并影响y”的证据：

from statsmodels.tsa.stattools import grangercausalitytests

# 整理成x和y的时序数据集
data = test[['x', 'y']]
# 检验x是否格兰杰引起y，maxlag可根据数据调整
grangercausalitytests(data[['y', 'x']], maxlag=3)

• 查看结果中的p值，如果p<0.05，说明x的滞后项能显著预测y的变化，支持x在时间维度上的先导性，但仍不能绝对证明因果。

3. 控制时间趋势（排除共同趋势干扰） #

从你的数据看，y和x都随行号（假设是时间）递增，可能是共同的时间趋势导致的高相关，而非x直接影响y。可以加入时间变量控制这种趋势：

# 添加时间索引作为控制变量
test['time'] = test.index
# 拟合包含时间的回归模型
X = sm.add_constant(test[['x', 'time']])
model = sm.OLS(test['y'], X).fit()
print(model.summary())

• 如果加入时间后x的系数仍然显著且数值变化不大，说明x对y的关联不受时间趋势干扰；如果系数变得不显著，说明之前的高相关是时间趋势导致的。

4. 倾向得分匹配（PSM，适用于有混淆变量的情况） #

如果你有其他混淆变量（比如影响x和y的第三方变量），可以用PSM匹配x不同取值下混淆变量相似的样本，再比较y的差异，以此估计因果效应：

from causalinference import CausalModel

# 先将连续的x转为二元处理变量（比如以0.5为分界）
test['x_treated'] = (test['x'] > 0.5).astype(int)
# 初始化因果模型，X为混淆变量（替换为你的实际变量）
cm = CausalModel(
    Y=test['y'].values,
    D=test['x_treated'].values,
    X=test[['confounder1', 'confounder2']].values  # 替换为你的混淆变量
)
# 估计倾向得分
cm.est_propensity_score()
# 通过匹配估计因果效应
cm.est_via_matching()
# 输出结果
print(cm.estimates)

• 输出中的ATE（平均处理效应）就是x对y的因果效应估计值，前提是你能识别并控制所有混淆变量。

5. 工具变量法（IV，解决内生性问题） #

如果存在双向因果（y也影响x）或未观测的混淆变量，可以找一个工具变量z：z只影响x，不直接影响y。通过两阶段最小二乘法估计因果效应：

from linearmodels import IV2SLS

# 用公式拟合IV模型，z是工具变量
model = IV2SLS.from_formula(
    'y ~ 1 + x + [x ~ z]',
    data=test
)
results = model.fit()
print(results.summary())

• 这里的x系数就是因果效应的估计值，但工具变量的选择必须严格满足“只影响x，不直接影响y”的条件，这需要领域知识支撑。

三、关键总结 #

1. 观测数据无法绝对证明因果，所有方法都依赖假设（如无混淆、工具变量有效等），这些假设需要领域知识验证。
2. 优先检查数据的时序性和共同趋势，排除虚假相关。
3. 量化因果效应的核心是尽可能控制混淆变量，或通过实验设计（如果可能）获取随机对照数据。

内容的提问来源于stack exchange，提问作者hbstha123