拓端tecdat：R語言用貝葉斯線性回歸、貝葉斯模型平均 (BMA)來預測工人工資

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

原文出處：拓端數據部落公眾號

背景

貝葉斯模型提供了變量選擇技術，確保變量選擇的可靠性。對社會經濟因素如何影響收入和工資的研究為應用這些技術提供了充分的機會，同時也為從性別歧視到高等教育的好處等主題提供了洞察力。下面，貝葉斯信息准則（BIC）和貝葉斯模型平均法被應用於構建一個簡明的收入預測模型。

 這些數據是從 935 名受訪者的隨機樣本中收集的。該數據集是計量經濟學數據集系列的一部分 。

加載包

數據將首先使用該dplyr 包進行探索 ，並使用該ggplot2 包進行可視化 。稍后，實現逐步貝葉斯線性回歸和貝葉斯模型平均 (BMA)。

數據

數據集網頁提供了以下變量描述表：

變量	描述
wage	每周收入（元）
hours	每周平均工作時間
IQ	智商分數
kww	對世界工作的了解得分
educ	受教育年數
exper	多年工作經驗
tenure	在現任雇主工作的年數
age	年齡
married	=1 如果已婚
black	=1 如果是黑人
south	=1 如果住在南方
urban	=1 如果居住在都市
sibs	兄弟姐妹的數量
brthord	出生順序
meduc	母親的教育（年）
feduc	父親的教育（年）
lwage	工資自然對數 wage

探索數據

與任何新數據集一樣，一個好的起點是標准的探索性數據分析。匯總表是簡單的第一步。

1.  
   # 數據集中所有變量的匯總表--包括連續變量和分類變量
2.  
   summary(wage)

因變量（工資）的直方圖給出了合理預測應該是什么樣子的。

1.  
   #工資數據的簡單柱狀圖
2.  
   hst(wge$wae, breks = 30)

直方圖還可用於大致了解哪些地方不太可能出現結果。

1.  
   # 檢查圖表 "尾部 "的點的數量
2.  
   sm(wage$ge < 300)

## [1] 6

sm(wae$wge > 2000)

## [1] 20

簡單線性回歸

由於周工資（'wage'）是該分析中的因變量，我們想探索其他變量作為預測變量的關系。我們在數據中看到的工資變化的一種可能的、簡單的解釋是更聰明的人賺更多的錢。下圖顯示了每周工資和 IQ 分數之間的散點圖。

gplot(wae, es(iq, wge)) + gom_oint() +gom_smoth() 

IQ 分數和工資之間似乎存在輕微的正線性關系，但僅靠 IQ 並不能可靠地預測工資。盡管如此，這種關系可以通過擬合一個簡單的線性回歸來量化，它給出：

工資 i = α + β⋅iqi + ϵiwagei = α + β⋅iqi + ϵi

1.  
   m_wg_iq = lm(wge ~ iq, dta = age)
2.  
   coefients

工資 i = 116.99 + 8.3 ⋅iqi + ϵiwagei = 116.99 + 8.3 ⋅iqi + ϵi

在轉向貝葉斯改進這個模型之前，請注意貝葉斯建模假設誤差 (ϵi) 以恆定方差正態分布。通過檢查模型的殘差分布來檢查該假設。如果殘差高度非正態或偏斜，則違反假設並且任何后續推論都無效。要檢查假設，請按如下方式繪制殘差：

1.  
   # 用散點圖和模型誤差殘差的直方圖來檢查正態性假設
2.  
    
3.  
   glot(dta = mwag_q, es(x = .ite, y = .rd)) +
4.  
   gemittr() +
5.  
    

1.  
   plot(dta = m_g_iq, aes(x = .reid)) +
2.  
   histgm(bnwth = 10)
3.  
    

變量變換

兩個圖都顯示殘差是右偏的。因此，IQ（因為它目前存在於數據集中）不應用作貝葉斯預測模型。但是，對 僅具有正值的偏斜因變量使用（自然）對數變換 通常可以解決問題。下面，該模型使用轉換后的工資變量進行了重新擬合。

1.  
   # 用IQ的自然對數擬合th模型
2.  
   lm(lage ~ iq, data = wae)

1.  
   # 殘差sctterplot和轉換后數據的柱狀圖
2.  
   plt(data = m_lag_iq, es(x = .fited, y = .reid))
3.  
   geiter() +

1.  
   ggpot(dta = m_lwgeiq, as(x = .resd)) +
2.  
   gostgam(binwth = .1) +
3.  
    

殘差確實大致呈正態分布。然而，由此產生的 IQ 系數非常小（只有 0.0088），這是可以預料的，因為 IQ 分數提高 1 分幾乎不會對工資產生太大影響。需要進一步細化。數據集包含更多信息。

多元線性回歸和 BIC

我們可以首先在回歸模型中包含所有潛在的解釋變量，來粗略地嘗試解釋盡可能多的工資變化。

1.  
   # 對數據集中的所有變量運行一個線性模型，使用'.'約定。
2.  
   full = lm(lwge ~ . - wage, dta = wge)

完整線性模型的上述總結表明，自變量的許多系數在統計上並不顯着（請參閱第 4 個數字列中的 p 值）。選擇模型變量的一種方法是使用貝葉斯信息准則 (BIC)。BIC 是模型擬合的數值評估，它也會按樣本大小的比例懲罰更多的參數。這是完整線性模型的 BIC：

BIC(full)

BIC 值越小表示擬合越好。因此，BIC 可以針對各種縮減模型進行計算，然后與完整模型 BIC 進行比較，以找到適合工資預測工作的最佳模型。當然，R 有一個功能可以系統地執行這些 BIC 調整。

1.  
   # 用step計算模型
2.  
   pIC(lwge ~ . - wge, dta = na.oi(wge))lg(lgth(na.mit(wge))))
3.  
    

1.  
   # 顯示逐步模型的BIC
2.  
   BIC(se_mol)

調用 step找到產生最低 BIC 的變量組合，並提供它們的系數。很不錯。

貝葉斯模型平均（BMA）

即使BIC處於最低值，我們能有多大把握確定所得到的模型是真正的 "最佳擬合"？答案很可能取決於基礎數據的規模和穩定性。在這些不確定的時候，貝葉斯模型平均化（BMA）是有幫助的。BMA對多個模型進行平均化，獲得系數的后驗值和新數據的預測值。下面，BMA被應用於工資數據（排除NA值后）。

 

1.  
   # 不包括NA
2.  
   a_ona = na.omt(wae)
3.  
    
4.  
   # 運行BMA，指定BIC作為判斷結果模型的標准
5.  
   BMA(wge ~ . -wge, daa= ae_o_a,
6.  
   pror = "BIC",
7.  
   moepor = ufom())
8.  
    
9.  
   # 顯示結果
10.  
    summary

 

結果表顯示了五個最有可能的模型，以及每個系數被包含在真實模型中的概率。我們看到，出生順序和是否有兄弟姐妹是最不可能被包含的變量，而教育和智商變量則被鎖定。BMA模型的排名也可以用圖像圖來顯示，它清楚地顯示哪些變量在所有模型中，哪些變量被排除在所有模型之外，以及那些介於兩者之間的變量。

 

ge(b_lge, tp.oels)

我們還可以提供模型系數的95%置信區間。下面的結果支持了關於包括或排除系數的決定。例如，在區間包含零，有大量證據支持排除該變量。

confint(ceflae)

進行預測

構建模型后，pediction 只是插入數據的問題：

1.  
   # 用一個虛構的工人的統計資料來預測數據的例子
2.  
    
3.  
   # 進行預測
4.  
   redict = pedct(e_odl, newdt = wrkr,eitr = "BMA")
5.  
    
6.  
   # 將結果轉換為元
7.  
   exp(wk_pedct)

預計這名化妝工作人員的周薪為 745 元。這到底有多准確？你得問她，但我們對我們的變量選擇很有信心，並對現有的數據盡了最大努力。應用的貝葉斯技術使我們對結果有信心。

---

最受歡迎的見解

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模型實現