原文鏈接：?http://tecdat.cn/?p=23652

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

本文為讀者提供了如何進(jìn)行貝葉斯回歸的基本教程。包括完成導(dǎo)入數(shù)據(jù)文件、探索匯總統(tǒng)計(jì)和回歸分析。

在本文中，我們首先使用軟件的默認(rèn)先驗(yàn)設(shè)置。在第二步中，我們將應(yīng)用用戶指定的先驗(yàn)，對(duì)自己的數(shù)據(jù)使用貝葉斯。

準(zhǔn)備工作

本教程要求:

• 已安裝的JAGS

• 安裝R軟件。

• 假設(shè)檢驗(yàn)的基本知識(shí)

• 相關(guān)性和回歸的基本知識(shí)

• 貝葉斯推理的基本知識(shí)

• R語言編碼的基本知識(shí)

數(shù)據(jù)實(shí)例

我們?cè)谶@個(gè)練習(xí)中使用的數(shù)據(jù)是基于一項(xiàng)關(guān)于預(yù)測博士生完成論文時(shí)間的研究（Van de Schoot, Yerkes, Mouw and Sonneveld 2013）。研究人員詢問了博士生完成他們的博士論文需要多長時(shí)間（n=333）。結(jié)果顯示，博士學(xué)位獲得者平均花了59.8個(gè)月（5年4個(gè)月）來完成他們的博士學(xué)位。變量B3衡量計(jì)劃和實(shí)際項(xiàng)目時(shí)間之間的差異，以月為單位（平均=9.97，最小=-31，最大=91，sd=14.43）。

對(duì)于目前的工作，我們感興趣的問題是，博士學(xué)位獲得者的年齡（M=31.7，SD=6.86）是否與他們項(xiàng)目的延期有關(guān)。

預(yù)計(jì)完成時(shí)間和年齡之間的關(guān)系是非線性的。這可能是由于在人生的某個(gè)階段（即三十多歲），家庭生活比你在二十多歲時(shí)或年長時(shí)占用了你更多的時(shí)間。

因此，在我們的模型中，差距（B3）是因變量，年齡和年齡平方是預(yù)測因素。

問題：請(qǐng)寫出零假設(shè)和備擇假設(shè)。?寫下代表這個(gè)問題的無效假設(shè)和備選假設(shè)。你認(rèn)為哪個(gè)假設(shè)更有可能？

回答：

H0:年齡與博士項(xiàng)目的延期無關(guān)。

H1:?年齡與博士項(xiàng)目的延期有關(guān)。

H0:age2與博士項(xiàng)目的延期無關(guān)。

H1:age2與博士項(xiàng)目的延期有關(guān)。

準(zhǔn)備--導(dǎo)入和探索數(shù)據(jù)

數(shù)據(jù)是一個(gè).csv文件，但你可以使用以下語法直接將其加載到R中。

一旦你加載了你的數(shù)據(jù)，建議你檢查一下你的數(shù)據(jù)導(dǎo)入是否順利。因此，首先看看你的數(shù)據(jù)的匯總統(tǒng)計(jì)。你可以使用describe()函數(shù)。

問題：?你所有的數(shù)據(jù)都被正確地載入了嗎？也就是說，所有的數(shù)據(jù)點(diǎn)都有實(shí)質(zhì)性的意義嗎？

回答：

描述性統(tǒng)計(jì)有意義。

差異。平均值（9.97），SE（0.79）。

年齡。平均值（31.68），SE（0.38）。

age2。平均值（1050.22），SE（35.97）。

繪圖

在繼續(xù)分析數(shù)據(jù)之前，我們還可以繪制期望的關(guān)系。

回歸

在這個(gè)練習(xí)中，你將研究博士生的年齡和age2對(duì)他們的項(xiàng)目時(shí)間延期的影響，這作為結(jié)果變量使用回歸分析。

如你所知，貝葉斯推理包括將先驗(yàn)分布與從數(shù)據(jù)中獲得的似然性相結(jié)合。指定先驗(yàn)分布是貝葉斯推斷中最關(guān)鍵的一點(diǎn)，應(yīng)該受到高度重視（例如Van de Schoot等人，2017）。在本教程中，我們將首先依賴默認(rèn)的先驗(yàn)設(shè)置。

要用運(yùn)行多元回歸，首先要指定模型，然后擬合模型，最后獲得總結(jié)。模型的指定方法如下。

1. 我們想要預(yù)測的因變量。

2. "~"，我們用它來表示我們現(xiàn)在給其他感興趣的變量。(相當(dāng)于回歸方程的"="）。

3. 用求和符號(hào)'+'分隔的不同自變量。

4. 最后，我們插入因變量有一個(gè)方差，有一個(gè)截距。

下面的代碼是如何指定回歸模型的。

然后，我們需要使用以下代碼來擬合模型。我們指定target = "jags "來使用Jags而不是Stan編譯器。

現(xiàn)在我們用summary(fit.bayes)來看看總結(jié)。

顯示輸出

頻率主義模型與貝葉斯分析模型所提供的結(jié)果確實(shí)不同。

貝葉斯統(tǒng)計(jì)推斷和頻率主義統(tǒng)計(jì)方法之間的關(guān)鍵區(qū)別在于估計(jì)的未知參數(shù)的性質(zhì)。在頻率主義框架中，一個(gè)感興趣的參數(shù)被假定為未知的，但卻是固定的。也就是說，假設(shè)在人口中只有一個(gè)真實(shí)的人口參數(shù)，例如，一個(gè)真實(shí)的平均值或一個(gè)真實(shí)的回歸系數(shù)。在貝葉斯的主觀概率觀中，所有的未知參數(shù)都被視為不確定的，因此要用一個(gè)概率分布來描述。每個(gè)參數(shù)都是未知的，而所有未知的東西都會(huì)得到一個(gè)分布。

這就是為什么在頻率推斷中，你主要得到的是一個(gè)未知但固定的群體參數(shù)的點(diǎn)估計(jì)。這是一個(gè)參數(shù)值，考慮到數(shù)據(jù)，它最有可能出現(xiàn)在人群中。附帶的置信區(qū)間試圖讓你進(jìn)一步了解這個(gè)估計(jì)值的不確定性。重要的是要認(rèn)識(shí)到，置信區(qū)間只是構(gòu)成一個(gè)模擬量。在從人口中抽取的無限多的樣本中，構(gòu)建（95%）置信區(qū)間的程序?qū)⑹蛊湓?5%的時(shí)間內(nèi)包含真實(shí)的人口值。這并沒有為你提供任何信息，即人口參數(shù)位于你所分析的非常具體和唯一的樣本中的置信區(qū)間邊界內(nèi)的可能性有多大。

在貝葉斯分析中，你推斷的關(guān)鍵是感興趣的參數(shù)的后驗(yàn)分布。它滿足了概率分布的每一個(gè)屬性，并量化了人口參數(shù)位于某些區(qū)域的概率。一方面，你可以通過它的模式來描述后驗(yàn)的特點(diǎn)。這是一個(gè)參數(shù)值，考慮到數(shù)據(jù)和它的先驗(yàn)概率，它在人群中是最有可能的。另外，你也可以使用后驗(yàn)的平均數(shù)或中位數(shù)。使用相同的分布，你可以構(gòu)建一個(gè)95%的置信區(qū)間，與頻率主義統(tǒng)計(jì)中的置信區(qū)間相對(duì)應(yīng)。除了置信區(qū)間之外，貝葉斯的對(duì)應(yīng)區(qū)間直接量化了人口值在一定范圍內(nèi)的概率。所關(guān)注的參數(shù)值有95%的概率位于95%置信區(qū)間的邊界內(nèi)。與置信區(qū)間不同，這不僅僅是一個(gè)模擬量，而是一個(gè)簡明直觀的概率聲明。

問題：解釋估計(jì)效果、其區(qū)間和后驗(yàn)分布

答案：年齡似乎是預(yù)測博士延期的一個(gè)相關(guān)因素，后驗(yàn)平均回歸系數(shù)為2.317，95%HPD（可信區(qū)間）[1.194 3.417]。另外，age2似乎也是預(yù)測博士延期的一個(gè)相關(guān)因素，后驗(yàn)平均值為-0.022，95%可信區(qū)間為[-0.033-0.01]。95%的HPD顯示，人口中的這些回歸系數(shù)有95%的概率位于相應(yīng)的區(qū)間內(nèi)，也請(qǐng)看下面的數(shù)字中的后驗(yàn)分布。由于0不包含在可信區(qū)間內(nèi)，我們可以相當(dāng)肯定存在影響。

問題：?每個(gè)貝葉斯模型都使用一個(gè)先驗(yàn)分布。描述一下回歸系數(shù)的先驗(yàn)分布的形狀。

回答：

檢查使用了哪些默認(rèn)的先驗(yàn)。

（Jags）利用一個(gè)非常寬的正態(tài)分布來得出這個(gè)無信息的先驗(yàn)。默認(rèn)情況下，平均值為0，標(biāo)準(zhǔn)差為10（精度為0.01）。

回歸--用戶指定的先驗(yàn)

你也可以手動(dòng)指定你的先驗(yàn)分布。理論上，你可以使用你喜歡的任何一種分布來指定你的先驗(yàn)知識(shí)。然而，如果你的先驗(yàn)分布不遵循與你的似然相同的參數(shù)形式，計(jì)算模型可能會(huì)很麻煩。?共軛先驗(yàn)避免了這個(gè)問題，因?yàn)樗鼈儾捎昧四銟?gòu)建的模型的函數(shù)形式。對(duì)于你的正態(tài)線性回歸模型，如果你的回歸參數(shù)的預(yù)設(shè)是用正態(tài)分布來指定的，就可以達(dá)到共軛性（殘差得到一個(gè)反伽馬分布，這里忽略不計(jì)）。你可以很靈活地指定信息性先驗(yàn)。

讓我們用共軛先驗(yàn)來重新指定上面練習(xí)的回歸模型。我們暫時(shí)不涉及截距和殘差的預(yù)設(shè)。關(guān)于你的回歸參數(shù)，你需要指定其正態(tài)分布的超參數(shù)，即均值和方差。平均值表示你認(rèn)為哪一個(gè)參數(shù)值最有可能。方差表示你對(duì)此的確定程度。我們?yōu)棣履挲g回歸系數(shù)和β年齡2系數(shù)嘗試了4種不同的先驗(yàn)規(guī)范。

首先，我們使用以下先驗(yàn)。

Age?~?N(3,0.4)

Age2 ~?N(0,0.1)

先驗(yàn)指標(biāo)是在模型制定步驟中設(shè)置的。請(qǐng)注意，精度而不是正態(tài)分布的方差。精度是方差的倒數(shù)，所以方差為0.1對(duì)應(yīng)的精度為10，方差為0.4對(duì)應(yīng)的精度為2.5。

先驗(yàn)參數(shù)在代碼中呈現(xiàn)如下。

現(xiàn)在擬合模型并提供匯總統(tǒng)計(jì)。

回答：

問題 在下表中填寫結(jié)果。

年齡默認(rèn)情況下 先驗(yàn)N(3,.4)N(3,1000)N(20,.4)N(20,1000)后驗(yàn)平均值2.317后驗(yàn)標(biāo)準(zhǔn)差0.568

年齡2默認(rèn)情況下 先驗(yàn)N(0,.1)N(0,1000)N(20,.1)N(20,1000)后驗(yàn)平均值-0.022后驗(yàn)標(biāo)準(zhǔn)差0.006

回答：

年齡默認(rèn)情況下 先驗(yàn)N(3,.4)N(3,1000)N(20,.4)N(20,1000)后驗(yàn)平均值2.3172.625后驗(yàn)標(biāo)準(zhǔn)差0.5680.408

年齡2默認(rèn)情況下 先驗(yàn)N(0,.1)N(0,1000)N(20,.1)N(20,1000)后驗(yàn)平均值-0.022-0.026后驗(yàn)標(biāo)準(zhǔn)差0.0060.004

下一步，嘗試改編代碼，使用其他列的先驗(yàn)規(guī)范，然后完成該表。

回答：

年齡默認(rèn)情況下 先驗(yàn)N(3,.4)N(3,1000)N(20,.4)N(20,1000)后驗(yàn)平均值2.3172.6252.42211.1432.457后驗(yàn)標(biāo)準(zhǔn)差0.5680.4080.5020.5360.576

年齡2默認(rèn)情況下 先驗(yàn)N(0,.1)N(0,1000)N(20,.1)N(20,1000)后驗(yàn)平均值-0.022-0.026-0.023-0.113-0.024后驗(yàn)標(biāo)準(zhǔn)差0.0060.0040.0050.0060.006

問題：比較不同先驗(yàn)的結(jié)果。結(jié)果是否與默認(rèn)模型有可比性？

問題：使用不同的先驗(yàn)，我們最終的結(jié)論是否相似？

要回答這些問題，請(qǐng)按以下步驟進(jìn)行。我們可以計(jì)算出相對(duì)偏差來表示這種差異。我們將只計(jì)算兩個(gè)回歸系數(shù)的偏差，比較默認(rèn)（無信息）模型和使用N(20,.4)和N(20,.1)先驗(yàn)的模型。

Beta[1,2,1]和Beta[1,3,1]指數(shù)分別代表了βage和βage2。Beta[x,x,x]是回歸系數(shù)（按照我們?cè)谀Ｐ椭兄付ǖ捻樞?，所以首先是age，然后是age2），alpha[x,x,x]是截距，psi[x,x,x]是方差，def[x,x,x]是間接效應(yīng)（如果你的模型中有這些）。它們的排列順序與summary()輸出中的順序相同。因此，首先是回歸系數(shù)，然后是截距，然后是協(xié)方差，然后是間接效應(yīng)。

我們還可以通過繪制我們運(yùn)行的五個(gè)不同模型的后驗(yàn)和先驗(yàn)來繪制這些差異。在這個(gè)例子中，我們只繪制年齡βage的回歸系數(shù)。

首先我們提取5個(gè)不同模型的MCMC鏈，只針對(duì)這一個(gè)參數(shù)（βage=Beta[1,2,1]）。

然后，我們可以通過使用以下代碼繪制不同的后驗(yàn)和前驗(yàn)。

?現(xiàn)在，根據(jù)表格中的信息、偏差估計(jì)和圖表，你可以回答關(guān)于先驗(yàn)因素對(duì)結(jié)果的影響的兩個(gè)問題。

答案

我們看到，這種高信息量先驗(yàn)的影響對(duì)兩個(gè)回歸系數(shù)的影響分別為375%和397%左右。這是一個(gè)很大的差異，因此我們肯定不會(huì)得出類似的結(jié)論。

不同的先驗(yàn)，結(jié)果會(huì)發(fā)生變化，但仍具有可比性。只有對(duì)年齡使用N(20,.4)，才會(huì)產(chǎn)生真正不同的系數(shù)，因?yàn)檫@個(gè)先驗(yàn)均值離數(shù)據(jù)的均值很遠(yuǎn)，而其方差卻相當(dāng)確定。然而，一般來說，其他的結(jié)果是可以比較的。因?yàn)槲覀兪褂昧艘粋€(gè)大的數(shù)據(jù)集，先驗(yàn)的影響相對(duì)較小。如果使用一個(gè)較小的數(shù)據(jù)集，先驗(yàn)的影響就會(huì)更大。為了檢查這一點(diǎn)，你可以所有案例的大約20%進(jìn)行抽樣，然后重新進(jìn)行同樣的分析。結(jié)果當(dāng)然會(huì)不同，因?yàn)槲覀兪褂玫陌咐倭撕芏?。使用這段代碼。

我們做了一個(gè)新的數(shù)據(jù)集，從原始數(shù)據(jù)集的333個(gè)觀測值中隨機(jī)選擇了60個(gè)。用同樣的代碼重復(fù)分析，只改變數(shù)據(jù)集的名稱，以觀察先驗(yàn)因素對(duì)較小數(shù)據(jù)集的影響。用這個(gè)新的數(shù)據(jù)集運(yùn)行priors2模型

參考文獻(xiàn)

Benjamin, D. J., Berger, J., Johannesson, M., Nosek, B. A., Wagenmakers, E.,… Johnson, V. (2017, July 22).?Redefine statistical significance. Retrieved from psyarxiv.com/mky9j

Greenland, S., Senn, S. J., Rothman, K. J., Carlin, J. B., Poole, C., Goodman, S. N. Altman, D. G. (2016).?Statistical tests, P values, confidence intervals, and power: a guide to misinterpretations.?European Journal of Epidemiology 31 (4).?https://doi.org/10.1007/s10654-016-0149-3?_ _

van de Schoot R, Yerkes MA, Mouw JM, Sonneveld H (2013)?What Took Them So Long? Explaining PhD Delays among Doctoral Candidates.?PLoS ONE 8(7): e68839.?https://doi.org/10.1371/journal.pone.0068839

最受歡迎的見解

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

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

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

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

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

6.Python用PyMC3實(shí)現(xiàn)貝葉斯線性回歸模型

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

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

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