本研究介绍了贝叶斯网络个体化风险评估模型的构建方法，使用实例与传统基于回归的Logistic模型进行对比，对该模型的性能进行评价，并使用R语言进行演示说明，为研究者理解与应用贝叶斯网络模型提供参考。
随着大量的医疗记录转变为电子格式，医疗领域积累了海量的临床数据，为更加精准地评估患者疾病风险及选择干预措施提供了机会。而从患者的角度出发，个体化风险评估具有广泛的应用前景，可以帮助临床医生有针对性的为每个患者提供个性化的预防和治疗方案，降低疾病的发生和死亡率。传统的预测评估模型通常基于患者的临床或非临床特征，用于预测其在特定时间内发生特定健康状况的风险和可能性[1]。然而，受限于所纳入变量数[2]、大规模数据处理能力、模型复杂度以及难以捕捉预后因素之间的非线性关系[3]等原因难以满足临床需求。因此，机器学习方法受到越来越多的关注，该类方法在处理复杂数据上的优势已使其广泛应用于医疗卫生领域[4,5]。其中贝叶斯网络（Bayesian networks，BN）模型通过构建有向无环图（direct acyclic graph，DAG）和条件概率表，以其擅长阐明变量之间的潜在因果关系等优势[6,7]，为个体化风险评估研究提供了新途径[8,9]。本研究旨在介绍BN的基本原理和模型构建步骤，并基于实例采用R语言实现个体化风险预测模型的构建，为帮助研究者更好理解、应用该方法提供参考。
1 BN基本原理
BN又称信念网络（belief network），由Judea Pear于1985年首次提出[10]，是一种利用DAG表示随机变量之间依赖关系的概率图形模型，作为一种图形化的建模工具，BN提供了一种表示变量之间因果关系的方法，用来发现隐藏在数据中的信息，BN将DAG与概率论有机结合，在不确定推理方面发挥了很大的优势。
BN由两部分组成，DAG和条件概率分布。DAG定性描述了变量之间的依赖和独立关系，而条件概率分布则对变量与其父节点之间的依赖关系进行了定量描述，是联合概率分布（joint probability distributions，JPD）的一种表现形式。JPD包含了所有的概率信息[11]。利用链式法则，JPD可以表示为：
 |
条件独立性是BN的核心。DAG中的节点为随机变量，表示为
，有向边表示变量之间的概率依赖关系，箭头的指向是对变量之间的定性描述。当观察到有向边从
指向
时，说明
对
存在影响，在
和
中，
是父节点，而
是子节点。每个节点都有一个条件概率分布表用来定量描述所有父节点对该节点的作用效果。假设
是
所有父节点的集合，即：
 |
 |
1.1 BN模型与传统回归模型优缺点比较
与传统基于回归的模型相比，BN模型是对JPD简洁直观的图形表示，将与节点相关的条件概率分布分解为直接影响节点的变量，采用网络的形式表示变量及其之间的复杂的相互作用关系，存在许多优势，但仍存在一些不足[14]。具体比较内容见表1。
表1
BN模型与传统回归模型优缺点比较
1.2 BN模型的学习
BN的学习过程主要分为结构学习和参数学习[15]。
结构学习是通过训练集学习到相应的DAG结构，主要采用基于约束和基于分数的方法[11]。基于约束的方法利用数据中的条件独立性推断模型结构，例如Grow-Shrink算法。而基于评分的方法则引入评分函数评估每个候选模型，以找到最适合数据的模型结构，如禁忌搜索（tabu search）算法。另外，常用的方法还包括爬山算法（hill climbing algorithm），这是一种局部搜索算法，从空白网络结构开始，通过添加、删除或反转变量之间的边来搜索最优网络结构，直到达到最优评分指标。同时，BN允许研究者将已知的时间关系、直接关系和直接因果关系纳入成为结构的一部分。
参数学习旨在确定BN网络模型各节点的条件概率分布，从而量化网络中变量之间的依赖性。常见的学习方法包括最大似然估计（maximum likelihood estimation，MLE）方法和贝叶斯估计（Bayesian estimation，BE）方法等。MLE方法通过最大化观测数据的似然函数来估计参数值，而BE方法则基于贝叶斯定理和先验概率来计算后验概率分布，从而得到参数的估计值。
2 BN个体化风险评估的R语言实现
2.1 数据来源
本研究采用来自Kaggle的卒中预测数据集[16]进行演示。该数据集的目标是根据患者临床指标分析预测患者是否会发生卒中。数据集包含5 110个样本，每个样本有11个变量和一个结果变量。具体变量介绍见表2。
表2
Kaggle卒中预测数据变量及内容
2.2 BN个体化风险评估建模思路
变量的选择关系到BN模型的质量，本研究使用单因素分析方法与多因素Logistic回归方法对变量进行筛选。数据按照7∶3的比例随机分为训练集和测试集，使用训练集构建BN模型，使用测试集对模型进行验证。为解决数据不平衡与防止数据泄露的问题，分别对训练集和测试集使用了合成少数过采样技术（synthetic minority oversampling technique，SMOTE）[17]。该方法不仅是对少数样本进行重复，而是通过在附近实例之间插值，在算法上生成新的少数类实例，从而增强数据集平衡[18,19]。本研究以爬山算法构建BN结构，MLE方法进行参数学习。使用受试者工作特征（receiver operating characteristic，ROC）曲线下面积（area under the curve，AUC）、灵敏度、特异度和准确度对模型进行评价。同时，使用Delong检验比较BN模型与Logistic模型之间性能是否存在差异。所有统计分析使用R 4.2.2软件进行，以双侧检验P<0.05为差异具有统计学意义。实现的R语言代码见框1。
2.3 主要结果
2.3.1 缺血性脑卒中患者影响因素分析
数据集共纳入5 110例患者，其中男性患者2 994例（58.59%），女性患者2 116例（41.41%）。对患者进行随访，IS患者249例（4.87%），非IS患者4 861例（95.13%）。Logistic回归结果显示年龄、高血压史、心脏病史及平均血糖与缺血性脑卒中患者预后存在统计学意义（表3）。
表3
缺血性脑卒中发病影响因素的多因素Logistic回归分析
2.3.2 BN模型构建
R软件中，常使用bnlearn包进行学习BN图形结构，并对参数进行估计。操作过程再在R Studio中进行。数据包加载完成后，对研究所需数据进行录入。
本研究构建离散型BN模型，将连续型变量平均血糖分为正常（<3.9 mmol/L）和高血糖（≥3.9 mmol/L）。研究使用了爬山算法构建DAG。在此基础上，结合相关文献和专家经验，对初始网络结构进行了人工调整，包括删除某些不合理的有向边、增加必要的有向边以及改变某些有向边的方向。通常使用白名单（whitelist）和黑名单（blacklist）来进一步完善BN的结构。白名单指定了所需的有向边存在，而黑名单指定了有向边不存在。使用graphviz.plot()函数绘制所构建的BN结构。
构建的BN模型结构如图1所示。年龄、心脏病、高血压和高血糖与IS直接相关。此外，年龄还可通过心脏病、高血压和高血糖与IS间接相关，高血糖也可通过影响心脏病和高血压，间接与IS相关。
图1
缺血性脑卒中发病的贝叶斯网络模型结构
Age：年龄；AG：平均血糖；HBP：高血压；Heart disease：心脏病；IS：缺血性脑卒中。
2.3.3 输出条件概率表
构建BN结构后，需要对模型拟合参数。当所有节点都是离散变量时，这些参数被列在通常被称为条件概率表的表格中[11]。本研究使用MLE对局部分布的参数进行拟合，使其形成条件概率表。bn.fit()函数默认使用MLE的方法，也可以根据method对方法进行更换，例如可以更换为method="bayes"使用BN。
fitted$IS可以输出每个变量的条件概率表，在这里我们的目标节点是缺血性脑卒中。条件概率表如表4所示。60岁以上70岁以下，有高血压、心脏病和高血糖的患者，缺血性脑卒中发病概率为0.870。而70岁以上，有高血压和高血糖的患者，缺血性脑卒中发病概率为0.764。
表4
缺血性脑卒中发病BN预测模型条件概率分布表
2.3.4 BN交叉验证
研究使用bnlearn包中的bn.cv()函数进行了10折交叉验证。该函数提供了三种交叉验证方法，包括k折交叉验证（默认）、自定义折叠交叉验证（"custom-fold"）和Hold-out交叉验证（"hold-out"），通过method对使用的交叉验证方法进行修改。我们使用对数似然损失（log-likelihood loss）对BN模型性能进行评价，一般而言，期望损失值越低，表示模型在训练集上的性能越好。本研究得到训练集的期望损失为2.858，表明该模型在训练集上表现良好。
2.3.5 模型效果评价
使用pROC包构建BN模型与Logistic回归模型的ROC曲线，并分别计算AUC值、灵敏度、特异度和准确度，对模型效果进行评价（表5、图2）。
表5
BN模型与Logistic回归模型性能评估
图2
贝叶斯网络模型与Logistic回归模型的ROC曲线
与Logistic回归模型相比，在测试集和训练集中BN模型均表现出更好的性能。在训练集和测试集中，BN模型的准确率分别为79.37%和83.32%，AUC分别为0.85［95%CI（0.83，0.86）］和0.88［95%CI（0.86，0.89）］，灵敏度分别为59.64%和72.80%，特异度分别为87.48%和87.59%。训练集与测试集的Delong检验证实了两个模型之间AUC之间存在统计学差异（P<0.05）。此外，BN模型的校准图显示，在训练集和测试集中预测概率和实际观测值之间较为一致，见图3。
图3
贝叶斯网络模型校准曲线
3 讨论
本研究基于实例采用R语言实现个体化风险预测模型的构建，研究结果显示，BN模型的预测性能优于传统Logistic回归模型。且通过DAG图可以清晰的展示出变量之间的相互作用关系，年龄、心脏病、高血压和高血糖与缺血性脑卒中发病直接相关。此外，年龄还可通过心脏病、高血压和高血糖与缺血性脑卒中间接相关，高血糖也可通过影响心脏病和高血压，间接与缺血性脑卒中相关。
基于BN的风险预测模型相较于统计回归的方法具有多个优势[20-22]。首先，BN模型通过其图形结构和条件概率分布清晰地阐明了各变量之间复杂的相互作用关系，以及变量在网络中的相对位置，描述了变量对结果的贡献程度[7,23]。这种清晰的关系表达提高了对结果风险的评估准确性，并增强了结果的可解释性[14]。其次，BN模型可以充分利用先验知识和现有数据，有效地评估风险，并且在加入新数据时无需重新训练模型[8,24]，这使得其在个体化风险预测方面更加灵活和实用。随着时间的推移，更多的病历信息被纳入，BN模型可以获得进一步的提升。此外，BN模型能够进行从原因到结果和从结果到原因的推理，以及处理变量间的线性和非线性关系，因而具有更广泛的适用性[14]。
相比之下，基于回归的方法局限于单一结果的预测，在加入新数据时需要重新估计，且在面对数据缺失或个体特异性较强的情况下，基于回归的模型复杂程度较高且准确性有限。此外，基于回归的方法在纳入变量和结果之间的非线性关系时，通常涉及变量的转换、加权等操作，往往使得结果难以解释[14]。因此，BN模型在风险预测中展现出了显著的优势，并在实践中得到了广泛应用。
在临床中，疾病的进展往往与时间密切相关，BN能够与时间信息相结合，构成动态BN[25]，从而更好地反映变量的发展变化规律。同时，BN也可以与多种机器学习技术和传统统计方法结合使用，更好地处理数据，如与深度学习结合的贝叶斯深度学习模型[26]，这为解决复杂问题提供了新的途径。
BN为处理复杂的医疗卫生数据提供了一种强大且灵活的分析方法。数据规模的扩大、数据复杂性的增加等挑战的不断出现，传统的数据分析方法面临诸多限制。而BN能够应对这些挑战，不断更新信息。随着“精准医学”时代的到来，新的风险评估与决策制定的方法逐渐引起人们的关注，而BN凭借其高度的灵活性、适应性以及强大的可解释性，为大数据的分析提供了一种新的思路，值得进一步探索。本文提供了清晰的代码和代码解释，以实例的形式进行演示，方便研究人员更好地理解和应用BN模型。
随着大量的医疗记录转变为电子格式，医疗领域积累了海量的临床数据，为更加精准地评估患者疾病风险及选择干预措施提供了机会。而从患者的角度出发，个体化风险评估具有广泛的应用前景，可以帮助临床医生有针对性的为每个患者提供个性化的预防和治疗方案，降低疾病的发生和死亡率。传统的预测评估模型通常基于患者的临床或非临床特征，用于预测其在特定时间内发生特定健康状况的风险和可能性[1]。然而，受限于所纳入变量数[2]、大规模数据处理能力、模型复杂度以及难以捕捉预后因素之间的非线性关系[3]等原因难以满足临床需求。因此，机器学习方法受到越来越多的关注，该类方法在处理复杂数据上的优势已使其广泛应用于医疗卫生领域[4,5]。其中贝叶斯网络（Bayesian networks，BN）模型通过构建有向无环图（direct acyclic graph，DAG）和条件概率表，以其擅长阐明变量之间的潜在因果关系等优势[6,7]，为个体化风险评估研究提供了新途径[8,9]。本研究旨在介绍BN的基本原理和模型构建步骤，并基于实例采用R语言实现个体化风险预测模型的构建，为帮助研究者更好理解、应用该方法提供参考。
1 BN基本原理
BN又称信念网络（belief network），由Judea Pear于1985年首次提出[10]，是一种利用DAG表示随机变量之间依赖关系的概率图形模型，作为一种图形化的建模工具，BN提供了一种表示变量之间因果关系的方法，用来发现隐藏在数据中的信息，BN将DAG与概率论有机结合，在不确定推理方面发挥了很大的优势。
BN由两部分组成，DAG和条件概率分布。DAG定性描述了变量之间的依赖和独立关系，而条件概率分布则对变量与其父节点之间的依赖关系进行了定量描述，是联合概率分布（joint probability distributions，JPD）的一种表现形式。JPD包含了所有的概率信息[11]。利用链式法则，JPD可以表示为：
|
条件独立性是BN的核心。DAG中的节点为随机变量，表示为，有向边表示变量之间的概率依赖关系，箭头的指向是对变量之间的定性描述。当观察到有向边从指向时，说明对存在影响，在和中，是父节点，而是子节点。每个节点都有一个条件概率分布表用来定量描述所有父节点对该节点的作用效果。假设是所有父节点的集合，即：
|
|
1.1 BN模型与传统回归模型优缺点比较
与传统基于回归的模型相比，BN模型是对JPD简洁直观的图形表示，将与节点相关的条件概率分布分解为直接影响节点的变量，采用网络的形式表示变量及其之间的复杂的相互作用关系，存在许多优势，但仍存在一些不足[14]。具体比较内容见表1。
表1
BN模型与传统回归模型优缺点比较
1.2 BN模型的学习
BN的学习过程主要分为结构学习和参数学习[15]。
结构学习是通过训练集学习到相应的DAG结构，主要采用基于约束和基于分数的方法[11]。基于约束的方法利用数据中的条件独立性推断模型结构，例如Grow-Shrink算法。而基于评分的方法则引入评分函数评估每个候选模型，以找到最适合数据的模型结构，如禁忌搜索（tabu search）算法。另外，常用的方法还包括爬山算法（hill climbing algorithm），这是一种局部搜索算法，从空白网络结构开始，通过添加、删除或反转变量之间的边来搜索最优网络结构，直到达到最优评分指标。同时，BN允许研究者将已知的时间关系、直接关系和直接因果关系纳入成为结构的一部分。
参数学习旨在确定BN网络模型各节点的条件概率分布，从而量化网络中变量之间的依赖性。常见的学习方法包括最大似然估计（maximum likelihood estimation，MLE）方法和贝叶斯估计（Bayesian estimation，BE）方法等。MLE方法通过最大化观测数据的似然函数来估计参数值，而BE方法则基于贝叶斯定理和先验概率来计算后验概率分布，从而得到参数的估计值。
2 BN个体化风险评估的R语言实现
2.1 数据来源
本研究采用来自Kaggle的卒中预测数据集[16]进行演示。该数据集的目标是根据患者临床指标分析预测患者是否会发生卒中。数据集包含5 110个样本，每个样本有11个变量和一个结果变量。具体变量介绍见表2。
表2
Kaggle卒中预测数据变量及内容
2.2 BN个体化风险评估建模思路
变量的选择关系到BN模型的质量，本研究使用单因素分析方法与多因素Logistic回归方法对变量进行筛选。数据按照7∶3的比例随机分为训练集和测试集，使用训练集构建BN模型，使用测试集对模型进行验证。为解决数据不平衡与防止数据泄露的问题，分别对训练集和测试集使用了合成少数过采样技术（synthetic minority oversampling technique，SMOTE）[17]。该方法不仅是对少数样本进行重复，而是通过在附近实例之间插值，在算法上生成新的少数类实例，从而增强数据集平衡[18,19]。本研究以爬山算法构建BN结构，MLE方法进行参数学习。使用受试者工作特征（receiver operating characteristic，ROC）曲线下面积（area under the curve，AUC）、灵敏度、特异度和准确度对模型进行评价。同时，使用Delong检验比较BN模型与Logistic模型之间性能是否存在差异。所有统计分析使用R 4.2.2软件进行，以双侧检验P<0.05为差异具有统计学意义。实现的R语言代码见框1。
2.3 主要结果
2.3.1 缺血性脑卒中患者影响因素分析
数据集共纳入5 110例患者，其中男性患者2 994例（58.59%），女性患者2 116例（41.41%）。对患者进行随访，IS患者249例（4.87%），非IS患者4 861例（95.13%）。Logistic回归结果显示年龄、高血压史、心脏病史及平均血糖与缺血性脑卒中患者预后存在统计学意义（表3）。
表3
缺血性脑卒中发病影响因素的多因素Logistic回归分析
2.3.2 BN模型构建
R软件中，常使用bnlearn包进行学习BN图形结构，并对参数进行估计。操作过程再在R Studio中进行。数据包加载完成后，对研究所需数据进行录入。
本研究构建离散型BN模型，将连续型变量平均血糖分为正常（<3.9 mmol/L）和高血糖（≥3.9 mmol/L）。研究使用了爬山算法构建DAG。在此基础上，结合相关文献和专家经验，对初始网络结构进行了人工调整，包括删除某些不合理的有向边、增加必要的有向边以及改变某些有向边的方向。通常使用白名单（whitelist）和黑名单（blacklist）来进一步完善BN的结构。白名单指定了所需的有向边存在，而黑名单指定了有向边不存在。使用graphviz.plot()函数绘制所构建的BN结构。
构建的BN模型结构如图1所示。年龄、心脏病、高血压和高血糖与IS直接相关。此外，年龄还可通过心脏病、高血压和高血糖与IS间接相关，高血糖也可通过影响心脏病和高血压，间接与IS相关。
图1
缺血性脑卒中发病的贝叶斯网络模型结构
Age：年龄；AG：平均血糖；HBP：高血压；Heart disease：心脏病；IS：缺血性脑卒中。
2.3.3 输出条件概率表
构建BN结构后，需要对模型拟合参数。当所有节点都是离散变量时，这些参数被列在通常被称为条件概率表的表格中[11]。本研究使用MLE对局部分布的参数进行拟合，使其形成条件概率表。bn.fit()函数默认使用MLE的方法，也可以根据method对方法进行更换，例如可以更换为method="bayes"使用BN。
fitted$IS可以输出每个变量的条件概率表，在这里我们的目标节点是缺血性脑卒中。条件概率表如表4所示。60岁以上70岁以下，有高血压、心脏病和高血糖的患者，缺血性脑卒中发病概率为0.870。而70岁以上，有高血压和高血糖的患者，缺血性脑卒中发病概率为0.764。
表4
缺血性脑卒中发病BN预测模型条件概率分布表
2.3.4 BN交叉验证
研究使用bnlearn包中的bn.cv()函数进行了10折交叉验证。该函数提供了三种交叉验证方法，包括k折交叉验证（默认）、自定义折叠交叉验证（"custom-fold"）和Hold-out交叉验证（"hold-out"），通过method对使用的交叉验证方法进行修改。我们使用对数似然损失（log-likelihood loss）对BN模型性能进行评价，一般而言，期望损失值越低，表示模型在训练集上的性能越好。本研究得到训练集的期望损失为2.858，表明该模型在训练集上表现良好。
2.3.5 模型效果评价
使用pROC包构建BN模型与Logistic回归模型的ROC曲线，并分别计算AUC值、灵敏度、特异度和准确度，对模型效果进行评价（表5、图2）。
表5
BN模型与Logistic回归模型性能评估
图2
贝叶斯网络模型与Logistic回归模型的ROC曲线
与Logistic回归模型相比，在测试集和训练集中BN模型均表现出更好的性能。在训练集和测试集中，BN模型的准确率分别为79.37%和83.32%，AUC分别为0.85［95%CI（0.83，0.86）］和0.88［95%CI（0.86，0.89）］，灵敏度分别为59.64%和72.80%，特异度分别为87.48%和87.59%。训练集与测试集的Delong检验证实了两个模型之间AUC之间存在统计学差异（P<0.05）。此外，BN模型的校准图显示，在训练集和测试集中预测概率和实际观测值之间较为一致，见图3。
图3
贝叶斯网络模型校准曲线
3 讨论
本研究基于实例采用R语言实现个体化风险预测模型的构建，研究结果显示，BN模型的预测性能优于传统Logistic回归模型。且通过DAG图可以清晰的展示出变量之间的相互作用关系，年龄、心脏病、高血压和高血糖与缺血性脑卒中发病直接相关。此外，年龄还可通过心脏病、高血压和高血糖与缺血性脑卒中间接相关，高血糖也可通过影响心脏病和高血压，间接与缺血性脑卒中相关。
基于BN的风险预测模型相较于统计回归的方法具有多个优势[20-22]。首先，BN模型通过其图形结构和条件概率分布清晰地阐明了各变量之间复杂的相互作用关系，以及变量在网络中的相对位置，描述了变量对结果的贡献程度[7,23]。这种清晰的关系表达提高了对结果风险的评估准确性，并增强了结果的可解释性[14]。其次，BN模型可以充分利用先验知识和现有数据，有效地评估风险，并且在加入新数据时无需重新训练模型[8,24]，这使得其在个体化风险预测方面更加灵活和实用。随着时间的推移，更多的病历信息被纳入，BN模型可以获得进一步的提升。此外，BN模型能够进行从原因到结果和从结果到原因的推理，以及处理变量间的线性和非线性关系，因而具有更广泛的适用性[14]。
相比之下，基于回归的方法局限于单一结果的预测，在加入新数据时需要重新估计，且在面对数据缺失或个体特异性较强的情况下，基于回归的模型复杂程度较高且准确性有限。此外，基于回归的方法在纳入变量和结果之间的非线性关系时，通常涉及变量的转换、加权等操作，往往使得结果难以解释[14]。因此，BN模型在风险预测中展现出了显著的优势，并在实践中得到了广泛应用。
在临床中，疾病的进展往往与时间密切相关，BN能够与时间信息相结合，构成动态BN[25]，从而更好地反映变量的发展变化规律。同时，BN也可以与多种机器学习技术和传统统计方法结合使用，更好地处理数据，如与深度学习结合的贝叶斯深度学习模型[26]，这为解决复杂问题提供了新的途径。
BN为处理复杂的医疗卫生数据提供了一种强大且灵活的分析方法。数据规模的扩大、数据复杂性的增加等挑战的不断出现，传统的数据分析方法面临诸多限制。而BN能够应对这些挑战，不断更新信息。随着“精准医学”时代的到来，新的风险评估与决策制定的方法逐渐引起人们的关注，而BN凭借其高度的灵活性、适应性以及强大的可解释性，为大数据的分析提供了一种新的思路，值得进一步探索。本文提供了清晰的代码和代码解释，以实例的形式进行演示，方便研究人员更好地理解和应用BN模型。
