课程视频|R语言bnlearn包：贝叶斯网络的构造及参数学习的原理和实例(上)：https://developer.aliyun.com/article/1496661
我们可以在热图中看到两个集群：第一个集群包括dCoGo、dGoPg和dCoA，第二个集群包括Treatment、dANB和dCoA。第一个聚类在临床上很有意思，因为它包括治疗和两个都与唐氏A点有关的变量，这为治疗的主要效果提供了一些线索。
plot(ug )
模型#1：作为差异模型的静态贝叶斯网络
在这里，我们使用保存在diff中的差异来为数据建模，而不是原始值；我们将使用GBN处理，因为所有变量都是数字。对差异进行建模会导致局部分布，其形式为回归模型
其中 
对于其他回归因子，以此类推。我们可以将这种回归改写为 
这是一组微分方程，对变化率进行建模，其关系被假定为很好地近似于线性关系。然而，这种表述仍然意味着原始值随时间线性变化，因为变化率取决于其他变量的变化率，但不取决于时间本身。要有一个非线性的趋势，我们需要
此外，包括增长变量意味着我们可以有以下形式的回归模型
从而允许不同的变化率，这取决于病人是否在畸形中表现出积极的发展，以及他是否正在接受治疗。
学习贝叶斯网络
学习结构
学习BN的第一步是学习其结构，即DAG 
. 我们可以使用数据（来自不同的数据框架）结合先验知识来做这件事；结合后者可以减少我们必须探索的模型空间，并生成更强大的BN。一个直接的方法是将那些编码我们知道不可能/真实的关系的弧列入黑名单; 并将那些编码我们知道存在的关系的弧列入白名单。
黑名单只是一个矩阵（或一个数据框），其中有from和to两列，列出了我们不希望在BN中出现的弧。
- 我们把任何指向正畸变量中的dT、治疗和生长的弧列入黑名单。
- 我们将从dT到Treatment的弧列入黑名单。这意味着一个病人是否被治疗不会随时间而改变。
- 我们将从生长到dT和治疗的弧线列入黑名单。这意味着病人是否接受治疗不会随时间变化，而且显然不会因预后而变化。
白名单的结构与黑名单相同。
- 我们将依赖结构dANB → dIMPA ← dPPPM列入白名单。
- 我们将从dT到Growth的弧线列入白名单，这使得预后可以随时间变化。
一个简单的学习 
方法是在整个数据上找到具有最佳拟合度的网络结构。例如，使用hc()与默认分数（BIC）和整个diff数据框架。
至于绘图，关键函数是plot()。
plot(dag, , highlight )
然而，dag的质量关键取决于变量是否是正态分布，以及连接它们的关系是否是线性的；从探索性分析来看，并不清楚所有的变量都是如此。我们也不知道哪些弧线代表强关系，也就是说，它们能抵抗数据的扰动。我们可以用boot来解决这两个问题。
- 使用bootstrap对数据重新取样。
- 从每个bootstrap样本中学习一个单独的网络。
- 检查每个可能的弧在网络中出现的频率。
- 用出现频率较高的弧构建一个共识网络。
booth(diff, R = 200)
boot.strength()的返回值包括，对于每一对节点，连接它们的弧的强度（例如，我们观察到dANB → dPPPM或dPPPM → dANB的频率）及其方向的强度（例如，当我们观察到dANB和dPPPM之间有弧时，我们观察到dANB → dPPPM的频率）。
attr( "threshold")
因此，averaged.network()取所有强度至少为0.585的弧，并返回一个平均的共识网络，除非指定不同的阈值。
> avg.diff = averaged.network(str.diff)
纳入我们现在拥有的关于弧线强度的信息。
> strength.plot(avg.diff, str.diff, shape = "ellipse", highlight = list(arcs = wl))
我们如何将平均的网络（avg.diff）与我们最初从所有数据中学习到的网络（dag）进行比较？最定性的方法是将两个网络并排绘制，节点位置相同，并突出显示一个网络中出现而另一个网络中没有的弧，或者出现的方向不同的弧。
> par(mfrow = c(1, 2)) > graphviz.compare(avg.diff, dag, shape = "ellipse", main = c("averaged DAG", "single DAG"))
我们可以看到，Treatment→dIMPa、dANB→dGoPg和dCoGo→dPPPM这些弧线只出现在平均网络中，而dPPPM→dANB只出现在我们从所有数据中学到的网络中。我们可以假设，前三个弧被数据的噪声加上小样本量和偏离常态的情况所隐藏。编程可以返回真阳性（出现在两个网络中的弧）和假阳性/阴性（只出现在两个网络中的一个的弧）的数量。
> compare
或弧=TRUE。
但是，考虑到网络是用BIC学习的，而BIC是等价的，那么所有的弧线方向是否都很确定？看一下dag和avg.diff的CPDAGs（并考虑到白名单和黑名单），我们看到没有无方向的弧。所有弧的方向都是唯一的。
最后，我们可以结合来进行原则性的比较，如果两个弧被唯一确定为不同，我们就说它们是不同的。
也可以看一下关于弧长分布的阈值：平均的网络是相当密集的（9个节点有17个弧），很难阅读。
> plot(str.diff) > abline(v = 0.75, col = "tomato", lty = 2, lwd = 2) > abline(v = 0.85, col = "steelblue", lty = 2, lwd = 2)
因此，把阈值提高一点，多剔除几个弧就好了。看一下上面的图，由于弧长分布的差距，较高的阈值的两个自然选择是0.75（红色虚线）和0.85（蓝色虚线）。
> nrow( strength > "threshold" direction > 0.5, ]) [1] 18 trength > 0.75 & direction > 0.5 [1] 15 strength > 0.85 & direction > 0.5 [1] 12
我们通过在 network()中设置阈值=0.85得到的更简单的网络如下所示；从定性的角度来看，它当然更容易推理。
> avg.simpler = averaged.network(str.diff, threshold = 0.85) > strength.plot(avg.simpler, str.diff, shape = "ellipse", highlight = list(arcs = wl))
学习参数
在学习了结构之后，我们现在可以学习参数。由于我们正在处理连续变量，我们选择用GBN来建模。因此，如果我们使用最大似然估计来拟合网络的参数，我们就会发现每个局部分布都是一个典型的线性回归。
fit(avg, diff)
我们可以通过比较bn.fit()和lm()产生的模型，例如dANB，很容易确认这是事实。
> summary(lm(dANB ~ Growth + Treatment, data = diff))
我们会不会有共线性的问题？理论上是可能的，但在实践中，从数据中学习的网络结构大多不是问题。原因是，如果两个变量 
和 
是共线性的，在增加（比如说）Xi←Xj之后，那么Xj←Xk将不再显著提高BIC，因为Xj和Xk（在某种程度上）提供了关于Xi的相同信息。
> # 逐渐增加解释变量之间的关联性。 > for (rho 5)) { + # 更新相关矩阵并生成数据。 + R = R = rho + data = as.data.frame(mvrnorm(1000)) + # 比较线性模型 + cat( " BIC:", + }
课程视频|R语言bnlearn包：贝叶斯网络的构造及参数学习的原理和实例（下）：https://developer.aliyun.com/article/1496666
