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

最近我們被客戶要求撰寫關(guān)于貝葉斯回歸模型的研究報(bào)告，包括一些圖形和統(tǒng)計(jì)輸出。

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

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

準(zhǔn)備工作

本教程要求:

• 已安裝的JAGS

• 安裝R軟件。

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

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

• 貝葉斯推理的基本知識

• R語言編碼的基本知識

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

我們在這個(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）。

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

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

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

問題：請寫出零假設(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)。

向下滑動查看結(jié)果▼

• 
  

• 
  

準(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ì)性的意義嗎？

describe(data)

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

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

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

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

向下滑動查看結(jié)果▼

**

繪圖

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

plot(aes(x?=?age, ?????????????y?=?diff))

回歸

在這個(gè)練習(xí)中，你將研究博士生的年齡和age2對他們的項(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è)置。

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

R語言實(shí)現(xiàn)貝葉斯分位數(shù)回歸、lasso和自適應(yīng)lasso貝葉斯分位數(shù)回歸分析

左右滑動查看更多

01

02

03

04

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

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

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

• 用求和符號'+'分隔的不同自變量。

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

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

#?1)?指定模型?'#回歸模型????????????????????diff?~?age?+?age2 ????????????????????#顯示因變量有方差????????????????????diff?~~?diff ????????????????????#有一個(gè)截距????????????????????diff?~~?1'

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

fitbayes(model,?data,?target?=?"jags",?test?=?"none",?seed?=?c(12,34,56)?)# test="none "的輸入停止了一些后驗(yàn)的計(jì)算，我們現(xiàn)在不需要，加快了計(jì)算過程。#?種子命令只是為了保證在多次運(yùn)行采樣器時(shí)有相同的準(zhǔn)確結(jié)果。你不需要設(shè)置這個(gè)。當(dāng)使用Jags時(shí)，你需要設(shè)置盡可能多的種子鏈（默認(rèn)）。

現(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ù)都是未知的，而所有未知的東西都會得到一個(gè)分布。

這就是為什么在_頻率_推斷中，你主要得到的是一個(gè)未知但固定的群體參數(shù)的點(diǎn)估計(jì)。這是一個(gè)參數(shù)值，考慮到數(shù)據(jù)，它最有可能出現(xiàn)在人群中。附帶的置信區(qū)間試圖讓你進(jìn)一步了解這個(gè)估計(jì)值的不確定性。重要的是要認(rèn)識到，置信區(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ū)間相對應(yīng)。除了置信區(qū)間之外，貝葉斯的對應(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)，也請看下面的數(shù)字中的后驗(yàn)分布。由于0不包含在可信區(qū)間內(nèi)，我們可以相當(dāng)肯定存在影響。

向下滑動查看結(jié)果▼

問題：?每個(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）。

向下滑動查看結(jié)果▼

• 
  

• 
  

**

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

你也可以手動指定你的先驗(yàn)分布。理論上，你可以使用你喜歡的任何一種分布來指定你的先驗(yàn)知識。然而，如果你的先驗(yàn)分布不遵循與你的似然相同的參數(shù)形式，計(jì)算模型可能會很麻煩。?共軛先驗(yàn)避免了這個(gè)問題，因?yàn)樗鼈儾捎昧四銟?gòu)建的模型的函數(shù)形式。對于你的正態(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ù)值最有可能。方差表示你對此的確定程度。我們?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è)置的。請注意，精度而不是正態(tài)分布的方差。精度是方差的倒數(shù)，所以方差為0.1對應(yīng)的精度為10，方差為0.4對應(yīng)的精度為2.5。

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

??'#帶有priors的回歸模型prior("dnorm(3,2.5)")*age?+?prior("dnorm(0,10)")*age2

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

回答：

summary(fit)

向下滑動查看結(jié)果▼

問題：在下表中填寫結(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ī)范，然后完成該表。

?'#有priors的回歸模型 ??~?prior("dnorm(3,0.001)")*age?+?prior("dnorm(0,0.001)")*age2summary(prior2)

summary(prior3)

summary(prior4)

年齡默認(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

向下滑動查看結(jié)果▼

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

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

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

#1）減去MCMC鏈的內(nèi)容fitbayes(?what?=?"mcmc")#2)?綁定不同的鏈，計(jì)算回歸系數(shù)的平均值（估計(jì)值）。?colMeans(as.matrix(mcmc.list)#3)?計(jì)算偏差100*((estimat-estimat)/estimat)

beta[1,2,1]和beta[1,3,1]指數(shù)分別代表了βage和βage2。Beta[x,x,x]是回歸系數(shù)（按照我們在模型中指定的順序，所以首先是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鏈，只針對這一個(gè)參數(shù)（βage=beta[1,2,1]）。

?binrows(posterior1.5,?prior1.5)

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

qplot(data?=?post,)+ ??density(size?=?1)+ ??scale_x_continuous(limit?)

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

回答：

#1）減去MCMC鏈fit.bayes(what?=?"mcmc")#2)?綁定不同的鏈，計(jì)算回歸系數(shù)的平均值（估計(jì)值）。colMeans(as.matrix(mcmc.list)#3)?計(jì)算偏差100*((imstim-estima)/estimat)##?beta[1,2,1]?beta[1,3,1].?##?374.55?397.31

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

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

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

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

fit.bayes(model?=?priors2,smalldata) summary(fit)

向下滑動查看結(jié)果▼

參考文獻(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

點(diǎn)擊文末?“閱讀原文”

獲取全文完整代碼數(shù)據(jù)資料。

本文選自《R語言JAGS貝葉斯回歸模型分析博士生延期畢業(yè)完成論文時(shí)間》。

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

R語言Gibbs抽樣的貝葉斯簡單線性回歸仿真分析
python貝葉斯隨機(jī)過程：馬爾可夫鏈Markov-Chain，MC和Metropolis-Hastings，MH采樣算法可視化
Python貝葉斯推斷Metropolis-Hastings（M-H）MCMC采樣算法的實(shí)現(xiàn)
Metropolis Hastings采樣和貝葉斯泊松回歸Poisson模型
Matlab用BUGS馬爾可夫區(qū)制轉(zhuǎn)換Markov switching隨機(jī)波動率模型、序列蒙特卡羅SMC、M H采樣分析時(shí)間序列R語言RSTAN MCMC：NUTS采樣算法用LASSO 構(gòu)建貝葉斯線性回歸模型分析職業(yè)聲望數(shù)據(jù)
R語言BUGS序列蒙特卡羅SMC、馬爾可夫轉(zhuǎn)換隨機(jī)波動率SV模型、粒子濾波、Metropolis Hasting采樣時(shí)間序列分析
R語言Metropolis Hastings采樣和貝葉斯泊松回歸Poisson模型
R語言貝葉斯MCMC：用rstan建立線性回歸模型分析汽車數(shù)據(jù)和可視化診斷
R語言貝葉斯MCMC：GLM邏輯回歸、Rstan線性回歸、Metropolis Hastings與Gibbs采樣算法實(shí)例
R語言貝葉斯Poisson泊松-正態(tài)分布模型分析職業(yè)足球比賽進(jìn)球數(shù)
R語言用Rcpp加速M(fèi)etropolis-Hastings抽樣估計(jì)貝葉斯邏輯回歸模型的參數(shù)
R語言邏輯回歸、Naive Bayes貝葉斯、決策樹、隨機(jī)森林算法預(yù)測心臟病
R語言中貝葉斯網(wǎng)絡(luò)（BN）、動態(tài)貝葉斯網(wǎng)絡(luò)、線性模型分析錯(cuò)頜畸形數(shù)據(jù)
R語言中的block Gibbs吉布斯采樣貝葉斯多元線性回歸
Python貝葉斯回歸分析住房負(fù)擔(dān)能力數(shù)據(jù)集
R語言實(shí)現(xiàn)貝葉斯分位數(shù)回歸、lasso和自適應(yīng)lasso貝葉斯分位數(shù)回歸分析
Python用PyMC3實(shí)現(xiàn)貝葉斯線性回歸模型
R語言用WinBUGS 軟件對學(xué)術(shù)能力測驗(yàn)建立層次（分層）貝葉斯模型
R語言Gibbs抽樣的貝葉斯簡單線性回歸仿真分析
R語言和STAN,JAGS：用RSTAN,RJAG建立貝葉斯多元線性回歸預(yù)測選舉數(shù)據(jù)
R語言基于copula的貝葉斯分層混合模型的診斷準(zhǔn)確性研究
R語言貝葉斯線性回歸和多元線性回歸構(gòu)建工資預(yù)測模型
R語言貝葉斯推斷與MCMC：實(shí)現(xiàn)Metropolis-Hastings 采樣算法示例
R語言stan進(jìn)行基于貝葉斯推斷的回歸模型
R語言中RStan貝葉斯層次模型分析示例
R語言使用Metropolis-Hastings采樣算法自適應(yīng)貝葉斯估計(jì)與可視化
R語言隨機(jī)搜索變量選擇SSVS估計(jì)貝葉斯向量自回歸（BVAR）模型
WinBUGS對多元隨機(jī)波動率模型：貝葉斯估計(jì)與模型比較
R語言實(shí)現(xiàn)MCMC中的Metropolis–Hastings算法與吉布斯采樣
R語言貝葉斯推斷與MCMC：實(shí)現(xiàn)Metropolis-Hastings 采樣算法示例
R語言使用Metropolis-Hastings采樣算法自適應(yīng)貝葉斯估計(jì)與可視化
視頻：R語言中的Stan概率編程MCMC采樣的貝葉斯模型
R語言MCMC:Metropolis-Hastings采樣用于回歸的貝葉斯估計(jì)