原文連結：http://tecdat.cn/?p=23652

原文出處：拓端資料部落公眾號

本文為讀者提供瞭如何進行貝葉斯迴歸的基本教程。包括完成匯入資料檔案、探索彙總統計和迴歸分析。

在本文中，我們首先使用軟體的預設先驗設定。在第二步中，我們將應用使用者指定的先驗，對自己的資料使用貝葉斯。

準備工作

本教程要求:

• 已安裝的JAGS
• 安裝R軟體。
• 假設檢驗的基本知識
• 相關性和迴歸的基本知識
• 貝葉斯推理的基本知識
• R語言編碼的基本知識

資料例項

我們在這個練習中使用的資料是基於一項關於預測博士生完成論文時間的研究（Van de Schoot, Yerkes, Mouw and Sonneveld 2013）。研究人員詢問了博士生完成他們的博士論文需要多長時間（n=333）。結果顯示，博士學位獲得者平均花了59.8個月（5年4個月）來完成他們的博士學位。變數B3衡量計劃和實際專案時間之間的差異，以月為單位（平均=9.97，最小=-31，最大=91，sd=14.43）。

對於目前的工作，我們感興趣的問題是，博士學位獲得者的年齡（M=31.7，SD=6.86）是否與他們專案的延期有關。

預計完成時間和年齡之間的關係是非線性的。這可能是由於在人生的某個階段（即三十多歲），家庭生活比你在二十多歲時或年長時佔用了你更多的時間。

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

問題：請寫出零假設和備擇假設。寫下代表這個問題的無效假設和備選假設。你認為哪個假設更有可能？

回答：

H0:年齡與博士專案的延期無關。

H1:年齡與博士專案的延期有關。

H0:age2與博士專案的延期無關。

H1:age2與博士專案的延期有關。

準備--匯入和探索資料

資料是一個.csv檔案，但你可以使用以下語法直接將其載入到R中。

一旦你載入了你的資料，建議你檢查一下你的資料匯入是否順利。因此，首先看看你的資料的彙總統計。你可以使用describe()函式。

問題：你所有的資料都被正確地載入了嗎？也就是說，所有的資料點都有實質性的意義嗎？

回答：

describe(data)

描述性統計有意義。

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

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

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

繪圖

在繼續分析資料之前，我們還可以繪製期望的關係。

1. plot(aes(x = age,
2. y = diff))

迴歸

在這個練習中，你將研究博士生的年齡和age2對他們的專案時間延期的影響，這作為結果變數使用迴歸分析。

如你所知，貝葉斯推理包括將先驗分佈與從資料中獲得的似然性相結合。指定先驗分佈是貝葉斯推斷中最關鍵的一點，應該受到高度重視（例如Van de Schoot等人，2017）。在本教程中，我們將首先依賴預設的先驗設定。

要用執行多元迴歸，首先要指定模型，然後擬合模型，最後獲得總結。模型的指定方法如下。

1. 我們想要預測的因變數。
2. "~"，我們用它來表示我們現在給其他感興趣的變數。(相當於迴歸方程的"="）。
3. 用求和符號'+'分隔的不同自變數。
4. 最後，我們插入因變數有一個方差，有一個截距。

下面的程式碼是如何指定迴歸模型的。

1. # 1) 指定模型
2. '#迴歸模型
3. diff ~ age + age2
5. #顯示因變數有方差
6. diff ~~ diff
8. #有一個截距
9. diff ~~ 1'

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

1. fitbayes(model, data, target = "jags", test = "none", seed = c(12,34,56) )
3. # test="none "的輸入停止了一些後驗的計算，我們現在不需要，加快了計算過程。
4. # 種子命令只是為了保證在多次執行取樣器時有相同的準確結果。你不需要設定這個。當使用Jags時，你需要設定儘可能多的種子鏈（預設）。

現在我們用summary(fit.bayes)來看看總結。

顯示輸出

頻率主義模型與貝葉斯分析模型所提供的結果確實不同。

貝葉斯統計推斷和頻率主義統計方法之間的關鍵區別在於估計的未知引數的性質。在頻率主義框架中，一個感興趣的引數被假定為未知的，但卻是固定的。也就是說，假設在人口中只有一個真實的人口引數，例如，一個真實的平均值或一個真實的迴歸係數。在貝葉斯的主觀概率觀中，所有的未知引數都被視為不確定的，因此要用一個概率分佈來描述。每個引數都是未知的，而所有未知的東西都會得到一個分佈。

這就是為什麼在頻率推斷中，你主要得到的是一個未知但固定的群體引數的點估計。這是一個引數值，考慮到資料，它最有可能出現在人群中。附帶的置信區間試圖讓你進一步瞭解這個估計值的不確定性。重要的是要認識到，置信區間只是構成一個模擬量。在從人口中抽取的無限多的樣本中，構建（95%）置信區間的程式將使其在95%的時間內包含真實的人口值。這並沒有為你提供任何資訊，即人口引數位於你所分析的非常具體和唯一的樣本中的置信區間邊界內的可能性有多大。

在貝葉斯分析中，你推斷的關鍵是感興趣的引數的後驗分佈。它滿足了概率分佈的每一個屬性，並量化了人口引數位於某些區域的概率。一方面，你可以通過它的模式來描述後驗的特點。這是一個引數值，考慮到資料和它的先驗概率，它在人群中是最有可能的。另外，你也可以使用後驗的平均數或中位數。使用相同的分佈，你可以構建一個95%的置信區間，與頻率主義統計中的置信區間相對應。除了置信區間之外，貝葉斯的對應區間直接量化了人口值在一定範圍內的概率。所關注的引數值有95%的概率位於95%置信區間的邊界內。與置信區間不同，這不僅僅是一個模擬量，而是一個簡明直觀的概率宣告。

問題：解釋估計效果、其區間和後驗分佈

答案：年齡似乎是預測博士延期的一個相關因素，後驗平均迴歸係數為2.317，95%HPD（可信區間）[1.194 3.417]。另外，age2似乎也是預測博士延期的一個相關因素，後驗平均值為-0.022，95%可信區間為[-0.033-0.01]。95%的HPD顯示，人口中的這些迴歸係數有95%的概率位於相應的區間內，也請看下面的數字中的後驗分佈。由於0不包含在可信區間內，我們可以相當肯定存在影響。

問題：每個貝葉斯模型都使用一個先驗分佈。描述一下回歸係數的先驗分佈的形狀。

回答：

檢查使用了哪些預設的先驗。

（Jags）利用一個非常寬的正態分佈來得出這個無資訊的先驗。預設情況下，平均值為0，標準差為10（精度為0.01）。

迴歸--使用者指定的先驗

你也可以手動指定你的先驗分佈。理論上，你可以使用你喜歡的任何一種分佈來指定你的先驗知識。然而，如果你的先驗分佈不遵循與你的似然相同的引數形式，計算模型可能會很麻煩。共軛先驗避免了這個問題，因為它們採用了你構建的模型的函式形式。對於你的正態線性迴歸模型，如果你的迴歸引數的預設是用正態分佈來指定的，就可以達到共軛性（殘差得到一個反伽馬分佈，這裡忽略不計）。你可以很靈活地指定資訊性先驗。

讓我們用共軛先驗來重新指定上面練習的迴歸模型。我們暫時不涉及截距和殘差的預設。關於你的迴歸引數，你需要指定其正態分佈的超引數，即均值和方差。平均值表示你認為哪一個引數值最有可能。方差表示你對此的確定程度。我們為β年齡迴歸係數和β年齡2係數嘗試了4種不同的先驗規範。

首先，我們使用以下先驗。

Age~N(3,0.4)

Age2 ~N(0,0.1)

先驗指標是在模型制定步驟中設定的。請注意，精度而不是正態分佈的方差。精度是方差的倒數，所以方差為0.1對應的精度為10，方差為0.4對應的精度為2.5。

先驗引數在程式碼中呈現如下。

1. '#帶有priors的迴歸模型
2. prior("dnorm(3,2.5)")*age + prior("dnorm(0,10)")*age2

現在擬合模型並提供彙總統計。

回答：

summary(fit)

問題 在下表中填寫結果。

回答：

下一步，嘗試改編程式碼，使用其他列的先驗規範，然後完成該表。

回答：

1. '#有priors的迴歸模型
2. ~ prior("dnorm(3,0.001)")*age + prior("dnorm(0,0.001)")*age2

summary(prior2)

summary(prior3)

summary(prior4)

問題：比較不同先驗的結果。結果是否與預設模型有可比性？

問題：使用不同的先驗，我們最終的結論是否相似？

要回答這些問題，請按以下步驟進行。我們可以計算出相對偏差來表示這種差異。我們將只計算兩個迴歸係數的偏差，比較預設（無資訊）模型和使用N(20,.4)和N(20,.1)先驗的模型。

1. #1）減去MCMC鏈的內容
2. fitbayes( what = "mcmc")
4. #2) 繫結不同的鏈，計算迴歸係數的平均值（估計值）。
5. colMeans(as.matrix(mcmc.list)
7. #3) 計算偏差
8. 100*((estimat-estimat)/estimat)

beta[1,2,1]和beta[1,3,1]指數分別代表了βage和βage2。Beta[x,x,x]是迴歸係數（按照我們在模型中指定的順序，所以首先是age，然後是age2），alpha[x,x,x]是截距，psi[x,x,x]是方差，def[x,x,x]是間接效應（如果你的模型中有這些）。它們的排列順序與summary()輸出中的順序相同。因此，首先是迴歸係數，然後是截距，然後是協方差，然後是間接效應。

我們還可以通過繪製我們執行的五個不同模型的後驗和先驗來繪製這些差異。在這個例子中，我們只繪製年齡βage的迴歸係數。

首先我們提取5個不同模型的MCMC鏈，只針對這一個引數（βage=beta[1,2,1]）。

 binrows(posterior1.5, prior1.5)

然後，我們可以通過使用以下程式碼繪製不同的後驗和前驗。

1. qplot(data = post,)+
2. density(size = 1)+
3. scale_x_continuous(limit )

現在，根據表格中的資訊、偏差估計和圖表，你可以回答關於先驗因素對結果的影響的兩個問題。

答案

1. #1）減去MCMC鏈
2. fit.bayes(what = "mcmc")
4. #2) 繫結不同的鏈，計算迴歸係數的平均值（估計值）。
5. colMeans(as.matrix(mcmc.list)
8. #3) 計算偏差
9. 100*((imstim-estima)/estimat)
11. ## beta[1,2,1] beta[1,3,1].
12. ## 374.55 397.31

我們看到，這種高資訊量先驗的影響對兩個迴歸係數的影響分別為375%和397%左右。這是一個很大的差異，因此我們肯定不會得出類似的結論。

不同的先驗，結果會發生變化，但仍具有可比性。只有對年齡使用N(20,.4)，才會產生真正不同的係數，因為這個先驗均值離資料的均值很遠，而其方差卻相當確定。然而，一般來說，其他的結果是可以比較的。因為我們使用了一個大的資料集，先驗的影響相對較小。如果使用一個較小的資料集，先驗的影響就會更大。為了檢查這一點，你可以所有案例的大約20%進行抽樣，然後重新進行同樣的分析。結果當然會不同，因為我們使用的案例少了很多。使用這段程式碼。

1. indices <- sample.int(333, 60)
2. smalldata <- data[indices,]

我們做了一個新的資料集，從原始資料集的333個觀測值中隨機選擇了60個。用同樣的程式碼重複分析，只改變資料集的名稱，以觀察先驗因素對較小資料集的影響。用這個新的資料集執行priors2模型

1. fit.bayes(model = priors2,smalldata)
2. summary(fit)

參考文獻

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使用貝葉斯優化的深度學習

2.matlab貝葉斯隱馬爾可夫hmm模型實現

3.R語言Gibbs抽樣的貝葉斯簡單線性迴歸模擬

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

5.R語言中的Stan概率程式設計MCMC取樣的貝葉斯模型

6.Python用PyMC3實現貝葉斯線性迴歸模型

7.R語言使用貝葉斯 層次模型進行空間資料分析

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

9.matlab貝葉斯隱馬爾可夫hmm模型實現