R語言中貝葉斯網(wǎng)絡(luò)（BN）、動態(tài)貝葉斯網(wǎng)絡(luò)、線性模型分析錯頜畸形數(shù)據(jù)|附代碼數(shù)據(jù)
全文鏈接：http://tecdat.cn/?p=22956
最近我們被客戶要求撰寫關(guān)于貝葉斯網(wǎng)絡(luò)的研究報告，包括一些圖形和統(tǒng)計輸出。
貝葉斯網(wǎng)絡(luò)（BN）是一種基于有向無環(huán)圖的概率模型，它描述了一組變量及其相互之間的條件依賴性。它是一個圖形模型，我們可以很容易地檢查變量的條件依賴性和它們在圖中的方向
在這篇文章中，我將簡要地學(xué)習(xí)如何用R來使用貝葉斯網(wǎng)絡(luò)。
本教程旨在介紹貝葉斯網(wǎng)絡(luò)學(xué)習(xí)和推理的基礎(chǔ)知識，使用真實世界的數(shù)據(jù)來探索圖形建模的典型數(shù)據(jù)分析工作流程。關(guān)鍵點將包括：
預(yù)處理數(shù)據(jù);
學(xué)習(xí)貝葉斯網(wǎng)絡(luò)的結(jié)構(gòu)和參數(shù)。
使用網(wǎng)絡(luò)作為預(yù)測模型。
使用網(wǎng)絡(luò)進(jìn)行推理。
通過與外部信息的對比來驗證網(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)的主要作用是通過圖形分離來表達(dá)模型中各變量之間的條件獨立性關(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ù)的矩陣。
點擊標(biāo)題查閱往期內(nèi)容
R語言BUGS/JAGS貝葉斯分析: 馬爾科夫鏈蒙特卡洛方法（MCMC）采樣
左右滑動查看更多
01
02
03
04
如上弧線顯示，在我們的數(shù)據(jù)中存在'類型'到'顏色'，以及'狀態(tài)'到'接受'和'元素'的關(guān)系。'類型'和'狀態(tài)'是兩個獨立的組，它們之間不存在相互依賴關(guān)系。
接下來，我們將用數(shù)據(jù)來擬合模型。
simd_fitted
基于上述訓(xùn)練數(shù)據(jù)，我們可以進(jìn)行條件概率查詢。
我們檢查 "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ù):
我們學(xué)習(xí)一個BN，并使用它來確定和可視化在成長和治療過程中各種III類錯位頜面特征之間的相互作用。
我們通過驗證一些普遍接受的關(guān)于這些骨骼不平衡演變的假說來檢驗其一致性。
我們表明，與接受快速上頜擴張和面罩治療的正畸患者相比，未經(jīng)治療的受試者形成了不同的III類顱面生長模式。
在接受治療的患者中，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。同樣地，我們對成長進(jìn)行重新編碼，0表示壞，1表示好。
由于我們將使用高斯BN進(jìn)行分析，檢查這些變量是否是正態(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ù)字。對差異進(jìn)行建模會導(dǎo)致局部分布，其形式為回歸模型
其中
?對于其他回歸因子，以此類推。我們可以將這種回歸改寫為
這是一組微分方程，對變化率進(jìn)行建模，其關(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()與默認(rèn)分?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來解決這兩個問題。
使用bootstrap對數(shù)據(jù)重新取樣。
從每個bootstrap樣本中學(xué)習(xí)一個單獨的網(wǎng)絡(luò)。
檢查每個可能的弧在網(wǎng)絡(luò)中出現(xiàn)的頻率。
用出現(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）進(jìn)行比較？最定性的方法是將兩個網(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é)合來進(jìn)行原則性的比較，如果兩個弧被唯一確定為不同，我們就說它們是不同的。
也可以看一下關(guān)于弧長分布的閾值：平均的網(wǎng)絡(luò)是相當(dāng)密集的（9個節(jié)點有17個?。茈y閱讀。
>?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（藍(lá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
我們通過在 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，很容易確認(rèn)這是事實。
?>?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。
只看網(wǎng)絡(luò)結(jié)構(gòu)：如果學(xué)習(xí)BN的主要目標(biāo)是識別弧和路徑，當(dāng)BN被解釋為因果模型時，通常是這種情況，我們可以進(jìn)行本質(zhì)上的路徑分析和研究弧的強度。
將BN視為一個整體，包括參數(shù)：如果學(xué)習(xí)BN的主要目標(biāo)是將其作為一個專家模型，那么我們可能想。
根據(jù)其他一些變量的值，預(yù)測新個體的一個或多個變量的值；以及
將CP查詢的結(jié)果與專家知識進(jìn)行比較，以確認(rèn)BN反映了關(guān)于我們正在建模的現(xiàn)象的最佳知識。
預(yù)測準(zhǔn)確性
我們可以用通常的方法來衡量我們所選擇的學(xué)習(xí)策略的預(yù)測準(zhǔn)確性，即交叉驗證。實現(xiàn)了：
k-fold交叉驗證;
指定的k進(jìn)行交叉驗證;
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ù)的變體，它使用從所有其他變量計算的后驗期望值進(jìn)行預(yù)測?；镜膿p失函數(shù)（cor, mse, pred）僅僅從它們的父代來預(yù)測每個節(jié)點的值，這在處理很少或沒有父代的節(jié)點時是沒有意義的。
用專家知識進(jìn)行確認(rèn)
確認(rèn)BN是否有意義的另一種方法是把它當(dāng)作的工作模型，看看它是否表達(dá)了關(guān)于關(guān)鍵事實，這些事實在學(xué)習(xí)過程中沒有作為先驗知識使用。否則，我們將只是拿回我們放在先驗中的信息）。一些例子。
"CoGo的過度增長應(yīng)該會引起PPPM的減少"。
我們通過為存儲在 fitted.simpler中的BN生成dCoGo和dPPPM的樣本，并假設(shè)沒有發(fā)生任何處理，來測試這個假設(shè)。隨著dCoGo的增加（這表明增長越來越快），DPPPM變得越來越負(fù)（這表明假設(shè)角度最初是正的，則角度會減少。>?sim?=?dist(fitted.simpler) >?plot(sim?) >?abline(v?=?0,?col?=?2,?lty?=?2,?lwd?=?2)
"CoGo的小幅增長應(yīng)該會引起PPPM的增長。"
從上圖來看，CoGo的負(fù)增長或空增長（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
"如果ANB減少，IMPA就會減少以進(jìn)行補償。"
像以前一樣通過模擬測試，我們正在尋找與IMPA（相同）的負(fù)值相關(guān)的dANB的負(fù)值（這表明假設(shè)角度最初是正的，就會減少）。從下圖中可以看出，dANB與dIMPA成正比，所以其中一個的減少表明另一個的減少；兩者的平均趨勢（黑線）同時為負(fù)。
>?plot(sim?) >?abline(coef(lm(dIMPA?~?dANB?))
"如果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
_"治療試圖阻止ANB的減少。如果我們固定ANB，治療過的病人和未治療過的病人是否有區(qū)別？"
首先，我們可以檢查在沒有任何干預(yù)的情況下，dANB≈0的病人的治療和增長之間的關(guān)系（即使用我們從數(shù)據(jù)中得知的BN）。_``` dist(fitted?) table(TREATMENT?=?Treatment?<?0.5,?GOOD.GROWTH?=??Growth?>?0.5) ```
估計的P（GOOD.GROWTH ∣ TREATMENT）對于接受治療和未接受治療的病人是不同的（≈0.65對≈0.52）。
如果我們模擬一個正式的干預(yù)措施（如Judea Pearl），并從外部設(shè)置dANB=0（從而使其獨立于其父母，并刪除相應(yīng)的?。覀兙蜁l(fā)現(xiàn)GOOD.GROWTH對于接受治療和未接受治療的病人來說實際上具有相同的分布，從而變得與TREATMENT無關(guān)。這表明，有利的預(yù)后確實是由防止ANB的變化決定的，而治療的其他成分（如果有的話）就變得不重要了。 ``` table(TREATMENT?=??Treatment?<?0.5,?GOOD.GROWTH?=??Growth?>?0.5) ```
_"治療試圖阻止ANB的減少。如果我們固定ANB，治療和未治療的病人之間是否有區(qū)別？"
評估的方法之一是檢查在保持GoPg固定的情況下，A點和B點之間的角度（ANB）是否在治療和未治療的病人之間發(fā)生變化。_
假設(shè)GoPg不發(fā)生變化，對于接受治療的病人來說，A點和B點之間的角度會增加（強烈的負(fù)值表示水平不平衡，所以正的變化率表示不平衡的減少），而對于未接受治療的病人來說則會減少（不平衡會隨著時間慢慢惡化）。
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ù)。然而，我們選擇使用治療變量而不是生長變量作為變量來表達(dá)受試者可能正在接受醫(yī)療的事實。原因是生長變量是一個衡量第二次測量時的預(yù)后的變量，它的值在第一次測量時是未知的；而治療變量在兩次測量時都是相同的。
學(xué)習(xí)結(jié)構(gòu)
首先，我們將變量分為三組：時間為t2的變量，時間為t1=t2-1的變量，以及與時間無關(guān)的變量，因為它們在t1和t1取值相同。
>?t2.variables
然后我們引入一個黑名單，其中。
我們將所有從臨床變量到T1、T2和治療的弧線列入黑名單，因為我們知道，年齡和治療不是由臨床測量決定的。
我們將所有進(jìn)入Treatment和t1時間段的所有變量的弧列入黑名單，因為我們假設(shè)t1時間段的變量之間的弧與t2時間段的相應(yīng)變量是一樣的，兩次學(xué)習(xí)它們是沒有意義的。
我們將所有從t2到t1的弧列入黑名單。
grid(from?=?setdiff(names(ortho),?c("T1",?"T2")),
?to?=?c("T1",?"T2"))
相比之下，我們只對T1→T2的弧線進(jìn)行白名單，因為第二次測量的年齡顯然取決于第一次測量的年齡。
>??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)
然后，我們對變量進(jìn)行分組，以方便區(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?fn26??1??1
學(xué)習(xí)參數(shù)
由于Treatment是一個離散變量，BN是一個CLGBN。這意味著以Treatment為父節(jié)點的連續(xù)節(jié)點的參數(shù)化與其他節(jié)點不同。
fit(dynavg)
我們可以看到，ANB2取決于ANB（所以，在前一個時間點的同一變量）和治療。ANB是連續(xù)的，所以它被用作ANB2的回歸因子。?治療變量是離散的，決定了線性回歸的成分。
模型驗證和推理
我們可以對這個新模型提出另一組問題
"在不同的治療下，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)
```
density(~?diff,?groups?=?treatment)
```我們知道，治療試圖阻止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))?{
+???#?進(jìn)行聯(lián)合預(yù)測，目前用predict()無法實現(xiàn)。
+???dist(dyn.fitted,?nodes?=?c(),
+???intervals[i,]?=?weighted.mean(ANB2,?weights)
+???intervals[i,]?=?weighted.mean(CoA2,?weights)
本文摘選?《?R語言中貝葉斯網(wǎng)絡(luò)（BN）、動態(tài)貝葉斯網(wǎng)絡(luò)、線性模型分析錯頜畸形數(shù)據(jù)?》?，點擊“閱讀原文”獲取全文完整資料。
點擊標(biāo)題查閱往期內(nèi)容
使用貝葉斯層次模型進(jìn)行空間數(shù)據(jù)分析MCMC的rstan貝葉斯回歸模型和標(biāo)準(zhǔn)線性回歸模型比較
python貝葉斯隨機過程：馬爾可夫鏈Markov-Chain，MC和Metropolis-Hastings，MH采樣算法可視化
Python貝葉斯推斷Metropolis-Hastings（M-H）MCMC采樣算法的實現(xiàn)
matlab貝葉斯隱馬爾可夫hmm模型實現(xiàn)
貝葉斯線性回歸和多元線性回歸構(gòu)建工資預(yù)測模型
Metropolis Hastings采樣和貝葉斯泊松回歸Poisson模型
貝葉斯分位數(shù)回歸、lasso和自適應(yīng)lasso貝葉斯分位數(shù)回歸分析免疫球蛋白、前列腺癌數(shù)據(jù)
R語言RSTAN MCMC：NUTS采樣算法用LASSO 構(gòu)建貝葉斯線性回歸模型分析職業(yè)聲望數(shù)據(jù)
R語言STAN貝葉斯線性回歸模型分析氣候變化影響北半球海冰范圍和可視化檢查模型收斂性
PYTHON用戶流失數(shù)據(jù)挖掘：建立邏輯回歸、XGBOOST、隨機森林、決策樹、支持向量機、樸素貝葉斯和KMEANS聚類用戶畫像
貝葉斯分位數(shù)回歸、lasso和自適應(yīng)lasso貝葉斯分位數(shù)回歸分析免疫球蛋白、前列腺癌數(shù)據(jù)R語言JAGS貝葉斯回歸模型分析博士生延期畢業(yè)完成論文時間
R語言Metropolis Hastings采樣和貝葉斯泊松回歸Poisson模型
Python決策樹、隨機森林、樸素貝葉斯、KNN（K-最近鄰居）分類分析銀行拉新活動挖掘潛在貸款客戶
R語言貝葉斯MCMC：用rstan建立線性回歸模型分析汽車數(shù)據(jù)和可視化診斷
R語言貝葉斯MCMC：GLM邏輯回歸、Rstan線性回歸、Metropolis Hastings與Gibbs采樣算法實例
R語言貝葉斯Poisson泊松-正態(tài)分布模型分析職業(yè)足球比賽進(jìn)球數(shù)
隨機森林優(yōu)化貝葉斯預(yù)測分析汽車燃油經(jīng)濟性
R語言邏輯回歸、Naive Bayes貝葉斯、決策樹、隨機森林算法預(yù)測心臟病
R語言用Rcpp加速Metropolis-Hastings抽樣估計貝葉斯邏輯回歸模型的參數(shù)
R語言中的block Gibbs吉布斯采樣貝葉斯多元線性回歸
Python貝葉斯回歸分析住房負(fù)擔(dān)能力數(shù)據(jù)集
R語言實現(xiàn)貝葉斯分位數(shù)回歸、lasso和自適應(yīng)lasso貝葉斯分位數(shù)回歸分析
Python用PyMC3實現(xiàn)貝葉斯線性回歸模型
R語言用WinBUGS 軟件對學(xué)術(shù)能力測驗建立層次（分層）貝葉斯模型
R語言Gibbs抽樣的貝葉斯簡單線性回歸仿真分析
R語言和STAN,JAGS：用RSTAN,RJAG建立貝葉斯多元線性回歸預(yù)測選舉數(shù)據(jù)
R語言基于copula的貝葉斯分層混合模型的診斷準(zhǔn)確性研究
R語言貝葉斯線性回歸和多元線性回歸構(gòu)建工資預(yù)測模型
R語言貝葉斯推斷與MCMC：實現(xiàn)Metropolis-Hastings 采樣算法示例
R語言stan進(jìn)行基于貝葉斯推斷的回歸模型
R語言中RStan貝葉斯層次模型分析示例
R語言使用Metropolis-Hastings采樣算法自適應(yīng)貝葉斯估計與可視化
R語言隨機搜索變量選擇SSVS估計貝葉斯向量自回歸（BVAR）模型
WinBUGS對多元隨機波動率模型：貝葉斯估計與模型比較
R語言實現(xiàn)MCMC中的Metropolis–Hastings算法與吉布斯采樣
R語言貝葉斯推斷與MCMC：實現(xiàn)Metropolis-Hastings 采樣算法示例
R語言使用Metropolis-Hastings采樣算法自適應(yīng)貝葉斯估計與可視化
視頻：R語言中的Stan概率編程MCMC采樣的貝葉斯模型
R語言MCMC:Metropolis-Hastings采樣用于回歸的貝葉斯估計
