【視頻】線性回歸中的貝葉斯推斷與R語言預測工人工資數(shù)據(jù)|數(shù)據(jù)分享
原文鏈接：http://tecdat.cn/?p=24141
原文出處：拓端數(shù)據(jù)部落公眾號
在這個視頻中，我們轉向簡單線性回歸中的貝葉斯推斷。?我們將使用一個參照先驗分布，它提供了頻率主義解決方案和貝葉斯答案之間的聯(lián)系。?然后在R語言中用貝葉斯線性回歸、貝葉斯模型平均 (BMA)來預測工人工資數(shù)據(jù)?（查看文末了解數(shù)據(jù)獲取方式）。
視頻：線性回歸中的貝葉斯推斷與R語言預測工人工資數(shù)據(jù)案例
貝葉斯推斷線性回歸與R語言預測工人工資數(shù)據(jù)
，時長09:58
為了說明這些想法，我們將使用一個例子來預測身體脂肪。?獲得準確的體脂測量值是很昂貴的，而且在家里也不容易做到。?
相反，可以使用現(xiàn)成的測量值（如腹圍）來預測體脂百分比的預測模型是很容易使用的，而且價格低廉。
圖中顯示了從水下稱重得到的體脂百分比和多名男性的腹圍。?
為了預測體脂，我們可以從一個線性回歸開始，在散點圖中加入了回歸線。
這條線的估計斜率為0.63，截距約為-39%。?每增加一厘米，我們預計體脂會增加0.63%。?負的截距過程作為一個物理模型沒有意義，但預測一個腰圍為0厘米的男性也沒有意義。?盡管如此，這個線性回歸可能是一個準確的近似值，用于預測目的。?隨著測量結果在這個人群的觀察范圍內(nèi)。
我們的最佳估計線使用α和β的通常平方估計值來獲得擬合值和預測值。
提供擬合誤差估計的殘差是觀察值和預測值之間的差異，用于診斷以及估計sigma平方。?這是通過平均平方誤差，即平方誤差之和除以自由度。
記住，剩余自由度是樣本量減去模型中的回歸系數(shù)。
模型和先驗
讓我們從基準的角度來看看估計的情況。?我們將從與OLS相同的模型開始，但有一個額外的假設，即誤差是正態(tài)分布的常數(shù)變體。
我們經(jīng)常加上這個假設，以獲得具有常數(shù)平方的置信區(qū)間。?因此，在頻率主義回歸中經(jīng)常使用的任何診斷圖都可以在這里用來檢查這個假設。
對于共軛分析，先驗是以西格瑪為條件的回歸系數(shù)的雙變量正態(tài)分布。Sigma在這里提供了響應的單位比例。
從邊際上看，α是正態(tài)分布，給定西格瑪?shù)钠椒剑宰⒁獠辉谶@個無單位規(guī)模參數(shù)作為子α的方差控制。?同樣地，β是正態(tài)分布，給定西格瑪?shù)钠椒?，其平均值為b0，方差控制由參數(shù)S次β。
阿爾法和貝塔之間的協(xié)方差，是西格瑪平方乘以一個參數(shù)，這個參數(shù)描述了我們對阿爾法和貝塔如何共同變化的先驗信念。?如果s alpha beta是0，那么先驗的信念是，alpha將獨立于beta，條件是sigma平方。
為了完成規(guī)范，我們使用共軛先驗sigma平方，其中1動方差有一個自由度為N的gamma分布，sigma平方不是對sigma平方的先驗估計。?這與我們在以前的視頻中使用的先驗相似。?由于共軛的關系，后驗分布將是正態(tài)伽馬分布，有簡單的規(guī)則來更新參數(shù)。?
參照先驗和后驗分布
提供一個貝葉斯分析作為起點是很有用的。?參照先驗分布是作為那些方差參數(shù)到無窮大時的正態(tài)分布的極限得到的，并且在α和β中是平坦或均勻的。?而共軛伽馬先驗的極限，先驗自由度為0，提供了我們在以前視頻中使用的西格瑪平方的參照先驗。
與之相關的β的參照后驗是一個以OLS估計值為中心的學生t分布，其尺度與OLS標準誤差相同。?自由度是n減2，就像在頻率分析中一樣。
同樣，α的邊際后驗分布也是一個學生t分布，其中心和尺度由OLS估計值給出。
利用這些參數(shù)的聯(lián)合分布，我們可以得到x的任何值的預期體脂分布，這是一個學生t分布，有n減2的自由度。?估計值由OLS估計值給出。
使用統(tǒng)計軟件，我們可以得到參數(shù)估計值、標準差和95%置信區(qū)間。?你應該能夠確認，當我們使用這個參照先驗時，貝葉斯估計值與頻率主義估計值相同。
參數(shù)估計
參數(shù)、擬合值或預測的置信區(qū)間都是以后驗均值加或減適當?shù)膖-quantile乘以標準差得到的。?主要區(qū)別在于對區(qū)間的解釋。?例如，根據(jù)數(shù)據(jù)，我們現(xiàn)在認為，腰圍每增加10厘米，體脂將增加5.8%至6.9%的可能性是95%。?當然，這個模型是一個近似值，所以要小心任何因果關系的解釋，然而，它仍然可以對預測有幫助。
預測體脂
對于預測身體脂肪，我們將使用后驗預測分布。
如果我們知道參數(shù)，一個新的觀察將只是通過取平均值來獲得，基于我們的人口回歸方程加上相關的不確定性，描述個體偏離人口平均值和x的程度。
鑒于手頭的數(shù)據(jù)，我們的后驗預測分布是一個具有n-2個自由度的t分布。?預測新值的最佳估計是人口線的后驗平均值。?這與我們之前計算的擬合值是一樣的，但是基于預測標準差的比例參數(shù)。
預測標準差包含了x處的回歸線的后驗不確定性。從最后一個項中，我們可以看到，當我們在數(shù)據(jù)的平均值附近進行預測時，變異性將是最小的。?隨著處于x范圍兩端的x的值的變異性增加。
還有一個額外的估計值sigma squared，它來自于誤差epsilon的不確定性，或者說我們期望觀測值偏離回歸線的程度。?
圖中顯示了數(shù)據(jù)，預測方程是后驗平均數(shù)，用紅色表示，估計總體平均數(shù)的95%的點狀區(qū)間用灰色虛線表示。?預測體脂的95%可信區(qū)間是外側的虛線。?大多數(shù)數(shù)據(jù)都在預測區(qū)間內(nèi)，正如人們所期望的那樣。?然而，橙色圈出的點的體脂率遠遠低于模型的預期。
?
R語言用貝葉斯線性回歸、貝葉斯模型平均 (BMA)來預測工人工資案例
下面，貝葉斯信息準則（BIC）和貝葉斯模型平均法被應用于構建一個簡明的收入預測模型。
這些數(shù)據(jù)是從 935 名受訪者的隨機樣本中收集的。該數(shù)據(jù)集是_計量經(jīng)濟學數(shù)據(jù)集_系列的一部分 。
加載包
數(shù)據(jù)將首先使用該dplyr?包進行探索?，并使用該ggplot2?包進行可視化?。稍后，實現(xiàn)逐步貝葉斯線性回歸和貝葉斯模型平均 (BMA)。
數(shù)據(jù)
數(shù)據(jù)集網(wǎng)頁提供了以下變量描述表：
變量描述wage每周收入（元）hours每周平均工作時間IQ智商分數(shù)kww對世界工作的了解得分educ受教育年數(shù)exper多年工作經(jīng)驗tenure在現(xiàn)任雇主工作的年數(shù)age年齡married=1 如果已婚black=1 如果是黑人south=1 如果住在南方urban=1 如果居住在都市sibs兄弟姐妹的數(shù)量brthord出生順序meduc母親的教育（年）feduc父親的教育（年）lwage工資自然對數(shù)?wage
探索數(shù)據(jù)
與任何新數(shù)據(jù)集一樣，一個好的起點是標準的探索性數(shù)據(jù)分析。匯總表是簡單的第一步。
#?數(shù)據(jù)集中所有變量的匯總表--包括連續(xù)變量和分類變量summary(wage)
因變量（工資）的直方圖給出了合理預測應該是什么樣子的。
#工資數(shù)據(jù)的簡單柱狀圖hst(wge$wae,?breks?=?30)
直方圖還可用于大致了解哪些地方不太可能出現(xiàn)結果。
#?檢查圖表?"尾部?"的點的數(shù)量sm(wage$ge?<?300)##?\[1\]?6sm(wae$wge?>?2000)##?\[1\]?20
簡單線性回歸
由于周工資（'wage'）是該分析中的因變量，我們想探索其他變量作為預測變量的關系。我們在數(shù)據(jù)中看到的工資變化的一種可能的、簡單的解釋是更聰明的人賺更多的錢。下圖顯示了每周工資和 IQ 分數(shù)之間的散點圖。
gplot(wae,?es(iq,?wge))?+?gom\_oint()?+gom\_smoth()
IQ 分數(shù)和工資之間似乎存在輕微的正線性關系，但僅靠 IQ 并不能可靠地預測工資。盡管如此，這種關系可以通過擬合一個簡單的線性回歸來量化，它給出：
工資 i = α + β?iqi + ?iwagei = α + β?iqi + ?i
m\_wg\_iq?=?lm(wge?~?iq,?dta?=?age)coefients
工資 i = 116.99 + 8.3 ?iqi + ?iwagei = 116.99 + 8.3 ?iqi + ?i
在轉向貝葉斯改進這個模型之前，請注意貝葉斯建模假設誤差 (?i) 以恒定方差正態(tài)分布。通過檢查模型的殘差分布來檢查該假設。如果殘差高度非正態(tài)或偏斜，則違反假設并且任何后續(xù)推論都無效。要檢查假設，請按如下方式繪制殘差：
#?用散點圖和模型誤差殘差的直方圖來檢查正態(tài)性假設glot(dta?=?mwag_q,?es(x?=?.ite,?y?=?.rd))?+??gemittr()?+
plot(dta?=?m\_g\_iq,?aes(x?=?.reid))?+??histgm(bnwth?=?10)
變量變換
兩個圖都顯示殘差是右偏的。因此，IQ（因為它目前存在于數(shù)據(jù)集中）不應用作貝葉斯預測模型。但是，對?僅具有正值的偏斜_因_變量使用（自然）對數(shù)變換?通?？梢越鉀Q問題。下面，該模型使用轉換后的工資變量進行了重新擬合。
#?用IQ的自然對數(shù)擬合th模型lm(lage?~?iq,?data?=?wae)
#?殘差sctterplot和轉換后數(shù)據(jù)的柱狀圖plt(data?=?m\_lag\_iq,?es(x?=?.fited,?y?=?.reid))???geiter()?+
ggpot(dta?=?m_lwgeiq,?as(x?=?.resd))?+??gostgam(binwth?=?.1)?+
殘差確實大致呈正態(tài)分布。然而，由此產(chǎn)生的 IQ 系數(shù)非常?。ㄖ挥?0.0088），這是可以預料的，因為 IQ 分數(shù)提高 1 分幾乎不會對工資產(chǎn)生太大影響。需要進一步細化。數(shù)據(jù)集包含更多信息。
多元線性回歸和 BIC
我們可以首先在回歸模型中包含所有潛在的解釋變量，來粗略地嘗試解釋盡可能多的工資變化。
#?對數(shù)據(jù)集中的所有變量運行一個線性模型，使用'.'約定。full?=?lm(lwge?~?.?-?wage,?dta?=?wge)
完整線性模型的上述總結表明，自變量的許多系數(shù)在統(tǒng)計上并不顯著（請參閱第 4 個數(shù)字列中的 p 值）。選擇模型變量的一種方法是使用貝葉斯信息準則 (BIC)。BIC 是模型擬合的數(shù)值評估，它也會按樣本大小的比例懲罰更多的參數(shù)。這是完整線性模型的 BIC：
BIC(full)
BIC 值越小表示擬合越好。因此，BIC 可以針對各種縮減模型進行計算，然后與完整模型 BIC 進行比較，以找到適合工資預測工作的最佳模型。當然，R 有一個功能可以系統(tǒng)地執(zhí)行這些 BIC 調(diào)整。
#?用step計算模型pIC(lwge?~?.?-?wge,?dta?=?na.oi(wge))lg(lgth(na.mit(wge))))
#?顯示逐步模型的BICBIC(se_mol)
調(diào)用 step找到產(chǎn)生最低 BIC 的變量組合，并提供它們的系數(shù)。很不錯。
貝葉斯模型平均（BMA）
即使BIC處于最低值，我們能有多大把握確定所得到的模型是真正的 "最佳擬合"？答案很可能取決于基礎數(shù)據(jù)的規(guī)模和穩(wěn)定性。在這些不確定的時候，貝葉斯模型平均化（BMA）是有幫助的。BMA對多個模型進行平均化，獲得系數(shù)的后驗值和新數(shù)據(jù)的預測值。下面，BMA被應用于工資數(shù)據(jù)（排除NA值后）。
#?不包括NAa_ona?=?na.omt(wae)#?運行BMA，指定BIC作為判斷結果模型的標準BMA(wge?~?.?-wge,?daa=?ae\_o\_a,???????????????????pror?=?"BIC",????????????????????moepor?=?ufom())#?顯示結果summary
結果表顯示了五個最有可能的模型，以及每個系數(shù)被包含在真實模型中的概率。我們看到，出生順序和是否有兄弟姐妹是最不可能被包含的變量，而教育和智商變量則被鎖定。BMA模型的排名也可以用圖像圖來顯示，它清楚地顯示哪些變量在所有模型中，哪些變量被排除在所有模型之外，以及那些介于兩者之間的變量。
ge(b_lge,?tp.oels)
我們還可以提供模型系數(shù)的95%置信區(qū)間。下面的結果支持了關于包括或排除系數(shù)的決定。例如，在區(qū)間包含零，有大量證據(jù)支持排除該變量。
confint(ceflae)
進行預測
構建模型后，pediction 只是插入數(shù)據(jù)的問題：
#?用一個虛構的工人的統(tǒng)計資料來預測數(shù)據(jù)的例子#?進行預測redict?=?pedct(e_odl,?newdt?=?wrkr,eitr?=?"BMA")#?將結果轉換為元exp(wk_pedct)
預計這名化妝工作人員的周薪為 745 元。這到底有多準確？你得問她，但我們對我們的變量選擇很有信心，并對現(xiàn)有的數(shù)據(jù)盡了最大努力。應用的貝葉斯技術使我們對結果有信心。
數(shù)據(jù)獲取
在下面公眾號后臺回復“工資數(shù)據(jù)”，可獲取完整數(shù)據(jù)。
本文摘選《R語言用貝葉斯線性回歸、貝葉斯模型平均 (BMA)來預測工人工資》，點擊“閱讀原文”獲取全文完整資料。
點擊標題查閱往期內(nèi)容
R語言邏輯回歸、Naive Bayes貝葉斯、決策樹、隨機森林算法預測心臟病
PYTHON銀行機器學習:回歸、隨機森林、KNN近鄰、決策樹、高斯樸素貝葉斯、支持向量機SVM分析營銷活動數(shù)據(jù)|數(shù)據(jù)分享
PYTHON用戶流失數(shù)據(jù)挖掘：建立邏輯回歸、XGBOOST、隨機森林、決策樹、支持向量機、樸素貝葉斯和KMEANS聚類用戶畫像
MATLAB隨機森林優(yōu)化貝葉斯預測分析汽車燃油經(jīng)濟性
R語言中貝葉斯網(wǎng)絡（BN）、動態(tài)貝葉斯網(wǎng)絡、線性模型分析錯頜畸形數(shù)據(jù)
使用貝葉斯層次模型進行空間數(shù)據(jù)分析
MCMC的rstan貝葉斯回歸模型和標準線性回歸模型比較
python貝葉斯隨機過程：馬爾可夫鏈Markov-Chain，MC和Metropolis-Hastings，MH采樣算法可視化
Python貝葉斯推斷Metropolis-Hastings（M-H）MCMC采樣算法的實現(xiàn)
matlab貝葉斯隱馬爾可夫hmm模型實現(xiàn)
貝葉斯線性回歸和多元線性回歸構建工資預測模型
Metropolis Hastings采樣和貝葉斯泊松回歸Poisson模型
貝葉斯分位數(shù)回歸、lasso和自適應lasso貝葉斯分位數(shù)回歸分析免疫球蛋白、前列腺癌數(shù)據(jù)
R語言RSTAN MCMC：NUTS采樣算法用LASSO 構建貝葉斯線性回歸模型分析職業(yè)聲望數(shù)據(jù)
R語言STAN貝葉斯線性回歸模型分析氣候變化影響北半球海冰范圍和可視化檢查模型收斂性
PYTHON用戶流失數(shù)據(jù)挖掘：建立邏輯回歸、XGBOOST、隨機森林、決策樹、支持向量機、樸素貝葉斯和KMEANS聚類用戶畫像
貝葉斯分位數(shù)回歸、lasso和自適應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è)足球比賽進球數(shù)
隨機森林優(yōu)化貝葉斯預測分析汽車燃油經(jīng)濟性
R語言邏輯回歸、Naive Bayes貝葉斯、決策樹、隨機森林算法預測心臟病
R語言用Rcpp加速Metropolis-Hastings抽樣估計貝葉斯邏輯回歸模型的參數(shù)
R語言中的block Gibbs吉布斯采樣貝葉斯多元線性回歸
Python貝葉斯回歸分析住房負擔能力數(shù)據(jù)集
R語言實現(xiàn)貝葉斯分位數(shù)回歸、lasso和自適應lasso貝葉斯分位數(shù)回歸分析
Python用PyMC3實現(xiàn)貝葉斯線性回歸模型
R語言用WinBUGS 軟件對學術能力測驗建立層次（分層）貝葉斯模型
R語言Gibbs抽樣的貝葉斯簡單線性回歸仿真分析
R語言和STAN,JAGS：用RSTAN,RJAG建立貝葉斯多元線性回歸預測選舉數(shù)據(jù)
R語言基于copula的貝葉斯分層混合模型的診斷準確性研究
R語言貝葉斯線性回歸和多元線性回歸構建工資預測模型
R語言貝葉斯推斷與MCMC：實現(xiàn)Metropolis-Hastings 采樣算法示例
R語言stan進行基于貝葉斯推斷的回歸模型
R語言中RStan貝葉斯層次模型分析示例
R語言使用Metropolis-Hastings采樣算法自適應貝葉斯估計與可視化
R語言隨機搜索變量選擇SSVS估計貝葉斯向量自回歸（BVAR）模型
WinBUGS對多元隨機波動率模型：貝葉斯估計與模型比較
R語言實現(xiàn)MCMC中的Metropolis–Hastings算法與吉布斯采樣
R語言貝葉斯推斷與MCMC：實現(xiàn)Metropolis-Hastings 采樣算法示例
R語言使用Metropolis-Hastings采樣算法自適應貝葉斯估計與可視化
視頻：R語言中的Stan概率編程MCMC采樣的貝葉斯模型
R語言MCMC:Metropolis-Hastings采樣用于回歸的貝葉斯估計
