简介Python-Causality 是面向数据分析与因果推断学习者的开源工具包适合希望在真实数据集中从相关性进一步探索变量因果结构的研究人员、数据科学从业者及进阶学习者可用于社会科学、经济、生物医学与营销等场景的因果发现与效应评估。压缩包共29个文件约74KB以18个Python源码文件为主体辅以3张PNG图示、2份CSV示例数据、2份Markdown说明文档及少量配置与测试文件其中代码按因果推断、估计与工具函数等模块组织便于按目录理解实现思路。已有3293人学习说明其在因果分析入门与实践中具有一定参考价值。内容涉及结构方程模型、贝叶斯网络、格兰杰因果检验、倾向得分匹配、工具变量法和Do-Calculus等常见方法并配有示例数据与图示读者可据此梳理从数据预处理、模型选择到因果发现与结果可视化的完整流程也可作为课程实验或项目原型的代码参考。1. 相关性跑在前面、因果图落在后面这份工具要解决的就是这个错位有条业务线做留存分析回归一跑功能使用次数的系数 0.4p 值小到打印不出来产品同学当场决定全量推。三个月后留存几乎没动。原因不复杂使用次数和留存互为因果二者还共同受活跃度影响回归把三条路径揉成了一个数。Python-Causality 这类数据集因果分析工具要做的就是把「先定图、再估效应」这套流程产品化——结构学习给出候选因果图领域知识做约束用后门准则算出可识别的调整集最后接一个估计量给出数值和区间。它适合手里只有观测数据、做不了对照实验、又必须交出一个能解释的效应量的数据分析同学会用 pandas 做过 groupby、看得懂 p 值和置信区间后面每一步都能跟下来。2. 资源包拆开看六个目录怎么串成一条从 CSV 到效应值的链路2.1 目录分工谁负责算图谁负责算数解压之后主体是六个平级目录职责切得比较干净基本按「数据进、图出来、数出来」的顺序排。很多人上来就翻estimate/找算法其实真正决定结论对错的是前面两个目录。目录职责输入输出data/示例数据与字段字典csv / parquet标准 DataFramepreprocess/类型推断、缺失处理、编码、时间切分原始表建模矩阵 元数据discovery/结构学习PC / GES / LiNGAM建模矩阵 先验约束邻接矩阵、CPDAG、DAGidentify/后门、前门调整集求解DAG 处理变量 结果变量调整集列表与识别性说明estimate/匹配、IPW、AIPW、bootstrap 区间调整集 数据效应值、标准误、区间report/图可视化、平衡性诊断、报告导出全部产物markdown / html模块之间不靠全局变量传递而是靠一个配置对象加一份元数据对象串起来preprocess/把列类型、缺失率、是否处理后变量写进元数据discovery/和identify/都读这份元数据做校验。这个设计在实操里最省事的地方是——你把某个列标成「处理后变量」后面任何一个估计量都不会偷偷把它塞进调整集省掉了一类最难查的沉默错误。2.2 先定数据契约再动代码因果分析跟普通建模最大的差别是普通建模少一列顶多掉点精度因果分析多一列就直接改结论。所以配置要写在跑代码之前而且写成文件别写成临时变量——三个月后回看结论时你需要知道当时到底调了哪些。# config/dataset_card.yaml treatment: promo_expose # 处理变量二值 0/1连续处理另走 dose 模式 outcome: retention_14d # 结果变量0/1 或连续 time_index: dt # 时间列用于样本外切分和安慰剂检验 covariates: include: [age, tenure_days, city_tier, session_cnt_7d] exclude: [post_click_cnt, order_cnt_after] # 处理后变量绝不能进调整集 min_unique: 2 # 取值数少于一律丢弃避免常数陷阱列 missing: drop_threshold: 0.3 # 单列缺失率超过 30% 直接丢 fill: median # 数值列中位数填充类别列补 __NA__参数里最值得盯的是exclude和drop_threshold。前者是因果分析的安全带凡是发生在处理动作之后、或者由处理动作引起的列都要显式列进去后者不要设得太激进缺失率高的列往往正好是信息量大的列30% 是个保守起点真要放宽得先看缺失是不是随机产生。min_unique放在预处理阶段做比等到结构学习报「奇异矩阵」再回头查要省事得多。提示exclude里的列不会从表里删掉只是不参与调整集计算后面做敏感性分析还用得上。2.3 最小链路四行调用跑出第一个效应值配置写完之后先跑最小链路不要一上来就调算法参数。这条链路的价值在于验证「数据读得进、图算得出、调整集非空、估计量收敛」这四件事任何一步断了后面调参都是白调。from causality import Dataset, Discovery, Identify, Estimate # 1) 读数据并联接字段契约元数据随数据集对象走 ds Dataset.from_csv(data/demo.csv, cardconfig/dataset_card.yaml) ds.clean() # 按契约里的 missing 段做缺失与类型处理 # 2) 结构学习先用最保守的 PCalpha 卡 0.05 dag Discovery(methodpc, alpha0.05, max_cond_size3, seed42).fit(ds) # 3) 识别给定处理变量和结果变量求后门调整集 adj Identify(dag).adjustment_set(treatmentpromo_expose, outcomeretention_14d) print(adj) # 形如 [age, tenure_days, city_tier] # 4) 估计AIPW bootstrap 区间 est Estimate(methodaipw, n_boot500, seed42) res est.fit(ds, adj, promo_expose, retention_14d) print(res.ate, res.ci_lower, res.ci_upper)clean()不带参数时全部读契约想临时改缺失阈值就传ds.clean(drop_threshold0.2)但改完记得同步回契约文件否则下次复现就是两套逻辑。max_cond_size3是 PC 里限制条件集大小的开关样本一两千行的时候 3 已经偏激进后面第 3 章会细说。n_boot500是 bootstrap 次数跑一次大概几十秒调参阶段可以先设n_boot0只出点估计定稿前再补区间。seed一定要写死PC 的条件独立性检验在某些实现里会受样本顺序影响bootstrap 更是必须固定随机源。2.4 跑通之后的第一件事是看元数据不是看结果链路跑通的瞬间绝大多数人第一反应是去看效应值。血泪经验是先打印元数据。元数据里记录了每列的推断类型、缺失填充方式、是否被排除这几项任何一个出问题效应值都会看起来「挺合理」但完全错。最典型的两种情形一是某个整数列被判别成类别变量进了结构学习就变成几十个哑变量图直接稠密到没法看二是日期列没被识别成时间列time_index那一栏空着后面所有基于时间的切分和安慰剂检验都会静默失效——不报错只是不做。所以跑通后的第一件事是把ds.meta.summary()打出来扫一眼把类型不对的列写进契约的dtypes覆盖段再重新跑一遍链路。这一步花两分钟能省掉后面半天的排查。3. 因果发现PC、GES、LiNGAM 的假设差异与四个必调参数3.1 三种算法的假设边界选错后面全白做结构学习没有「哪个更强」这回事只有「假设成不成立」。资源包里三个算法的差异主要在函数形式、噪声分布和能否处理隐变量三件事上选之前先拿自己的数据对一遍假设。算法核心假设输出形式适合的数据明显不适用的场景PC因果马尔可夫性 忠实性靠条件独立性检验逐条剪边CPDAG含无向边样本量中等、关系可能非线性样本量少、变量多检验功效不足GES用 BIC 之类的评分函数做加边减边搜索CPDAG稠密图、样本量较大变量数上百搜索空间爆掉LiNGAM线性 噪声非高斯靠独立性做定向完整定向的 DAG连续变量、噪声明显非高斯变量离散、或关系强非线性实操里的判断顺序我一般是这样先看变量是不是都是连续或近似连续的是的话给 LiNGAM 一次机会它出的是全定向图可解释性最好如果有一半是类别变量退到 PC如果 PC 出来的图稀得像断线的风筝再换 GES 试。三者都跑一遍对比边集不是为了挑一个「最好看的」而是看哪些边在三种方法下都稳定出现——那部分才敢往调整集里放。3.2 alpha 与 max_cond_size最容易被随手写死的两个参数alpha是条件独立性检验的显著性水平max_cond_size是条件集的最大阶数。这两个参数在教程里通常写成 0.05 和 3然后一辈子不动。实际项目里它们的取值直接决定图的密度alpha调小独立性更难被拒绝边更多、图更稠调整集随之膨胀max_cond_size调小高阶条件独立性检验被跳过本该被剪掉的边留下来。rows [] for alpha in (0.01, 0.05, 0.10): for mcs in (1, 2, 3): d Discovery(methodpc, alphaalpha, max_cond_sizemcs, seed42).fit(ds) rows.append({ alpha: alpha, max_cond_size: mcs, n_edges: d.n_edges(), # 子采样 20 次看边集重合度衡量稳定性 jaccard: d.stability(ds, n_split20, seed42), }) print(sorted(rows, keylambda r: -r[jaccard])[:5])stability()的做法是把数据随机切一半各跑一次统计两次边集的 Jaccard 相似度重复n_split轮取均值。参数上alpha不要超过 0.1否则假阳性边会明显增多max_cond_size一般在 2 到 4 之间变量数五十以内选 3上百个变量的时候果断压到 2不然单次检验的样本被切得太碎功效低到没有意义。选参数的标准不是「边最少」也不是「p 值最好看」而是找 Jaccard 的平台区随着alpha变小、max_cond_size变大Jaccard 会先升后平落在平台区上的那组参数就是稳定的如果整片区域 Jaccard 都低于 0.6说明数据量根本撑不起这个变量集该做的是减变量不是继续调参。3.3 把领域知识写成约束白名单、黑名单与时间分层纯靠数据学出来的图边方向经常是玄学。好在结构学习这一层留了约束入口能把「结果不可能导致处理」「入职天数一定在处理之前」这类常识硬编码进去代价小、收益大。priors { # 黑名单禁止出现的边方向结果变量不可能反向导致处理变量 forbidden: [(retention_14d, promo_expose)], # 白名单必须保留的边业务上确定存在 required: [(tenure_days, promo_expose)], # 时间分层层号小的变量不能由层号大的变量导致 tiers: {tenure_days: 0, promo_expose: 1, session_cnt_7d: 2, retention_14d: 3}, } dag Discovery(methodpc, alpha0.05, max_cond_size3).fit(ds, priorspriors) print(dag.violations()) # 返回被约束改写的边逐条看有没有改错forbidden和required是硬约束搜索过程中直接排除tiers是软约束只禁止逆时间的边方向同层之间的边不受影响。三个参数里tiers最好用因为它把业务常识里的时序关系一次性表达完了不用一条条列边。风险在于约束本身写错——把本该是结果的变量写进上层图会朝着错误方向收敛而且不报错。所以每次加完约束都要看一遍violations()确认被改写的边和你的预期一致。3.4 图不稳的时候先看抖动来源再换算法同一批数据跑两次边集差一截第一反应往往是换个算法。但在换之前值得花十分钟定位抖动从哪来因为抖动本身是信息。常见的三个来源一是分类变量的稀有类别某个取值只有十几条样本条件独立性检验在这种格子上的统计量完全不可信二是缺失填充引入的伪相关中位数填充会把分布压出一个尖峰跟处理变量一交叉就出假边三是时间维度的分布漂移把两年数据混在一起跑等于把两个不同的因果结构平均成一个。这三种情况换算法都救不了只能回到预处理阶段改稀有类别合并、缺失改成加一列缺失指示、按时间窗分段各跑一遍。分段跑出来如果图结构差异很大那结论本身就该分群报告而不是硬凑一个全局效应。这一步看着麻烦但它决定了后面估计出来的那个数到底能不能拿去汇报。4. 效应估计从调整集到置信区间中间隔着三道坎4.1 拿到调整集先做重叠性检查调整集算出来只是拿到一张候选名单能不能估取决于处理组和对照组在协变量上有没有可比的人。这一步叫重叠性检查跳过它直接上估计量是最常见的翻车姿势。from causality.diagnostics import propensity_overlap ps propensity_overlap(ds, treatmentpromo_expose, covariatesadj) print(ps.summary()) # 输出包含倾向得分模型 AUC、共同支撑区间、越界样本占比 ps.plot(report/overlap.png)summary()里三个数要连起来看。AUC 衡量协变量对处理的解释力0.6 到 0.8 属于正常超过 0.9 意味着两组人几乎不重叠任何估计量都只能靠模型外推撑着。共同支撑区间是两组倾向得分都覆盖的范围区间越窄、越界占比越高估计结果越是「向中间那一小撮人」的效应外推到全量就会失真。plot()出的图我一般会连着summary()一起存档汇报的时候比一个孤零零的数字有说服力。4.2 倾向得分匹配caliper 与平衡性表怎么看匹配是最直观的估计方式也是最容易被误用的——很多人匹配完只看效应值不看平衡性。匹配的目的本来就是让两组协变量分布对齐没对齐的话效应值算得再精细也是错的。from causality import Match m Match(caliper0.05, ratio1, replaceFalse, seed42) m.fit(ds, promo_expose, adj) print(m.balance_table()) # 每列匹配前后的标准化均值差 SMD print(m.att()) # 处理组的平均处理效应caliper0.05是在倾向得分尺度上设卡尺经验值取倾向得分标准差的 0.2 倍左右0.05 对应中等严格程度数据量大可以收紧到 0.02样本紧张就放到 0.1但超过 0.1 之后匹配就名存实亡了。ratio1是一对一replaceFalse表示不放回放回能提高匹配质量但会重复使用少数对照样本让区间偏窄。balance_table()是这一步的核心产物看的是每列匹配前后的 SMD判据一般取 0.1匹配后所有列的 SMD 都在 0.1 以内说明平衡性达标有列超过 0.1 就回去改匹配参数或者把这一列从调整集里拆出来单独分层匹配。有些列怎么调都平衡不了通常是倾向得分模型设定有问题比如该加交互项没加或者非线性关系被当成线性的了。4.3 AIPW 与 bootstrap点估计一样、区间差一倍的原因双重稳健估计AIPW的好处是倾向得分模型和结果模型只要有一个设定对估计就一致。但它的区间比点估计更容易出问题因为区间依赖于重抽样的做法。est Estimate(methodaipw, n_boot500, stratifyTrue, seed42) res est.fit(ds, adj, promo_expose, retention_14d) print(res.ate, res.ci_lower, res.ci_upper) print(res.diagnostics()) # 结果模型残差、倾向得分的极端值个数stratifyTrue是关键参数意思是重抽样时按处理标签分层保证每个 bootstrap 样本里处理组和对照组的比例不变。不分层的话处理组占比本来就低的数据里某些重抽样样本可能只剩个位数处理样本倾向得分模型直接崩区间要么异常宽、要么异常窄——这就是同一个点估计配出两个完全不同区间的常见原因。n_boot500用来出报告够用要发论文级别可以上到 2000但边际收益递减得很快。diagnostics()里重点看倾向得分的极端值个数接近 0 或 1 的样本如果超过总量 5%说明重叠性本身就有问题这时候该做的是截断把倾向得分卡在 0.05 到 0.95 之间并报告截断前后的差异而不是硬报一个区间。4.4 ATE、ATT 与目标人群别把两个数混着报ATE 是全体人群的平均效应ATT 是处理组人群的平均效应两者的目标人群不一样数值可以差出一倍以上。业务上最常见的情形是处理只覆盖了一部分人比如只对高活跃用户推了活动这时候能估准的是 ATTATE 得靠模型外推把它报出去容易误导决策。资源包里通常用target参数切换targetatt会走匹配路径targetate走 IPW 或 AIPW。切换之后结论方向都变了的话说明重叠性有问题——同一批数据在两种目标人群下效应符号相反不是发现了什么新规律而是外推区间里根本没数据。这种情况我一般两个数都报配上重叠性诊断和人群占比让看报告的人自己判断。5. 避坑记录五条翻车过的问题按现象、原因、解决写5.1 变量层的两个坑时间关系和稀有类别坑一处理后变量混进调整集。现象是效应值被压得极小有时甚至反号换估计量也没用。原因是某个字段看着像协变量实际发生在处理动作之后——比如活动曝光之后的点击次数、下单后的订单数。这类变量是处理变量的下游放进调整集等于把处理到结果的通路截断了一部分。解决办法是在契约里显式写进exclude更稳的做法是给每列标注发生时间结构学习阶段用tiers做时序约束兜底。这件事没有后悔药只能靠数据字典和时间戳逐个核。坑二稀有类别造出假边。现象是同一批数据重跑两次某几条边时有时无Jaccard 相似度死活上不去。原因是某个分类变量的某个取值只有十几条样本条件独立性检验在这种小格子上统计量不可信检验结果随样本组成左右横跳。解决是把占比低于 1% 的类别合并成「其他」或者干脆把这一列从发现阶段排除、留到估计阶段做分层。判断标准很直接跑一次groupby看取值分布长尾明显就动手合并。5.2 结构学习层的两个坑小样本与离散化坑三小样本 PC 输出近似全连接图。现象是画出来的图密密麻麻调整集列出十几个变量估计阶段直接因为共线性崩掉。原因是 PC 的剪边全靠条件独立性检验样本量不足时检验功效低本该被剪掉的边全都活下来。解决是三步先把max_cond_size压到 2再把alpha提到 0.1 放宽拒绝标准还不行就加forbidden/required约束把没争议的边先定死。变量数超过样本量十分之一的时候别指望纯数据学出可用的图减变量才是正解。坑四连续变量分箱后再做结构学习。现象是明摆着的强关联被判成独立。原因是把连续变量切成五档之后原本的单调关系被抹成了阶梯条件独立性检验在分箱后的列上抓不到信号更糟的是分箱边界由人拍换个切法结论就变。解决是发现阶段用原始连续变量非参数检验撑得住确实需要离散化的放到估计阶段做分层并且把分箱规则固化下来单独汇报。5.3 估计与推断层的一个坑结论随估计量翻转坑五换估计量结论就反。现象是匹配出来的效应对照组更高AIPW 出来的反过来两个数都在置信区间内但方向不同汇报时没法开口。原因通常有两个一是重叠性差不同估计量对外推的处理方式不一样差的那部分样本主导了结果二是结果模型严重误设比如真实关系是阈值型的模型却按线性拟合。解决顺序是先查重叠性再查模型倾向得分 AUC 高于 0.9 或者越界样本超过 10%就把结论限定在共同支撑区间内报同时给出截断前后的对比重叠性没问题的话回来检查结果模型有没有加交互项和非线性项可以拿简单分层均值差做一次交叉验证——三种方法方向一致才敢下结论。6. 让结论扛得住追问三个检验和一份能复现的报告效应值算出来只是半成品真正让结论立得住的是后面三个检验。第一个是安慰剂处理把处理标签随机打散重跑整套估计流程重复两百次真实效应应该落在安慰剂分布的尾部之外如果真实值稳稳待在安慰剂分布的正中间说明整条链路在某处泄漏了处理信息。import numpy as np placebo [] for i in range(200): d ds.shuffle_treatment(promo_expose, seedi) # 打散处理标签保留协变量 r Estimate(methodaipw, n_boot0, seedi).fit(d, adj, promo_expose, retention_14d) placebo.append(r.ate) placebo np.array(placebo) p (np.abs(placebo) abs(res.ate)).mean() # 经验 p 值 print(fplacebo mean{placebo.mean():.4f}, p{p:.3f})shuffle_treatment只打散处理列协变量和结果列不动这样协变量与结果的关系保留下来检验的是「处理效应是否真的来自处理变量」。n_boot0是因为安慰剂检验只看点估计的分布跑区间纯属浪费时间。经验 p 值小于 0.05 才算过得去接近 0 更好。第二个是负对照结果找一个理论上不该受处理影响的变量比如用户注册渠道把它当结果重跑一遍效应应该接近 0。第三个是时间样本外验证用time_index切出后一段数据在前一段上定参数、在后一段上复算效应方向一致才算稳。检验做法期望结果异常时的排查方向安慰剂处理打散处理标签重跑 200 次真实值落在分布尾部经验 p 0.05检查是否有处理后变量泄漏负对照结果换成理论上无关的结果变量效应接近 0 且区间跨 0检查共同原因是否被漏掉时间样本外前段定参、后段复算方向一致数值差距可解释检查分布漂移与结构变化这三项我现在的习惯是写进一个report/validation.md跟效应值、平衡性表、重叠性图放在一起形成一份别人能照着复现的产物。参数、随机种子、契约文件路径全部记进去谁想复核就按文件重跑一遍。有次项目赶进度我把重叠性检查跳了直接报了一个漂亮的效应值两周后被追问「对照组里根本没有类似的人」——那批数据里处理组的倾向得分集中在 0.8 以上对照组集中在 0.2 以下估计量全程在外推。从那以后我每次跑完估计都强制走一遍「重叠性 → 平衡性 → 安慰剂」这三关哪怕结论已经很明显也不省这一步。资源包里的东西够用但顺序不能乱。希望帮到你。本文还有配套的精品资源点击获取