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

最近我們被客戶(hù)要求撰寫(xiě)關(guān)于貝葉斯網(wǎng)絡(luò)的研究報(bào)告，包括一些圖形和統(tǒng)計(jì)輸出。

貝葉斯網(wǎng)絡(luò)（BN）是一種基于有向無(wú)環(huán)圖的概率模型，它描述了一組變量及其相互之間的條件依賴(lài)性。它是一個(gè)圖形模型，我們可以很容易地檢查變量的條件依賴(lài)性和它們?cè)趫D中的方向

在這篇文章中，我將簡(jiǎn)要地學(xué)習(xí)如何用R來(lái)使用貝葉斯網(wǎng)絡(luò)。

本教程旨在介紹貝葉斯網(wǎng)絡(luò)學(xué)習(xí)和推理的基礎(chǔ)知識(shí)，使用真實(shí)世界的數(shù)據(jù)來(lái)探索圖形建模的典型數(shù)據(jù)分析工作流程。關(guān)鍵點(diǎn)將包括：

• 預(yù)處理數(shù)據(jù);

• 學(xué)習(xí)貝葉斯網(wǎng)絡(luò)的結(jié)構(gòu)和參數(shù)。

• 使用網(wǎng)絡(luò)作為預(yù)測(cè)模型。

• 使用網(wǎng)絡(luò)進(jìn)行推理。

• 通過(guò)與外部信息的對(duì)比來(lái)驗(yàn)證網(wǎng)絡(luò)的有效性。

快速介紹

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

定義

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

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

• , 其中每個(gè)節(jié)點(diǎn)?

• ?對(duì)應(yīng)于一個(gè)隨機(jī)變量?

• ;

• 一個(gè)全局概率分布?

• ?(帶參數(shù)?

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

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

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

; 而??

?要比

小得多，因?yàn)樵S多參數(shù)是固定的，因?yàn)樗鼈兯鶎俚淖兞渴仟?dú)立的。

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

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

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

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

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

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

>?plot(hc_simd)

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

點(diǎn)擊標(biāo)題查閱往期內(nèi)容

R語(yǔ)言BUGS/JAGS貝葉斯分析: 馬爾科夫鏈蒙特卡洛方法（MCMC）采樣

左右滑動(dòng)查看更多

01

02

03

04

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

simd_fitted

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

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

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

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

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

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

任務(wù):

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

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

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

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

數(shù)據(jù)

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

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

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

• ANB：唐氏點(diǎn)A和B之間的角度（度）。

• IMPA：門(mén)牙-下頜平面角（度）。

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

• CoA：上頜骨從髁狀突到唐氏點(diǎn)A的總長(zhǎng)度（mm）。

• GoPg：下頜體從齒齦到齒齦的長(zhǎng)度（mm）。

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

所有的測(cè)量都是通過(guò)X射線(xiàn)掃描得出的，使用一套參考點(diǎn)建立的圖，如以下。

>?str(data)

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

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

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

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

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

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

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

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

最后，我們可以看看這些變量是否以任何方式聚在一起，因?yàn)榫墼谝黄鸬淖兞扛锌赡茉贐N中發(fā)生聯(lián)系。

>?heatmap(rho)

我們可以在熱圖中看到兩個(gè)集群：第一個(gè)集群包括dCoGo、dGoPg和dCoA，第二個(gè)集群包括Treatment、dANB和dCoA。第一個(gè)聚類(lèi)在臨床上很有意思，因?yàn)樗ㄖ委熀蛢蓚€(gè)都與唐氏A點(diǎn)有關(guān)的變量，這為治療的主要效果提供了一些線(xiàn)索。

plot(ug?)

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

在這里，我們使用保存在diff中的差異來(lái)為數(shù)據(jù)建模，而不是原始值；我們將使用GBN處理，因?yàn)樗凶兞慷际菙?shù)字。對(duì)差異進(jìn)行建模會(huì)導(dǎo)致局部分布，其形式為回歸模型

其中

?對(duì)于其他回歸因子，以此類(lèi)推。我們可以將這種回歸改寫(xiě)為

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

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

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

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

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

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

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

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

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

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

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

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

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

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

一個(gè)簡(jiǎn)單的學(xué)習(xí)?

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

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

plot(dag,?,?highlight?)

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

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

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

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

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

booth(diff,?R?=?200)

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

attr(?"threshold")

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

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

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

>?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）進(jìn)行比較？最定性的方法是將兩個(gè)網(wǎng)絡(luò)并排繪制，節(jié)點(diǎn)位置相同，并突出顯示一個(gè)網(wǎng)絡(luò)中出現(xiàn)而另一個(gè)網(wǎng)絡(luò)中沒(méi)有的弧，或者出現(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)只出現(xiàn)在平均網(wǎng)絡(luò)中，而dPPPM→dANB只出現(xiàn)在我們從所有數(shù)據(jù)中學(xué)到的網(wǎng)絡(luò)中。我們可以假設(shè)，前三個(gè)弧被數(shù)據(jù)的噪聲加上小樣本量和偏離常態(tài)的情況所隱藏。編程可以返回真陽(yáng)性（出現(xiàn)在兩個(gè)網(wǎng)絡(luò)中的?。┖图訇?yáng)性/陰性（只出現(xiàn)在兩個(gè)網(wǎng)絡(luò)中的一個(gè)的弧）的數(shù)量。

>?compare

或弧=TRUE。

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

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

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

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

因此，把閾值提高一點(diǎn)，多剔除幾個(gè)弧就好了。看一下上面的圖，由于弧長(zhǎng)分布的差距，較高的閾值的兩個(gè)自然選擇是0.75（紅色虛線(xiàn)）和0.85（藍(lán)色虛線(xiàn)）。

>?nrow(?strength?>??"threshold"?direction?>?0.5,?])[1]?18?trength?>?0.75?&??direction?>?0.5?[1]?15?strength?>?0.85?&??direction?>?0.5?[1]?12

我們通過(guò)在 network()中設(shè)置閾值=0.85得到的更簡(jiǎn)單的網(wǎng)絡(luò)如下所示；從定性的角度來(lái)看，它當(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來(lái)建模。因此，如果我們使用最大似然估計(jì)來(lái)擬合網(wǎng)絡(luò)的參數(shù)，我們就會(huì)發(fā)現(xiàn)每個(gè)局部分布都是一個(gè)典型的線(xiàn)性回歸。

fit(avg,?diff)

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

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

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

?和

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

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

比較線(xiàn)性模型?

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

?dANB ?dANB?=?penalized(?dANB) ?dANB

模型驗(yàn)證

有兩種主要的方法來(lái)驗(yàn)證一個(gè)BN。

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

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

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

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

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

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

• k-fold交叉驗(yàn)證;

• 指定的k進(jìn)行交叉驗(yàn)證;

• hold-out?交叉驗(yàn)證

對(duì)于:

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

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

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

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

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

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

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

用專(zhuān)家知識(shí)進(jìn)行確認(rèn)

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

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

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

1. "CoGo的小幅增長(zhǎng)應(yīng)該會(huì)引起PPPM的增長(zhǎng)。"

   從上圖來(lái)看，CoGo的負(fù)增長(zhǎng)或空增長(zhǎng)（dCoGo ? 0）對(duì)應(yīng)于PPPM的正增長(zhǎng)，概率為≈0.60。對(duì)于CoGo的小幅增長(zhǎng)（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

1. "如果ANB減少，IMPA就會(huì)減少以進(jìn)行補(bǔ)償。"

   像以前一樣通過(guò)模擬測(cè)試，我們正在尋找與IMPA（相同）的負(fù)值相關(guān)的dANB的負(fù)值（這表明假設(shè)角度最初是正的，就會(huì)減少）。從下圖中可以看出，dANB與dIMPA成正比，所以其中一個(gè)的減少表明另一個(gè)的減少；兩者的平均趨勢(shì)（黑線(xiàn)）同時(shí)為負(fù)。

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

1. "如果GoPg強(qiáng)烈增加，那么ANB和IMPA都會(huì)減少。" 如果我們從BN中模擬dGoPg、dANB和dIMPA，假設(shè)dGoPg>5（即GoPg在增加），我們估計(jì)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

1. _"治療試圖阻止ANB的減少。如果我們固定ANB，治療過(guò)的病人和未治療過(guò)的病人是否有區(qū)別？"

   首先，我們可以檢查在沒(méi)有任何干預(yù)的情況下，dANB≈0的病人的治療和增長(zhǎng)之間的關(guān)系（即使用我們從數(shù)據(jù)中得知的BN）。_``` dist(fitted?) table(TREATMENT?=?Treatment?<?0.5,?GOOD.GROWTH?=??Growth?>?0.5) ```

1. 估計(jì)的P（GOOD.GROWTH ∣ TREATMENT）對(duì)于接受治療和未接受治療的病人是不同的（≈0.65對(duì)≈0.52）。

   如果我們模擬一個(gè)正式的干預(yù)措施（如Judea Pearl），并從外部設(shè)置dANB=0（從而使其獨(dú)立于其父母，并刪除相應(yīng)的?。?，我們就會(huì)發(fā)現(xiàn)GOOD.GROWTH對(duì)于接受治療和未接受治療的病人來(lái)說(shuō)實(shí)際上具有相同的分布，從而變得與TREATMENT無(wú)關(guān)。這表明，有利的預(yù)后確實(shí)是由防止ANB的變化決定的，而治療的其他成分（如果有的話(huà)）就變得不重要了。 ``` table(TREATMENT?=??Treatment?<?0.5,?GOOD.GROWTH?=??Growth?>?0.5) ```

1. _"治療試圖阻止ANB的減少。如果我們固定ANB，治療和未治療的病人之間是否有區(qū)別？"

   評(píng)估的方法之一是檢查在保持GoPg固定的情況下，A點(diǎn)和B點(diǎn)之間的角度（ANB）是否在治療和未治療的病人之間發(fā)生變化。_

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

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

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

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

我們用于這個(gè)模型的數(shù)據(jù)是我們?cè)诜治鲩_(kāi)始時(shí)存儲(chǔ)到正交的原始數(shù)據(jù)。然而，我們選擇使用治療變量而不是生長(zhǎng)變量作為變量來(lái)表達(dá)受試者可能正在接受醫(yī)療的事實(shí)。原因是生長(zhǎng)變量是一個(gè)衡量第二次測(cè)量時(shí)的預(yù)后的變量，它的值在第一次測(cè)量時(shí)是未知的；而治療變量在兩次測(cè)量時(shí)都是相同的。

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

首先，我們將變量分為三組：時(shí)間為t2的變量，時(shí)間為t1=t2-1的變量，以及與時(shí)間無(wú)關(guān)的變量，因?yàn)樗鼈冊(cè)趖1和t1取值相同。

>?t2.variables

然后我們引入一個(gè)黑名單，其中。

1. 我們將所有從臨床變量到T1、T2和治療的弧線(xiàn)列入黑名單，因?yàn)槲覀冎溃挲g和治療不是由臨床測(cè)量決定的。

2. 我們將所有進(jìn)入Treatment和t1時(shí)間段的所有變量的弧列入黑名單，因?yàn)槲覀兗僭O(shè)t1時(shí)間段的變量之間的弧與t2時(shí)間段的相應(yīng)變量是一樣的，兩次學(xué)習(xí)它們是沒(méi)有意義的。

3. 我們將所有從t2到t1的弧列入黑名單。

grid(from?=?setdiff(names(ortho),?c("T1",?"T2")), ?to?=?c("T1",?"T2"))

相比之下，我們只對(duì)T1→T2的弧線(xiàn)進(jìn)行白名單，因?yàn)榈诙螠y(cè)量的年齡顯然取決于第一次測(cè)量的年齡。

>??data.frame(from?=?c("T1"),?to?=?c("T2"))

最后我們可以用bl和wl來(lái)學(xué)習(xí)BN的結(jié)構(gòu)。

>?dyn.dag

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

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

plot(dyn,?render?=?FALSE)

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

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

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

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

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

>?boot?(ortho?) >?plot(dyn)

avernet(dyn.str)

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

compare(dag,?avg) tp?fp?fn26??1??1

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

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

fit(dynavg)

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

模型驗(yàn)證和推理

我們可以對(duì)這個(gè)新模型提出另一組問(wèn)題

1. "在不同的治療下，ANB從第一次測(cè)量到第二次測(cè)量的轉(zhuǎn)變程度如何？"
   我們可以用cpdist()生成一對(duì)(ANB, ANB2)，條件是治療方法等于NT、TB和TG，并觀(guān)察其分布。

   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的左邊這一事實(shí)相一致。未經(jīng)治療的病人病情繼續(xù)惡化；治療無(wú)效的病人沒(méi)有真正改善，但他們的病情也沒(méi)有惡化；而治療有效的病人則有改善。

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

對(duì)CoA的模擬軌跡是比較現(xiàn)實(shí)的：它隨著年齡的增長(zhǎng)而減慢。這與ANB不同，它的發(fā)生是因?yàn)镃oA2同時(shí)取決于T1和T2。(ANB2則兩者都不依賴(lài))。

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

本文摘選?《?R語(yǔ)言中貝葉斯網(wǎng)絡(luò)（BN）、動(dòng)態(tài)貝葉斯網(wǎng)絡(luò)、線(xiàn)性模型分析錯(cuò)頜畸形數(shù)據(jù)?》?，點(diǎn)擊“閱讀原文”獲取全文完整資料。

點(diǎn)擊標(biāo)題查閱往期內(nèi)容

使用貝葉斯層次模型進(jìn)行空間數(shù)據(jù)分析MCMC的rstan貝葉斯回歸模型和標(biāo)準(zhǔn)線(xiàn)性回歸模型比較
python貝葉斯隨機(jī)過(guò)程：馬爾可夫鏈Markov-Chain，MC和Metropolis-Hastings，MH采樣算法可視化
Python貝葉斯推斷Metropolis-Hastings（M-H）MCMC采樣算法的實(shí)現(xiàn)
matlab貝葉斯隱馬爾可夫hmm模型實(shí)現(xiàn)
貝葉斯線(xiàn)性回歸和多元線(xiàn)性回歸構(gòu)建工資預(yù)測(cè)模型
Metropolis Hastings采樣和貝葉斯泊松回歸Poisson模型
貝葉斯分位數(shù)回歸、lasso和自適應(yīng)lasso貝葉斯分位數(shù)回歸分析免疫球蛋白、前列腺癌數(shù)據(jù)
R語(yǔ)言RSTAN MCMC：NUTS采樣算法用LASSO 構(gòu)建貝葉斯線(xiàn)性回歸模型分析職業(yè)聲望數(shù)據(jù)
R語(yǔ)言STAN貝葉斯線(xiàn)性回歸模型分析氣候變化影響北半球海冰范圍和可視化檢查模型收斂性
PYTHON用戶(hù)流失數(shù)據(jù)挖掘：建立邏輯回歸、XGBOOST、隨機(jī)森林、決策樹(shù)、支持向量機(jī)、樸素貝葉斯和KMEANS聚類(lèi)用戶(hù)畫(huà)像
貝葉斯分位數(shù)回歸、lasso和自適應(yīng)lasso貝葉斯分位數(shù)回歸分析免疫球蛋白、前列腺癌數(shù)據(jù)R語(yǔ)言JAGS貝葉斯回歸模型分析博士生延期畢業(yè)完成論文時(shí)間
R語(yǔ)言Metropolis Hastings采樣和貝葉斯泊松回歸Poisson模型
Python決策樹(shù)、隨機(jī)森林、樸素貝葉斯、KNN（K-最近鄰居）分類(lèi)分析銀行拉新活動(dòng)挖掘潛在貸款客戶(hù)
R語(yǔ)言貝葉斯MCMC：用rstan建立線(xiàn)性回歸模型分析汽車(chē)數(shù)據(jù)和可視化診斷
R語(yǔ)言貝葉斯MCMC：GLM邏輯回歸、Rstan線(xiàn)性回歸、Metropolis Hastings與Gibbs采樣算法實(shí)例
R語(yǔ)言貝葉斯Poisson泊松-正態(tài)分布模型分析職業(yè)足球比賽進(jìn)球數(shù)
隨機(jī)森林優(yōu)化貝葉斯預(yù)測(cè)分析汽車(chē)燃油經(jīng)濟(jì)性
R語(yǔ)言邏輯回歸、Naive Bayes貝葉斯、決策樹(shù)、隨機(jī)森林算法預(yù)測(cè)心臟病
R語(yǔ)言用Rcpp加速M(fèi)etropolis-Hastings抽樣估計(jì)貝葉斯邏輯回歸模型的參數(shù)
R語(yǔ)言中的block Gibbs吉布斯采樣貝葉斯多元線(xiàn)性回歸
Python貝葉斯回歸分析住房負(fù)擔(dān)能力數(shù)據(jù)集
R語(yǔ)言實(shí)現(xiàn)貝葉斯分位數(shù)回歸、lasso和自適應(yīng)lasso貝葉斯分位數(shù)回歸分析
Python用PyMC3實(shí)現(xiàn)貝葉斯線(xiàn)性回歸模型
R語(yǔ)言用WinBUGS 軟件對(duì)學(xué)術(shù)能力測(cè)驗(yàn)建立層次（分層）貝葉斯模型
R語(yǔ)言Gibbs抽樣的貝葉斯簡(jiǎn)單線(xiàn)性回歸仿真分析
R語(yǔ)言和STAN,JAGS：用RSTAN,RJAG建立貝葉斯多元線(xiàn)性回歸預(yù)測(cè)選舉數(shù)據(jù)
R語(yǔ)言基于copula的貝葉斯分層混合模型的診斷準(zhǔn)確性研究
R語(yǔ)言貝葉斯線(xiàn)性回歸和多元線(xiàn)性回歸構(gòu)建工資預(yù)測(cè)模型
R語(yǔ)言貝葉斯推斷與MCMC：實(shí)現(xiàn)Metropolis-Hastings 采樣算法示例
R語(yǔ)言stan進(jìn)行基于貝葉斯推斷的回歸模型
R語(yǔ)言中RStan貝葉斯層次模型分析示例
R語(yǔ)言使用Metropolis-Hastings采樣算法自適應(yīng)貝葉斯估計(jì)與可視化
R語(yǔ)言隨機(jī)搜索變量選擇SSVS估計(jì)貝葉斯向量自回歸（BVAR）模型
WinBUGS對(duì)多元隨機(jī)波動(dòng)率模型：貝葉斯估計(jì)與模型比較
R語(yǔ)言實(shí)現(xiàn)MCMC中的Metropolis–Hastings算法與吉布斯采樣
R語(yǔ)言貝葉斯推斷與MCMC：實(shí)現(xiàn)Metropolis-Hastings 采樣算法示例
R語(yǔ)言使用Metropolis-Hastings采樣算法自適應(yīng)貝葉斯估計(jì)與可視化
視頻：R語(yǔ)言中的Stan概率編程MCMC采樣的貝葉斯模型
R語(yǔ)言MCMC:Metropolis-Hastings采樣用于回歸的貝葉斯估計(jì)