全文鏈接：http://tecdat.cn/?p=32462

原文出處：拓端數(shù)據(jù)部落公眾號

分析師：Chang Zhang

貝葉斯網(wǎng)絡(luò)（BN）是一種基于有向無環(huán)圖的概率模型，它描述了一組變量及其相互之間的條件依賴性。貝葉斯網(wǎng)絡(luò)在信息不完備的情況下通過可以觀察隨機變量推斷不可觀察的隨機變量，對于解決復(fù)雜的不確定性和關(guān)聯(lián)性問題有很強的優(yōu)勢。

基本概念

有向無環(huán)圖（DAG）：DAG的結(jié)點V包括隨機變量（類別和特征），有向連接E（A->B）表示結(jié)點A是結(jié)點B的parent，且B與A是有依賴關(guān)系的（不獨立）。 條件概率表集合：引入了一個條件性獨立（conditional independence）概念，即圖中任意結(jié)點v在給定v的parent結(jié)點的情況下，與圖中其它結(jié)點都是獨立的，也就是說P(v|par(v),x1,x2...,xn) = P(v|par(v))。這里par(v)表示v的parent結(jié)點集，x1,x2,...,xn表示圖中其它結(jié)點。

貝葉斯網(wǎng)絡(luò)實施步驟

1、確定隨機變量間的拓撲關(guān)系，形成DAG。 2、訓(xùn)練貝葉斯網(wǎng)絡(luò)參數(shù)——估計出各節(jié)點的條件概率表。這一步也就是要完成條件概率表的構(gòu)造。

結(jié)構(gòu)學(xué)習(xí)

1、基于評分搜索的方法：將結(jié)構(gòu)學(xué)習(xí)視為組合優(yōu)化問題，首先通過定義評分函數(shù)對BN結(jié)構(gòu)空間中的不同元素與樣本數(shù)據(jù)的擬合程度進行度量，然后利用搜索算法確定評分最高的網(wǎng)絡(luò)結(jié)構(gòu)，即與數(shù)據(jù)擬合最好的網(wǎng)絡(luò)結(jié)構(gòu)。

2、基于約束的方法：利用統(tǒng)計或信息論的方法定量的分析變量間的依賴關(guān)系以獲取最優(yōu)地表達這些關(guān)系的網(wǎng)絡(luò)結(jié)構(gòu)。首先對訓(xùn)練數(shù)據(jù)集進行統(tǒng)計測試，尤其是條件獨立性測試，確定出變量之間的條件獨立性，然后利用變量之間的條件獨立性構(gòu)造一個有向無環(huán)圖，以盡可能多的涵蓋這些條件獨立性。

3、基于評分搜索和約束相混合的方法：首先采用條件獨立性檢驗縮減搜索空間，然后將此搜索空間作為輸入以限制基于評分的搜索。

結(jié)構(gòu)學(xué)習(xí)——基于評分搜索的方法

將結(jié)構(gòu)學(xué)習(xí)視為組合優(yōu)化問題，首先通過定義評分函數(shù)對BN結(jié)構(gòu)空間中的不同元素與樣本數(shù)據(jù)的擬合程度進行度量，然后利用搜索算法確定評分最高的網(wǎng)絡(luò)結(jié)構(gòu)，即與數(shù)據(jù)擬合最好的網(wǎng)絡(luò)結(jié)構(gòu)。

1、爬山法（hc）：從一個初始的網(wǎng)絡(luò)結(jié)構(gòu)開始，通過加邊、減邊和轉(zhuǎn)邊三個搜索算子對當(dāng)前網(wǎng)絡(luò)結(jié)構(gòu)進行修改，并對修改后的網(wǎng)絡(luò)結(jié)構(gòu)進行評分，最后找出評分最高的結(jié)構(gòu)

2、禁忌搜索（tabu）:是一種亞啟發(fā)式隨機搜索算法，通過對局部領(lǐng)域逐步搜索得到全局最優(yōu)解。為了避免陷入局部最優(yōu)和重復(fù)迭代，禁忌搜索算法采用了一種禁忌表，利用該表對解的移動搜索過程進行記錄和選擇，當(dāng)搜索的最優(yōu)解已經(jīng)在禁忌表中時，則會自動放棄該解，避免重復(fù)搜索。

結(jié)構(gòu)學(xué)習(xí)——基于約束的方法（Grow-Shrink）

結(jié)構(gòu)學(xué)習(xí)——基于混合的方法

MMHC算法結(jié)合了條件獨立性測試、搜索和評分等方法, 該算法分為2個階段, 在第1階段利用MMPC (max-min parents and children) 算法構(gòu)建貝葉斯網(wǎng)絡(luò)結(jié)構(gòu)的框架, 在第2階段執(zhí)行評分搜索確定網(wǎng)絡(luò)結(jié)構(gòu)的邊以及邊的方向.

參數(shù)學(xué)習(xí)

R語言貝葉斯網(wǎng)絡(luò)快速介紹

貝葉斯網(wǎng)絡(luò)

定義

貝葉斯網(wǎng)絡(luò)（BNs）的定義是：

• 一個網(wǎng)絡(luò)結(jié)構(gòu)，一個有向無環(huán)圖?

• , 其中每個節(jié)點?

• ?對應(yīng)于一個隨機變量?

• ;

• 一個全局概率分布?

• ?(帶參數(shù)?

• ), 它可以根據(jù)圖中存在的弧被分解成更小的局部概率分布。

網(wǎng)絡(luò)結(jié)構(gòu)的主要作用是通過圖形分離來表達模型中各變量之間的條件獨立性關(guān)系，從而指定全局分布的因子化。

每個局部分布都有自己的參數(shù)集?

; 而??

?要比

小得多，因為許多參數(shù)是固定的，因為它們所屬的變量是獨立的。

R實現(xiàn)了以下學(xué)習(xí)算法。

基于約束的：PC, GS, IAMB, MMPC, Hilton-PC
基于分?jǐn)?shù)的：爬山算法、Tabu Search
配對的：ARACNE, Chow-Liu
混合：MMHC, RSMAX2

我們使用基于分?jǐn)?shù)的學(xué)習(xí)算法，希爾算法。首先，我們將先為本教程生成簡單的數(shù)據(jù)集。

在這個數(shù)據(jù)集中，'狀態(tài)'與'元素'和'接受'列有關(guān)系。而'類型'與'顏色'列有關(guān)系。當(dāng)你創(chuàng)建一個帶有分類數(shù)據(jù)的數(shù)據(jù)框時，列應(yīng)該是一個因子類型。否則，該數(shù)據(jù)框不能用于BN結(jié)構(gòu)的創(chuàng)建。

?
接下來，我們將創(chuàng)建學(xué)習(xí)結(jié)構(gòu)。

我們可以在一個圖中看到結(jié)構(gòu)。

> plot(hc_simd)

在這個圖中，狀態(tài)、元素、接受、類型和顏色被稱為節(jié)點。節(jié)點之間的方向用弧線描述，弧線是一個包含從元素到元素方向數(shù)據(jù)的矩陣。

如上弧線顯示，在我們的數(shù)據(jù)中存在'類型'到'顏色'，以及'狀態(tài)'到'接受'和'元素'的關(guān)系。'類型'和'狀態(tài)'是兩個獨立的組，它們之間不存在相互依賴關(guān)系。
接下來，我們將用數(shù)據(jù)來擬合模型。

simd_fitted

基于上述訓(xùn)練數(shù)據(jù)，我們可以進行條件概率查詢。

我們檢查 "Outlier "和 "Target "的狀態(tài)概率。

該樣本成為 "離群 "的概率為51%。

狀態(tài)成為 "目標(biāo) "的概率是0%。

錯頜畸形數(shù)據(jù)的貝葉斯網(wǎng)絡(luò)分析

問題：?受第三類錯牙合畸形影響的患者（以下牙弓突出為特征），其骨骼不平衡在生命早期就產(chǎn)生，在青春期和骨骼成熟前會變得更加明顯。在單個III類患者中早期預(yù)測治療的成功或失敗，使其更容易矯正，但僅從少量的形態(tài)決定因素中預(yù)測是很難做到的。原因是III類錯頜畸形很少是單一顱面部件異常的結(jié)果，所以單個的臨床和放射學(xué)測量值可能不如測量值本身的相互作用具有指示性。

任務(wù):

1. 我們學(xué)習(xí)一個BN，并使用它來確定和可視化在成長和治療過程中各種III類錯位頜面特征之間的相互作用。

2. 我們通過驗證一些普遍接受的關(guān)于這些骨骼不平衡演變的假說來檢驗其一致性。

3. 我們表明，與接受快速上頜擴張和面罩治療的正畸患者相比，未經(jīng)治療的受試者形成了不同的III類顱面生長模式。

4. 在接受治療的患者中，CoA段（上頜骨長度）和ANB角（上頜骨與下頜骨的前后關(guān)系）似乎是接受治療的主要影響的骨骼亞空間。

數(shù)據(jù)

我們將使用的數(shù)據(jù)集包含143名患者，在T1和T2年齡段有兩組測量數(shù)據(jù)（以年為單位），用于以下變量。

• 治療：未經(jīng)治療（NT），治療后效果不好（TB），治療后效果好（TG）。

• 生長：一個二元變量，數(shù)值為好或壞。

• ANB：唐氏點A和B之間的角度（度）。

• IMPA：門牙-下頜平面角（度）。

• PPPM：腭平面-下頜平面的角度（度）。

• CoA：上頜骨從髁狀突到唐氏點A的總長度（mm）。

• GoPg：下頜體從齒齦到齒齦的長度（mm）。

• CoGo：下頜骨的長度，從髁狀突到齒狀突（mm）。

所有的測量都是通過X射線掃描得出的，使用一套參考點建立的圖，如以下。

> str(data)

預(yù)處理和探索性數(shù)據(jù)分析

首先，我們創(chuàng)建一個數(shù)據(jù)框架，其中包括所有變量的差異以及增長和治療。

生長和治療變量帶有關(guān)于病人預(yù)后的冗余信息，這一點從TB和TG之間生長良好的病人比例的差異中可以看出。

為了避免在模型中包括這兩個變量所導(dǎo)致的混雜，我們將治療重新編碼為一個二元變量，0表示NT，1表示TB或TG。同樣地，我們對成長進行重新編碼，0表示壞，1表示好。

由于我們將使用高斯BN進行分析，檢查這些變量是否是正態(tài)分布；從下面的圖來看，似乎并非所有的變量都是如此。

+ ? hist(x, prob = TRUE ) + ? lines(density(x), lwd = 2 )

這些變量是通過線性關(guān)系聯(lián)系起來的嗎？其中一些是，但不是全部。

> pairs(diff[, setdiff(names(diff)

最后，我們可以看看這些變量是否以任何方式聚在一起，因為聚在一起的變量更有可能在BN中發(fā)生聯(lián)系。

> heatmap(rho)

我們可以在熱圖中看到兩個集群：第一個集群包括dCoGo、dGoPg和dCoA，第二個集群包括Treatment、dANB和dCoA。第一個聚類在臨床上很有意思，因為它包括治療和兩個都與唐氏A點有關(guān)的變量，這為治療的主要效果提供了一些線索。

plot(ug )

模型#1：作為差異模型的靜態(tài)貝葉斯網(wǎng)絡(luò)

在這里，我們使用保存在diff中的差異來為數(shù)據(jù)建模，而不是原始值；我們將使用GBN處理，因為所有變量都是數(shù)字。對差異進行建模會導(dǎo)致局部分布，其形式為回歸模型

其中

?對于其他回歸因子，以此類推。我們可以將這種回歸改寫為

這是一組微分方程，對變化率進行建模，其關(guān)系被假定為很好地近似于線性關(guān)系。然而，這種表述仍然意味著原始值隨時間線性變化，因為變化率取決于其他變量的變化率，但不取決于時間本身。要有一個非線性的趨勢，我們需要

此外，包括增長變量意味著我們可以有以下形式的回歸模型?

從而允許不同的變化率，這取決于病人是否在畸形中表現(xiàn)出積極的發(fā)展，以及他是否正在接受治療。

學(xué)習(xí)貝葉斯網(wǎng)絡(luò)

學(xué)習(xí)結(jié)構(gòu)

學(xué)習(xí)BN的第一步是學(xué)習(xí)其結(jié)構(gòu)，即DAG?

. 我們可以使用數(shù)據(jù)（來自不同的數(shù)據(jù)框架）結(jié)合先驗知識來做這件事；結(jié)合后者可以減少我們必須探索的模型空間，并生成更強大的BN。一個直接的方法是將那些編碼我們知道不可能/真實的關(guān)系的弧列入黑名單; 并將那些編碼我們知道存在的關(guān)系的弧列入白名單。

黑名單只是一個矩陣（或一個數(shù)據(jù)框），其中有from和to兩列，列出了我們不希望在BN中出現(xiàn)的弧。

• 我們把任何指向正畸變量中的dT、治療和生長的弧列入黑名單。

• 我們將從dT到Treatment的弧列入黑名單。這意味著一個病人是否被治療不會隨時間而改變。

• 我們將從生長到dT和治療的弧線列入黑名單。這意味著病人是否接受治療不會隨時間變化，而且顯然不會因預(yù)后而變化。

白名單的結(jié)構(gòu)與黑名單相同。

• 我們將依賴結(jié)構(gòu)dANB → dIMPA ← dPPPM列入白名單。

• 我們將從dT到Growth的弧線列入白名單，這使得預(yù)后可以隨時間變化。

一個簡單的學(xué)習(xí)?

?方法是在整個數(shù)據(jù)上找到具有最佳擬合度的網(wǎng)絡(luò)結(jié)構(gòu)。例如，使用hc()與默認分?jǐn)?shù)（BIC）和整個diff數(shù)據(jù)框架。

至于繪圖，關(guān)鍵函數(shù)是plot()。

plot(dag, , highlight )

然而，dag的質(zhì)量關(guān)鍵取決于變量是否是正態(tài)分布，以及連接它們的關(guān)系是否是線性的；從探索性分析來看，并不清楚所有的變量都是如此。我們也不知道哪些弧線代表強關(guān)系，也就是說，它們能抵抗數(shù)據(jù)的擾動。我們可以用boot來解決這兩個問題。

1. 使用bootstrap對數(shù)據(jù)重新取樣。

2. 從每個bootstrap樣本中學(xué)習(xí)一個單獨的網(wǎng)絡(luò)。

3. 檢查每個可能的弧在網(wǎng)絡(luò)中出現(xiàn)的頻率。

4. 用出現(xiàn)頻率較高的弧構(gòu)建一個共識網(wǎng)絡(luò)。

booth(diff, R = 200)

boot.strength()的返回值包括，對于每一對節(jié)點，連接它們的弧的強度（例如，我們觀察到dANB → dPPPM或dPPPM → dANB的頻率）及其方向的強度（例如，當(dāng)我們觀察到dANB和dPPPM之間有弧時，我們觀察到dANB → dPPPM的頻率）。?

attr( "threshold")

因此，averaged.network()取所有強度至少為0.585的弧，并返回一個平均的共識網(wǎng)絡(luò)，除非指定不同的閾值。?

> avg.diff = averaged.network(str.diff)

納入我們現(xiàn)在擁有的關(guān)于弧線強度的信息。?

> strength.plot(avg.diff, str.diff, shape = "ellipse", highlight = list(arcs = wl))

?

我們?nèi)绾螌⑵骄木W(wǎng)絡(luò)（avg.diff）與我們最初從所有數(shù)據(jù)中學(xué)習(xí)到的網(wǎng)絡(luò)（dag）進行比較？最定性的方法是將兩個網(wǎng)絡(luò)并排繪制，節(jié)點位置相同，并突出顯示一個網(wǎng)絡(luò)中出現(xiàn)而另一個網(wǎng)絡(luò)中沒有的弧，或者出現(xiàn)的方向不同的弧。?

> par(mfrow = c(1, 2)) > graphviz.compare(avg.diff, dag, shape = "ellipse", main = c("averaged DAG", "single DAG"))

我們可以看到，Treatment→dIMPa、dANB→dGoPg和dCoGo→dPPPM這些弧線只出現(xiàn)在平均網(wǎng)絡(luò)中，而dPPPM→dANB只出現(xiàn)在我們從所有數(shù)據(jù)中學(xué)到的網(wǎng)絡(luò)中。我們可以假設(shè)，前三個弧被數(shù)據(jù)的噪聲加上小樣本量和偏離常態(tài)的情況所隱藏。編程可以返回真陽性（出現(xiàn)在兩個網(wǎng)絡(luò)中的弧）和假陽性/陰性（只出現(xiàn)在兩個網(wǎng)絡(luò)中的一個的?。┑臄?shù)量。

> compare

或弧=TRUE。

但是，考慮到網(wǎng)絡(luò)是用BIC學(xué)習(xí)的，而BIC是等價的，那么所有的弧線方向是否都很確定？看一下dag和avg.diff的CPDAGs（并考慮到白名單和黑名單），我們看到?jīng)]有無方向的弧。所有弧的方向都是唯一的。?

最后，我們可以結(jié)合來進行原則性的比較，如果兩個弧被唯一確定為不同，我們就說它們是不同的。

也可以看一下關(guān)于弧長分布的閾值：平均的網(wǎng)絡(luò)是相當(dāng)密集的（9個節(jié)點有17個弧），很難閱讀。

> plot(str.diff) > abline(v = 0.75, col = "tomato", lty = 2, lwd = 2) > abline(v = 0.85, col = "steelblue", lty = 2, lwd = 2)

因此，把閾值提高一點，多剔除幾個弧就好了?？匆幌律厦娴膱D，由于弧長分布的差距，較高的閾值的兩個自然選擇是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()中設(shè)置閾值=0.85得到的更簡單的網(wǎng)絡(luò)如下所示；從定性的角度來看，它當(dāng)然更容易推理。

> avg.simpler = averaged.network(str.diff, threshold = 0.85) > strength.plot(avg.simpler, str.diff, shape = "ellipse", highlight = list(arcs = wl))

學(xué)習(xí)參數(shù)

在學(xué)習(xí)了結(jié)構(gòu)之后，我們現(xiàn)在可以學(xué)習(xí)參數(shù)。由于我們正在處理連續(xù)變量，我們選擇用GBN來建模。因此，如果我們使用最大似然估計來擬合網(wǎng)絡(luò)的參數(shù)，我們就會發(fā)現(xiàn)每個局部分布都是一個典型的線性回歸。

fit(avg, diff)

我們可以通過比較bn.fit()和lm()產(chǎn)生的模型，例如dANB，很容易確認這是事實。

> summary(lm(dANB ~ Growth + Treatment, data = diff))

我們會不會有共線性的問題？理論上是可能的，但在實踐中，從數(shù)據(jù)中學(xué)習(xí)的網(wǎng)絡(luò)結(jié)構(gòu)大多不是問題。原因是，如果兩個變量?

?和

是共線性的，在增加（比如說）Xi←Xj之后，那么Xj←Xk將不再顯著提高BIC，因為Xj和Xk（在某種程度上）提供了關(guān)于Xi的相同信息。

> # 逐漸增加解釋變量之間的關(guān)聯(lián)性。 > for (rho 5)) { + ? # 更新相關(guān)矩陣并生成數(shù)據(jù)。 + ? R ?= R = rho + ? data = as.data.frame(mvrnorm(1000)) + ? # 比較線性模型 + ? cat( " BIC:", + }

?

比較線性模型?

如果參數(shù)估計因任何原因出現(xiàn)問題，我們可以用一組新的、來自不同方法的估計值來取代它們。

dANB dANB = penalized( dANB) dANB

模型驗證

有兩種主要的方法來驗證一個BN。

1. 只看網(wǎng)絡(luò)結(jié)構(gòu)：如果學(xué)習(xí)BN的主要目標(biāo)是識別弧和路徑，當(dāng)BN被解釋為因果模型時，通常是這種情況，我們可以進行本質(zhì)上的路徑分析和研究弧的強度。

2. 將BN視為一個整體，包括參數(shù)：如果學(xué)習(xí)BN的主要目標(biāo)是將其作為一個專家模型，那么我們可能想。

• 根據(jù)其他一些變量的值，預(yù)測新個體的一個或多個變量的值；以及
  

• 將CP查詢的結(jié)果與專家知識進行比較，以確認BN反映了關(guān)于我們正在建模的現(xiàn)象的最佳知識。

預(yù)測準(zhǔn)確性

我們可以用通常的方法來衡量我們所選擇的學(xué)習(xí)策略的預(yù)測準(zhǔn)確性，即交叉驗證。實現(xiàn)了：

• k-fold交叉驗證;

• 指定的k進行交叉驗證;

• hold-out?交叉驗證

對于:

• 結(jié)構(gòu)學(xué)習(xí)算法（結(jié)構(gòu)和參數(shù)都是從數(shù)據(jù)中學(xué)習(xí)的）。

• 參數(shù)學(xué)習(xí)算法（結(jié)構(gòu)由用戶提供，參數(shù)從數(shù)據(jù)中學(xué)習(xí)）。

首先，我們檢查Growth，它編碼了錯牙合畸形的演變（0表示壞，1表示好）。我們檢查它，把它轉(zhuǎn)回離散變量并計算預(yù)測誤差。

cv(diff) > for (i in 1:10) { + ? err[i] = (sum(tt) - sum(diag(tt))) / sum(tt) + } >

其他變量是連續(xù)的，所以我們可以估計它們的預(yù)測相關(guān)性來代替。

> for (var in names(predcor)) { + ? xval = cv(diff) + ? ? predcor[var] = mean(sapply(xval, function(x) attr(x, "mean"))) + }

在這兩種情況下，我們使用損失函數(shù)的變體，它使用從所有其他變量計算的后驗期望值進行預(yù)測?；镜膿p失函數(shù)（cor, mse, pred）僅僅從它們的父代來預(yù)測每個節(jié)點的值，這在處理很少或沒有父代的節(jié)點時是沒有意義的。

用專家知識進行確認

確認BN是否有意義的另一種方法是把它當(dāng)作的工作模型，看看它是否表達了關(guān)于關(guān)鍵事實，這些事實在學(xué)習(xí)過程中沒有作為先驗知識使用。否則，我們將只是拿回我們放在先驗中的信息）。一些例子。

1. "CoGo的過度增長應(yīng)該會引起PPPM的減少"。
   我們通過為存儲在 fitted.simpler中的BN生成dCoGo和dPPPM的樣本，并假設(shè)沒有發(fā)生任何處理，來測試這個假設(shè)。隨著dCoGo的增加（這表明增長越來越快），DPPPM變得越來越負（這表明假設(shè)角度最初是正的，則角度會減少。

   > sim = dist(fitted.simpler) > plot(sim ) > abline(v = 0, col = 2, lty = 2, lwd = 2)

1. "CoGo的小幅增長應(yīng)該會引起PPPM的增長。"
   
   從上圖來看，CoGo的負增長或空增長（dCoGo ? 0）對應(yīng)于PPPM的正增長，概率為≈0.60。對于CoGo的小幅增長（dCoGo∈[0, 2]），不幸的是，dPPPM ?0，概率≈0.50，所以BN不支持這一假設(shè)。

   > 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

2. "如果ANB減少，IMPA就會減少以進行補償。"
   
   像以前一樣通過模擬測試，我們正在尋找與IMPA（相同）的負值相關(guān)的dANB的負值（這表明假設(shè)角度最初是正的，就會減少）。從下圖中可以看出，dANB與dIMPA成正比，所以其中一個的減少表明另一個的減少；兩者的平均趨勢（黑線）同時為負。

   > plot(sim ) > abline(coef(lm(dIMPA ~ dANB ))

1. "如果GoPg強烈增加，那么ANB和IMPA都會減少。" 如果我們從BN中模擬dGoPg、dANB和dIMPA，假設(shè)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

2. "治療試圖阻止ANB的減少。如果我們固定ANB，治療過的病人和未治療過的病人是否有區(qū)別？"
   首先，我們可以檢查在沒有任何干預(yù)的情況下，dANB≈0的病人的治療和增長之間的關(guān)系（即使用我們從數(shù)據(jù)中得知的BN）。

   dist(fitted ) table(TREATMENT = Treatment < 0.5, GOOD.GROWTH = ?Growth > 0.5)

1. 估計的P（GOOD.GROWTH ∣ TREATMENT）對于接受治療和未接受治療的病人是不同的（≈0.65對≈0.52）。
   如果我們模擬一個正式的干預(yù)措施（如Judea Pearl），并從外部設(shè)置dANB=0（從而使其獨立于其父母，并刪除相應(yīng)的?。?，我們就會發(fā)現(xiàn)GOOD.GROWTH對于接受治療和未接受治療的病人來說實際上具有相同的分布，從而變得與TREATMENT無關(guān)。這表明，有利的預(yù)后確實是由防止ANB的變化決定的，而治療的其他成分（如果有的話）就變得不重要了。

   table(TREATMENT = ?Treatment < 0.5, GOOD.GROWTH = ?Growth > 0.5)

1. "治療試圖阻止ANB的減少。如果我們固定ANB，治療和未治療的病人之間是否有區(qū)別？"
   評估的方法之一是檢查在保持GoPg固定的情況下，A點和B點之間的角度（ANB）是否在治療和未治療的病人之間發(fā)生變化。

假設(shè)GoPg不發(fā)生變化，對于接受治療的病人來說，A點和B點之間的角度會增加（強烈的負值表示水平不平衡，所以正的變化率表示不平衡的減少），而對于未接受治療的病人來說則會減少（不平衡會隨著時間慢慢惡化）。

Treatment = c("UNTREATED", "TREATED")[(Treatment > 0.5) + 1L] boxplot(dANB ~ Treatment)

模型#2：動態(tài)貝葉斯網(wǎng)絡(luò)

動態(tài)貝葉斯網(wǎng)絡(luò)在預(yù)測方面的效果不如1號模型好，同時更加復(fù)雜。這是動態(tài)貝葉斯網(wǎng)絡(luò)所固有的，即模擬隨機過程的貝葉斯網(wǎng)絡(luò)：每個變量都與被模擬的每個時間點的不同節(jié)點相關(guān)。(通常情況下，我們假設(shè)過程是一階馬爾可夫，所以我們在BN中有兩個時間點：t和t-1。)然而，我們探索它的目的是為了說明這樣一個BN可以被學(xué)習(xí)并用于bnlearn。

我們用于這個模型的數(shù)據(jù)是我們在分析開始時存儲到正交的原始數(shù)據(jù)。然而，我們選擇使用治療變量而不是生長變量作為變量來表達受試者可能正在接受醫(yī)療的事實。原因是生長變量是一個衡量第二次測量時的預(yù)后的變量，它的值在第一次測量時是未知的；而治療變量在兩次測量時都是相同的。

學(xué)習(xí)結(jié)構(gòu)

首先，我們將變量分為三組：時間為t2的變量，時間為t1=t2-1的變量，以及與時間無關(guān)的變量，因為它們在t1和t1取值相同。

> t2.variables

然后我們引入一個黑名單，其中。

1. 我們將所有從臨床變量到T1、T2和治療的弧線列入黑名單，因為我們知道，年齡和治療不是由臨床測量決定的。

2. 我們將所有進入Treatment和t1時間段的所有變量的弧列入黑名單，因為我們假設(shè)t1時間段的變量之間的弧與t2時間段的相應(yīng)變量是一樣的，兩次學(xué)習(xí)它們是沒有意義的。

3. 我們將所有從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來學(xué)習(xí)BN的結(jié)構(gòu)。

> dyn.dag

很明顯，這個BN比前一個更復(fù)雜：它有更多的節(jié)點（16對9），更多的?。?7對19），因此有更多的參數(shù)（218對37）。

繪制這個新模型的最好方法是用plot()開始。

plot(dyn, render = FALSE)

然后，我們對變量進行分組，以方便區(qū)分const、t1.variables和t2.variables；我們選擇從左到右而不是從上到下繪制網(wǎng)絡(luò)。

+ ? ? ? ?attrs = list(graph = list(rankdir = "LR"))) > Graph(gR)

與前一個模型一樣，治療作用于ANB：從治療出去的唯一弧是治療→ANB2和治療→CoA2。同樣，這兩個子節(jié)點都與Down的A點有關(guān)。

結(jié)構(gòu)學(xué)習(xí)中的模型平均化

我們想評估這個動態(tài)BN結(jié)構(gòu)的穩(wěn)定性，就像我們之前對靜態(tài)BN所做的那樣，我們可以再次做到這一點。

> boot (ortho ) > plot(dyn)

avernet(dyn.str)

?

平均下來的avg和dag幾乎是一樣的：它們只相差兩道弧。這表明結(jié)構(gòu)學(xué)習(xí)產(chǎn)生了一個穩(wěn)定的輸出。

compare(dag, avg) tp fp fn 26 ?1 ?1

學(xué)習(xí)參數(shù)

由于Treatment是一個離散變量，BN是一個CLGBN。這意味著以Treatment為父節(jié)點的連續(xù)節(jié)點的參數(shù)化與其他節(jié)點不同。?

fit(dynavg)

?

我們可以看到，ANB2取決于ANB（所以，在前一個時間點的同一變量）和治療。?ANB是連續(xù)的，所以它被用作ANB2的回歸因子。?治療變量是離散的，決定了線性回歸的成分。

模型驗證和推理

我們可以對這個新模型提出另一組問題

1. "在不同的治療下，ANB從第一次測量到第二次測量的轉(zhuǎn)變程度如何？"
   我們可以用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)

   ?

1. density(~ diff, groups = treatment)

1. 我們知道，治療試圖阻止ANB的下降；這與NT的分布是在TB的左邊，而TB是在TG的左邊這一事實相一致。未經(jīng)治療的病人病情繼續(xù)惡化；治療無效的病人沒有真正改善，但他們的病情也沒有惡化；而治療有效的病人則有改善。

相比之下，這是一個未經(jīng)治療的病人在相同初始條件下的模擬軌跡。?

對CoA的模擬軌跡是比較現(xiàn)實的：它隨著年齡的增長而減慢。這與ANB不同，它的發(fā)生是因為CoA2同時取決于T1和T2。(ANB2則兩者都不依賴)。?

> for (i in seq(nrow(interv)) { + ? # 進行聯(lián)合預(yù)測，目前用predict()無法實現(xiàn)。 + ? dist(dyn.fitted, nodes = c(), + ? intervals[i,] = weighted.mean(ANB2, weights) + ? intervals[i,] = weighted.mean(CoA2, weights)

?

關(guān)于分析師

在此對Chang Zhang對本文所作的貢獻表示誠摯感謝，她在河海大學(xué)完成了信息管理專業(yè)學(xué)位，專注數(shù)據(jù)挖掘、數(shù)據(jù)分析、自然語言處理。

若您是個人學(xué)習(xí)者或培訓(xùn)負責(zé)人，希望學(xué)習(xí)拓端的線上/線下課程， 欲了解更多拓端學(xué)堂課程及信息，請點擊文末“閱讀原文”咨詢。

?最受歡迎的見解

1.matlab使用貝葉斯優(yōu)化的深度學(xué)習(xí)

2.matlab貝葉斯隱馬爾可夫hmm模型實現(xiàn)

3.R語言Gibbs抽樣的貝葉斯簡單線性回歸仿真

4.R語言中的block Gibbs吉布斯采樣貝葉斯多元線性回歸

5.R語言中的Stan概率編程MCMC采樣的貝葉斯模型

6.R語言貝葉斯Poisson泊松-正態(tài)分布模型分析職業(yè)足球比賽進球數(shù)

7.R語言使用貝葉斯 層次模型進行空間數(shù)據(jù)分析

8.R語言隨機搜索變量選擇SSVS估計貝葉斯向量自回歸（BVAR）模型

9.matlab貝葉斯隱馬爾可夫hmm模型實現(xiàn)