六种因果库对比：哪种贝叶斯方法能发掘数据中的隐藏因果关系？

告别相关关系，直击因果本质

在数据科学的世界里，我们常常满足于构建一个精准的预测模型。但很多时候，我们真正想知道的是：“为什么？”

是哪个变量真正驱动了结果的变化？哪些变量只是搭便车的乘客，仅仅与结果相关却无直接因果？

理解变量间的因果关系至关重要。无论是进行预测性维护、优化营销策略，还是制定关键业务决策，识别出真正的“驱动变量”都能让我们事半功倍。

然而，面对众多功能各异的贝叶斯和因果推断库，初学者往往感到无所适从。别担心，本文将对Python中六大流行的因果推断库进行横向对比，通过同一数据集上的实战示例，帮你一次性理清它们的优缺点。

在因果推断中，我们超越了预测，去识别事件发生的机制和途径。

参与评测的六位选手是：

• Bnlearn
• Pgmpy
• CausalNex
• DoWhy
• Pyagrum
• CausalImpact

文章结尾有详细的总结对比表格，帮你快速选出最适合你项目的工具库！

因果推断：从“是什么”到“为什么”

简单来说，因果推断就是确定变量之间因果关系的过程。它让我们超越预测，去识别事件发生的机制和路径。

想象一下，我们发现发动机故障与地理位置强相关。因果推断能帮助我们找到背后的驱动变量——也许是该地区特定的空气质量或操作规范，而地理位置本身只是一个乘客变量。

评测环境与数据

为了公平起见，我们使用同一个数据集进行评测：美国人口普查收入数据集。我们将围绕一个核心问题展开分析：“拥有研究生学历是否会增加年收入超过5万美元的概率？”

以下是数据加载和清洗的代码，为后续所有分析打下基础：

# 安装并导入必要库
pip install datazets
import datazets as dz 
import pandas as pd 
import numpy as np
import matplotlib.pyplot as plt

# 加载数据集并清洗
df = dz.import_example(data='census_income')
# 剔除连续型变量和敏感特征
drop_cols = ['age', 'fnlwgt', 'education-num', 'capital-gain', 'capital-loss', 'hours-per-week', 'race', 'sex']
df.drop(labels=drop_cols, axis=1, inplace=True)

# 查看数据
print(df.head())

1. Bnlearn

一站式的贝叶斯网络工具箱

一句话点评： 功能全面、开箱即用，是初学者和专家的理想选择。

Bnlearn 提供了一套完整的工具，用于创建、分析和可视化贝叶斯网络。它支持离散、连续和混合数据集，封装了从结构学习、参数估计到推理验证的核心流程。

bnlearn库提供了一系列用于在 Python 中创建、分析和可视化贝叶斯网络的工具。它支持离散、混合和连续数据集，设计简洁易用，同时包含了因果学习中最基本的贝叶斯流程。借助 bnlearn，可以执行结构学习、参数估计、推理，并创建合成数据，使其既适用于探索性分析，也适用于高级因果发现。只需在初始化期间指定参数，即可应用一系列统计检验。此外，它还提供了一些辅助函数，用于转换数据集、推导整个图的拓扑顺序、比较两个网络以及生成各种富有洞察力的图表。

核心功能：

• 结构学习：从数据中学习网络结构或整合专家知识。
• 参数学习：根据观测数据估计条件概率分布。
• 因果推断：查询网络以计算概率、模拟干预措施并检验因果效应。
• 合成数据：从学习到的贝叶斯网络生成新的数据集，用于模拟、基准测试或数据增强。
• 离散化数据：使用内置的离散化方法将连续变量转换为离散状态，以便更轻松地进行建模和推理。
• 模型评估：使用评分函数、统计检验和性能指标来比较网络。
• 可视化：交互式图形绘制、概率热图和结构叠加，直观展示因果关系。

实战演示：

我们让 Bnlearn 无监督地学习数据中的因果结构，不预设目标变量。

Bnlearn for Python的一大优势在于它能够以无监督的方式学习因果结构。这意味着你无需指定任何目标变量，它就能自行找出答案。为此，它实现了多种搜索和评分策略：

• 搜索策略：爬山搜索、穷举搜索、约束搜索、Chow-Liu 算法、朴素贝叶斯算法和 TAN 算法
• 评分策略： BIC、K2、BDEU。

某些策略需要设置根节点，例如树增强朴素贝叶斯（TAN），当已知结果（或目标）值时建议采用此策略。这将显著降低计算负担，在特征数量庞大的情况下尤为推荐。此外，通过独立性检验可轻松从模型中剪枝掉虚假边。在下例中，我将采用基于BIC评分类型的爬坡搜索法进行结构学习。本例中不预先定义目标值，而是让Bnlearn完全基于数据本身决定因果结构。

import bnlearn as bn

# 结构学习
model = bn.structure_learning.fit(df, methodtype='hillclimbsearch', scoretype='bic')
# 独立性检验，修剪不显著的边
model = bn.independence_test(model, df, test="chi_square", alpha=0.05, prune=True)
# 参数学习
model = bn.parameter_learning.fit(model, df)
# 绘制交互式网络图
bn.plot(model, interactive=True) #生成的因果网络图会清晰展示`教育`、`职业`、`婚姻状况`等变量如何影响`薪资`
print(model['model_edges'])
"""
[('education', 'occupation'),
 ('marital-status', 'relationship'),
 ('occupation', 'workclass'),
 ('relationship', 'salary'),
 ('relationship', 'education'),
 ('salary', 'education')]
"""

要确定有向无环图（DAG），我们需要按上文代码段所示指定输入数据框。模型拟合完成后，结果将存储在模型字典中，可用于后续分析。下图展示了因果结构的静态图和交互式图。

通过学习到的有向无环图（图1），我们可以估计条件概率分布（CPDs，参见下文代码部分），并使用do-calculus进行推断。换言之，我们可以开始向数据提出问题。

# 使用估计的边学习CPD。
# 注意我们也可以自定义边或手动提供有向无环图。
# model = bn.make_DAG(model['model_edges'])

# 通过提供模型和数据框学习CPD
model = bn.parameter_learning.fit(model, df)

# 打印CPD
CPD = bn.print_CPD(model)

开始提问：

问题一：拥有博士学位的人，薪资超过5万的概率有多大？

直观上，我们可能预期概率较高，因为受教育程度为博士。让我们通过贝叶斯模型推导后验概率。在下面的代码段中，我们得出概率P=0.7093。这证实了当受教育程度为博士时，获得>50K年薪的概率确实高于非博士学历者。

query = bn.inference.fit(model, variables=['salary'], evidence={'education':'Doctorate'})
print(query)

"""
+---------------+---------------+
| salary        |   phi(salary) |
+===============+===============+
| salary(<=50K) |        0.2907 |
+---------------+---------------+
| salary(>50K)  |        0.7093 |
+---------------+---------------+

Summary for variables: ['salary']
Given evidence: education=Doctorate

salary outcomes:
- salary: <=50K (29.1%)
- salary: >50K (70.9%)
"""

结果：P(>50K | 博士) = 70.9%。果然，高学历带来高收入！

问题二：那么只有高中学历的人呢？

由此得出后验概率P=0.1615。根据该数据集，学习对获得更高薪资具有显著益处。但需注意，我们未采用任何其他可能影响结果的约束条件作为证据。

query = bn.inference.fit(model, variables=['salary'], evidence={'education':'HS-grad'})
print(query)

"""
+---------------+---------------+
| salary        |   phi(salary) |
+===============+===============+
| salary(<=50K) |        0.8385 |
+---------------+---------------+
| salary(>50K)  |        0.1615 |
+---------------+---------------+

Summary for variables: ['salary']
Given evidence: education=HS-grad

salary outcomes:
- salary: <=50K (83.8%)
- salary: >50K (16.2%)
"""

结果：P(>50K | 高中) = 16.2%。对比非常鲜明。

总结：

• 优势： pipeline 完整，易于上手，可视化效果佳，支持多种数据类型。
• 输入数据： 离散、连续、混合数据均可。

2. Pgmpy

底层构建模块的乐高大师

一句话点评： 灵活度高，但需要你亲手搭建一切，适合资深玩家。

Pgmpy 提供了概率图模型的底层构建块。它的优势在于灵活性，你可以自由组合各种学习算法和推理方法。但这也意味着你需要更多的专业知识来搭建完整流程。

实战演示：

实现与 Bnlearn 类似的结构学习和推理，但代码量明显更多。

!pip install pgmpy
from pgmpy.estimators import HillClimbSearch, BicScore, BayesianEstimator
from pgmpy.models import BayesianNetwork
from pgmpy.inference import VariableElimination

# 结构学习
est = HillClimbSearch(df)
scoring_method = BicScore(df)
model = est.estimate(scoring_method=scoring_method)
print("发现的边：", model.edges())

# ... (需要手动构建模型、拟合参数，过程更复杂)
# [('education', 'salary'),
# ('marital-status', 'relationship'),
# ('occupation', 'workclass'),
# ('occupation', 'education'),
# ('relationship', 'salary'),
# ('relationship', 'occupation')]
# 推理查询

要使用do-calculus对数据进行推断，我们首先需要估计条件概率分布（CPDs），如下面的代码段所示。

vec = {
'source': ['education', 'marital-status', 'occupation', 'relationship', 'relationship', 'salary'],
'target': ['occupation', 'relationship', 'workclass', 'education', 'salary', 'education'],
'weight': [True, True, True, True, True, True]
}
vec = pd.DataFrame(vec)

# Create Bayesian model
bayesianmodel = BayesianNetwork(vec)

# Fit the model
bayesianmodel.fit(df, estimator=BayesianEstimator, prior_type='bdeu', equivalent_sample_size=1000)

# Create model for variable elimination
model_infer = VariableElimination(bayesianmodel)

# Query
query = model_infer.query(variables=['salary'], evidence={'education':'Doctorate'})
print(query) # 同样得到 P=0.709

"""
+---------------+---------------+
| salary        |   phi(salary) |
+===============+===============+
| salary(<=50K) |        0.2907 |
+---------------+---------------+
| salary(>50K)  |        0.7093 |
+---------------+---------------+
"""

总结：

• 优势： 极其灵活，适合自定义复杂模型和研究算法。
• 劣势： 学习曲线陡峭，需要手动处理许多步骤。
• 输入数据： 必须是离散数据。

3. CausalNex

基于NOTEARS的结构探索者

一句话点评： 功能强大，但环境配置和数据预处理稍显繁琐。

CausalNex 基于著名的 NOTEARS 算法，擅长从数据中学习因果结构。但它仅支持离散数据，且对类别变量多的数据拟合效果不佳，需要额外的预处理来降低基数。

实战演示：

!pip install causalnex
from causalnex.structure.notears import from_pandas
from causalnex.network import BayesianNetwork
import networkx as nx

# 数据必须转换为数值型
from sklearn.preprocessing import LabelEncoder
le = LabelEncoder()
df_num = df.apply(le.fit_transform)

# 结构学习与阈值处理
sm = from_pandas(df_num)
sm.remove_edges_below_threshold(0.8) # 过滤弱边

# 可视化
nx.draw_networkx(sm, ...)

总结：

• 优势： 集成了先进的NOTEARS算法。
• 劣势： 预处理复杂，对Python和新版库的兼容性较差。
• 输入数据： 必须是数值型离散数据。

4. DoWhy

专注于因果效应估计的专家

一句话点评： 不关心完整网络，只专注于回答“这个干预有多大效果？”

DoWhy 的思维方式与前几个库不同。它要求你明确定义处理变量和结果变量，然后专注于估计处理变量对结果的平均因果效应。它不擅长从数据中学习因果结构，但非常擅长在给定因果图的情况下进行稳健的效应估计。

实战演示：

我们研究“拥有博士学位”对“高薪”的因果效应。

！pip install dowhy
from dowhy import CausalModel
import dowhy.datasets
import datazets as dz
from sklearn.preprocessing import LabelEncoder
import numpy as np
le = LabelEncoder()
# 导入数据集并删除连续型和敏感特征
df = dz.get(data='census_income')

drop_cols = ['age', 'fnlwgt', 'education-num', 'capital-gain', 'capital-loss', 'hours-per-week', 'race', 'sex']
df.drop(labels=drop_cols, axis=1, inplace=True)

# 处理变量必须为二元
df['education'] = df['education']=='Doctorate'

# 接下来需将数据转换为数值型
df_num = df.copy()
for col in df_num.columns:
    df_num[col] = le.fit_transform(df_num[col])

# 指定处理变量、结果变量及潜在混杂变量
treatment = “education”
outcome = “salary”

# 步骤1. 创建因果图
model = CausalModel(
        data=df_num,
        treatment=treatment,
        outcome=outcome,
        common_causes=list(df.columns[~np.isin(df.columns, [treatment, outcome])]),
        graph_builder='ges',
        alpha=0.05,
        )

# 展示模型
model.view_model()

如上文代码部分所见，处理变量必须为二元变量（设为“博士学位”），所有分类变量均需编码为数值。在下文代码部分，我们将利用因果图的特性来识别因果效应。结果或许不足为奇：拥有博士学位确实能提高获得优厚薪资的概率。

# 步骤2：识别因果效应并返回目标估计量
identified_estimand = model.identify_effect(proceed_when_unidentifiable=True)

# 结果
print(identified_estimand)
"""
[16-09-2025 21:06:21] [dowhy.causal_identifier.auto_identifier] [INFO] Causal effect can be identified.
[16-09-2025 21:06:21] [dowhy.causal_identifier.auto_identifier] [INFO] Instrumental variables for treatment and outcome:[]
[16-09-2025 21:06:21] [dowhy.causal_identifier.auto_identifier] [INFO] Frontdoor variables for treatment and outcome:[]
[16-09-2025 21:06:21] [dowhy.causal_identifier.auto_identifier] [INFO] Number of general adjustment sets found: 1
[16-09-2025 21:06:21] [dowhy.causal_identifier.auto_identifier] [INFO] Causal effect can be identified.
"""

print(identified_estimand)

"""
Estimand type: EstimandType.NONPARAMETRIC_ATE

### Estimand : 1
Estimand name: backdoor
Estimand expression:
     d                                                                         ↪
────────────(E[salary|native-country,occupation,marital-status,workclass,relat ↪
d[education]                                                                   ↪

↪          
↪ ionship])
↪          
Estimand assumption 1, Unconfoundedness: If U→{education} and U→salary then P(salary|education,native-country,occupation,marital-status,workclass,relationship,U) = P(salary|education,native-country,occupation,marital-status,workclass,relationship)

### Estimand : 2
Estimand name: iv
No such variable(s) found!

### Estimand : 3
Estimand name: frontdoor
No such variable(s) found!

### Estimand : 4
Estimand name: general_adjustment
Estimand expression:
     d                                                                         ↪
────────────(E[salary|native-country,marital-status,occupation,workclass,relat ↪
d[education]                                                                   ↪

↪          
↪ ionship])
↪          
Estimand assumption 1, Unconfoundedness: If U→{education} and U→salary then P(salary|education,native-country,marital-status,occupation,workclass,relationship,U) = P(salary|education,native-country,marital-status,occupation,workclass,relationship)
"""

# 步骤3：使用统计方法估计目标估计量。
estimate = model.estimate_effect(identified_estimand, method_name="backdoor.propensity_score_stratification"

# 结果
print(estimate)
"""
*** Causal Estimate ***

## Identified estimand
Estimand type: EstimandType.NONPARAMETRIC_ATE

### Estimand : 1
Estimand name: backdoor
Estimand expression:
     d                                                                         ↪
────────────(E[salary|native-country,occupation,marital-status,workclass,relat ↪
d[education]                                                                   ↪

↪          
↪ ionship])
↪          
Estimand assumption 1, Unconfoundedness: If U→{education} and U→salary then P(salary|education,native-country,occupation,marital-status,workclass,relationship,U) = P(salary|education,native-country,occupation,marital-status,workclass,relationship)

## Realized estimand
b: salary~education+native-country+occupation+marital-status+workclass+relationship
Target units: ate

## Estimate
Mean value: 0.46973242081553496
## Estimate
# Mean value: 0.4697157228651772

"""
# 步骤4：通过多重稳健性检验来反驳所得估计值。
refute_results = model.refute_estimate(identified_estimand, estimate, method_name="random_common_cause")

总结：

• 优势： 因果效应估计的框架非常严谨，内置多种鲁棒性检验。
• 劣势： 无法学习因果图，输出结果较复杂，需要统计知识解读。
• 输入数据： 需要明确的处理变量和结果变量。

5. PyAgrum

多才多艺的图模型学者

一句话点评： 来自学术界的全能选手，功能丰富但文档相对学院派。

PyAgrum 是一个功能强大的概率图模型库，不仅支持贝叶斯网络，还支持马尔可夫网络等。它提供了多种学习算法和推理工具，但安装和入门可能对新手有些挑战。

实战演示：

# 安装 pyagrum
pip install pyagrum

# 安装 graphviz（可视化所需）
pip install setgraphviz
import datazets as dz
import pandas as pd
import pyagrum as gum
import pyagrum.lib.notebook as gnb
import pyagrum.lib.explain as explain
import pyagrum.lib.bn_vs_bn as bnvsbn

# 导入库以显示点图
from setgraphviz import setgraphviz
# 设置路径
setgraphviz()

# 导入数据集并剔除连续型和敏感特征
df = dz.get(data='census_income')
# 数据清洗
drop_cols = ['age', 'fnlwgt', 'education-num', 'capital-gain', 'capital-loss', 'hours-per-week', 'race', 'sex']
df.drop(labels=drop_cols, axis=1, inplace=True)
# 删除存在任何缺失值的行
df2 = df.dropna().copy()
# 将所有缺失值替换为占位符字符串
df2 = df2.fillna(“missing”).replace(“?”, “missing”)

# 确保分类列为类别型数据类型（pyAgrum要求离散变量）
for col in df2.columns:
    df2[col] = df2[col].astype(“category”)

# 从pandas数据框创建学习器
learner = gum.BNLearner(df2)

# 配置评分与搜索策略
learner.useScoreBIC()
learner.useGreedyHillClimbing()

# 训练网络
bn = learner.learnBN()

# 学习参数
bn2 = learner.learnParameters(bn.dag())

# 绘制网络图
gnb.showBN(bn2)

总结：

• 优势： 功能全面，模型种类多。
• 劣势： 对数据预处理要求严格，可视化依赖Graphviz，社区资源相对较少。
• 输入数据： 必须是完整的离散数据集。

6. CausalImpact

时间序列因果推断的利器

一句话点评： 专为时间序列设计，评估政策、活动等干预的短期影响。

CausalImpact 与其他库完全不同。它使用贝叶斯结构时间序列模型，来评估一次干预（如新政策上线、营销活动）对时间序列数据的影响。它的核心是比较干预发生后的实际值与“假设未发生干预”的预测值之间的差异。

实战演示：

假设我们有一段时间内的网站流量数据，并在第70天进行了一次改版。

from causalimpact import CausalImpact
import pandas as pd

# 模拟数据：y是流量，x1是控制变量，第71天后y有一个提升
data = pd.DataFrame({'y': y_data, 'x1': x1_data})

# 定义干预前后时间段
impact = CausalImpact(data, pre_period=[0, 69], post_period=[70, 99])
impact.run()
impact.plot() # 生成包含三个面板的详细效果图
impact.summary()
#                              Average    Cumulative
# Actual                           130          3773
# Predicted                        120          3501
# 95% CI                    [118, 122]  [3447, 3555]

# Absolute Effect                    9           272
# 95% CI                       [11, 7]    [326, 218]

# Relative Effect                 7.8%          7.8%
# 95% CI                  [9.3%, 6.2%]  [9.3%, 6.2%]

# P-value                         0.0%              
# Prob. of Causal Effect        100.0%

# Summary report
impact.summary(output="report")

平均值列描述了干预后时段的平均时间。累计值列汇总了各个时间点的数据，当响应变量代表流量指标（如查询量、点击量、访问量、安装量、销售额或收入）时，这种视角尤为有用。本例中可见因果效应概率为100%，P值为0——这符合预期，因为数据中已包含该信息。

在图4中，我们可以看到三个子图，其描述如下：在子图1，我们既能看到原始数据，也能看到针对治疗后时期的反事实预测。子图2展示了观测数据与反事实预测之间的差异，该差异即模型所确定的点因果效应估计值。子图3通过聚合前一面板的点效应贡献，描绘了干预措施的累积效应。需特别注意的是，这些推论的有效性依赖于协变量未受干预本身影响的假设。此外，模型还假定协变量与处理时间序列在前期建立的关系，在后期始终保持一致。

总结：

• 优势： 解决时间序列因果问题的专用工具，结果解释直观。
• 劣势： 仅适用于时间序列数据，不适用于构建通用因果图。
• 输入数据： 时间序列数据，需明确干预点。

写在最后

如何选择？看这里：

1. 如果你想快速入门，从数据中自动发现因果关系并进行推理，Bnlearn 是你的不二之选。
2. 如果你是高级用户，希望完全控制建模的每一个细节，请选择 Pgmpy。
3. 如果你的核心任务是评估某个处理变量对结果变量的因果效应，DoWhy 提供了最严谨的框架。
4. 如果你的数据是时间序列，并且想评估一次干预的因果影响，CausalImpact 是唯一答案。

希望这篇详尽的对比能帮助你在因果推断的旅程中，找到那盏最明亮的指路明灯！

互动时间：你在工作中遇到过哪些因果推断的挑战？或者用过哪些好用的因果推断库？欢迎在评论区留言分享！

（关注我们，获取更多数据科学干货！）