课程视频|R语言bnlearn包：贝叶斯网络的构造及参数学习的原理和实例(中)：https://developer.aliyun.com/article/1496663
比较线性模型
如果参数估计因任何原因出现问题，我们可以用一组新的、来自不同方法的估计值来取代它们。
dANB dANB = penalized( dANB) dANB
模型验证
有两种主要的方法来验证一个BN。
- 只看网络结构：如果学习BN的主要目标是识别弧和路径，当BN被解释为因果模型时，通常是这种情况，我们可以进行本质上的路径分析和研究弧的强度。
- 将BN视为一个整体，包括参数：如果学习BN的主要目标是将其作为一个专家模型，那么我们可能想。
- 根据其他一些变量的值，预测新个体的一个或多个变量的值；以及
- 将CP查询的结果与专家知识进行比较，以确认BN反映了关于我们正在建模的现象的最佳知识。
预测准确性
我们可以用通常的方法来衡量我们所选择的学习策略的预测准确性，即交叉验证。实现了：
- k-fold交叉验证;
- 指定的k进行交叉验证;
- hold-out 交叉验证
对于:
- 结构学习算法（结构和参数都是从数据中学习的）。
- 参数学习算法（结构由用户提供，参数从数据中学习）。
首先，我们检查Growth，它编码了错牙合畸形的演变（0表示坏，1表示好）。我们检查它，把它转回离散变量并计算预测误差。
cv(diff) > for (i in 1:10) { + err[i] = (sum(tt) - sum(diag(tt))) / sum(tt) + } >
其他变量是连续的，所以我们可以估计它们的预测相关性来代替。
> for (var in names(predcor)) { + xval = cv(diff) + predcor[var] = mean(sapply(xval, function(x) attr(x, "mean"))) + }
在这两种情况下，我们使用损失函数的变体，它使用从所有其他变量计算的后验期望值进行预测。基本的损失函数（cor, mse, pred）仅仅从它们的父代来预测每个节点的值，这在处理很少或没有父代的节点时是没有意义的。
用专家知识进行确认
确认BN是否有意义的另一种方法是把它当作的工作模型，看看它是否表达了关于关键事实，这些事实在学习过程中没有作为先验知识使用。否则，我们将只是拿回我们放在先验中的信息）。一些例子。
- "CoGo的过度增长应该会引起PPPM的减少"。
我们通过为存储在 fitted.simpler中的BN生成dCoGo和dPPPM的样本，并假设没有发生任何处理，来测试这个假设。随着dCoGo的增加（这表明增长越来越快），DPPPM变得越来越负（这表明假设角度最初是正的，则角度会减少。
> sim = dist(fitted.simpler) > plot(sim ) > abline(v = 0, col = 2, lty = 2, lwd = 2)
-
- "CoGo的小幅增长应该会引起PPPM的增长。"
从上图来看，CoGo的负增长或空增长（dCoGo ⋜ 0）对应于PPPM的正增长，概率为≈0.60。对于CoGo的小幅增长（dCoGo∈[0, 2]），不幸的是，dPPPM ⋜0，概率≈0.50，所以BN不支持这一假设。
> nrow(sim[( dCoGo <= 0) & ( PPPM > 0), ]) / nrow(sim[( dCoGo <= 0), ]) [1] 0.6112532 > nrow(sim[( dCoGo > 0) & ( dCoGo < 2) & ( dPPPM > 0), ]) / + nrow(sim[( CoGo) > 0 & ( dCoGo < 2), ]) [1] 0.4781784
"如果ANB减少，IMPA就会减少以进行补偿。"
像以前一样通过模拟测试，我们正在寻找与IMPA（相同）的负值相关的dANB的负值（这表明假设角度最初是正的，就会减少）。从下图中可以看出，dANB与dIMPA成正比，所以其中一个的减少表明另一个的减少；两者的平均趋势（黑线）同时为负。
> plot(sim ) > abline(coef(lm(dIMPA ~ dANB ))
"如果GoPg强烈增加，那么ANB和IMPA都会减少。" 如果我们从BN中模拟dGoPg、dANB和dIMPA，假设dGoPg>5（即GoPg在增加），我们估计dANB>0（即ANB在增加）的概率为≈0.70，dIMPA<0的概率仅为≈0.58。
> nrow(sim[(dGoPg > 5) & (dANB < 0), ]) / nrow(sim[(dGoPg > 5), ]) [1] 0.695416 > nrow(sim[(dGoPg > 5) & (dIMPA < 0), ]) / nrow(sim[(dGoPg > 5), ]) [1] 0.5756936
"治疗试图阻止ANB的减少。如果我们固定ANB，治疗过的病人和未治疗过的病人是否有区别？"
首先，我们可以检查在没有任何干预的情况下，dANB≈0的病人的治疗和增长之间的关系（即使用我们从数据中得知的BN）。
dist(fitted ) table(TREATMENT = Treatment < 0.5, GOOD.GROWTH = Growth > 0.5)
估计的P（GOOD.GROWTH ∣ TREATMENT）对于接受治疗和未接受治疗的病人是不同的（≈0.65对≈0.52）。
如果我们模拟一个正式的干预措施（如Judea Pearl），并从外部设置dANB=0（从而使其独立于其父母，并删除相应的弧），我们就会发现GOOD.GROWTH对于接受治疗和未接受治疗的病人来说实际上具有相同的分布，从而变得与TREATMENT无关。这表明，有利的预后确实是由防止ANB的变化决定的，而治疗的其他成分（如果有的话）就变得不重要了。
table(TREATMENT = Treatment < 0.5, GOOD.GROWTH = Growth > 0.5)
"治疗试图阻止ANB的减少。如果我们固定ANB，治疗和未治疗的病人之间是否有区别？"
评估的方法之一是检查在保持GoPg固定的情况下，A点和B点之间的角度（ANB）是否在治疗和未治疗的病人之间发生变化。
假设GoPg不发生变化，对于接受治疗的病人来说，A点和B点之间的角度会增加（强烈的负值表示水平不平衡，所以正的变化率表示不平衡的减少），而对于未接受治疗的病人来说则会减少（不平衡会随着时间慢慢恶化）。
Treatment = c("UNTREATED", "TREATED")[(Treatment > 0.5) + 1L] boxplot(dANB ~ Treatment)
模型#2：动态贝叶斯网络
动态贝叶斯网络在预测方面的效果不如1号模型好，同时更加复杂。这是动态贝叶斯网络所固有的，即模拟随机过程的贝叶斯网络：每个变量都与被模拟的每个时间点的不同节点相关。(通常情况下，我们假设过程是一阶马尔可夫，所以我们在BN中有两个时间点：t和t-1。)然而，我们探索它的目的是为了说明这样一个BN可以被学习并用于bnlearn。
我们用于这个模型的数据是我们在分析开始时存储到正交的原始数据。然而，我们选择使用治疗变量而不是生长变量作为变量来表达受试者可能正在接受医疗的事实。原因是生长变量是一个衡量第二次测量时的预后的变量，它的值在第一次测量时是未知的；而治疗变量在两次测量时都是相同的。
学习结构
首先，我们将变量分为三组：时间为t2的变量，时间为t1=t2-1的变量，以及与时间无关的变量，因为它们在t1和t1取值相同。
> t2.variables
然后我们引入一个黑名单，其中。
- 我们将所有从临床变量到T1、T2和治疗的弧线列入黑名单，因为我们知道，年龄和治疗不是由临床测量决定的。
- 我们将所有进入Treatment和t1时间段的所有变量的弧列入黑名单，因为我们假设t1时间段的变量之间的弧与t2时间段的相应变量是一样的，两次学习它们是没有意义的。
- 我们将所有从t2到t1的弧列入黑名单。
grid(from = setdiff(names(ortho), c("T1", "T2")), to = c("T1", "T2"))
相比之下，我们只对T1→T2的弧线进行白名单，因为第二次测量的年龄显然取决于第一次测量的年龄。
> data.frame(from = c("T1"), to = c("T2"))
最后我们可以用bl和wl来学习BN的结构。
> dyn.dag
很明显，这个BN比前一个更复杂：它有更多的节点（16对9），更多的弧（27对19），因此有更多的参数（218对37）。
绘制这个新模型的最好方法是用plot()开始。
plot(dyn, render = FALSE)
然后，我们对变量进行分组，以方便区分const、t1.variables和t2.variables；我们选择从左到右而不是从上到下绘制网络。
+ attrs = list(graph = list(rankdir = "LR"))) > Graph(gR)
与前一个模型一样，治疗作用于ANB：从治疗出去的唯一弧是治疗→ANB2和治疗→CoA2。同样，这两个子节点都与Down的A点有关。
结构学习中的模型平均化
我们想评估这个动态BN结构的稳定性，就像我们之前对静态BN所做的那样，我们可以再次做到这一点。
> boot (ortho ) > plot(dyn)
avernet(dyn.str)
平均下来的avg和dag几乎是一样的：它们只相差两道弧。这表明结构学习产生了一个稳定的输出。
compare(dag, avg) tp fp fn 26 1 1
学习参数
由于Treatment是一个离散变量，BN是一个CLGBN。这意味着以Treatment为父节点的连续节点的参数化与其他节点不同。
fit(dynavg)
我们可以看到，ANB2取决于ANB（所以，在前一个时间点的同一变量）和治疗。ANB是连续的，所以它被用作ANB2的回归因子。 治疗变量是离散的，决定了线性回归的成分。
模型验证和推理
我们可以对这个新模型提出另一组问题
"在不同的治疗下，ANB从第一次测量到第二次测量的转变程度如何？"
我们可以用cpdist()生成一对(ANB, ANB2)，条件是治疗方法等于NT、TB和TG，并观察其分布。
data.frame( diff = c(nt[, 2] - nt[, 1], tb[, 2] - tb[, 1], tg[, 2] - tg[, 1]), > by(effect$diff, effect$treatment, FUN = mean)
density(~ diff, groups = treatment)
我们知道，治疗试图阻止ANB的下降；这与NT的分布是在TB的左边，而TB是在TG的左边这一事实相一致。未经治疗的病人病情继续恶化；治疗无效的病人没有真正改善，但他们的病情也没有恶化；而治疗有效的病人则有改善。
相比之下，这是一个未经治疗的病人在相同初始条件下的模拟轨迹。
对CoA的模拟轨迹是比较现实的：它随着年龄的增长而减慢。这与ANB不同，它的发生是因为CoA2同时取决于T1和T2。(ANB2则两者都不依赖)。
> for (i in seq(nrow(interv)) { + # 进行联合预测，目前用predict()无法实现。 + dist(dyn.fitted, nodes = c(), + intervals[i,] = weighted.mean(ANB2, weights) + intervals[i,] = weighted.mean(CoA2, weights)
